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
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
|
// spherical_surface.cc -- interface with SphericalSurface thorn
// $Header$
//
// <<<access to persistent data>>>
// AHFinderDirect_store_SS_info - ... SphericalSurface information
/// store_diagnostic_info - store BH_diagnostics info in SphericalSurface vars
/// store_horizon_shape - store horizon shape in SphericalSurface vars
//
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "config.h"
#include "stdc.h"
#include "../jtutil/util.hh"
#include "../jtutil/array.hh"
#include "../jtutil/cpm_map.hh"
#include "../jtutil/linear_map.hh"
#include "../patch/coords.hh"
#include "../patch/grid.hh"
#include "../patch/fd_grid.hh"
#include "../patch/patch.hh"
#include "../patch/patch_edge.hh"
#include "../patch/patch_interp.hh"
#include "../patch/ghost_zone.hh"
#include "../patch/patch_system.hh"
#include "../elliptic/Jacobian.hh"
#include "../gr/gfns.hh"
#include "../gr/gr.hh"
#include "horizon_sequence.hh"
#include "BH_diagnostics.hh"
#include "driver.hh"
// all the code in this file is inside this namespace
namespace AHFinderDirect
{
using jtutil::error_exit;
//******************************************************************************
//
// ***** prototypes for functions local to this file *****
//
namespace {
void store_diagnostic_info(CCTK_ARGUMENTS,
const patch_system& ps,
const struct BH_diagnostics& BH_diagnostics,
const int sn);
void store_horizon_shape(CCTK_ARGUMENTS,
const patch_system& ps,
const int sn);
}
//******************************************************************************
//
// ***** access to persistent data *****
//
extern struct state state;
//******************************************************************************
//
// This function is called by the Cactus scheduler, to store any desired
// apparent horizon info to the SphericalSurface variables.
//
// Cactus parameters used:
// SphericalSurface::nsurfaces = (in)
//
extern "C"
void AHFinderDirect_store_SS_info(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
const int N_surfaces = /* SphericalSurface:: */ nsurfaces;
const struct verbose_info& verbose_info = state.verbose_info;
for (int hn = 1; hn <= N_horizons; ++ hn)
{
assert(state.AH_data_array[hn] != NULL);
const struct AH_data& AH_data = *state.AH_data_array[hn];
assert(AH_data.ps_ptr != NULL);
const patch_system& ps = *AH_data.ps_ptr;
if (! AH_data.store_info_in_SS_vars)
then continue; // *** LOOP CONTROL ***
const int sn = AH_data.SS_surface_number;
assert(sn >= 0);
assert(sn < N_surfaces);
if (! state.find_now(cctk_iteration))
then {
// we didn't try to find this (or any!) horizon
// at this time step
/* SphericalSurface:: */ sf_valid[sn] = 0;
continue; // *** LOOP CONTROL ***
}
if (! AH_data.found_flag)
then {
// we tried to find this horizon, but failed
/* SphericalSurface:: */ sf_valid[sn] = -1;
continue; // *** LOOP CONTROL ***
}
// get to here ==> we found this horizon
/* SphericalSurface:: */ sf_valid[sn] = 1;
if (verbose_info.print_algorithm_highlights)
then CCTK_VInfo(CCTK_THORNSTRING,
"storing horizon %d info in SphericalSurface surface %d",
hn, sn);
const struct BH_diagnostics& BH_diagnostics = AH_data.BH_diagnostics;
store_diagnostic_info(CCTK_PASS_CTOC, ps, BH_diagnostics, sn);
store_horizon_shape(CCTK_PASS_CTOC, ps, sn);
}
}
//******************************************************************************
//
// Assuming that we found this horizon, this function stores various
// diagnostic info about it in the corresponding SphericalSurface variables.
//
// Arguments:
// sn = (in) The SphericalSurface surface number to store the information in.
//
// Cactus variables:
// sf_* = (out) The SphericalSurface variables.
//
namespace {
void store_diagnostic_info(CCTK_ARGUMENTS,
const patch_system& ps,
const struct BH_diagnostics& BH_diagnostics,
const int sn)
{
DECLARE_CCTK_ARGUMENTS
/* SphericalSurface:: */ sf_origin_x[sn] = ps.origin_x();
/* SphericalSurface:: */ sf_origin_y[sn] = ps.origin_y();
/* SphericalSurface:: */ sf_origin_z[sn] = ps.origin_z();
/* SphericalSurface:: */ sf_mean_radius [sn] = BH_diagnostics.mean_radius;
/* SphericalSurface:: */ sf_min_radius [sn] = BH_diagnostics.min_radius;
/* SphericalSurface:: */ sf_max_radius [sn] = BH_diagnostics.max_radius;
/* SphericalSurface:: */ sf_area [sn] = BH_diagnostics.area;
/* SphericalSurface:: */ sf_centroid_x [sn] = BH_diagnostics.centroid_x;
/* SphericalSurface:: */ sf_centroid_y [sn] = BH_diagnostics.centroid_y;
/* SphericalSurface:: */ sf_centroid_z [sn] = BH_diagnostics.centroid_z;
/* SphericalSurface:: */ sf_quadrupole_xx[sn] = BH_diagnostics.quadrupole_xx;
/* SphericalSurface:: */ sf_quadrupole_xy[sn] = BH_diagnostics.quadrupole_xy;
/* SphericalSurface:: */ sf_quadrupole_xz[sn] = BH_diagnostics.quadrupole_xz;
/* SphericalSurface:: */ sf_quadrupole_yy[sn] = BH_diagnostics.quadrupole_yy;
/* SphericalSurface:: */ sf_quadrupole_yz[sn] = BH_diagnostics.quadrupole_yz;
/* SphericalSurface:: */ sf_quadrupole_zz[sn] = BH_diagnostics.quadrupole_zz;
/* SphericalSurface:: */ sf_min_x [sn] = BH_diagnostics.min_x;
/* SphericalSurface:: */ sf_max_x [sn] = BH_diagnostics.max_x;
/* SphericalSurface:: */ sf_min_y [sn] = BH_diagnostics.min_y;
/* SphericalSurface:: */ sf_max_y [sn] = BH_diagnostics.max_y;
/* SphericalSurface:: */ sf_min_z [sn] = BH_diagnostics.min_z;
/* SphericalSurface:: */ sf_max_z [sn] = BH_diagnostics.max_z;
}
}
//******************************************************************************
//
// This function stores our information about a specified horizon to the
// SphericalSurface variables.
//
// Arguments:
// sn = (in) The SphericalSurface surface number to store the information in.
//
// Cactus variables:
// sf_maxn{theta,phi} = (in) The SphericalSurface array dimensions for
// the 2-D array indexing.
// sf_n{theta,phi} = (in) The SphericalSurface array dimensions for the
// part of the 2-D array that's actually used.
// sf_radius = (out) The SphericalSurface radius.
//
// FIXME:
// * The present implementation is quite inefficient, as it calls
// patch_system::radius_in_local_xyz_direction() (which does a
// separate interpolator call) for each point. It would be more
// efficient to have a new patch_system:: API which operated
// on a whole array of points at once, to amortize the interpolator
// overhead.
// * It would be cleaner to abstract out the (complicated) indexing
// calculation for the SphericalSurface shape array into a separate
// access structure/function, the way we do with struct cactus_grid_info
// for grid functions.
//
namespace {
void store_horizon_shape(CCTK_ARGUMENTS,
const patch_system& ps,
const int sn)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
const double origin_theta = /* SphericalSurface:: */ sf_origin_theta[sn];
const double origin_phi = /* SphericalSurface:: */ sf_origin_phi [sn];
const double delta_theta = /* SphericalSurface:: */ sf_delta_theta [sn];
const double delta_phi = /* SphericalSurface:: */ sf_delta_phi [sn];
const int N_theta = /* SphericalSurface:: */ ntheta[sn];
const int N_phi = /* SphericalSurface:: */ nphi [sn];
const int max_N_theta = /* SphericalSurface:: */ maxntheta;
const int max_N_phi = /* SphericalSurface:: */ maxnphi;
// we want Fortran loop nesting here for cache efficiency in storing
// to the SphericalSurface::sf_radius array (see comment below)
for (int i_phi = 0 ; i_phi < N_phi ; ++i_phi)
{
const double phi = origin_phi + i_phi*delta_phi;
const double sin_phi = sin(phi);
const double cos_phi = cos(phi);
for (int i_theta = 0 ; i_theta < N_theta ; ++i_theta)
{
const double theta = origin_theta + i_theta*delta_theta;
const double sin_theta = sin(theta);
const double cos_theta = cos(theta);
const double local_x = sin_theta * cos_phi;
const double local_y = sin_theta * sin_phi;
const double local_z = cos_theta;
const double local_r = ps.radius_in_local_xyz_direction
(gfns::gfn__h,
local_x, local_y, local_z);
// SphericalSurface::sf_radius is actually stored as
// a 3-D contiguous array, with indices
// theta (contiguous, i.e. stride=1)
// phi
// surface (largest stride)
const int sub = i_theta + max_N_theta * (i_phi + max_N_phi*sn);
/* SphericalSurface:: */ sf_radius[sub] = local_r;
}
}
}
}
//******************************************************************************
} // namespace AHFinderDirect
|
