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
|
// 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 "gfn.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_interpolator_info& gii,
const struct Jacobian_info& Jac_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, gii, true, &H_norms);
CCTK_VInfo(CCTK_THORNSTRING,
" H rms-norm=%.2e, infinity-norm=%.2e",
H_norms.rms_norm(), H_norms.infinity_norm());
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 ***
}
// take a Newton step
horizon_Jacobian(ps, cgi, gii, Jac_info, Jac);
CCTK_VInfo(CCTK_THORNSTRING,
"solving linear system for %d unknowns", Jac.NN());
Jac.solve_linear_system(nominal_gfns::gfn__H, // rhs gridfn
nominal_gfns::gfn__Delta_h); // soln gridfn
CCTK_VInfo(CCTK_THORNSTRING,
"done solving linear system");
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(nominal_gfns::gfn__Delta_h,
irho,isigma);
Delta_h_norms.data(Delta_h);
p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) -= Delta_h;
}
}
}
CCTK_VInfo(CCTK_THORNSTRING,
"moved h by Delta_h (rhs-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));
return true; // *** NORMAL RETURN ***
}
}
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"Newton_solve: no convergence in %d iterations!\n",
solver_info.max_Newton_iterations);
return false; // *** ERROR RETURN ***
}
|
