aboutsummaryrefslogtreecommitdiff
path: root/src/patch
diff options
context:
space:
mode:
Diffstat (limited to 'src/patch')
-rw-r--r--src/patch/coords.cc80
-rw-r--r--src/patch/coords.hh5
-rw-r--r--src/patch/fd_grid.cc2
-rw-r--r--src/patch/fd_grid.hh2
-rw-r--r--src/patch/ghost_zone.cc2
-rw-r--r--src/patch/ghost_zone.hh3
-rw-r--r--src/patch/make.code.defn14
-rw-r--r--src/patch/patch.cc57
-rw-r--r--src/patch/patch.hh9
-rw-r--r--src/patch/patch_info.cc2
-rw-r--r--src/patch/patch_interp.cc8
-rw-r--r--src/patch/patch_system.cc702
-rw-r--r--src/patch/patch_system.hh119
-rw-r--r--src/patch/test_coords.cc6
-rw-r--r--src/patch/test_coords2.cc6
-rw-r--r--src/patch/test_patch_system.cc2
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;
//******************************************************************************
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-10 07:49:56 +0000