diff options
Diffstat (limited to 'src/patch')
-rw-r--r-- | src/patch/coords.cc | 80 | ||||
-rw-r--r-- | src/patch/coords.hh | 5 | ||||
-rw-r--r-- | src/patch/fd_grid.cc | 2 | ||||
-rw-r--r-- | src/patch/fd_grid.hh | 2 | ||||
-rw-r--r-- | src/patch/ghost_zone.cc | 2 | ||||
-rw-r--r-- | src/patch/ghost_zone.hh | 3 | ||||
-rw-r--r-- | src/patch/make.code.defn | 14 | ||||
-rw-r--r-- | src/patch/patch.cc | 57 | ||||
-rw-r--r-- | src/patch/patch.hh | 9 | ||||
-rw-r--r-- | src/patch/patch_info.cc | 2 | ||||
-rw-r--r-- | src/patch/patch_interp.cc | 8 | ||||
-rw-r--r-- | src/patch/patch_system.cc | 702 | ||||
-rw-r--r-- | src/patch/patch_system.hh | 119 | ||||
-rw-r--r-- | src/patch/test_coords.cc | 6 | ||||
-rw-r--r-- | src/patch/test_coords2.cc | 6 | ||||
-rw-r--r-- | src/patch/test_patch_system.cc | 2 |
16 files changed, 916 insertions, 103 deletions
diff --git a/src/patch/coords.cc b/src/patch/coords.cc index 828553a..1b229c2 100644 --- a/src/patch/coords.cc +++ b/src/patch/coords.cc @@ -33,22 +33,17 @@ #include "config.h" #include "stdc.h" #include "../jtutil/util.hh" +using jtutil::error_exit; +using jtutil::arctan_xy; +using jtutil::signum; +using jtutil::pow2; +using jtutil::hypot3; #include "coords.hh" -// FIXME: sometimes we get assertion failures in these checks, -// but I've never been able to track down just what's wrong -// ==> for now, just disable them :( :( -#undef XYZ_OF_RMNP_SIGN_CHECKS - // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; -using jtutil::arctan_xy; -using jtutil::signum; -using jtutil::pow2; -using jtutil::hypot3; //****************************************************************************** //****************************************************************************** @@ -134,15 +129,13 @@ z = sign_z * r * temp; x = x_over_z * z; y = y_over_z * z; -#ifdef XYZ_OF_RMNP_SIGN_CHECKS -if (jtutil::fuzzy<fp>::NE(r, 0.0)) - then { - if (jtutil::fuzzy<fp>::NE(x, 0.0)) - then assert( signum(x) == sign_x ); - if (jtutil::fuzzy<fp>::NE(y, 0.0)) - then assert( signum(y) == sign_y ); - } -#endif +// if (jtutil::fuzzy<fp>::NE(r, 0.0)) +// then { +// if (jtutil::fuzzy<fp>::NE(x, 0.0)) +// then assert( signum(x) == sign_x ); +// if (jtutil::fuzzy<fp>::NE(y, 0.0)) +// then assert( signum(y) == sign_y ); +// } } } @@ -179,15 +172,13 @@ y = sign_y * r * temp; z = z_over_y * y; x = x_over_y * y; -#ifdef XYZ_OF_RMNP_SIGN_CHECKS -if (jtutil::fuzzy<fp>::NE(r, 0.0)) - then { - if (jtutil::fuzzy<fp>::NE(z, 0.0)) - then assert( signum(z) == sign_z ); - if (jtutil::fuzzy<fp>::NE(x, 0.0)) - then assert( signum(x) == sign_x ); - } -#endif +// if (jtutil::fuzzy<fp>::NE(r, 0.0)) +// then { +// if (jtutil::fuzzy<fp>::NE(z, 0.0)) +// then assert( signum(z) == sign_z ); +// if (jtutil::fuzzy<fp>::NE(x, 0.0)) +// then assert( signum(x) == sign_x ); +// } } } @@ -223,15 +214,13 @@ x = sign_x * r * temp; z = z_over_x * x; y = y_over_x * x; -#ifdef XYZ_OF_RMNP_SIGN_CHECKS -if (jtutil::fuzzy<fp>::NE(r, 0.0)) - then { - if (jtutil::fuzzy<fp>::NE(z, 0.0)) - then assert( signum(z) == sign_z ); - if (jtutil::fuzzy<fp>::NE(y, 0.0)) - then assert( signum(y) == sign_y ); - } -#endif +// if (jtutil::fuzzy<fp>::NE(r, 0.0)) +// then { +// if (jtutil::fuzzy<fp>::NE(z, 0.0)) +// then assert( signum(z) == sign_z ); +// if (jtutil::fuzzy<fp>::NE(y, 0.0)) +// then assert( signum(y) == sign_y ); +// } } } @@ -321,7 +310,8 @@ const fp tan2_nu = pow2(tan_nu); fp x, y, z; xyz_of_r_mu_nu(r,mu,nu, x,y,z); -assert( jtutil::fuzzy<fp>::NE(r, 0.0) ); +// Comment this out, accept a nan instead +// assert( jtutil::fuzzy<fp>::NE(r, 0.0) ); const fp rinv = 1.0/r; partial_x_wrt_r = x*rinv; partial_y_wrt_r = y*rinv; @@ -363,7 +353,8 @@ const fp tan2_phi_bar = pow2(tan_phi_bar); fp x, y, z; xyz_of_r_mu_phi(r,mu,phi, x,y,z); -assert( jtutil::fuzzy<fp>::NE(r, 0.0) ); +// Comment this out, accept a nan instead +// assert( jtutil::fuzzy<fp>::NE(r, 0.0) ); const fp rinv = 1.0/r; partial_x_wrt_r = x*rinv; partial_y_wrt_r = y*rinv; @@ -404,7 +395,8 @@ const fp tan2_phi = pow2(tan_phi); fp x, y, z; xyz_of_r_nu_phi(r,nu,phi, x,y,z); -assert( jtutil::fuzzy<fp>::NE(r, 0.0) ); +// Comment this out, accept a nan instead +// assert( jtutil::fuzzy<fp>::NE(r, 0.0) ); const fp rinv = 1.0/r; partial_x_wrt_r = x*rinv; partial_y_wrt_r = y*rinv; @@ -659,13 +651,13 @@ else if (S == coords_set_nu) then return "nu"; else if (S == coords_set_phi) then return "phi"; -else if (S == coords_set_mu|coords_set_nu) +else if (S == (coords_set_mu|coords_set_nu)) then return "{mu,nu}"; -else if (S == coords_set_mu|coords_set_phi) +else if (S == (coords_set_mu|coords_set_phi)) then return "{mu,phi}"; -else if (S == coords_set_nu|coords_set_phi) +else if (S == (coords_set_nu|coords_set_phi)) then return "{nu,phi}"; -else if (S == coords_set_mu|coords_set_nu|coords_set_phi) +else if (S == (coords_set_mu|coords_set_nu|coords_set_phi)) then return "{mu,nu,phi}"; else error_exit(PANIC_EXIT, "***** local_coords::mu_nu_phi::name_of_coords_set:\n" diff --git a/src/patch/coords.hh b/src/patch/coords.hh index e6ca9bd..f362474 100644 --- a/src/patch/coords.hh +++ b/src/patch/coords.hh @@ -278,6 +278,11 @@ public: fp origin_y() const { return origin_y_; } fp origin_z() const { return origin_z_; } + // set global (x,y,z) coordinates of local origin point + void origin_x(const fp x) { origin_x_=x; } + void origin_y(const fp y) { origin_y_=y; } + void origin_z(const fp z) { origin_z_=z; } + // radius of given (x,y,z) with respect to local origin point #ifdef NOT_USED fp radius_of_local_xyz(fp local_x, fp local_y, fp local_z) const diff --git a/src/patch/fd_grid.cc b/src/patch/fd_grid.cc index 80d1e1b..48e490e 100644 --- a/src/patch/fd_grid.cc +++ b/src/patch/fd_grid.cc @@ -20,6 +20,7 @@ #include "../jtutil/util.hh" #include "../jtutil/array.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; #include "coords.hh" #include "grid.hh" @@ -28,7 +29,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; //***************************************************************************** diff --git a/src/patch/fd_grid.hh b/src/patch/fd_grid.hh index c750a75..9ae24c4 100644 --- a/src/patch/fd_grid.hh +++ b/src/patch/fd_grid.hh @@ -309,7 +309,7 @@ namespace AHFinderDirect #define FD_GRID__ORDER4__DXX__KPM2 ( 1.0/12.0) #define FD_GRID__ORDER4__DXX__KPM1 (16.0/12.0) -#define FD_GRID__ORDER4__DXX__K0 (30.0/12.0) +#define FD_GRID__ORDER4__DXX__K0 (30.0/12.0) #define FD_GRID__ORDER4__DXX(inv_delta_x_, data_, \ irho_plus_m_, isigma_plus_m_) \ const fp data_p2 = data_(ghosted_gfn, \ diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc index 32c073f..985c1ca 100644 --- a/src/patch/ghost_zone.cc +++ b/src/patch/ghost_zone.cc @@ -34,6 +34,7 @@ #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; #include "coords.hh" #include "grid.hh" @@ -46,7 +47,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; //****************************************************************************** //****************************************************************************** diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh index a35aa6f..22916ad 100644 --- a/src/patch/ghost_zone.hh +++ b/src/patch/ghost_zone.hh @@ -793,7 +793,8 @@ private: // ------------------- // partial gridfn at y // - mutable int Jacobian_min_y_ipar_m_, Jacobian_max_y_ipar_m_; + mutable int Jacobian_min_y_ipar_m_; + mutable int Jacobian_max_y_ipar_m_; // other patch's y ipar posn for a Jacobian row // ... subscripts are (oiperp, ipar) diff --git a/src/patch/make.code.defn b/src/patch/make.code.defn index 15194f2..a24223f 100644 --- a/src/patch/make.code.defn +++ b/src/patch/make.code.defn @@ -15,17 +15,3 @@ SRCS = coords.cc \ # Subdirectories containing source files SUBDIRS = - -# disable automatic template instantiation on DEC Alphas cxx compiler -ifeq ($(shell uname), OSF1) - ifeq ($(CXX), cxx) - CXXFLAGS += -nopt - endif -endif - -# disable automagic template instantiation on SGI Irix CC compiler -ifneq (,$(findstring IRIX,$(shell uname))) - ifeq ($(notdir $(CXX)), CC) - CXXFLAGS += -no_auto_include - endif -endif diff --git a/src/patch/patch.cc b/src/patch/patch.cc index f4305f6..6b355fd 100644 --- a/src/patch/patch.cc +++ b/src/patch/patch.cc @@ -43,6 +43,7 @@ #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; #include "coords.hh" #include "grid.hh" @@ -55,7 +56,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; //****************************************************************************** //****************************************************************************** @@ -300,7 +300,7 @@ else error_exit(PANIC_EXIT, " which is a multiple of 90 degrees!\n" " %s patch: [min,max]_dsigma()=[%g,%g]\n" , - name(), min_dsigma(), max_dsigma()); + name(), double(min_dsigma()), double(max_dsigma())); const fp sigma = sigma_of_dsigma(dsigma); const int isigma = isigma_of_sigma(sigma); @@ -371,7 +371,7 @@ else error_exit(PANIC_EXIT, " which is a multiple of 90 degrees!\n" " %s patch: [min,max]_drho()=[%g,%g]\n" , - name(), min_drho(), max_drho()); + name(), double(min_drho()), double(max_drho())); const fp rho = rho_of_drho(drho); const int irho = irho_of_rho(rho); @@ -611,6 +611,57 @@ fp sum = 0.0; return delta_rho() * delta_sigma() * sum; } +fp patch::integrate_gridpoint(int unknown_src_gfn, + int ghosted_radius_gfn, + int g_xx_gfn, int g_xy_gfn, int g_xz_gfn, + int g_yy_gfn, int g_yz_gfn, + int g_zz_gfn, + enum integration_method method, + int irho, int isigma) + const +{ +const bool src_is_ghosted = is_valid_ghosted_gfn(unknown_src_gfn); + +const fp fn = unknown_gridfn(src_is_ghosted, + unknown_src_gfn, irho,isigma); + +const fp rho = rho_of_irho (irho); +const fp sigma = sigma_of_isigma(isigma); +const fp r = ghosted_gridfn(ghosted_radius_gfn, irho,isigma); +const fp partial_surface_r_wrt_rho + = partial_rho (ghosted_radius_gfn, irho,isigma); +const fp partial_surface_r_wrt_sigma + = partial_sigma(ghosted_radius_gfn, irho,isigma); + +const fp g_xx = gridfn(g_xx_gfn, irho,isigma); +const fp g_xy = gridfn(g_xy_gfn, irho,isigma); +const fp g_xz = gridfn(g_xz_gfn, irho,isigma); +const fp g_yy = gridfn(g_yy_gfn, irho,isigma); +const fp g_yz = gridfn(g_yz_gfn, irho,isigma); +const fp g_zz = gridfn(g_zz_gfn, irho,isigma); + +fp g_rho_rho, g_rho_sigma, g_sigma_sigma; +const fp Jac = rho_sigma_metric(r, rho, sigma, + partial_surface_r_wrt_rho, + partial_surface_r_wrt_sigma, + g_xx, g_xy, g_xz, + g_yy, g_yz, + g_zz, + g_rho_rho, g_rho_sigma, + g_sigma_sigma); + +const fp coeff_rho = integration_coeff(method, + max_irho()-min_irho(), + irho -min_irho()); +const fp coeff_sigma = integration_coeff(method, + max_isigma()-min_isigma(), + isigma -min_isigma()); + +const fp val = coeff_rho*coeff_sigma * fn * sqrt(jtutil::abs(Jac)); + +return delta_rho() * delta_sigma() * val; +} + //****************************************************************************** // diff --git a/src/patch/patch.hh b/src/patch/patch.hh index d18c53c..6482154 100644 --- a/src/patch/patch.hh +++ b/src/patch/patch.hh @@ -416,6 +416,15 @@ public: enum integration_method method) const; + fp integrate_gridpoint(int unknown_src_gfn, + int ghosted_radius_gfn, + int g_xx_gfn, int g_xy_gfn, int g_xz_gfn, + int g_yy_gfn, int g_yz_gfn, + int g_zz_gfn, + enum integration_method method, + int irho, int isigma) + const; + // compute integration coefficient $c_i$ where // $\int_{x_0}^{x_N} f(x) \, dx // \approx \Delta x \, \sum_{i=0}^N c_i f(x_i)$ diff --git a/src/patch/patch_info.cc b/src/patch/patch_info.cc index 002b63b..c293cf6 100644 --- a/src/patch/patch_info.cc +++ b/src/patch/patch_info.cc @@ -19,6 +19,7 @@ #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; #include "coords.hh" #include "grid.hh" @@ -27,7 +28,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; //****************************************************************************** //****************************************************************************** diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc index 580ad3f..a7ac765 100644 --- a/src/patch/patch_interp.cc +++ b/src/patch/patch_interp.cc @@ -26,6 +26,7 @@ #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; #include "coords.hh" #include "grid.hh" @@ -38,7 +39,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; //***************************************************************************** //***************************************************************************** @@ -306,9 +306,9 @@ if ( MSS_is_fn_of_interp_coords || MSS_is_fn_of_input_array_values " MSS_is_fn_of_input_array_values=(int)%d (we only grok 0)\n" " Jacobian_is_fn_of_input_array_values=(int)%d (we only grok 0)\n" , - MSS_is_fn_of_interp_coords, - MSS_is_fn_of_input_array_values, - Jacobian_is_fn_of_input_array_values); /*NOTREACHED*/ + int(MSS_is_fn_of_interp_coords), + int(MSS_is_fn_of_input_array_values), + int(Jacobian_is_fn_of_input_array_values)); /*NOTREACHED*/ } //****************************************************************************** diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc index 98974e6..bbd29f7 100644 --- a/src/patch/patch_system.cc +++ b/src/patch/patch_system.cc @@ -30,6 +30,8 @@ // // patch_system::set_gridfn_to_constant // patch_system::gridfn_copy +// patch_system::ghosted_gridfn_copy +// patch_system::add_to_gridfn // patch_system::add_to_ghosted_gridfn // patch_system::gridfn_norms // patch_system::ghosted_gridfn_norms @@ -39,9 +41,11 @@ // // patch_system::patch_containing_local_xyz // patch_system::radius_in_local_xyz_direction +// patch_system::radii_in_local_xyz_directions // // patch_system::print_unknown_gridfn // patch_system::read_unknown_gridfn +// patch_system::output_unknown_gridfn // // patch_system::synchronize // patch_system::compute_synchronize_Jacobian @@ -51,20 +55,37 @@ /// patch_system::ghost_zone_Jacobian // +#include <assert.h> #include <stdio.h> #include <math.h> #include <string.h> #include <assert.h> #include <limits.h> +#include <sys/stat.h> +#include <sys/types.h> + +#include <vector> +#include <sstream> +#include <string> #include "cctk.h" +#include "CactusBase/IOUtil/src/ioGH.h" #include "config.h" #include "stdc.h" +#include "../gr/gfns.hh" #include "../jtutil/util.hh" #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; + +#ifdef CCTK_HDF5 +// Some macros to fix compatibility issues as long as 1.8.0 is in beta +// phase +# define H5_USE_16_API 1 +# include <hdf5.h> +#endif #include "coords.hh" #include "grid.hh" @@ -80,7 +101,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; //****************************************************************************** //****************************************************************************** @@ -325,12 +345,13 @@ delete[] ghosted_gridfn_storage_; delete[] gridfn_storage_; delete[] ghosted_starting_gpn_; delete[] starting_gpn_; -delete[] all_patches_; for (int pn = N_patches()-1 ; pn >= 0 ; --pn) { delete &ith_patch(pn); } + +delete[] all_patches_; } //****************************************************************************** @@ -342,6 +363,8 @@ delete[] all_patches_; // of this patch system. This function also correctly sets // N_grid_points_ // N_ghosted_grid_points_ +// max_N_additional_points_ +// N_additional_points_ // all_patches_[] // starting_gpn_[] // ghosted_starting_gpn_[] @@ -353,6 +376,9 @@ void patch_system::create_patches(const struct patch_info patch_info_in[], int N_zones_per_right_angle, bool print_msg_flag) { +max_N_additional_points_ = 1; +// N_additional_points_ = 0; +N_additional_points_ = 1; N_grid_points_ = 0; ghosted_N_grid_points_ = 0; for (int pn = 0 ; pn < N_patches() ; ++pn) @@ -461,8 +487,8 @@ const int N_gridfns_in = jtutil::how_many_in_range(min_gfn_in, max_gfn_in); const int ghosted_N_gridfns_in = jtutil::how_many_in_range(ghosted_min_gfn_in, ghosted_max_gfn_in); -const int gfn_stride = N_grid_points(); -const int ghosted_gfn_stride = ghosted_N_grid_points(); +const int gfn_stride = N_grid_points() + max_N_additional_points(); +const int ghosted_gfn_stride = ghosted_N_grid_points() + max_N_additional_points(); const int N_storage = gfn_stride * N_gridfns_in; const int ghosted_N_storage = ghosted_gfn_stride * ghosted_N_gridfns_in; @@ -531,7 +557,7 @@ const patch& plast = ith_patch(N_patches()-1); const fp* const gfn1_first_ptr = & pfirst.gridfn(gfn+1, pfirst.min_irho(), pfirst.min_isigma()); - if (! (gfn1_first_ptr == gfn_last_ptr+1)) + if (! (gfn1_first_ptr == gfn_last_ptr+max_N_additional_points()+1)) then error_exit(PANIC_EXIT, "***** patch_system::setup_gridfn_storage():\n" " nominal-grid gridfns don't partition overall storage array!" @@ -561,7 +587,7 @@ const patch& plast = ith_patch(N_patches()-1); = & pfirst.ghosted_gridfn(ghosted_gfn+1, pfirst.ghosted_min_irho(), pfirst.ghosted_min_isigma()); - if (! (ghosted_gfn1_first_ptr == ghosted_gfn_last_ptr+1)) + if (! (ghosted_gfn1_first_ptr == ghosted_gfn_last_ptr+max_N_additional_points()+1)) then error_exit(PANIC_EXIT, "***** patch_system::setup_gridfn_storage():\n" " ghosted-grid gridfns don't partition overall storage array!" @@ -635,7 +661,7 @@ const patch& plast = ith_patch(N_patches()-1); if (! (p1_first_ptr == p_last_ptr+1)) then error_exit(PANIC_EXIT, "***** patch_system::setup_gridfn_storage():\n" -" nominal-grid patches gridfns don't partition storage for gfn=%d!\n" +" ghosted-grid patches gridfns don't partition storage for gfn=%d!\n" " (this should never happen!)\n" " %s patch (pn=%d) last point at p_last_ptr=%p\n" " %s patch (pn=%d) first point at p1_first_ptr=%p\n" @@ -1489,6 +1515,29 @@ void patch_system::set_gridfn_to_constant(fp a, int dst_gfn) //****************************************************************************** // +// This function sets a ghosted gridfn to a constant value. +// +void patch_system::set_ghosted_gridfn_to_constant(fp a, int ghosted_dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = ith_patch(pn); + for (int irho = p.ghosted_min_irho() ; + irho <= p.ghosted_max_irho() ; ++irho) + { + for (int isigma = p.ghosted_min_isigma() ; + isigma <= p.ghosted_max_isigma() ; + ++isigma) + { + p.ghosted_gridfn(ghosted_dst_gfn, irho,isigma) = a; + } + } + } +} + +//****************************************************************************** + +// // This function copies one (nominal-grid) gridfn to another. // void patch_system::gridfn_copy(int src_gfn, int dst_gfn) @@ -1511,6 +1560,78 @@ void patch_system::gridfn_copy(int src_gfn, int dst_gfn) //****************************************************************************** // +// This function copies one ghosted gridfn to another. +// +void patch_system::ghosted_gridfn_copy(int ghosted_src_gfn, + int ghosted_dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = ith_patch(pn); + for (int irho = p.ghosted_min_irho() ; + irho <= p.ghosted_max_irho() ; ++irho) + { + for (int isigma = p.ghosted_min_isigma() ; + isigma <= p.ghosted_max_isigma() ; + ++isigma) + { + p.ghosted_gridfn(ghosted_dst_gfn, irho,isigma) + = p.ghosted_gridfn(ghosted_src_gfn, irho,isigma); + } + } + } +} + +//****************************************************************************** + +// +// This function copies a ghosted gridfn to a nominal gridfn. +// +void patch_system::ghosted_gridfn_copy_to_nominal(int ghosted_src_gfn, + int dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = 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) + { + p.gridfn(dst_gfn, irho,isigma) + = p.ghosted_gridfn(ghosted_src_gfn, irho,isigma); + } + } + } +} + +//****************************************************************************** + +// +// This function adds a scalar to a (nominal-grid) gridfn. +// +void patch_system::add_to_gridfn(fp delta, int dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = 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) + { + p.gridfn(dst_gfn, irho,isigma) += delta; + } + } + } +} + +//****************************************************************************** + +// // This function adds a scalar to a ghosted gridfn. // void patch_system::add_to_ghosted_gridfn(fp delta, int ghosted_dst_gfn) @@ -1535,6 +1656,30 @@ void patch_system::add_to_ghosted_gridfn(fp delta, int ghosted_dst_gfn) //****************************************************************************** // +// This function scales a ghosted gridfn with a scalar. +// +void patch_system::scale_ghosted_gridfn(fp factor, int ghosted_dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = ith_patch(pn); + for (int irho = p.ghosted_min_irho() ; + irho <= p.ghosted_max_irho() ; + ++irho) + { + for (int isigma = p.ghosted_min_isigma() ; + isigma <= p.ghosted_max_isigma() ; + ++isigma) + { + p.ghosted_gridfn(ghosted_dst_gfn, irho,isigma) *= factor; + } + } + } +} + +//****************************************************************************** + +// // This function computes norms of a nominal-grid gridfn. // void patch_system::gridfn_norms(int src_gfn, jtutil::norm<fp>& norms) @@ -1757,6 +1902,69 @@ default: return integral; } +fp patch_system::integrate_gridpoint(int unknown_src_gfn, + int ghosted_radius_gfn, + int g_xx_gfn, int g_xy_gfn, int g_xz_gfn, + int g_yy_gfn, int g_yz_gfn, + int g_zz_gfn, + enum patch::integration_method method, + int pn, int irho, int isigma) + const +{ +const patch& p = ith_patch(pn); +const fp val = p.integrate_gridpoint(unknown_src_gfn, + ghosted_radius_gfn, + g_xx_gfn, g_xy_gfn, g_xz_gfn, + g_yy_gfn, g_yz_gfn, + g_zz_gfn, + method, + irho, isigma); + +return val; +} +fp patch_system::integrate_correction(bool src_gfn_is_even_across_xy_plane, + bool src_gfn_is_even_across_xz_plane, + bool src_gfn_is_even_across_yz_plane) + const +{ +fp integral = 1.0; +// +// correct the integral +// for the fact that the patch system may not cover the full 2-sphere +// +switch (type()) + { +case patch_system__full_sphere: + break; +case patch_system__plus_z_hemisphere: + integral *= src_gfn_is_even_across_xy_plane ? 2.0 : 0.0; + break; +case patch_system__plus_xy_quadrant_mirrored: +case patch_system__plus_xy_quadrant_rotating: + integral *= src_gfn_is_even_across_xz_plane ? 2.0 : 0.0; + integral *= src_gfn_is_even_across_yz_plane ? 2.0 : 0.0; + break; +case patch_system__plus_xz_quadrant_mirrored: +case patch_system__plus_xz_quadrant_rotating: + integral *= src_gfn_is_even_across_xy_plane ? 2.0 : 0.0; + integral *= src_gfn_is_even_across_yz_plane ? 2.0 : 0.0; + break; +case patch_system__plus_xyz_octant_mirrored: +case patch_system__plus_xyz_octant_rotating: + integral *= src_gfn_is_even_across_xy_plane ? 2.0 : 0.0; + integral *= src_gfn_is_even_across_xz_plane ? 2.0 : 0.0; + integral *= src_gfn_is_even_across_yz_plane ? 2.0 : 0.0; + break; +default: + error_exit(PANIC_EXIT, +"***** patch_system::integrate_gridfn(): bad patch system type()=(int)%d!\n" +" (this should never happen!)\n", + int(type())); /*NOTREACHED*/ + } + +return integral; +} + //****************************************************************************** //****************************************************************************** //****************************************************************************** @@ -1814,7 +2022,6 @@ else error_exit(ERROR_EXIT, // patch-system symmetries. If the point coincides with the local origin, // we return the dummy value 1.0. // -// Bugs: // Due to the surface-interpolator overhead, repeatedly calling this // function is rather inefficient. // @@ -1940,13 +2147,177 @@ if (status < 0) , status, double(x), double(y), double(z), - p.name(), rho, sigma, + p.name(), double(rho), double(sigma), double(p.drho_of_rho(rho)), double(p.dsigma_of_sigma(sigma))); /*NOTREACHED*/ return xyz_radius; } +struct perpatch { + std::vector<int> n; + std::vector<fp> rho, sigma; + std::vector<fp> radii; + void reserve (int const npoints) + { + n.reserve (npoints); + rho.reserve (npoints); + sigma.reserve (npoints); + radii.reserve (npoints); + } + void addpoint (int const n_, fp const rho_, fp const sigma_) + { + n .push_back (n_); + rho .push_back (rho_); + sigma.push_back (sigma_); + radii.push_back (fp (0.0)); + } +}; + +void patch_system::radii_in_local_xyz_directions(int const ghosted_radius_gfn, + int const npoints, + fp const * const xp, + fp const * const yp, + fp const * const zp, + fp * const radii) + const +{ + std::vector<perpatch> patchpoints (N_patches()); + for (int pn=0; pn<N_patches(); ++pn) { + patchpoints.at(pn).reserve (npoints); + } + + for (int n=0; n<npoints; ++n) { + fp x = xp[n]; + fp y = yp[n]; + fp z = zp[n]; + + if ((x == 0.0) && (y == 0.0) && (z == 0.0)) { + radii[n] = 1.0; + } else { + // + // apply symmetries to map (x,y,z) into that part of the + // 2-sphere which is covered by the patch system + // + switch (type()) + { + case patch_system__full_sphere: + break; + case patch_system__plus_z_hemisphere: + z = fabs(z); + break; + case patch_system__plus_xy_quadrant_mirrored: + case patch_system__plus_xy_quadrant_rotating: + x = fabs(x); + y = fabs(y); + break; + case patch_system__plus_xz_quadrant_mirrored: + case patch_system__plus_xz_quadrant_rotating: + x = fabs(x); + z = fabs(z); + break; + case patch_system__plus_xyz_octant_mirrored: + case patch_system__plus_xyz_octant_rotating: + x = fabs(x); + y = fabs(y); + z = fabs(z); + break; + default: + error_exit(PANIC_EXIT, +"***** patch_system::radii_in_local_xyz_directions():\n" +" unknown patch system type()=(int)%d!\n" +" (this should never happen!)\n", + int(type())); /*NOTREACHED*/ + } + + const patch* p_ptr = patch_containing_local_xyz(x, y, z); + if (p_ptr == NULL) + then error_exit(ERROR_EXIT, +"***** patch_system::radius_in_local_xyz_direction():\n" +" can't find containing patch!\n" +" (this should never happen!)\n" +" [local] (x,y,z)=(%g,%g,%g)\n" + , + double(x), double(y), double(z)); /*NOTREACHED*/ + + const patch& p = *p_ptr; + const fp rho = p. rho_of_xyz(x,y,z); + const fp sigma = p.sigma_of_xyz(x,y,z); + + patchpoints.at(p.patch_number()).addpoint (n, rho, sigma); + } + } // for n + + for (int pn=0; pn<N_patches(); ++pn) { + const patch& p = ith_patch(pn); + perpatch& localpatchpoints = patchpoints.at(pn); + + // Set up the surface interpolator to interpolate the surface + // radius gridfn to the (rho,sigma) coordinates: + // + // Notes on the interpolator setup: + // * The interpolator assumes Fortran subscripting, so we take the + // coordinates in the order (sigma,rho) to match our C + // subscripting in the patch system. + // * To avoid having to set up min/max array subscripts in the + // parameter table, we treat the patch as using 0-origin + // (integer) array subscripts (irho - ghosted_min_irho(), isigma + // - ghosted_min_isigma()). However, we use the usual + // floating-point coordinates. + // + + const int N_dims = 2; + const CCTK_REAL coord_origin[N_dims] + = { p.ghosted_min_sigma(), p.ghosted_min_rho() }; + const CCTK_REAL coord_delta[N_dims] + = { p.delta_sigma(), p.delta_rho() }; + + const int N_interp_points = localpatchpoints.radii.size(); + const int interp_coords_type_code = CCTK_VARIABLE_REAL; + const void* const interp_coords[N_dims] + = { static_cast<const void*>(& localpatchpoints.sigma.front()), + static_cast<const void*>(& localpatchpoints.rho .front()) }; + + const int N_input_arrays = 1; + const CCTK_INT input_array_dims[N_dims] + = { p.ghosted_N_isigma(), p.ghosted_N_irho() }; + const CCTK_INT input_array_type_codes[N_input_arrays] + = { CCTK_VARIABLE_REAL }; + const void* const input_arrays[N_input_arrays] + = { static_cast<const void*> + (p.ghosted_gridfn_data_array(ghosted_radius_gfn)) }; + + const int N_output_arrays = 1; + const CCTK_INT output_array_type_codes[N_output_arrays] + = { CCTK_VARIABLE_REAL }; + void* const output_arrays[N_output_arrays] + = { static_cast<void*>(& localpatchpoints.radii.front()) }; + + const int status = CCTK_InterpLocalUniform + (N_dims, + surface_interp_handle_, + surface_interp_par_table_handle_, + coord_origin, coord_delta, + N_interp_points, + interp_coords_type_code, interp_coords, + N_input_arrays, + input_array_dims, input_array_type_codes, input_arrays, + N_output_arrays, + output_array_type_codes, output_arrays); + if (status < 0) + then error_exit(ERROR_EXIT, +"***** patch_system::radii_in_local_xyz_directions():\n" +" error return (status=%d) from surface interpolator!\n" + , + status); /*NOTREACHED*/ + + for (size_t nn=0; nn<localpatchpoints.radii.size(); ++nn) { + radii[localpatchpoints.n.at(nn)] = localpatchpoints.radii.at(nn); + } + + } // for p +} + //****************************************************************************** //****************************************************************************** //****************************************************************************** @@ -2034,8 +2405,9 @@ fprintf(output_fp, "\n"); const fp dpy = p.dpy_of_rho_sigma(rho, sigma); fprintf(output_fp, "%g\t%g\t%#.15g", - dpx, dpy, p.unknown_gridfn(ghosted_flag, - unknown_gfn, irho,isigma)); + double(dpx), double(dpy), + double(p.unknown_gridfn(ghosted_flag, + unknown_gfn, irho,isigma))); if (print_xyz_flag) then { const fp r = p.unknown_gridfn(radius_is_ghosted_flag, @@ -2049,7 +2421,8 @@ fprintf(output_fp, "\n"); const fp global_z = origin_z() + local_z; fprintf(output_fp, "\t%#.10g\t%#.10g\t%#.10g", - global_x, global_y, global_z); + double(global_x), double(global_y), + double(global_z)); } fprintf(output_fp, "\n"); } @@ -2122,12 +2495,13 @@ int line_number = 1; " dpx=%g dpy=%g\n" , p.name(), unknown_gfn, - irho, p.effective_min_irho(want_ghost_zones), - p.effective_max_irho(want_ghost_zones), - isigma, + int(irho), + p.effective_min_irho(want_ghost_zones), + p.effective_max_irho(want_ghost_zones), + int(isigma), p.effective_min_isigma(want_ghost_zones), p.effective_max_isigma(want_ghost_zones), - dpx, dpy); /*NOTREACHED*/ + double(dpx), double(dpy)); /*NOTREACHED*/ ++line_number; } while ((buffer[0] == '#') || (buffer[0] == '\n')); @@ -2150,8 +2524,8 @@ int line_number = 1; , p.name(), unknown_gfn, line_number, - dpx, dpy, - read_dpx, read_dpy); /*NOTREACHED*/ + double(dpx), double(dpy), + double(read_dpx), double(read_dpy)); /*NOTREACHED*/ p.unknown_gridfn(ghosted_flag, unknown_gfn, irho,isigma) = read_gridfn_value; @@ -2164,6 +2538,298 @@ fclose(input_fp); } //****************************************************************************** + +#ifdef CCTK_HDF5 + +static void WriteAttribute (const hid_t dataset, const char* name, char value); +static void WriteAttribute (const hid_t dataset, const char* name, const char* values); +static void WriteAttribute (const hid_t dataset, const char* name, const char* values, int nvalues); + +static void WriteAttribute (const hid_t dataset, const char* const name, const char value) +{ + WriteAttribute (dataset, name, &value, 1); +} + +static void WriteAttribute (const hid_t dataset, const char* const name, const char* const values) +{ + WriteAttribute (dataset, name, values, strlen(values)); +} + +static void WriteAttribute (const hid_t dataset, const char* const name, const char* const values, const int nvalues) +{ + assert (dataset>=0); + assert (name); + assert (values); + assert (nvalues>=0); + + herr_t herr; + + const hid_t dataspace = H5Screate (H5S_SCALAR); + assert (dataspace>=0); + + const hid_t datatype = H5Tcopy (H5T_C_S1); + assert (datatype>=0); + herr = H5Tset_size (datatype, nvalues); + assert (!herr); + + const hid_t attribute = H5Acreate (dataset, name, datatype, dataspace, H5P_DEFAULT); + assert (attribute>=0); + herr = H5Awrite (attribute, datatype, values); + assert (!herr); + herr = H5Aclose (attribute); + assert (!herr); + + herr = H5Tclose (datatype); + assert (!herr); + + herr = H5Sclose (dataspace); + assert (!herr); +} + +#endif + +// +// This function output an unknown-grid gridfn in HDF5 format to a +// named output file. +// +void patch_system::output_unknown_gridfn + (bool ghosted_flag, int unknown_gfn, + const char gfn_name[], const cGH *cctkGH, + bool output_xyz_flag, bool radius_is_ghosted_flag, + int unknown_radius_gfn, + const char output_file_name[], bool want_ghost_zones) + const +{ +if (want_ghost_zones && !ghosted_flag) + then error_exit(PANIC_EXIT, +"***** patch_system::output_unknown_gridfn(unknown_gfn=%d):\n" +" can't have want_ghost_zones && !ghosted_flag !\n" +, + unknown_gfn); /*NOTREACHED*/ +if (want_ghost_zones && output_xyz_flag && !radius_is_ghosted_flag) + then error_exit(PANIC_EXIT, +"***** patch_system::output_unknown_gridfn(unknown_gfn=%d):\n" +" can't have want_ghost_zones && output_xyz_flag\n" +" && !radius_is_ghosted_flag!\n" +" unknown_radius_gfn=%d\n" +, + unknown_gfn, + unknown_radius_gfn); /*NOTREACHED*/ + +#ifdef CCTK_HDF5 + + using std::ostringstream; + using std::string; + using std::vector; + + herr_t herr; + + hid_t writer; + + // Get grid hierarchy extentsion from IOUtil + const ioGH * const iogh = (const ioGH *)CCTK_GHExtension (cctkGH, "IO"); + assert (iogh); + + struct stat fileinfo; + if (! iogh->recovered || stat(output_file_name, &fileinfo) != 0) { + writer = H5Fcreate (output_file_name, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + } else { + writer = H5Fopen (output_file_name, H5F_ACC_RDWR, H5P_DEFAULT); + } + assert (writer >= 0); + + hid_t datatype; + if (sizeof(fp) == sizeof(float)) { + datatype = H5T_NATIVE_FLOAT; + } else if (sizeof(fp) == sizeof(double)) { + datatype = H5T_NATIVE_DOUBLE; + } else if (sizeof(fp) == sizeof(long double)) { + datatype = H5T_NATIVE_LDOUBLE; + } else { + assert (0); + } + assert (datatype >= 0); + + for (int pn=0; pn<N_patches(); ++pn) { + const patch& p = ith_patch(pn); + + hsize_t shape[2]; + shape[0] = p.effective_N_isigma(want_ghost_zones); + shape[1] = p.effective_N_irho(want_ghost_zones); + + hid_t dataspace = H5Screate_simple (2, shape, NULL); + assert (dataspace >= 0); + + + + ostringstream datasetnamebuf; + datasetnamebuf << gfn_name + << " it=" << cctkGH->cctk_iteration + << " patch=" << pn; + string datasetnamestr = datasetnamebuf.str(); + assert (datasetnamestr.size() <= 256); // limit dataset name length + const char * const datasetname = datasetnamestr.c_str(); + assert (datasetname); + + ostringstream datasetnamebufx; + datasetnamebufx << "x" + << " it=" << cctkGH->cctk_iteration + << " patch=" << pn; + string datasetnamestrx = datasetnamebufx.str(); + assert (datasetnamestrx.size() <= 256); // limit dataset name length + const char * const datasetnamex = datasetnamestrx.c_str(); + assert (datasetnamex); + + ostringstream datasetnamebufy; + datasetnamebufy << "y" + << " it=" << cctkGH->cctk_iteration + << " patch=" << pn; + string datasetnamestry = datasetnamebufy.str(); + assert (datasetnamestry.size() <= 256); // limit dataset name length + const char * const datasetnamey = datasetnamestry.c_str(); + assert (datasetnamey); + + ostringstream datasetnamebufz; + datasetnamebufz << "z" + << " it=" << cctkGH->cctk_iteration + << " patch=" << pn; + string datasetnamestrz = datasetnamebufz.str(); + assert (datasetnamestrz.size() <= 256); // limit dataset name length + const char * const datasetnamez = datasetnamestrz.c_str(); + assert (datasetnamez); + + + + hid_t dataset = H5Dcreate (writer, datasetname, datatype, dataspace, H5P_DEFAULT); + assert (dataset >= 0); + + hid_t datasetx; + hid_t datasety; + hid_t datasetz; + if (output_xyz_flag) { + + datasetx = H5Dcreate (writer, datasetnamex, datatype, dataspace, H5P_DEFAULT); + assert (datasetx >= 0); + + datasety = H5Dcreate (writer, datasetnamey, datatype, dataspace, H5P_DEFAULT); + assert (datasety >= 0); + + datasetz = H5Dcreate (writer, datasetnamez, datatype, dataspace, H5P_DEFAULT); + assert (datasetz >= 0); + + } // if output_xyf_flag + + + + vector<fp> tmpdata; + vector<fp> tmpx, tmpy, tmpz; + tmpdata.resize (shape[0] * shape[1]); + if (output_xyz_flag) { + tmpx.resize (shape[0] * shape[1]); + tmpy.resize (shape[0] * shape[1]); + tmpz.resize (shape[0] * shape[1]); + } + + int elt = 0; + for (int irho = p.effective_min_irho(want_ghost_zones); + irho <= p.effective_max_irho(want_ghost_zones); + ++irho) + { + for (int isigma = p.effective_min_isigma(want_ghost_zones); + isigma <= p.effective_max_isigma(want_ghost_zones); + ++isigma) + { + + tmpdata.at(elt) = p.unknown_gridfn (ghosted_flag, unknown_gfn, irho,isigma); + if (output_xyz_flag) { + const fp r = p.unknown_gridfn(radius_is_ghosted_flag, + unknown_radius_gfn, + irho,isigma); + const fp rho = p.rho_of_irho(irho); + const fp sigma = p.sigma_of_isigma(isigma); + fp local_x, local_y, local_z; + p.xyz_of_r_rho_sigma(r,rho,sigma, + local_x,local_y,local_z); + const fp global_x = origin_x() + local_x; + const fp global_y = origin_y() + local_y; + const fp global_z = origin_z() + local_z; + tmpx.at(elt) = global_x; + tmpy.at(elt) = global_y; + tmpz.at(elt) = global_z; + } + ++elt; + + } + } + + + + herr = H5Dwrite (dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &tmpdata.front()); + assert (! herr); + + WriteAttribute (dataset, "name", gfn_name); + + herr = H5Dclose (dataset); + assert (! herr); + + + + if (output_xyz_flag) { + + herr = H5Dwrite (datasetx, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &tmpx.front()); + assert (! herr); + + WriteAttribute (datasetx, "name", "x"); + + herr = H5Dclose (datasetx); + assert (! herr); + + + + herr = H5Dwrite (datasety, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &tmpy.front()); + assert (! herr); + + WriteAttribute (datasety, "name", "y"); + + herr = H5Dclose (datasety); + assert (! herr); + + + + herr = H5Dwrite (datasetz, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &tmpz.front()); + assert (! herr); + + WriteAttribute (datasetz, "name", "z"); + + herr = H5Dclose (datasetz); + assert (! herr); + + } // if output_xyz_flag + + + + herr = H5Sclose (dataspace); + assert (! herr); + + } // for pn + + + + herr = H5Fclose (writer); + assert (! herr); + +#else + error_exit(ERROR_EXIT, +"***** patch_system::output_unknown_gridfn(unknown_gfn=%d):\n" +" no HDF5 support compiled in. Cannot write output file \"%s\"\n!" +, + unknown_gfn, + output_file_name); /*NOTREACHED*/ + +#endif +} + +//****************************************************************************** //****************************************************************************** //****************************************************************************** diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh index 38e5ea3..11cf665 100644 --- a/src/patch/patch_system.hh +++ b/src/patch/patch_system.hh @@ -104,6 +104,11 @@ public: fp origin_y() const { return global_coords_.origin_y(); } fp origin_z() const { return global_coords_.origin_z(); } + // set global (x,y,z) coordinates of local origin point + void origin_x(const fp x) { global_coords_.origin_x(x); } + void origin_y(const fp y) { global_coords_.origin_y(y); } + void origin_z(const fp z) { global_coords_.origin_z(z); } + // // ***** meta-info about the entire patch system ***** @@ -132,6 +137,9 @@ public: // total number of grid points int N_grid_points() const { return N_grid_points_; } int ghosted_N_grid_points() const { return ghosted_N_grid_points_; } + int max_N_additional_points() const { return max_N_additional_points_; } + int N_additional_points() const { return N_additional_points_; } + int N_additional_points(const int n) { return N_additional_points_=n; } // @@ -263,14 +271,21 @@ private: // public: // dst = a - void set_gridfn_to_constant(fp a, int dst_gfn); + void set_gridfn_to_constant(fp a, int dst_gfn); + void set_ghosted_gridfn_to_constant(fp a, int ghosted_dst_gfn); // dst = src - void gridfn_copy(int src_gfn, int dst_gfn); + void gridfn_copy(int src_gfn, int dst_gfn); + void ghosted_gridfn_copy(int ghosted_src_gfn, int ghosted_dst_gfn); + void ghosted_gridfn_copy_to_nominal(int ghosted_src_gfn, int dst_gfn); // dst += delta + void add_to_gridfn(fp delta, int dst_gfn); void add_to_ghosted_gridfn(fp delta, int ghosted_dst_gfn); + // dst *= factor + void scale_ghosted_gridfn(fp factor, int ghosted_dst_gfn); + // compute norms of gridfn (only over nominal grid) void gridfn_norms(int src_gfn, jtutil::norm<fp>& norms) const; @@ -296,15 +311,18 @@ public: // radius of surface in direction of an (x,y,z) point, // taking into account any patch-system symmetries; // or dummy value 1.0 if point is identical to local origin - // - // FIXME: - // We should provide another API to compute this for a whole - // batch of points at once, since this would be more efficient - // (the interpolator overhead would be amortized over the whole batch) fp radius_in_local_xyz_direction(int ghosted_radius_gfn, fp x, fp y, fp z) const; + void radii_in_local_xyz_directions(int ghosted_radius_gfn, + int npoints, + fp const * xp, + fp const * yp, + fp const * zp, + fp * radii) + const; + // // ***** line/surface operations ***** @@ -349,6 +367,20 @@ public: enum patch::integration_method method) const; + fp integrate_gridpoint(int unknown_src_gfn, + int ghosted_radius_gfn, + int g_xx_gfn, int g_xy_gfn, int g_xz_gfn, + int g_yy_gfn, int g_yz_gfn, + int g_zz_gfn, + enum patch::integration_method method, + int pn, int irho, int isigma) + const; + + fp integrate_correction(bool src_gfn_is_even_across_xy_plane, + bool src_gfn_is_even_across_xz_plane, + bool src_gfn_is_even_across_yz_plane) + const; + // @@ -427,6 +459,74 @@ private: const char input_file_name[], bool want_ghost_zones); +public: + // output to a named file + // output format is HDF5 + void output_gridfn(int gfn, + const char gfn_name[], const cGH *cctkGH, + const char output_file_name[]) const + { + output_unknown_gridfn(false, gfn, + gfn_name, cctkGH, + false, false, 0, + output_file_name, false); + } + void output_ghosted_gridfn(int ghosted_gfn, + const char gfn_name[], const cGH *cctkGH, + const char output_file_name[], + bool want_ghost_zones = true) + const + { + output_unknown_gridfn(true, ghosted_gfn, + gfn_name, cctkGH, + false, false, 0, + output_file_name, want_ghost_zones); + } + + // output to a named file (newly (re)created) + // output format is + // dpx dpy gridfn global_x global_y global_z + // where global_[xyz} are derived from the angular position + // and a specified (unknown-grid) radius gridfn + void output_gridfn_with_xyz + (int gfn, + const char gfn_name[], const cGH *cctkGH, + bool radius_is_ghosted_flag, int unknown_radius_gfn, + const char output_file_name[]) + const + { + output_unknown_gridfn(false, gfn, + gfn_name, cctkGH, + true, radius_is_ghosted_flag, + unknown_radius_gfn, + output_file_name, false); + } + void output_ghosted_gridfn_with_xyz + (int ghosted_gfn, + const char gfn_name[], const cGH *cctkGH, + bool radius_is_ghosted_flag, int unknown_radius_gfn, + const char output_file_name[], + bool want_ghost_zones = true) + const + { + output_unknown_gridfn(true, ghosted_gfn, + gfn_name, cctkGH, + true, radius_is_ghosted_flag, + unknown_radius_gfn, + output_file_name, want_ghost_zones); + } + +private: + // ... internal worker functions + void output_unknown_gridfn + (bool ghosted_flag, int unknown_gfn, + const char gfn_name[], const cGH *cctkGH, + bool print_xyz_flag, bool radius_is_ghosted_flag, + int unknown_radius_gfn, + const char output_file_name[], bool want_ghost_zones) + const; + + // // ***** access to gridfns as 1-D arrays ***** @@ -587,6 +687,8 @@ private: enum patch_system_type type_; int N_patches_; int N_grid_points_, ghosted_N_grid_points_; + int max_N_additional_points_; + int N_additional_points_; // [pn] = --> individual patches // *** constructor initialization list ordering: @@ -607,7 +709,8 @@ private: fp* ghosted_gridfn_storage_; // min/max m over all ghost zone points - mutable int global_min_ym_, global_max_ym_; + mutable int global_min_ym_; + mutable int global_max_ym_; // info about the surface interpolator // ... used only by radius_in_local_xyz_direction() diff --git a/src/patch/test_coords.cc b/src/patch/test_coords.cc index bf40c74..0abe7d8 100644 --- a/src/patch/test_coords.cc +++ b/src/patch/test_coords.cc @@ -15,14 +15,14 @@ #include "config.h" #include "stdc.h" #include "../jtutil/util.hh" +using jtutil::error_exit; +using jtutil::radians_of_degrees; +using jtutil::degrees_of_radians; #include "coords.hh" using namespace AHFinderDirect; using namespace local_coords; -using jtutil::error_exit; -using jtutil::radians_of_degrees; -using jtutil::degrees_of_radians; //****************************************************************************** diff --git a/src/patch/test_coords2.cc b/src/patch/test_coords2.cc index cbf7a0a..bc53633 100644 --- a/src/patch/test_coords2.cc +++ b/src/patch/test_coords2.cc @@ -18,14 +18,14 @@ #include "coords.hh" -using namespace AHFinderDirect; -using namespace local_coords; - using jtutil::fuzzy; using jtutil::error_exit; using jtutil::radians_of_degrees; using jtutil::degrees_of_radians; +using namespace AHFinderDirect; +using namespace local_coords; + // prototypes namespace { void test_r_mu_nu_phi(fp x, fp y, fp z); diff --git a/src/patch/test_patch_system.cc b/src/patch/test_patch_system.cc index cad9711..64ff79e 100644 --- a/src/patch/test_patch_system.cc +++ b/src/patch/test_patch_system.cc @@ -41,6 +41,7 @@ #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; #include "coords.hh" #include "grid.hh" @@ -52,7 +53,6 @@ #include "patch_system.hh" using namespace AHFinderDirect; -using jtutil::error_exit; //****************************************************************************** |