diff options
Diffstat (limited to 'src/driver/spherical_surface.cc')
| -rw-r--r-- | src/driver/spherical_surface.cc | 264 |
1 files changed, 0 insertions, 264 deletions
diff --git a/src/driver/spherical_surface.cc b/src/driver/spherical_surface.cc deleted file mode 100644 index 96d3bc8..0000000 --- a/src/driver/spherical_surface.cc +++ /dev/null @@ -1,264 +0,0 @@ -// 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 |