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
|
// 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__interpolate_from_Cactus_grid,
geometry__Schwarzschild_EF // no comma
};
//
// this enum holds the (a) decoded Jacobian_method parameter,
// i.e. it specifies how we compute the (a) Jacobian matrix
//
enum Jacobian_method
{
Jacobian_method__numerical_perturb,
Jacobian_method__symbolic_diff_with_FD_dr,
Jacobian_method__symbolic_diff // no comma
};
//******************************************************************************
//
// This structure holds all the information we need about the Cactus grid
// and gridfns. At present we use the CCTK_InterpLocalUniform() local
// interpolator to access the Cactus gridfns. Alas, much of the information
// in this structure is specific to that API, and (alas) will need changing
// when we eventually switch to the CCTK_InterpGridArrays() API.
//
struct cactus_grid_info
{
cGH *GH; // --> Cactus grid hierarchy
// 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];
// 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;
// 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
// --> 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_method Jacobian_method;
enum Jacobian_type Jacobian_storage_method;
fp perturbation_amplitude;
};
//******************************************************************************
//
// prototypes for functions visible outside their source files
//
// horizon_function.cc
// ... returns true for successful computation,
// ... returns false for failure (eg horizon outside grid);
// if (print_msg_flag) then also reports CCTK_VWarn() at the
// appropriate warning level
namespace horizon_function_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 horizon_function(patch_system& ps,
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);
enum geometry_method
decode_geometry_method(const char geometry_method_string[]);
// horizon_Jacobian.cc
// ... returns true for successful computation, false for failure
bool horizon_Jacobian(patch_system& ps,
Jacobian& Jac,
const struct cactus_grid_info& cgi,
const struct geometry_info& gi,
const struct Jacobian_info& Jacobian_info,
bool print_msg_flag);
enum Jacobian_method
decode_Jacobian_method(const char Jacobian_method_string[]);
// Schwarzschild_EF.cc
void Schwarzschild_EF_geometry(patch_system& ps,
const struct geometry_info& gi,
bool msg_flag);
|
