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
|
// 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 "jt/stdc.h"
#include "jt/util.hh"
#include "jt/array.hh"
#include "jt/cpm_map.hh"
#include "jt/linear_map.hh"
using jtutil::error_exit;
#include "../config.hh"
#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 {
};
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// 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");
Jac.zero(II);
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 row index for this point
const int II = ps.gpn_of_patch_irho_isigma(xp,x_irho,x_isigma);
// 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 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);
const fp Jac_sum = Jac_rho + Jac_rho_rho;
if (xp.is_in_nominal_grid(x_irho+m_irho, x_isigma))
then {
const int JJ
= ps.gpn_of_patch_irho_isigma(xp, x_irho+m_irho,
x_isigma);
Jac.add_to_element(II,JJ, Jac_sum);
}
else {
const patch_edge& xe = (m_irho < 0)
? xp.min_rho_patch_edge()
: xp.max_rho_patch_edge()
const ghost_zone& gz = xp.ghost_zone_on_edge(xe);
const int x_iperp = xe.iperp_of_irho_isigma(x_irho,
x_isigma);
const int x_ipar = xe.ipar_of_irho_isigma(x_irho, x_isigma);
gz.compute_Jacobian(gfn__h, gfn__h);
const patch& yp = gz.Jacobian_y_patch();
const patch_edge& ye = gz.Jacobian_y_edge();
const int min_y_ipar_m = gz.Jacobian_min_y_ipar_m();
const int max_y_ipar_m = gz.Jacobian_max_y_ipar_m();
const int y_iperp = gz.Jacobian_y_iperp(x_iperp);
const int y_ipar_posn
= gz.Jacobian_y_ipar_posn(x_iperp, x_ipar);
for (int y_ipar_m = min_y_ipar_m ;
y_ipar_m <= max_y_ipar_m ;
++y_ipar_m)
{
const int y_ipar = y_ipar_posn + y_ipar_m;
const int y_irho
= ye.irho_of_iperp_ipar(y_iperp, y_ipar);
const int y_isigma
= ye.isigma_of_iperp_ipar(y_iperp, y_ipar);
const int JJ
= ps.gpn_of_patch_irho_isigma(yp, y_irho, y_isigma);
const fp gz_Jac = gz.Jacobian(x_iperp,x_ipar, y_ipar_m);
Jac.add_to_element(II,JJ, Jac_sum*gz_Jac);
}
}
}
}
}
}
}
|
