diff options
Diffstat (limited to 'src/driver/spherical_surface.cc')
| -rw-r--r-- | src/driver/spherical_surface.cc | 257 |
1 files changed, 257 insertions, 0 deletions
diff --git a/src/driver/spherical_surface.cc b/src/driver/spherical_surface.cc new file mode 100644 index 0000000..20617bd --- /dev/null +++ b/src/driver/spherical_surface.cc @@ -0,0 +1,257 @@ +// 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" +using jtutil::error_exit; + +#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 + { + +//****************************************************************************** + +// +// ***** 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.areal_radius; +/* SphericalSurface:: */ sf_min_radius [sn] = BH_diagnostics.min_radius; +/* SphericalSurface:: */ sf_max_radius [sn] = BH_diagnostics.max_radius; + +/* 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. +// +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 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] = r; + } + } +} + } + +//****************************************************************************** + + } // namespace AHFinderDirect |