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
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
|
// horizon_Jacobian.cc -- evaluate Jacobian matrix of LHS function H(h)
// $Id$
//
// <<<prototypes for functions local to this file>>>
// horizon_Jacobian - top-level driver
//
#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 "gfn.hh"
#include "Jacobian.hh"
#include "AHFinderDirect.hh"
//******************************************************************************
//
// ***** prototypes for functions local to this file *****
//
namespace {
void add_molecule_point_to_Jacobian(const patch_system& ps, const patch& xp,
int x_irho, int x_isigma,
int xm_irho, int xm_isigma,
fp mol);
}
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// This function computes the Jacobian matrix of the LHS function H(h)
// from the Jacobian coefficient (angular) gridfns. The computation is
// done a Jacobian row at a time, using equation (25) of my 1996 apparent
// horizon finding paper.
//
// Inputs (angular gridfns, on ghosted grid):
// h # shape of trial surface
// partial_H_wrt_partial_d_h, # Jacobian coefficients
// partial_H_wrt_partial_dd_h,
//
// Outputs:
// The Jacobian matrix is stored in the Jacobian object Jac.
void horizon_Jacobian(const patch_system& ps,
Jacobian& Jac);
{
CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian");
ps.compute_synchronize_Jacobian();
Jac.zero();
for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
{
patch& xp = ps.ith_patch(xpn);
for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho)
{
for (int x_isigma = xp.min_isigma() ;
x_isigma <= xp.max_isigma() ;
++x_isigma)
{
//
// compute the Jacobian row for this grid point, i.e.
// partial H(this point x, Jacobian row II)
// ---------------------------------------------
// partial h(other points y, Jacobian column JJ)
//
// FIXME FIXME: we still have to take into account the
// position-dependence of the coefficients,
// cf the difference between J[3H(h)] and J[2H(h)];
// here we're sort of computing the former, but
// Newton's method really wants the latter
//
// Jacobian coefficients for this point
const fp Jacobian_coeff_rho
= xp.gridfn(gfn__partial_H_wrt_partial_d_h_1, x_irho,x_isigma);
const fp Jacobian_coeff_sigma
= xp.gridfn(gfn__partial_H_wrt_partial_d_h_2, x_irho,x_isigma);
const fp Jacobian_coeff_rho_rho
= xp.gridfn(gfn__partial_H_wrt_partial_dd_h_11, x_irho,x_isigma);
const fp Jacobian_coeff_rho_sigma
= xp.gridfn(gfn__partial_H_wrt_partial_dd_h_12, x_irho,x_isigma);
const fp Jacobian_coeff_sigma_sigma
= xp.gridfn(gfn__partial_H_wrt_partial_dd_h_22, x_irho,x_isigma);
// partial_rho, partial_rho_rho
for (int m_irho = xp.molecule_min_m() ;
m_irho <= xp.molecule_max_m() ;
++m_irho)
{
const int xm_irho = x_irho + m_irho;
const fp Jac_rho = Jacobian_coeff_rho
* xp.partial_rho_coeff(m_rho);
const fp Jac_rho_rho = Jacobian_coeff_rho_rho
* xp.partial_rho_rho_coeff(m_rho);
add_molecule_point_to_Jacobian(ps, xp,
x_irho, x_isigma,
xm_irho, x_isigma,
Jac_rho + Jac_rho_rho);
}
// partial_sigma, partial_sigma_sigma
for (int m_isigma = xp.molecule_min_m() ;
m_isigma <= xp.molecule_max_m() ;
++m_isigma)
{
const int xm_isigma = x_isigma + m_isigma;
const fp Jac_sigma = Jacobian_coeff_sigma
* xp.partial_sigma_coeff(m_sigma);
const fp Jac_sigma_sigma = Jacobian_coeff_sigma_sigma
* xp.partial_sigma_sigma_coeff(m_sigma);
add_molecule_point_to_Jacobian(ps, xp,
x_irho, x_isigma,
x_irho, xm_isigma,
Jac_sigma + Jac_sigma_sigma);
}
// partial_rho_sigma
for (int m_irho = xp.molecule_min_m() ;
m_irho <= xp.molecule_max_m() ;
++m_irho)
{
for (int m_isigma = xp.molecule_min_m() ;
m_isigma <= xp.molecule_max_m() ;
++m_isigma)
{
const int xm_irho = x_irho + m_irho;
const int xm_isigma = x_isigma + m_isigma;
const fp Jac_rho_sigma
= Jacobian_coeff_rho_sigma
* xp.partial_rho_sigma_coeff(m_rho, m_sigma);
add_molecule_point_to_Jacobian(ps, xp,
x_irho, x_isigma,
xm_irho, xm_isigma,
Jac_rho_sigma);
}
}
}
}
}
}
//******************************************************************************
//
// This function adds all the elements for a single molecule point to
// a Jacobian matrix, including any contributions from other patches if
// the specified point lies in a ghost zone.
//
// Arguments:
// ps = The patch system.
// xp = The patch containing the center point of the molecule.
// x_(irho,isigma) = The coordinates of the center point of the molecule.
// xm_(irho,isigma) = The coordinates of the x+m point of the molecule.
// mol = The molecule coefficient.
//
namespace {
void add_molecule_point_to_Jacobian(const patch_system& ps, const patch& xp,
int x_irho, int x_isigma,
int xm_irho, int xm_isigma,
fp mol)
{
// Jacobian row
const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho, x_isigma);
if (xp.is_in_nominal_grid(xm_irho, xm_isigma))
then {
const int JJ = ps.gpn_of_patch_irho_isigma(xp, xm_irho, xm_isigma);
Jac.add_to_element(II,JJ, Jac_coeff);
}
else {
const ghost_zone& xmgz
= xp.ghost_zone_containing_point(xm_irho, xm_isigma);
const patch_edge& xme = xmgz.my_edge();
Const int xm_iperp = xme.iperp_of_irho_isigma(xm_irho, xm_isigma);
const int xm_ipar = xme. ipar_of_irho_isigma(xm_irho, xm_isigma);
// on what other points ym does this molecule point xm depend
// via the patch_system::synchronize() operation?
const patch& ymp = ps.synchronize_Jacobian_y_patch(xmgz);
const patch_edge& yme = ps.synchronize_Jacobian_y_edge (xmgz);
const int min_ym_ipar_m = ps.synchronize_Jacobian_min_y_ipar_m(xmgz);
const int max_ym_ipar_m = ps.synchronize_Jacobian_max_y_ipar_m(xmgz);
const int ym_iperp = ps.synchronize_Jacobian_y_iperp(xmgz, xm_iperp);
const int ym_ipar_posn
= ps.synchronize_Jacobian_y_ipar_posn(xmgz, xm_iperp, xm_ipar);
for (int ym_ipar_m = min_ym_ipar_m ;
ym_ipar_m <= max_ym_ipar_m ;
++ym_ipar_m)
{
const int ym_ipar = ym_ipar_posn + ym_ipar_m;
const int ym_irho = yme. irho_of_iperp_ipar(ym_iperp,ym_ipar);
const int ym_isigma = yme.isigma_of_iperp_ipar(ym_iperp,ym_ipar);
const int JJ = ps.gpn_of_patch_irho_isigma(ymp, ym_irho, ym_isigma);
const fp sync_Jac = ps.synchronize_Jacobian(xmgz, xm_iperp, xm_ipar,
ym_ipar_m);
Jac.add_to_element(II,JJ, mol*sync_Jac);
}
}
}
}
|
