diff options
-rw-r--r-- | src/driver/driver.hh | 4 | ||||
-rw-r--r-- | src/driver/find_horizons.cc | 20 | ||||
-rw-r--r-- | src/driver/setup.cc | 74 | ||||
-rw-r--r-- | src/gr/gfns.hh | 31 | ||||
-rw-r--r-- | src/gr/gr.hh | 8 | ||||
-rw-r--r-- | src/gr/horizon_function.cc | 318 |
6 files changed, 372 insertions, 83 deletions
diff --git a/src/driver/driver.hh b/src/driver/driver.hh index 554ce33..c71df91 100644 --- a/src/driver/driver.hh +++ b/src/driver/driver.hh @@ -176,6 +176,10 @@ extern "C" //************************************** +// setup.cc +const CCTK_REAL* Cactus_gridfn_data_ptr(const cGH *GH, + const char gridfn_name[]); + // Newton.cc // returns true for success, false for failure to converge bool Newton_solve(patch_system& ps, diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc index 9fc2545..a228915 100644 --- a/src/driver/find_horizons.cc +++ b/src/driver/find_horizons.cc @@ -84,11 +84,27 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS const struct verbose_info& verbose_info = state.verbose_info; -state.IO_info.time_iteration = cctk_iteration; - if (state.timer_handle >= 0) then CCTK_TimerResetI(state.timer_handle); +// what are the semantics of the Cactus gxx variables? +if (CCTK_Equals(metric_type, "physical")) + then state.cgi.Cactus_conformal_metric = false; +else if (CCTK_Equals(metric_type, "static conformal")) + then state.cgi.Cactus_conformal_metric = (conformal_state > 0); +else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, +"AHFinderDirect_find_horizons(): unknown metric_type=\"%s\"!", + metric_type); /*NOTREACHED*/ + +// if Cactus is using a conformal metric, +// we need to re-fetch the data pointer for the conformal factor, +// since the data array may have moved in memory since the last +// time we were called due to storage being switched off/on +if (state.cgi.Cactus_conformal_metric) + then state.cgi.psi_data = Cactus_gridfn_data_ptr(cctkGH, + "StaticConformal::psi"); + +state.IO_info.time_iteration = cctk_iteration; for (int hn = 1 ; hn <= state.N_horizons ; ++hn) { struct AH_info& AH_info = * state.AH_info_ptrs[hn]; diff --git a/src/driver/setup.cc b/src/driver/setup.cc index cd9efd7..375e458 100644 --- a/src/driver/setup.cc +++ b/src/driver/setup.cc @@ -5,6 +5,7 @@ // <<<access to persistent data>>> // // AHFinderDirect_driver - top-level driver to setup persistent data structures +// Cactus_gridfn_data_ptr - get Cactus data pointer for CCTK_REAL gridfn /// /// decode_method - decode the method parameter /// decode_verbose_level - decode the verbose_level parameter @@ -82,10 +83,8 @@ extern struct state state; //****************************************************************************** // -// This function is called by the Cactus scheduler to set up the -// following data structures: -// state.gi -// state.cgi +// This function is called by the Cactus scheduler to set up all our +// persistent data structures. (These are stored in struct state .) // extern "C" void AHFinderDirect_setup(CCTK_ARGUMENTS) @@ -153,6 +152,7 @@ solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence; solver_info.final_H_update_if_exit_x_H_small = (final_H_update_if_exit_x_H_small != 0); + // // set up the Cactus grid info // @@ -169,23 +169,23 @@ cgi.coord_delta[2] = cctk_delta_space[2]; cgi.gridfn_dims[0] = cctk_lsh[0]; cgi.gridfn_dims[1] = cctk_lsh[1]; cgi.gridfn_dims[2] = cctk_lsh[2]; -// CCTK_VarDataPtr() returns a void * , but we need a const fp *; -// since static_cast<> won't change const-ness, we need a 2-stage cast -// for this: -#define DATA_PTR_CAST(void_ptr_) \ - static_cast<const fp*>( const_cast<const void *>(void_ptr_) ) -cgi.g_dd_11_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxx")); -cgi.g_dd_12_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxy")); -cgi.g_dd_13_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxz")); -cgi.g_dd_22_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gyy")); -cgi.g_dd_23_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gyz")); -cgi.g_dd_33_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gzz")); -cgi.K_dd_11_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxx")); -cgi.K_dd_12_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxy")); -cgi.K_dd_13_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxz")); -cgi.K_dd_22_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kyy")); -cgi.K_dd_23_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kyz")); -cgi.K_dd_33_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kzz")); + +// dummy values -- real ones filled in in AHFinderDirect_find_horizons() +cgi.Cactus_conformal_metric = false; +cgi.psi_data = NULL; + +cgi.g_dd_11_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::gxx"); +cgi.g_dd_12_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::gxy"); +cgi.g_dd_13_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::gxz"); +cgi.g_dd_22_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::gyy"); +cgi.g_dd_23_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::gyz"); +cgi.g_dd_33_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::gzz"); +cgi.K_dd_11_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::kxx"); +cgi.K_dd_12_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::kxy"); +cgi.K_dd_13_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::kxz"); +cgi.K_dd_22_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::kyy"); +cgi.K_dd_23_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::kyz"); +cgi.K_dd_33_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::kzz"); // @@ -301,6 +301,38 @@ state.AH_info_ptrs.push_back(NULL); //****************************************************************************** // +// This function gets the Cactus data pointer for a gridfn, and checks +// to make sure this is non-NULL. +// +const CCTK_REAL* Cactus_gridfn_data_ptr(const cGH *GH, const char gridfn_name[]) +{ +const int time_level = 0; + +// +// CCTK_VarDataPtr() returns a void * , but we need a const CCTK_REAL*; +// since static_cast<> won't change const-ness, we need a 2-stage cast +// for this: +// +const CCTK_REAL *const data_ptr + = static_cast<const fp*>( + const_cast<const void *>( + CCTK_VarDataPtr(GH, time_level, gridfn_name) + ) + ); + +if (data_ptr == NULL) + then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, +"***** Cactus_gridfn_data_ptr(): can't find gridfn \"%s\"!", + gridfn_name); /*NOTREACHED*/ + +return data_ptr; +} + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// // This function decodes the method parameter (string) into an // internal enum for future use. // diff --git a/src/gr/gfns.hh b/src/gr/gfns.hh index 5404f46..f984f4f 100644 --- a/src/gr/gfns.hh +++ b/src/gr/gfns.hh @@ -15,22 +15,21 @@ enum { // nominal gridfns enum { nominal_min_gfn = 1, - // n.b. no macros for next 3 gridfns in "cg.hh" - gfn__global_x = nominal_min_gfn, - gfn__global_y, - gfn__global_z, + + // + // most of these gridfns have access macros in "cg.hh"; + // the ones that don't are marked explicitly + // + gfn__global_x = nominal_min_gfn, // no access macro + gfn__global_y, // no access macro + gfn__global_z, // no access macro + gfn__g_dd_11, gfn__g_dd_12, gfn__g_dd_13, gfn__g_dd_22, gfn__g_dd_23, gfn__g_dd_33, - gfn__K_dd_11, - gfn__K_dd_12, - gfn__K_dd_13, - gfn__K_dd_22, - gfn__K_dd_23, - gfn__K_dd_33, gfn__partial_d_g_dd_111, gfn__partial_d_g_dd_112, gfn__partial_d_g_dd_113, @@ -49,6 +48,18 @@ enum { gfn__partial_d_g_dd_322, gfn__partial_d_g_dd_323, gfn__partial_d_g_dd_333, + gfn__K_dd_11, + gfn__K_dd_12, + gfn__K_dd_13, + gfn__K_dd_22, + gfn__K_dd_23, + gfn__K_dd_33, + + gfn__psi, // no access macro + gfn__partial_d_psi_1, // no access macro + gfn__partial_d_psi_2, // no access macro + gfn__partial_d_psi_3, // no access macro + gfn__H, gfn__partial_H_wrt_partial_d_h_1, gfn__partial_H_wrt_partial_d_h_2, diff --git a/src/gr/gr.hh b/src/gr/gr.hh index 53f88c4..beb1866 100644 --- a/src/gr/gr.hh +++ b/src/gr/gr.hh @@ -41,6 +41,12 @@ struct cactus_grid_info // 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 gridfn data for grid posn (0,0,0) const fp* g_dd_11_data; const fp* g_dd_12_data; @@ -54,6 +60,8 @@ struct cactus_grid_info 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 }; // diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc index 3804be2..9b1e0a6 100644 --- a/src/gr/horizon_function.cc +++ b/src/gr/horizon_function.cc @@ -5,6 +5,7 @@ // horizon_function - top-level driver /// setup_xyz_posns - setup global xyz posns of grid points /// interpolate_geometry - interpolate g_ij and K_ij from Cactus 3-D grid +/// convert_conformal_to_physical - convert conformal gij to physical /// compute_H - compute H(h) given earlier setup /// @@ -53,6 +54,7 @@ bool interpolate_geometry(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_info& gi, bool print_msg_flag); +void convert_conformal_to_physical(patch_system& ps, bool print_msg_flag); void Schwarzschild_EF_geometry(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_info& gi, @@ -135,6 +137,8 @@ switch (gi.geometry_method) case geometry__interpolate_from_Cactus_grid: if (! interpolate_geometry(ps, cgi, gi, print_msg_flag)) then return false; // *** ERROR RETURN *** + if (cgi.Cactus_conformal_metric) + then convert_conformal_to_physical(ps, print_msg_flag); break; case geometry__Schwarzschild_EF: Schwarzschild_EF_geometry(ps, gi, print_msg_flag); @@ -200,34 +204,45 @@ if (print_msg_flag) //****************************************************************************** // -// This function interpolates the Cactus gridfns $g_{ij}$ and $K_{ij}$ -// to determine $g_{ij}$, $K_{ij}$, and the spatial derivatives -// $\partial_k g_{ij}$, on the trial horizon surface position given by h. -// -// Inputs (angular gridfns, on ghosted grid): -// ... defined on ghosted grid -// ... only values on nominal grid are actually used as input -// h # shape of trial surface +// This function interpolates the Cactus gridfns +// gxx...gzz +// kxx...kzz +// psi # optional +// to determine the nominal-grid angular gridfns +// g_dd_ij +// partial_d_g_dd_kij +// K_dd_ij +// psi # optional +// partial_d_psi_k # optional +// at the nominal-grid trial horizon surface positions given by the +// global_(x,y,z) angular gridfns. The psi interpolation is done if and +// only if the cgi.Cactus_conformal_metric flag is set. Note that this +// function ignores the physical-vs-conformal semantics of the gridfns; +// it just interpolates and takes derivatives of the stored gridfn values. // // Inputs (angular gridfns, all on the nominal grid): // global_[xyz] # xyz positions of grid points // // Inputs (Cactus 3-D gridfns): // gxx,gxy,gxz,gyy,gyz,gzz # 3-metric $g_{ij}$ +// # (may be either physical or conformal) // kxx,kxy,kxz,kyy,kyz,kzz # extrinsic curvature $K_{ij}$ +// psi # optional conformal factor $\psi$ // // Outputs (angular gridfns, all on the nominal grid): -// g_dd_{11,12,13,22,23,33} # $g_{ij}$ -// K_dd_{11,12,13,22,23,33} # $K_{ij}$ -// partial_d_g_dd_{1,2,3}{11,12,13,22,23,33} # $\partial_k g_{ij}$ +// g_dd_{11,12,13,22,23,33} # $\stored{g}_{ij}$ +// K_dd_{11,12,13,22,23,33} # $K_{ij}$ +// partial_d_g_dd_[123]{11,12,13,22,23,33} # $\partial_k \stored{g}_{ij}$ +// psi # (optional) $\psi$ +// partial_d_psi_[123] # (optional) $\partial_k \psi$ // -// On the first call this function also modifies the interpolator -// parameter table. +// This function may also modify the interpolator parameter table. // // Results: // This function returns true for a successful computation, or false // if the computation failed because the geometry interpolation would // need data outside the Cactus grid, or data from an excised region. +// // FIXME: excision isn't implemented yet :( // namespace { @@ -236,12 +251,32 @@ bool interpolate_geometry(patch_system& ps, const struct geometry_info& gi, bool print_msg_flag) { +// +// Implementation Notes: +// +// To handle the optional interpolation of psi, we set up all the data +// type and pointer arrays to include psi, but with the psi entries at +// the end. We then choose the array sizes passed to the interpolator +// to either include or exclude the psi entries as appropriate. +// +// We remember whether or not psi was interpolated on the previous call, +// and only modify the interpolator parameter table if this changes (and +// on our first call). +// + +const bool psi_flag = cgi.Cactus_conformal_metric; + if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, - " interpolating g_ij and Kij from Cactus 3-D grid"); + " interpolating %s from Cactus 3-D grid", + psi_flag ? "{g_ij, K_ij, psi}" : "{g_ij, K_ij}"); int status; + +// +// ***** interpolation points ***** +// const int N_interp_points = ps.N_grid_points(); const int interp_coords_type_code = CCTK_VARIABLE_REAL; const void* const interp_coords[N_GRID_DIMS] @@ -251,8 +286,12 @@ const void* const interp_coords[N_GRID_DIMS] static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_z)), }; -const int N_INPUT_ARRAYS = 12; -const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS] + +// +// ***** input arrays ***** +// + +const CCTK_INT input_array_type_codes[] = { // g_ij CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, @@ -262,8 +301,10 @@ const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS] CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + // psi + CCTK_VARIABLE_REAL, }; -const void* const input_arrays[N_INPUT_ARRAYS] +const void* const input_arrays[] = { static_cast<const void*>(cgi.g_dd_11_data), static_cast<const void*>(cgi.g_dd_12_data), @@ -277,11 +318,21 @@ const void* const input_arrays[N_INPUT_ARRAYS] static_cast<const void*>(cgi.K_dd_22_data), static_cast<const void*>(cgi.K_dd_23_data), static_cast<const void*>(cgi.K_dd_33_data), + static_cast<const void*>(cgi.psi_data), }; -const int N_OUTPUT_ARRAYS = 30; +const int N_input_arrays_for_psi = 1; +const int N_input_arrays_dim = sizeof(input_arrays) / sizeof(input_arrays[0]); +const int N_input_arrays_use + = psi_flag ? N_input_arrays_dim + : N_input_arrays_dim - N_input_arrays_for_psi; -const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS] + +// +// ***** output arrays ***** +// + +const CCTK_INT output_array_type_codes[] = { // g_ij partial_x g_ij partial_y g_ij partial_z g_ij CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, @@ -294,33 +345,36 @@ const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS] CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + // psi partial_x psi partial_y psi partial_z psi + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, }; -const CCTK_INT operand_indices[N_OUTPUT_ARRAYS] + +const CCTK_INT operand_indices[] = { - 0, 0, 0, 0, // g_dd_11 - 1, 1, 1, 1, // g_dd_12 - 2, 2, 2, 2, // g_dd_13 - 3, 3, 3, 3, // g_dd_22 - 4, 4, 4, 4, // g_dd_23 - 5, 5, 5, 5, // g_dd_33 - 6, 7, 8, // K_dd_{11,12,13} - 9, 10, // K_dd_{22,23} - 11, // K_dd_33 + 0, 0, 0, 0, // g_dd_11, partial_[xyz] g_dd_11 + 1, 1, 1, 1, // g_dd_12, partial_[xyz] g_dd_12 + 2, 2, 2, 2, // g_dd_13, partial_[xyz] g_dd_13 + 3, 3, 3, 3, // g_dd_22, partial_[xyz] g_dd_22 + 4, 4, 4, 4, // g_dd_23, partial_[xyz] g_dd_23 + 5, 5, 5, 5, // g_dd_33, partial_[xyz] g_dd_33 + 6, 7, 8, 9, 10, 11, // K_dd_{11,12,13,22,23,33} + 12, 12, 12, 12, // psi, partial_[xyz] psi }; #define DERIV(x) x -const CCTK_INT operation_codes[N_OUTPUT_ARRAYS] - = { - DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_11 - DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_12 - DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_13 - DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_22 - DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_23 - DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_33 - DERIV(0), DERIV(0), DERIV(0), // K_dd_{11,12,13} - DERIV(0), DERIV(0), // K_dd_{22,23} - DERIV(0), // K_dd_{33} - }; -void* const output_arrays[N_OUTPUT_ARRAYS] +const CCTK_INT operation_codes[] + = { + DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_11, partial_[xyz] g_dd_11 + DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_12, partial_[xyz] g_dd_12 + DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_13, partial_[xyz] g_dd_13 + DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_22, partial_[xyz] g_dd_22 + DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_23, partial_[xyz] g_dd_23 + DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_33, partial_[xyz] g_dd_33 + DERIV(0), DERIV(0), DERIV(0), DERIV(0), DERIV(0), DERIV(0), + // K_dd_{11,12,13,22,23,33} + DERIV(0), DERIV(1), DERIV(2), DERIV(3), // psi, partial_[xyz] psi + }; + +void* const output_arrays[] = { static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_11)), static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_111)), @@ -352,20 +406,48 @@ void* const output_arrays[N_OUTPUT_ARRAYS] static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_22)), static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_23)), static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_33)), + static_cast<void*>(ps.gridfn_data(gfns::gfn__psi)), + static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_psi_1)), + static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_psi_2)), + static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_psi_3)), }; -bool first_time = true; -if (first_time) - then { - first_time = false; +const int N_output_arrays_for_psi = 4; +const int N_output_arrays_dim = sizeof(output_arrays) + / sizeof(output_arrays[0]); +const int N_output_arrays_use + = psi_flag ? N_output_arrays_dim + : N_output_arrays_dim - N_output_arrays_for_psi; + +// +// ***** parameter table ***** +// + +// this flag is true if and only if the parameter table already has the +// operand_indices +// operand_codes +// entries for psi_flag == par_table_psi_flag . +static bool par_table_setup = false; + +// if par_table_setup, +// this flag is the value of psi_flag for the parameter table entries +// otherwise this flag is ignored +static bool par_table_psi_flag = false; + +if (par_table_setup && (psi_flag == par_table_psi_flag)) + then { + // parameter table is already set to just what we need + // ==> no-op here + } + else { // store derivative info in interpolator parameter table if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " setting up interpolator derivative info"); status = Util_TableSetIntArray(gi.param_table_handle, - N_OUTPUT_ARRAYS, operand_indices, + N_output_arrays_use, operand_indices, "operand_indices"); if (status < 0) then error_exit(ERROR_EXIT, @@ -376,7 +458,7 @@ if (first_time) status); /*NOTREACHED*/ status = Util_TableSetIntArray(gi.param_table_handle, - N_OUTPUT_ARRAYS, operation_codes, + N_output_arrays_use, operation_codes, "operation_codes"); if (status < 0) then error_exit(ERROR_EXIT, @@ -385,8 +467,15 @@ if (first_time) " Util_TableSetIntArray() status=%d\n" , status); /*NOTREACHED*/ + + par_table_setup = true; + par_table_psi_flag = psi_flag; } + +// +// ***** the actual interpolation ***** +// if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " calling interpolator (%d points)", @@ -397,13 +486,18 @@ status = CCTK_InterpLocalUniform(N_GRID_DIMS, N_interp_points, interp_coords_type_code, interp_coords, - N_INPUT_ARRAYS, + N_input_arrays_use, cgi.gridfn_dims, input_array_type_codes, input_arrays, - N_OUTPUT_ARRAYS, + N_output_arrays_use, output_array_type_codes, output_arrays); + + +// +// ***** handle any interpolation errors ***** +// if (status == CCTK_ERROR_INTERP_POINT_X_RANGE) then { // look in interpolator output table entries @@ -461,6 +555,7 @@ if (status < 0) , status); /*NOTREACHED*/ + return true; // *** NORMAL RETURN *** } } @@ -468,6 +563,129 @@ return true; // *** NORMAL RETURN *** //****************************************************************************** // +// This function converts the g_dd_ij and partial_d_g_dd_kij gridfns +// from the Cactus conformal semantics to the physical semantics. As +// documented in the CactusEinstein/ConformalState thorn guide, the +// Cactus conformal semantics are +// physical_gij = psi^4 * stored_gij +// From this it's trivial to derive the transformation for the derivatives, +// partial_k physical_gij = 4 psi^3 (partial_k psi) stored_gij +// + psi^4 partial_k stored_gij +// +// Inputs (angular gridfns, all on the nominal grid): +// psi # $\psi$ +// partial_d_psi_{1,2,3} # $\partial_k \psi$ +// g_dd_{11,12,13,22,23,33} # $\stored{g}_{ij}$ +// partial_d_g_dd_{1,2,3}{11,12,13,22,23,33} # $\partial_k \stored{g}_{ij}$ +// +// +// Outputs (angular gridfns, all on the nominal grid): +// g_dd_{11,12,13,22,23,33} # $g_{ij}$ +// partial_d_g_dd_{1,2,3}{11,12,13,22,23,33} # $\partial_k g_{ij}$ +// +namespace { +void convert_conformal_to_physical(patch_system& ps, bool print_msg_flag) +{ +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + " converting Cactus conformal gij --> physical gij"); + + for (int pn = 0 ; pn < ps.N_patches() ; ++pn) + { + patch& p = ps.ith_patch(pn); + + for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho) + { + for (int isigma = p.min_isigma() ; isigma <= p.max_isigma() ; ++isigma) + { + const fp psi = p.gridfn(gfns::gfn__psi, irho,isigma); + const fp psi3 = jtutil::pow3(psi); + const fp psi4 = jtutil::pow4(psi); + + const fp partial_d_psi_1 + = p.gridfn(gfns::gfn__partial_d_psi_1, irho,isigma); + const fp partial_d_psi_2 + = p.gridfn(gfns::gfn__partial_d_psi_2, irho,isigma); + const fp partial_d_psi_3 + = p.gridfn(gfns::gfn__partial_d_psi_3, irho,isigma); + + const fp stored_g_dd_11 = p.gridfn(gfns::gfn__g_dd_11, irho, isigma); + const fp stored_g_dd_12 = p.gridfn(gfns::gfn__g_dd_12, irho, isigma); + const fp stored_g_dd_13 = p.gridfn(gfns::gfn__g_dd_13, irho, isigma); + const fp stored_g_dd_22 = p.gridfn(gfns::gfn__g_dd_22, irho, isigma); + const fp stored_g_dd_23 = p.gridfn(gfns::gfn__g_dd_23, irho, isigma); + const fp stored_g_dd_33 = p.gridfn(gfns::gfn__g_dd_33, irho, isigma); + + p.gridfn(gfns::gfn__g_dd_11, irho, isigma) *= psi4; + p.gridfn(gfns::gfn__g_dd_12, irho, isigma) *= psi4; + p.gridfn(gfns::gfn__g_dd_13, irho, isigma) *= psi4; + p.gridfn(gfns::gfn__g_dd_22, irho, isigma) *= psi4; + p.gridfn(gfns::gfn__g_dd_23, irho, isigma) *= psi4; + p.gridfn(gfns::gfn__g_dd_33, irho, isigma) *= psi4; + + p.gridfn(gfns::gfn__partial_d_g_dd_111, irho,isigma) + = 4.0*psi3*partial_d_psi_1*stored_g_dd_11 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_111, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_112, irho,isigma) + = 4.0*psi3*partial_d_psi_1*stored_g_dd_12 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_112, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_113, irho,isigma) + = 4.0*psi3*partial_d_psi_1*stored_g_dd_13 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_113, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_122, irho,isigma) + = 4.0*psi3*partial_d_psi_1*stored_g_dd_22 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_122, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_123, irho,isigma) + = 4.0*psi3*partial_d_psi_1*stored_g_dd_23 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_123, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_133, irho,isigma) + = 4.0*psi3*partial_d_psi_1*stored_g_dd_33 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_133, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_211, irho,isigma) + = 4.0*psi3*partial_d_psi_2*stored_g_dd_11 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_211, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_212, irho,isigma) + = 4.0*psi3*partial_d_psi_2*stored_g_dd_12 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_212, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_213, irho,isigma) + = 4.0*psi3*partial_d_psi_2*stored_g_dd_13 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_213, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_222, irho,isigma) + = 4.0*psi3*partial_d_psi_2*stored_g_dd_22 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_222, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_223, irho,isigma) + = 4.0*psi3*partial_d_psi_2*stored_g_dd_23 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_223, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_233, irho,isigma) + = 4.0*psi3*partial_d_psi_2*stored_g_dd_33 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_233, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_311, irho,isigma) + = 4.0*psi3*partial_d_psi_3*stored_g_dd_11 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_311, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_312, irho,isigma) + = 4.0*psi3*partial_d_psi_3*stored_g_dd_12 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_312, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_313, irho,isigma) + = 4.0*psi3*partial_d_psi_3*stored_g_dd_13 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_313, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_322, irho,isigma) + = 4.0*psi3*partial_d_psi_3*stored_g_dd_22 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_322, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_323, irho,isigma) + = 4.0*psi3*partial_d_psi_3*stored_g_dd_23 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_323, irho,isigma); + p.gridfn(gfns::gfn__partial_d_g_dd_333, irho,isigma) + = 4.0*psi3*partial_d_psi_3*stored_g_dd_33 + + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_333, irho,isigma); + } + } + } +} + } + +//****************************************************************************** + +// // This function computes H(h), and optionally its Jacobian coefficients, // (from which the Jacobian matrix may be computed later). This function // uses a mixture of algebraic operations and (rho,sigma) finite differencing. |