diff options
author | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2002-08-12 15:18:56 +0000 |
---|---|---|
committer | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2002-08-12 15:18:56 +0000 |
commit | d291c9597a25dc140f5a926ba4466f3407709dc8 (patch) | |
tree | 3c57f66b2b22796fa3ff3d94603038d85d5727cf /src/gr | |
parent | 211d94c18235f7e95ad2f0dc70047777c6b2d48a (diff) |
add an option to do the NP Jacobian only for testing
(useful if SD Jacobian core-dumps)
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@694 f88db872-0e4f-0410-b76b-b9085cfa78c5
Diffstat (limited to 'src/gr')
-rw-r--r-- | src/gr/driver.cc | 148 |
1 files changed, 90 insertions, 58 deletions
diff --git a/src/gr/driver.cc b/src/gr/driver.cc index cf26d0c..5da7b3f 100644 --- a/src/gr/driver.cc +++ b/src/gr/driver.cc @@ -65,7 +65,8 @@ 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_NP, const Jacobian& Jac_SD_FDdr, + const Jacobian* Jac_NP, + bool do_non_NP_Jac, const Jacobian* Jac_SD_FDdr, const char file_name[]); }; @@ -267,24 +268,36 @@ if (STRING_EQUAL(method, "horizon function")) true, gfns::gfn__h, H_of_h_file_name); } -else if (STRING_EQUAL(method, "Jacobian test")) +else if ( STRING_EQUAL(method, "Jacobian test") + || STRING_EQUAL(method, "Jacobian test (NP only)") ) then { - // 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 = Jacobian__symbolic_diff_with_FD_dr; - horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac_SD_FDdr); + const bool NP_only = STRING_EQUAL(method, "Jacobian test (NP only)"); + const bool do_non_NP_Jac = ! NP_only; // numerical perturbation - Jacobian& Jac_NP = new_Jacobian(ps, Jacobian_storage_method); - horizon_function(ps, cgi, gi); + Jacobian* const Jac_NP = & new_Jacobian(ps, Jacobian_storage_method); + horizon_function(ps, cgi, gi, true); Jacobian_info.Jacobian_method = Jacobian__numerical_perturb; - horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac_NP); + horizon_Jacobian(ps, cgi, gi, Jacobian_info, *Jac_NP); + + // symbolic differentiation with finite diff d/dr + Jacobian* Jac_SD_FDdr = NULL; + if (do_non_NP_Jac) + then { + Jac_SD_FDdr = & new_Jacobian(ps, Jacobian_storage_method); + horizon_function(ps, cgi, gi, true); + Jacobian_info.Jacobian_method + = Jacobian__symbolic_diff_with_FD_dr; + horizon_Jacobian(ps, cgi, gi, Jacobian_info, *Jac_SD_FDdr); + } CCTK_VInfo(CCTK_THORNSTRING, - "writing Jacobians to \"%s\"", - Jacobian_file_name); - print_Jacobians(ps, Jac_NP, Jac_SD_FDdr, Jacobian_file_name); + "writing Jacobian%s to \"%s\"", + (NP_only ? "" : "s"), Jacobian_file_name); + print_Jacobians(ps, + Jac_NP, + do_non_NP_Jac, Jac_SD_FDdr, + Jacobian_file_name); } else if (STRING_EQUAL(method, "Newton solve")) then { @@ -505,12 +518,17 @@ CCTK_VInfo(CCTK_THORNSTRING, //****************************************************************************** // -// This function prints a pair of Jacobian matrices (and their difference) -// to a named output file. +// This function prints one or two Jacobian matrices (and their difference +// in the latter case) to a named output file. +// +// Arguments: +// do_non_NP_Jac = true if we should also print *Jac_SD_FDdr +// false to skip this (and print only *Jac_NP) // namespace { void print_Jacobians(const patch_system& ps, - const Jacobian& Jac_NP, const Jacobian& Jac_SD_FDdr, + const Jacobian* Jac_NP, + bool do_non_NP_Jac, const Jacobian* Jac_SD_FDdr, const char file_name[]) { FILE *fileptr = fopen(file_name, "w"); @@ -519,54 +537,67 @@ if (fileptr == NULL) "print_Jacobians(): can't open output file \"%s\"!", file_name); /*NOTREACHED*/ -fprintf(fileptr, "# column 1 = x II\n"); +fprintf(fileptr, "# column 1 = x patch name\n"); fprintf(fileptr, "# column 2 = x patch number\n"); fprintf(fileptr, "# column 3 = x irho\n"); fprintf(fileptr, "# column 4 = x isigma\n"); -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_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) +fprintf(fileptr, "# column 5 = x II\n"); +fprintf(fileptr, "# column 6 = y patch name\n"); +fprintf(fileptr, "# column 7 = y patch number\n"); +fprintf(fileptr, "# column 8 = y irho\n"); +fprintf(fileptr, "# column 9 = y isigma\n"); +fprintf(fileptr, "# column 10 = y JJ\n"); +fprintf(fileptr, "# column 11 = Jac_NP(II,JJ)\n"); +if (do_non_NP_Jac) + then { + fprintf(fileptr, "# column 12 = Jac_SD_FDdr(II,JJ)\n"); + fprintf(fileptr, "# column 13 = abs error in Jac_SD_FDdr\n"); + fprintf(fileptr, "# column 14 = rel error in Jac_SD_FDdr\n"); + } + + for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) + { + patch& xp = ps.ith_patch(xpn); + + for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho) + { + for (int x_isigma = xp.min_isigma() ; + x_isigma <= xp.max_isigma() ; + ++x_isigma) + { + const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma); + + for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) { - patch& xp = ps.ith_patch(xpn); + patch& yp = ps.ith_patch(ypn); - for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho) + for (int y_irho = yp.min_irho() ; + y_irho <= yp.max_irho() ; + ++y_irho) { - for (int x_isigma = xp.min_isigma() ; - x_isigma <= xp.max_isigma() ; - ++x_isigma) + for (int y_isigma = yp.min_isigma() ; + y_isigma <= yp.max_isigma() ; + ++y_isigma) { - const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma); + const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma); - for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) - { - patch& yp = ps.ith_patch(ypn); + if (! Jac_NP->is_explicitly_stored(II,JJ)) + then continue; // skip sparse points - for (int y_irho = yp.min_irho() ; - y_irho <= yp.max_irho() ; - ++y_irho) - { - for (int y_isigma = yp.min_isigma() ; - y_isigma <= yp.max_isigma() ; - ++y_isigma) - { - const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma); + const fp NP = (*Jac_NP)(II,JJ); + const fp SD_FDdr = do_non_NP_Jac ? (*Jac_SD_FDdr)(II,JJ) : 0.0; - if (! Jac_SD_FDdr.is_explicitly_stored(II,JJ)) - then continue; // skip sparse points + if ((NP == 0.0) && (SD_FDdr == 0.0)) + then continue; // skip zero values (if == ) - const fp NP = Jac_NP (II,JJ); - const fp SD_FDdr = Jac_SD_FDdr(II,JJ); - - if ((NP == 0.0) && (SD_FDdr == 0.0)) - then continue; // skip zero values (if == ) + fprintf(fileptr, + "%s %d %d %d %d\t%s %d %d %d %d\t%.10g", + xp.name(), xpn, x_irho, x_isigma, II, + yp.name(), ypn, y_irho, y_isigma, JJ, + double(NP)); + if (do_non_NP_Jac) + then { 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); @@ -574,17 +605,18 @@ fprintf(fileptr, "# column 12 = rel error in Jac_SD_FDdr\n"); 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(NP), double(SD_FDdr), + "\t%.10g\t%e\t%e", + double(SD_FDdr), double(SD_FDdr_abs_error), double(SD_FDdr_rel_error)); } - } - } + + fprintf(fileptr, "\n"); } } } + } + } + } fclose(fileptr); } |