aboutsummaryrefslogtreecommitdiff
path: root/src/gr/driver.cc
diff options
context:
space:
mode:
Diffstat (limited to 'src/gr/driver.cc')
-rw-r--r--src/gr/driver.cc180
1 files changed, 112 insertions, 68 deletions
diff --git a/src/gr/driver.cc b/src/gr/driver.cc
index 00695e9..a9b9a13 100644
--- a/src/gr/driver.cc
+++ b/src/gr/driver.cc
@@ -3,6 +3,7 @@
//
// <<<prototypes for functions local to this file>>>
// AHFinderDirect_driver - top-level driver
+/// decode_Jacobian_method - decode the Jacobian_method parameter
///
/// setup_Kerr_horizon - set up Kerr horizon in h (Kerr or Kerr-Schild coords)
/// setup_ellipsoid - setup up a coordiante ellipsoid in h
@@ -39,7 +40,7 @@ using jtutil::error_exit;
#include "../elliptic/Jacobian.hh"
-#include "gfn.hh"
+#include "gfns.hh"
#include "AHFinderDirect.hh"
//******************************************************************************
@@ -49,8 +50,10 @@ using jtutil::error_exit;
//
namespace {
+enum Jacobian_method
+ decode_Jacobian_method(const char Jacobian_method_string[]);
void setup_Kerr_horizon(patch_system& ps,
- fp x_center, fp y_center, fp z_center,
+ fp x_posn, fp y_posn, fp z_posn,
fp m, fp a,
bool Kerr_Schild_flag);
void setup_ellipsoid(patch_system& ps,
@@ -58,7 +61,7 @@ void setup_ellipsoid(patch_system& ps,
fp x_radius, fp y_radius, fp z_radius);
void print_Jacobians(const patch_system& ps,
- const Jacobian& Jac, const Jacobian& NP_Jac,
+ const Jacobian& Jac_NP, const Jacobian& Jac_SD_FDdr,
const char file_name[]);
};
@@ -77,18 +80,25 @@ CCTK_VInfo(CCTK_THORNSTRING, "initializing AHFinderDirect data structures");
//
-// set up the geometry interpolator
+// set up the geometry parameters and the geometry interpolator
//
-struct geometry_interpolator_info gii;
+struct geometry_info gi;
+gi.hardwire_Schwarzschild_EF = (hardwire_Schwarzschild_EF != 0);
+gi.hardwire_Schwarzschild_EF__x_posn = hardwire_Schwarzschild_EF__x_posn;
+gi.hardwire_Schwarzschild_EF__y_posn = hardwire_Schwarzschild_EF__y_posn;
+gi.hardwire_Schwarzschild_EF__z_posn = hardwire_Schwarzschild_EF__z_posn;
+gi.hardwire_Schwarzschild_EF__mass = hardwire_Schwarzschild_EF__mass;
+gi.hardwire_Schwarzschild_EF__epsilon = hardwire_Schwarzschild_EF__epsilon;
+gi.Delta_xyz = hardwire_Schwarzschild_EF__Delta_xyz;
+
CCTK_VInfo(CCTK_THORNSTRING, " setting up geometry interpolator");
-gii.operator_handle = CCTK_InterpHandle(geometry_interpolator_name);
-if (gii.operator_handle < 0)
+gi.operator_handle = CCTK_InterpHandle(geometry_interpolator_name);
+if (gi.operator_handle < 0)
then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"couldn't find interpolator \"%s\"!",
geometry_interpolator_name); /*NOTREACHED*/
-
-gii.param_table_handle = Util_TableCreateFromString(geometry_interpolator_pars);
-if (gii.param_table_handle < 0)
+gi.param_table_handle = Util_TableCreateFromString(geometry_interpolator_pars);
+if (gi.param_table_handle < 0)
then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"bad geometry-interpolator parameter(s) \"%s\"!",
geometry_interpolator_pars); /*NOTREACHED*/
@@ -150,8 +160,8 @@ cgi.K_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kzz"));
patch_system ps(origin_x, origin_y, origin_z,
patch_system::type_of_name(patch_system_type),
N_ghost_points, N_overlap_points, delta_drho_dsigma,
- nominal_gfns::min_gfn, nominal_gfns::max_gfn,
- ghosted_gfns::min_gfn, ghosted_gfns::max_gfn,
+ gfns::nominal_min_gfn, gfns::nominal_max_gfn,
+ gfns::ghosted_min_gfn, gfns::ghosted_max_gfn,
interp_handle, interp_param_table_handle);
@@ -163,7 +173,7 @@ if (STRING_EQUAL(initial_guess_method, "read from file"))
CCTK_VInfo(CCTK_THORNSTRING,
"reading initial guess h from \"%s\"",
h_file_name);
- ps.read_ghosted_gridfn(ghosted_gfns::gfn__h,
+ ps.read_ghosted_gridfn(gfns::gfn__h,
h_file_name,
false); // no ghost zones
}
@@ -178,19 +188,19 @@ if (STRING_EQUAL(initial_guess_method, "read from file"))
initial_guess__ellipsoid__z_radius);
else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr"))
then setup_Kerr_horizon(ps,
- initial_guess__Kerr_KerrSchild__x_center,
- initial_guess__Kerr_KerrSchild__y_center,
- initial_guess__Kerr_KerrSchild__z_center,
- initial_guess__Kerr_KerrSchild__mass,
- initial_guess__Kerr_KerrSchild__spin,
+ initial_guess__Kerr__x_posn,
+ initial_guess__Kerr__y_posn,
+ initial_guess__Kerr__z_posn,
+ initial_guess__Kerr__mass,
+ initial_guess__Kerr__spin,
false); // use Kerr coords
else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr-Schild"))
then setup_Kerr_horizon(ps,
- initial_guess__Kerr_KerrSchild__x_center,
- initial_guess__Kerr_KerrSchild__y_center,
- initial_guess__Kerr_KerrSchild__z_center,
- initial_guess__Kerr_KerrSchild__mass,
- initial_guess__Kerr_KerrSchild__spin,
+ initial_guess__Kerr__x_posn,
+ initial_guess__Kerr__y_posn,
+ initial_guess__Kerr__z_posn,
+ initial_guess__Kerr__mass,
+ initial_guess__Kerr__spin,
true); // use Kerr-Schild coords
else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"unknown initial_guess_method=\"%s\"!",
@@ -200,31 +210,36 @@ if (STRING_EQUAL(initial_guess_method, "read from file"))
CCTK_VInfo(CCTK_THORNSTRING,
"writing initial guess h to \"%s\"",
h_file_name);
- ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h,
- true, ghosted_gfns::gfn__h,
+ ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h,
+ true, gfns::gfn__h,
h_file_name,
false); // no ghost zones
}
//
-// find the apparent horizon
+// decode/copy parameters into structures
//
-struct Jacobian_info Jac_info;
-Jac_info.perturbation_amplitude = NP_Jacobian__perturbation_amplitude;
+
+struct Jacobian_info Jacobian_info;
+Jacobian_info.Jacobian_method = decode_Jacobian_method(Jacobian_method);
+Jacobian_info.perturbation_amplitude = Jacobian_perturbation_amplitude;
struct solver_info solver_info;
solver_info.max_Newton_iterations = max_Newton_iterations;
solver_info.H_norm_for_convergence = H_norm_for_convergence;
solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence;
+//
+// find the apparent horizon
+//
if (STRING_EQUAL(method, "horizon function"))
then {
jtutil::norm<fp> H_norms;
const int timer_handle = CCTK_TimerCreateI();
CCTK_TimerStartI(timer_handle);
- horizon_function(ps, cgi, gii, false, &H_norms);
+ horizon_function(ps, cgi, gi, false, true, &H_norms);
CCTK_TimerStopI(timer_handle);
printf("timer stats for evaluating H(h):\n");
CCTK_TimerPrintDataI(timer_handle, -1);
@@ -235,33 +250,36 @@ if (STRING_EQUAL(method, "horizon function"))
CCTK_VInfo(CCTK_THORNSTRING,
"writing H(h) to \"%s\"",
H_of_h_file_name);
- ps.print_gridfn_with_xyz(nominal_gfns::gfn__H,
- true, ghosted_gfns::gfn__h,
+ ps.print_gridfn_with_xyz(gfns::gfn__H,
+ true, gfns::gfn__h,
H_of_h_file_name);
}
else if (STRING_EQUAL(method, "Jacobian test"))
then {
- jtutil::norm<fp> H_norms;
- Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
- horizon_function(ps, cgi, gii, true, &H_norms);
- horizon_Jacobian(ps, cgi, gii, Jac_info, Jac);
-
- Jacobian& NP_Jac = create_Jacobian(ps, Jacobian_type);
- horizon_function(ps, cgi, gii, true, &H_norms);
- horizon_Jacobian_NP(ps, cgi, gii, Jac_info, NP_Jac);
+ // symbolic differentiation with finite diff d/dr
+ Jacobian& Jac_SD_FDdr = new_Jacobian(ps, Jacobian_storage_method);
+ horizon_function(ps, cgi, gi);
+ Jacobian_info.Jacobian_method = symbolic_differentiation_with_FD_dr;
+ horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac_SD_FDdr);
+
+ // numerical perturbation
+ Jacobian& Jac_NP = new_Jacobian(ps, Jacobian_storage_method);
+ horizon_function(ps, cgi, gi);
+ Jacobian_info.Jacobian_method = numerical_perturbation;
+ horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac_NP);
CCTK_VInfo(CCTK_THORNSTRING,
"writing Jacobians to \"%s\"",
Jacobian_file_name);
- print_Jacobians(ps, Jac, NP_Jac, Jacobian_file_name);
+ print_Jacobians(ps, Jac_NP, Jac_SD_FDdr, Jacobian_file_name);
}
else if (STRING_EQUAL(method, "Newton solve"))
then {
- Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
+ Jacobian& Jac = new_Jacobian(ps, Jacobian_storage_method);
const int timer_handle = CCTK_TimerCreateI();
CCTK_TimerStartI(timer_handle);
- Newton_solve(ps, cgi, gii, Jac_info, solver_info, Jac);
+ Newton_solve(ps, cgi, gi, Jacobian_info, solver_info, Jac);
CCTK_TimerStopI(timer_handle);
printf("timer stats for finding the apparent horizon:\n");
CCTK_TimerPrintDataI(timer_handle, -1);
@@ -269,15 +287,15 @@ else if (STRING_EQUAL(method, "Newton solve"))
CCTK_VInfo(CCTK_THORNSTRING,
"writing final horizon shape h to \"%s\"",
h_file_name);
- ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h,
- true, ghosted_gfns::gfn__h,
+ ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h,
+ true, gfns::gfn__h,
h_file_name,
false); // no ghost zones
CCTK_VInfo(CCTK_THORNSTRING,
"writing H(h) to \"%s\"",
H_of_h_file_name);
- ps.print_gridfn_with_xyz(nominal_gfns::gfn__H,
- true, ghosted_gfns::gfn__h,
+ ps.print_gridfn_with_xyz(gfns::gfn__H,
+ true, gfns::gfn__h,
H_of_h_file_name);
}
else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
@@ -286,6 +304,31 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
}
//******************************************************************************
+
+//
+// This function decodes the Jacobian_method parameter (string) into
+// an internal enum for future use.
+//
+namespace {
+enum Jacobian_method
+ decode_Jacobian_method(const char Jacobian_method_string[])
+{
+if (STRING_EQUAL(Jacobian_method_string,
+ "numerical perturbation"))
+ then return numerical_perturbation;
+else if (STRING_EQUAL(Jacobian_method_string,
+ "symbolic differentiation with finite diff d/dr"))
+ then return symbolic_differentiation_with_FD_dr;
+else if (STRING_EQUAL(Jacobian_method_string,
+ "symbolic differentiation"))
+ then return symbolic_differentiation;
+else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "unknown Jacobian_method_string=\"%s\"!",
+ Jacobian_method_string); /*NOTREACHED*/
+}
+ }
+
+//******************************************************************************
//******************************************************************************
//******************************************************************************
@@ -297,7 +340,7 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
// and the horizon is worked out on page 13.2 of my AHFinderDirect notes.
//
// Arguments:
-// [xyz]_center = The position of the Kerr black hole.
+// [xyz]_posn = The position of the Kerr black hole.
// (m,a) = Describe the Kerr black hole. Note that my convention has
// a=J/m^2 dimensionless, while MTW take a=J/m=m*(my a).
// Kerr_Schild_flag = false to use Kerr coordinates,
@@ -305,7 +348,7 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
//
namespace {
void setup_Kerr_horizon(patch_system& ps,
- fp x_center, fp y_center, fp z_center,
+ fp x_posn, fp y_posn, fp z_posn,
fp m, fp a,
bool Kerr_Schild_flag)
{
@@ -315,9 +358,9 @@ CCTK_VInfo(CCTK_THORNSTRING,
"setting up horizon for Kerr in %s coords",
name);
CCTK_VInfo(CCTK_THORNSTRING,
- " mass=%g, spin=J/m^2=%g, posn=(%g,%g,%g)",
- double(m), double(a),
- double(x_center), double(y_center), double(z_center));
+ " posn=(%g,%g,%g) mass=%g spin=J/m^2=%g",
+ double(x_posn), double(y_posn), double(z_posn),
+ double(m), double(a));
// horizon in Kerr coordinates is coordinate sphere
const fp r = m * (1.0 + sqrt(1.0 - a*a));
@@ -330,7 +373,7 @@ CCTK_VInfo(CCTK_THORNSTRING,
" horizon is coordinate %s",
Kerr_Schild_flag ? "ellipsoid" : "sphere");
setup_ellipsoid(ps,
- x_center, y_center, z_center,
+ x_posn, y_posn, z_posn,
xy_radius, xy_radius, z_radius);
}
}
@@ -417,7 +460,7 @@ CCTK_VInfo(CCTK_THORNSTRING,
/*NOTREACHED*/
// r = horizon radius at this grid point
- p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) = r;
+ p.ghosted_gridfn(gfns::gfn__h, irho,isigma) = r;
}
}
}
@@ -434,7 +477,7 @@ CCTK_VInfo(CCTK_THORNSTRING,
//
namespace {
void print_Jacobians(const patch_system& ps,
- const Jacobian& Jac, const Jacobian& NP_Jac,
+ const Jacobian& Jac_NP, const Jacobian& Jac_SD_FDdr,
const char file_name[])
{
FILE *fileptr = fopen(file_name, "w");
@@ -451,10 +494,10 @@ fprintf(fileptr, "# column 5 = y JJ\n");
fprintf(fileptr, "# column 6 = y patch number\n");
fprintf(fileptr, "# column 7 = y irho\n");
fprintf(fileptr, "# column 8 = y isigma\n");
-fprintf(fileptr, "# column 9 = Jac(II,JJ)\n");
-fprintf(fileptr, "# column 10 = NP_Jac(II,JJ)\n");
-fprintf(fileptr, "# column 11 = abs error = Jac - NP_Jac\n");
-fprintf(fileptr, "# column 12 = rel error = abs error / max(|Jac|,|NP_Jac|)\n");
+fprintf(fileptr, "# column 9 = Jac_NP(II,JJ)\n");
+fprintf(fileptr, "# column 10 = Jac_SD_FDdr(II,JJ)\n");
+fprintf(fileptr, "# column 11 = abs error in Jac_SD_FDdr\n");
+fprintf(fileptr, "# column 12 = rel error in Jac_SD_FDdr\n");
for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
{
@@ -482,26 +525,27 @@ fprintf(fileptr, "# column 12 = rel error = abs error / max(|Jac|,|NP_Jac|)\n");
{
const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma);
- if (! Jac.is_explicitly_stored(II,JJ))
+ if (! Jac_SD_FDdr.is_explicitly_stored(II,JJ))
then continue; // skip sparse points
- const fp SD = Jac(II,JJ);
- const fp NP = NP_Jac(II,JJ);
- const fp abs_error = SD - NP;
+ const fp NP = Jac_NP (II,JJ);
+ const fp SD_FDdr = Jac_SD_FDdr(II,JJ);
- if ((SD == 0.0) && (NP == 0.0))
+ if ((NP == 0.0) && (SD_FDdr == 0.0))
then continue; // skip zero values (if == )
- const fp abs_SD = jtutil::abs(SD);
- const fp abs_NP = jtutil::abs(NP);
- const fp rel_error = abs_error / jtutil::max(abs_SD, abs_NP);
+ const fp abs_NP = jtutil::abs(NP );
+ const fp abs_SD_FDdr = jtutil::abs(SD_FDdr);
+ const fp max_abs = jtutil::max(abs_SD_FDdr, abs_NP);
+ const fp SD_FDdr_abs_error = SD_FDdr - NP;
+ const fp SD_FDdr_rel_error = SD_FDdr_abs_error / max_abs;
fprintf(fileptr,
"%d %d %d %d\t%d %d %d %d\t%.10g\t%.10g\t%e\t%e\n",
II, xpn, x_irho, x_isigma,
JJ, ypn, y_irho, y_isigma,
- double(SD), double(NP),
- double(abs_error), double(rel_error));
+ double(NP), double(SD_FDdr),
+ double(SD_FDdr_abs_error), double(SD_FDdr_rel_error));
}
}
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-04 16:55:49 +0000