aboutsummaryrefslogtreecommitdiff
path: root/src/driver/aliased_functions.cc
diff options
context:
space:
mode:
Diffstat (limited to 'src/driver/aliased_functions.cc')
-rw-r--r--src/driver/aliased_functions.cc122
1 files changed, 60 insertions, 62 deletions
diff --git a/src/driver/aliased_functions.cc b/src/driver/aliased_functions.cc
index 43c15ed..8b1218a 100644
--- a/src/driver/aliased_functions.cc
+++ b/src/driver/aliased_functions.cc
@@ -12,10 +12,12 @@
#include <assert.h>
#include <math.h>
+#include <vector>
+
#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"
-#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
#include "config.h"
#include "stdc.h"
@@ -23,6 +25,7 @@
#include "../jtutil/array.hh"
#include "../jtutil/cpm_map.hh"
#include "../jtutil/linear_map.hh"
+using jtutil::error_exit;
#include "../patch/coords.hh"
#include "../patch/grid.hh"
@@ -45,7 +48,6 @@
// all the code in this file is inside this namespace
namespace AHFinderDirect
{
-using jtutil::error_exit;
//******************************************************************************
@@ -57,13 +59,8 @@ extern struct state state;
//******************************************************************************
//
-// This function is called (via the Cactus flesh function-aliasing mechanism)
-// by other thorns to find out our local coordinate origin for a given AH.
-//
-// Arguments:
-// horizon_number = (in) The horizon to inquire about; must be in the range
-// 1 to N_horizons inclusive.
-// *p_origin_[xyz] = (out) The local coordinate origin is stored here.
+// This function is called (via the magic of function aliasing) by
+// other thorns to find out our local coordinate origin for a given AH.
//
// Results:
// This function returns 0 for ok, or -1 if the horizon number is invalid.
@@ -77,7 +74,7 @@ const struct verbose_info& verbose_info = state.verbose_info;
if (! ((horizon_number >= 1) && (horizon_number <= state.N_horizons)) )
then {
- CCTK_VWarn(SERIOUS_WARNING, __LINE__, __FILE__, CCTK_THORNSTRING,
+ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"AHFinderDirect_local_coordinate_origin():\n"
" horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
,
@@ -108,10 +105,6 @@ return 0; // *** NORMAL RETURN ***
// by other thorns to query whether or not the specified horizon was found
// the last time we searched for it.
//
-// Arguments:
-// horizon_number = (in) The horizon to inquire about; must be in the range
-// 1 to N_horizons inclusive.
-//
// Results:
// This function returns
// 1 if the horizon was found
@@ -125,7 +118,7 @@ const struct verbose_info& verbose_info = state.verbose_info;
if (! ((horizon_number >= 1) && (horizon_number <= state.N_horizons)) )
then {
- CCTK_VWarn(SERIOUS_WARNING, __LINE__, __FILE__, CCTK_THORNSTRING,
+ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"AHFinderDirect_horizon_was_found():\n"
" horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
,
@@ -202,26 +195,16 @@ if (AH_data.found_flag)
// its local coordinate origin to a given (x,y,z) coordinate or coordinates.
//
// Arguments:
-// horizon_number = (in) The horizon to inquire about; must be in the range
-// 1 to N_horizons inclusive.
-// N_points = (in) Number of (x,y,z) coordinate tuples (should be >= 0).
-// x[], y[], z[] = (in) These arrays give the (x,y,z) coordinates.
-// radius[] = (out) This array is set to the horizon radius values
-// (Euclidean distances from the local coordinate origin),
-// or to all -1.0 if we didn't find this horizon the
-// last time we looked for it.
+// horizon_number = must be in the range 1 to N_horizons
+// N_points = should be >= 0
+// x[], y[], z[] = these give the (x,y,z) coordinates
+// radius[] = this is set to the horizon radius values (Euclidean distance
+// from the local coordinate origin), or to all -1.0 if we didn't
+// find this horizon the last time we looked for it
//
// Results:
// This function returns 0 for ok, or -1 if the horizon number is invalid.
//
-// FIXME:
-// * The current implementation is quite inefficient: it does a full 2-D
-// interpolation (to find the horizon radius) for each of the xyz points.
-// It would be more efficient to batch the interpolations.
-// * If an (x,y,z) coordinate tuple coincides with this horizon's local
-// coordinate origin, we return the dummy radius 1.0. It would be
-// better to return -1.0.
-//
extern "C"
CCTK_INT AHFinderDirect_radius_in_direction
(CCTK_INT horizon_number,
@@ -229,37 +212,52 @@ extern "C"
const CCTK_REAL* const x, const CCTK_REAL* const y, const CCTK_REAL* const z,
CCTK_REAL* const radius)
{
-const struct verbose_info& verbose_info = state.verbose_info;
-
-if (! ((horizon_number >= 1) && (horizon_number <= state.N_horizons)) )
- then {
- CCTK_VWarn(SERIOUS_WARNING, __LINE__, __FILE__, CCTK_THORNSTRING,
-"AHFinderDirect_distance_outside_thorn():\n"
-" horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
- ,
- int(horizon_number), state.N_horizons);
- return -1; // *** ERROR RETURN ***
- }
-
-assert(state.AH_data_array[horizon_number] != NULL);
-const struct AH_data& AH_data = *state.AH_data_array[horizon_number];
-
- for (int point = 0 ; point < N_points ; ++point)
- {
- if (AH_data.found_flag)
- then {
- assert(AH_data.ps_ptr != NULL);
- const patch_system& ps = *AH_data.ps_ptr;
-
- const fp local_x = x[point] - ps.origin_x();
- const fp local_y = y[point] - ps.origin_y();
- const fp local_z = z[point] - ps.origin_z();
- radius[point] = ps.radius_in_local_xyz_direction
- (gfns::gfn__h,
- local_x, local_y, local_z);
- }
- else radius[point] = -1.0;
- }
+ const struct verbose_info& verbose_info = state.verbose_info;
+
+ if (! ((horizon_number >= 1) && (horizon_number <= state.N_horizons)) ) {
+ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "AHFinderDirect_distance_outside_thorn():\n"
+ " horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
+ ,
+ int(horizon_number), state.N_horizons);
+ return -1; // *** ERROR RETURN ***
+ }
+
+ assert(state.AH_data_array[horizon_number] != NULL);
+ const struct AH_data& AH_data = *state.AH_data_array[horizon_number];
+
+ if (AH_data.found_flag) {
+
+ assert(AH_data.ps_ptr != NULL);
+ const patch_system& ps = *AH_data.ps_ptr;
+
+ std::vector<fp> local_xs(N_points);
+ std::vector<fp> local_ys(N_points);
+ std::vector<fp> local_zs(N_points);
+
+ for (int point = 0 ; point < N_points ; ++point) {
+
+ local_xs.at(point) = x[point] - ps.origin_x();
+ local_ys.at(point) = y[point] - ps.origin_y();
+ local_zs.at(point) = z[point] - ps.origin_z();
+
+ }
+
+ ps.radii_in_local_xyz_directions (gfns::gfn__h,
+ N_points,
+ & local_xs.front(),
+ & local_ys.front(),
+ & local_zs.front(),
+ radius);
+
+ } else {
+ // if not found
+
+ for (int point = 0 ; point < N_points ; ++point) {
+ radius[point] = -1.0;
+ }
+
+ } // if not found
return 0; // *** NORMAL RETURN ***
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 04:16:00 +0000