1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
|
// Newton.cc -- solve H(h) = 0 via Newton's method
// $Id$
//
// <<<prototypes for functions local to this file>>>
// Newton_solve - driver to solve H(h) = 0 via Newton's method
//
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <vector>
#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"
#include "stdc.h"
#include "config.hh"
#include "../jtutil/util.hh"
#include "../jtutil/array.hh"
#include "../jtutil/cpm_map.hh"
#include "../jtutil/linear_map.hh"
using jtutil::error_exit;
#include "../util/coords.hh"
#include "../util/grid.hh"
#include "../util/fd_grid.hh"
#include "../util/patch.hh"
#include "../util/patch_edge.hh"
#include "../util/patch_interp.hh"
#include "../util/ghost_zone.hh"
#include "../util/patch_system.hh"
#include "../elliptic/Jacobian.hh"
#include "gfns.hh"
#include "AHFinderDirect.hh"
//******************************************************************************
//
// ***** prototypes for functions local to this file *****
//
namespace {
}
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// This function solves H(h) = 0 for h via Newton's method.
//
bool Newton_solve(patch_system& ps,
const struct cactus_grid_info& cgi,
const struct geometry_info& gi,
const struct Jacobian_info& Jacobian_info,
const struct solver_info& solver_info,
Jacobian& Jac)
{
for (int iteration = 1 ;
iteration <= solver_info.max_Newton_iterations ;
++iteration)
{
CCTK_VInfo(CCTK_THORNSTRING,
"Newton iteration %d", iteration);
jtutil::norm<fp> H_norms;
horizon_function(ps, cgi, gi, true, true, &H_norms);
CCTK_VInfo(CCTK_THORNSTRING,
" H rms-norm=%.2e, infinity-norm=%.2e",
H_norms.rms_norm(), H_norms.infinity_norm());
if (solver_info.output_h_and_H_at_each_Newton_iteration)
then {
const char N_buffer = 100;
char file_name_buffer[100];
snprintf(file_name_buffer, N_buffer,
"%s.it%d",
solver_info.h_file_name, iteration);
CCTK_VInfo(CCTK_THORNSTRING,
" writing h to \"%s\"",
file_name_buffer);
ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h,
true, gfns::gfn__h,
file_name_buffer,
false); // no ghost zones
snprintf(file_name_buffer, N_buffer,
"%s.it%d",
solver_info.H_of_h_file_name, iteration);
CCTK_VInfo(CCTK_THORNSTRING,
" writing H(h) to \"%s\"",
file_name_buffer);
ps.print_gridfn_with_xyz(gfns::gfn__H,
true, gfns::gfn__h,
file_name_buffer);
}
if (H_norms.infinity_norm() <= solver_info.H_norm_for_convergence)
then {
CCTK_VInfo(CCTK_THORNSTRING,
"==> finished (||H|| < %g)",
double(solver_info.H_norm_for_convergence));
return true; // *** NORMAL RETURN ***
}
// compute the Newton step
horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac);
CCTK_VInfo(CCTK_THORNSTRING,
"solving linear system for %d unknowns Delta_h",
Jac.NN());
const fp rcond = Jac.solve_linear_system(gfns::gfn__H,
gfns::gfn__Delta_h);
if (rcond == 0.0)
then {
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"Newton_solve: Jacobian matrix is numerically singular!");
return false; // *** ERROR RETURN ***
}
if (rcond > 0.0)
then CCTK_VInfo(CCTK_THORNSTRING,
"done solving linear system (condition number %.2e)",
double(rcond));
else CCTK_VInfo(CCTK_THORNSTRING,
"done solving linear system");
// take the Newton step
jtutil::norm<fp> Delta_h_norms;
for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
{
patch& p = ps.ith_patch(pn);
for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
{
for (int isigma = p.min_isigma() ;
isigma <= p.max_isigma() ;
++isigma)
{
const fp Delta_h = p.gridfn(gfns::gfn__Delta_h, irho,isigma);
Delta_h_norms.data(Delta_h);
p.ghosted_gridfn(gfns::gfn__h, irho,isigma) -= Delta_h;
}
}
}
CCTK_VInfo(CCTK_THORNSTRING,
"h += Delta_h (rms-norm=%.2e, infinity-norm=%.2e)",
Delta_h_norms.rms_norm(), Delta_h_norms.infinity_norm());
if (Delta_h_norms.infinity_norm() <= solver_info
.Delta_h_norm_for_convergence)
then {
CCTK_VInfo(CCTK_THORNSTRING,
"==> finished (||Delta_h|| < %g)",
double(solver_info.Delta_h_norm_for_convergence));
if (solver_info.final_H_update_if_exit_x_H_small)
then {
CCTK_VInfo(CCTK_THORNSTRING,
" doing final H(h) evaluation");
horizon_function(ps, cgi, gi, false, true);
}
return true; // *** NORMAL RETURN ***
}
}
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"Newton_solve: no convergence in %d iterations!",
solver_info.max_Newton_iterations);
if (solver_info.final_H_update_if_exit_x_H_small)
then {
CCTK_VInfo(CCTK_THORNSTRING,
" doing final H(h) evaluation");
horizon_function(ps, cgi, gi, false, true);
}
return false; // *** ERROR RETURN ***
}
|
