From 720ce53a79ef5bdff53b3db18bc45c46101aa0fa Mon Sep 17 00:00:00 2001 From: jthorn Date: Fri, 26 Jul 2002 14:57:48 +0000 Subject: re-sync changes from laptop - parameter to control how Jacobian is computed - can hardwire geometry to Schwarzschild/EF git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@661 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/gr/driver.cc | 180 ++++++++++++++++++++++++++++++++++--------------------- 1 file changed, 112 insertions(+), 68 deletions(-) (limited to 'src/gr/driver.cc') 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 @@ // // <<>> // 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(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 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 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, @@ -285,6 +303,31 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, method); /*NOTREACHED*/ } +//****************************************************************************** + +// +// 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)); } } } -- cgit v1.2.3