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
|
// gr.hh -- header file for general relativity code
// $Header$
//******************************************************************************
//
// number of spatial dimensions in the main Cactus grid
// and in our trial-horizon-surface grid
//
enum { N_GRID_DIMS = 3, N_HORIZON_DIMS = 2 };
//
// this enum holds the (a) decoded geometry_method parameter,
// i.e. it specifies how we compute the (a) geometry
//
enum geometry_method
{
geometry__global_interp_from_Cactus_grid,
geometry__local_interp_from_Cactus_grid,
geometry__Schwarzschild_EF // no comma
};
//
// this enum holds the (a) decoded Jacobian_compute_method parameter,
// i.e. it specifies how we compute the (a) Jacobian matrix
//
enum Jacobian_compute_method
{
Jacobian__numerical_perturbation,
Jacobian__symbolic_diff_with_FD_dr,
Jacobian__symbolic_diff // no comma
};
//******************************************************************************
//
// This structure holds all the information we need about the Cactus grid
// and gridfns in order to interpolate the geometry information.
//
struct cactus_grid_info
{
cGH *GH; // --> Cactus grid hierarchy
// this describes the semantics of the Cactus gij gridfns:
// true ==> the Cactus gij are conformal values as per
// CactusEinstein/StaticConformal and the psi gridfn
// false ==> the Cactus gij are the physical metric
bool Cactus_conformal_metric;
//
// stuff for doing a global interpolation via CCTK_InterpGridArrays()
// i.e. geometry_info.geometry_method
// == geometry__global_interp_from_Cactus_grid
//
// Cactus coordinate system
int coord_system_handle;
// Cactus variable indices
int g_dd_11_varindex, g_dd_12_varindex, g_dd_13_varindex,
g_dd_22_varindex, g_dd_23_varindex,
g_dd_33_varindex;
int K_dd_11_varindex, K_dd_12_varindex, K_dd_13_varindex,
K_dd_22_varindex, K_dd_23_varindex,
K_dd_33_varindex;
int psi_varindex; // unused if !Cactus_conformal_metric
//
// stuff for doing a local interpolation via CCTK_InterpLocalUniform()
// i.e. geometry_info.geometry_method
// == geometry__local_interp_from_Cactus_grid
//
// Cactus coordinate system
fp coord_origin[N_GRID_DIMS]; // (x,y,z) of grid posn (0,0,0)
fp coord_delta[N_GRID_DIMS]; // (x,y,z) grid spacing
// dimensions of gridfn data, viewed as a 3-D array
// n.b. storage ordering is Fortran,
// i.e. i is contiguous, j has stride Ni, k has stride Ni*Nj
CCTK_INT gridfn_dims[N_GRID_DIMS];
// --> Cactus gridfn data for grid posn (0,0,0)
const fp* g_dd_11_data;
const fp* g_dd_12_data;
const fp* g_dd_13_data;
const fp* g_dd_22_data;
const fp* g_dd_23_data;
const fp* g_dd_33_data;
const fp* K_dd_11_data;
const fp* K_dd_12_data;
const fp* K_dd_13_data;
const fp* K_dd_22_data;
const fp* K_dd_23_data;
const fp* K_dd_33_data;
const fp* psi_data; // NULL if !Cactus_conformal_metric
};
//
// This structure holds information for computing the spacetime geometry.
// This is normally done by interpolating $g_{ij}$ and $K_{ij}$ from the
// usual Cactus grid, but can optionally instead by done by hard-wiring
// the Schwarzschild geometry in Eddington-Finkelstein coordinates.
//
struct geometry_info
{
enum geometry_method geometry_method;
// parameters for geometry_method = interpolate_from_Cactus_grid
int operator_handle; // Cactus handle to interpolation op
int param_table_handle; // Cactus handle to parameter table
// for the interpolator
// parameters for geometry_method = Schwarzschild_EF
fp geometry__Schwarzschild_EF__x_posn; // x posn of Schwarzschild BH
fp geometry__Schwarzschild_EF__y_posn; // y posn of Schwarzschild BH
fp geometry__Schwarzschild_EF__z_posn; // z posn of Schwarzschild BH
fp geometry__Schwarzschild_EF__mass; // mass of Schwarzschild BH
fp geometry__Schwarzschild_EF__epsilon; // use z axis limits if
// (x^2+y^2)/r^2 <= this
fp geometry__Schwarzschild_EF__Delta_xyz;// "grid spacing" for FD
// computation of partial_k g_ij
// should we check various gridfns for NaNs?
bool check_that_h_is_finite;
bool check_that_geometry_is_finite;
};
//
// This struct holds parameters for computing the Jacobian matrix.
//
struct Jacobian_info
{
enum Jacobian_compute_method Jacobian_compute_method;
enum Jacobian_store_solve_method Jacobian_store_solve_method;
fp perturbation_amplitude;
};
//******************************************************************************
//
// prototypes for functions visible outside their source files
//
// expansion.cc
// ... returns true for successful computation,
// ... returns true for dummy computation (ps_ptr == NULL)
// ... returns false for failure (eg horizon outside grid);
// if (print_msg_flag) then also reports CCTK_VWarn() at the
// appropriate warning level
namespace expansion_warning_levels
{
enum {
failure = 2, // point outside grid etc
nonfinite = 2, // nonfinite gridfn detected
skip_nonfinite = 3 // no finite() function available
// ==> skipping nonfinite check
};
}
bool expansion(patch_system* ps_ptr,
const struct cactus_grid_info& cgi,
const struct geometry_info& gi,
bool Jacobian_flag = false,
bool print_msg_flag = false,
jtutil::norm<fp>* H_norms_ptr = NULL);
// expansion_Jacobian.cc
// ... returns true for successful computation
// ... returns true for dummy computation (ps_ptr == NULL && Jac_ptr == NULL)
// ... returns false for failure
bool expansion_Jacobian(patch_system* ps_ptr,
Jacobian* Jac_ptr,
const struct cactus_grid_info& cgi,
const struct geometry_info& gi,
const struct Jacobian_info& Jacobian_info,
bool print_msg_flag = false);
// Schwarzschild_EF.cc
void Schwarzschild_EF_geometry(patch_system& ps,
const struct geometry_info& gi,
bool msg_flag);
// misc.cc
enum geometry_method
decode_geometry_method(const char geometry_method_string[]);
#ifdef NOT_USED
bool geometry_method_is_interp(enum geometry_method geometry_method);
#endif
enum Jacobian_compute_method
decode_Jacobian_compute_method(const char Jacobian_compute_method_string[]);
|
