diff options
Diffstat (limited to 'src/gr/expansion_Jacobian.cc')
| -rw-r--r-- | src/gr/expansion_Jacobian.cc | 415 |
1 files changed, 399 insertions, 16 deletions
diff --git a/src/gr/expansion_Jacobian.cc b/src/gr/expansion_Jacobian.cc index 372d114..04e00b5 100644 --- a/src/gr/expansion_Jacobian.cc +++ b/src/gr/expansion_Jacobian.cc @@ -17,6 +17,7 @@ #include "util_Table.h" #include "cctk.h" +#include "cctk_Arguments.h" #include "config.h" #include "stdc.h" @@ -24,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" @@ -42,7 +44,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; //****************************************************************************** @@ -53,7 +54,8 @@ using jtutil::error_exit; namespace { enum expansion_status expansion_Jacobian_NP - (patch_system& ps, Jacobian& Jac, fp add_to_expansion, + (patch_system& ps, Jacobian& Jac, + const struct what_to_compute& comput_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, @@ -73,7 +75,8 @@ void add_ghost_zone_Jacobian(const patch_system& ps, int xm_irho, int xm_isigma); enum expansion_status expansion_Jacobian_dr_FD - (patch_system* ps_ptr, Jacobian* Jac_ptr, fp add_to_expansion, + (patch_system* ps_ptr, Jacobian* Jac_ptr, + const struct what_to_compute& compute_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, @@ -109,7 +112,7 @@ enum expansion_status // enum expansion_status expansion_Jacobian(patch_system* ps_ptr, Jacobian* Jac_ptr, - fp add_to_expansion, + const struct what_to_compute& compute_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, @@ -125,7 +128,7 @@ case Jacobian__numerical_perturbation: if (active_flag) then { status = expansion_Jacobian_NP(*ps_ptr, *Jac_ptr, - add_to_expansion, + compute_info, cgi, gi, Jacobian_info, error_info, initial_flag, print_msg_flag); @@ -153,7 +156,8 @@ case Jacobian__symbolic_diff_with_FD_dr: // this function looks at ps_ptr and Jac_ptr (non-NULL vs NULL) // to choose a normal vs dummy computation { - status = expansion_Jacobian_dr_FD(ps_ptr, Jac_ptr, add_to_expansion, + status = expansion_Jacobian_dr_FD(ps_ptr, Jac_ptr, + compute_info, cgi, gi, Jacobian_info, error_info, initial_flag, print_msg_flag); @@ -171,6 +175,327 @@ default: int(Jacobian_info.Jacobian_compute_method)); /*NOTREACHED*/ } + if (active_flag) { + + switch (compute_info.surface_modification) { + + case modification_none: + // do nothing + break; + + case modification_radius: { + // multiply with the coordinate radius + // H_{(r)i} = H_i h_i + // J_{(r)ij} = J_{ij} h_i + H_i \delta_{ij} + const int np = ps_ptr->N_grid_points(); + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->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) { + const int i = ps_ptr->gpn_of_patch_irho_isigma(p, irho,isigma); + const fp radius = p.ghosted_gridfn(gfns::gfn__h, irho, isigma); + for (int j=0; j<np; ++j) { + if (Jac_ptr->is_explicitly_stored (i, j)) { + const fp val = Jac_ptr->element (i, j); + Jac_ptr->set_element (i, j, val * radius); + } + } + const fp Theta = (p.gridfn(gfns::gfn__Theta, irho, isigma) + + compute_info.desired_value) / radius; + Jac_ptr->sum_into_element (i, i, Theta); + } + } + } + break; + } + + case modification_radius2: { + // multiply with the square of the coordinate radius + // H_{(r2)i} = H_i h_i^2 + // J_{(r2)ij} = J_{ij} h_i^2 + 2 H_i h_i \delta_{ij} + const int np = ps_ptr->N_grid_points(); + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->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) { + const int i = ps_ptr->gpn_of_patch_irho_isigma(p, irho,isigma); + const fp radius = p.ghosted_gridfn(gfns::gfn__h, irho, isigma); + const fp radius2 = radius * radius; + for (int j=0; j<np; ++j) { + if (Jac_ptr->is_explicitly_stored (i, j)) { + const fp val = Jac_ptr->element (i, j); + Jac_ptr->set_element (i, j, val * radius2); + } + } + const fp Theta = (p.gridfn(gfns::gfn__Theta, irho, isigma) + + compute_info.desired_value) / radius2; + Jac_ptr->sum_into_element (i, i, 2 * Theta * radius); + } + } + } + break; + } + +#if 0 + case modification_mean_radius: { + // multiply with the average coordinate radius + // H_{(\bar r)i} = H_i \bar r + // J_{(\bar r)ij} = J_{ij} \bar r + H_i / N + // calculate average coordinate radius + const int np = ps_ptr->N_grid_points(); + fp sum_radius = 0; + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->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) { + sum_radius += p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + } + } + } + mean_radius = sum_radius / np; + // correct Jacobian + const int np = ps_ptr->N_grid_points(); + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->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) { + const int i = ps_ptr->gpn_of_patch_irho_isigma(p, irho,isigma); + for (int j=0; j<np; ++j) { + if (Jac_ptr->is_explicitly_stored (i, j)) { + const fp val = Jac_ptr->element (i, j); + Jac_ptr->set_element (i, j, val * mean_radius); + } + } +#error "unfinished" + const fp Theta = (p.gridfn(gfns::gfn__Theta, irho, isigma) + + compute_info.desired_value) / areal_radius; + const fp dRdh = 0.5 * areal_radius; + Jac_ptr->sum_into_element (i, i, Theta * dRdh); + } + } + } + break; + } + + case modification_areal_radius: { + // multiply with the areal radius + // H_{(R)i} = H_i R + // J_{(R)ij} = J_{ij} R + H_i dR/dh_j + // get surface area + const fp area = ps_ptr->integrate_gridfn + (gfns::gfn__one, true, true, true, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + patch::integration_method__automatic_choice); + const fp areal_radius = sqrt(area / (4.0*PI)); + // correct Jacobian + const int np = ps_ptr->N_grid_points(); + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->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) { + const int i = ps_ptr->gpn_of_patch_irho_isigma(p, irho,isigma); + for (int j=0; j<np; ++j) { + if (Jac_ptr->is_explicitly_stored (i, j)) { + const fp val = Jac_ptr->element (i, j); + Jac_ptr->set_element (i, j, val * areal_radius); + } + } + const fp Theta = (p.gridfn(gfns::gfn__Theta, irho, isigma) + + compute_info.desired_value) / areal_radius; + const fp dRdh = 0.5 * areal_radius; + Jac_ptr->sum_into_element (i, i, Theta * dRdh); + } + } + } + break; + } +#endif + + default: + assert (0); + } // switch surface_modification + + if (ps_ptr->N_additional_points()) { + switch (compute_info.surface_selection) { + + case selection_definition: { + // we want nothing special + const int np = ps_ptr->N_grid_points(); + for (int i=0; i<np; ++i) { + Jac_ptr->set_element (i, np, 0.0); + } + for (int j=0; j<np; ++j) { + Jac_ptr->set_element (np, j, 0.0); + } + Jac_ptr->set_element (np, np, 1.0); + break; + } + + case selection_mean_coordinate_radius: { + // Jac_ptr->set_element (II, JJ, x) == dTheta(II)/dh(JJ) + // \frac{\partial R}{\partial h_j} = 1 / N + const int np = ps_ptr->N_grid_points(); + for (int i=0; i<np; ++i) { + Jac_ptr->set_element (i, np, -1.0); + } + for (int j=0; j<np; ++j) { + fp val = 0; + for (int k=0; k<np; ++k) { + val += Jac_ptr->element (k, j) / np; + } + val -= 1.0 / np; + Jac_ptr->set_element (np, j, val); + } + Jac_ptr->set_element (np, np, -1.0); + break; + } + + case selection_areal_radius: { + // \frac{\partial R_a}{\partial h_j} + // = \sqrt{1 / 16 \pi A} \sum_k \sqrt{q_k} dS_k + // The "trapezoid" method is faster +// const enum patch::integration_method method +// = patch::integration_method__automatic_choice; + const enum patch::integration_method method + = patch::integration_method__trapezoid; + const fp area = ps_ptr->integrate_gridfn + (gfns::gfn__one, true, true, true, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + method); + const int np = ps_ptr->N_grid_points(); + for (int i=0; i<np; ++i) { + Jac_ptr->set_element (i, np, -1.0); + } + for (int j=0; j<np; ++j) { + fp val = 0; + for (int k=0; k<np; ++k) { + val += Jac_ptr->element (k, j) / np; + } + Jac_ptr->set_element (np, j, val); + } + for (int jpn = 0; jpn < ps_ptr->N_patches(); ++jpn) { + patch& jp = ps_ptr->ith_patch(jpn); + for (int jrho = jp.min_irho(); jrho <= jp.max_irho(); ++jrho) { + for (int jsigma = jp.min_isigma(); jsigma <= jp.max_isigma(); ++jsigma) { + const int j = ps_ptr->gpn_of_patch_irho_isigma(jp, jrho,jsigma); + // const fp radius = jp.ghosted_gridfn(gfns::gfn__h, jrho, jsigma); + const fp epsilon = Jacobian_info.perturbation_amplitude; + fp val1, val2; +#if 0 + // Re-calculate all points + // (this is slow, but it works) + val1 = area; + jp.ghosted_gridfn(gfns::gfn__h, jrho, jsigma) += epsilon; + val2 = ps_ptr->integrate_gridfn + (gfns::gfn__one, true, true, true, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + method); + jp.ghosted_gridfn(gfns::gfn__h, jrho, jsigma) -= epsilon; +#else + // Re-calculate all points with non-zero Jacobian entries + jp.ghosted_gridfn(gfns::gfn__h, jrho, jsigma) -= epsilon/2; + val1 = 0; + for (int ipn = 0; ipn < ps_ptr->N_patches(); ++ipn) { + patch& ip = ps_ptr->ith_patch(ipn); + for (int irho = ip.min_irho(); irho <= ip.max_irho(); ++irho) { + for (int isigma = ip.min_isigma(); isigma <= ip.max_isigma(); ++isigma) { + const int i = ps_ptr->gpn_of_patch_irho_isigma(ip, irho,isigma); + if (Jac_ptr->is_explicitly_stored (i, j)) { + val1 += ps_ptr->integrate_gridpoint + (gfns::gfn__one, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + method, + ipn, irho, isigma); + } + } + } + } + jp.ghosted_gridfn(gfns::gfn__h, jrho, jsigma) += epsilon; + val2 = 0; + for (int ipn = 0; ipn < ps_ptr->N_patches(); ++ipn) { + patch& ip = ps_ptr->ith_patch(ipn); + for (int irho = ip.min_irho(); irho <= ip.max_irho(); ++irho) { + for (int isigma = ip.min_isigma(); isigma <= ip.max_isigma(); ++isigma) { + const int i = ps_ptr->gpn_of_patch_irho_isigma(ip, irho,isigma); + if (Jac_ptr->is_explicitly_stored (i, j)) { + val2 += ps_ptr->integrate_gridpoint + (gfns::gfn__one, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + method, + ipn, irho, isigma); + } + } + } + } + jp.ghosted_gridfn(gfns::gfn__h, jrho, jsigma) -= epsilon/2; +#endif + const fp val = 1 / sqrt(16*PI*area) * ps_ptr->integrate_correction(true, true, true) * (val2 - val1) / epsilon; + Jac_ptr->sum_into_element (np, j, -val); + } + } + } + Jac_ptr->set_element (np, np, -1.0); + break; + } + + case selection_expansion_mean_coordinate_radius: { + // Jac_ptr->set_element (II, JJ, x) == dTheta(II)/dh(JJ) + const int np = ps_ptr->N_grid_points(); + fp sum_expansion = 0; + fp sum_radius = 0; + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->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) { + sum_expansion += p.gridfn(gfns::gfn__Theta, irho,isigma); + sum_radius += p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + } + } + } + for (int i=0; i<np; ++i) { + Jac_ptr->set_element (i, np, -1.0); + } + for (int j=0; j<np; ++j) { + fp val = 0; + for (int k=0; k<np; ++k) { + val += (Jac_ptr->element (k, j) / np) * (1.0 - sum_radius / np); + } + val -= (sum_expansion / np) / np; + Jac_ptr->set_element (np, j, val); + } + Jac_ptr->set_element (np, np, -1.0); + break; + } + + case selection_expansion_areal_radius: { + CCTK_WARN (0, "selection_expansion_areal_radius not implemented"); + break; + } + + default: + assert (0); + } // switch surface_selection + } else { + assert (compute_info.surface_selection == selection_definition); + } + + } // if active + return expansion_success; // *** NORMAL RETURN *** } @@ -214,7 +539,8 @@ return expansion_success; // *** NORMAL RETURN *** namespace { enum expansion_status expansion_Jacobian_NP - (patch_system& ps, Jacobian& Jac, fp add_to_expansion, + (patch_system& ps, Jacobian& Jac, + const struct what_to_compute& compute_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, @@ -227,6 +553,7 @@ if (print_msg_flag) const fp epsilon = Jacobian_info.perturbation_amplitude; ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta); +ps.gridfn_copy(gfns::gfn__mean_curvature, gfns::gfn__save_mean_curvature); for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) { @@ -247,7 +574,8 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta); const fp save_h_y = yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma); yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) += epsilon; const - enum expansion_status status = expansion(&ps, add_to_expansion, + enum expansion_status status = expansion(&ps, + compute_info, cgi, gi, error_info, initial_flag); if (status != expansion_success) @@ -275,12 +603,19 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta); } } + if (ps.N_additional_points()) + { + const int np = ps.N_grid_points(); + Jac.set_element(np,JJ, 0.0); // insert dummy value + } + yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) = save_h_y; } } } ps.gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta); +ps.gridfn_copy(gfns::gfn__save_mean_curvature, gfns::gfn__mean_curvature); return expansion_success; // *** NORMAL RETURN *** } } @@ -440,6 +775,12 @@ ps.compute_synchronize_Jacobian(); } } + if (ps.N_additional_points()) + { + const int np = ps.N_grid_points(); + Jac.set_element(II,np, 0.0); // insert dummy value + } + } } } @@ -519,11 +860,12 @@ patch& yp = ye.my_patch(); // It's illegal for one but not both of ps_ptr and Jac_ptr to be NULL. // // The basic algorithm is that -// Jac += diag[ (Theta(h+epsilon) - Theta(h)) / epsilon ] +// Jac += diag[ (Theta(h+epsilon/2) - Theta(h-epsilon/2)) / epsilon ] // // Inputs (angular gridfns, on ghosted grid): // h # shape of trial surface // Theta # Theta(h) assumed to already be computed +// # (saved and restored, but not used) // // Outputs: // Jac += d/dr terms @@ -535,7 +877,8 @@ patch& yp = ye.my_patch(); namespace { enum expansion_status expansion_Jacobian_dr_FD - (patch_system* ps_ptr, Jacobian* Jac_ptr, fp add_to_expansion, + (patch_system* ps_ptr, Jacobian* Jac_ptr, + const struct what_to_compute& compute_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, @@ -550,18 +893,52 @@ if (print_msg_flag) const fp epsilon = Jacobian_info.perturbation_amplitude; -// compute Theta(h+epsilon) +what_to_compute this_compute_info (compute_info); +this_compute_info.surface_modification = modification_none; +this_compute_info.surface_selection = selection_definition; +this_compute_info.desired_value = 0.0; + +fp additional_save_Theta; + +// compute Theta(h-epsilon/2) if (active_flag) then { ps_ptr->gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta); - ps_ptr->add_to_ghosted_gridfn(epsilon, gfns::gfn__h); + ps_ptr->gridfn_copy(gfns::gfn__mean_curvature, gfns::gfn__save_mean_curvature); + if (ps_ptr->N_additional_points()) + then { + const int np = ps_ptr->N_grid_points(); + additional_save_Theta = ps_ptr->gridfn_data(gfns::gfn__Theta)[np]; + } + ps_ptr->add_to_ghosted_gridfn(-epsilon/2, gfns::gfn__h); } const - enum expansion_status status = expansion(ps_ptr, add_to_expansion, + enum expansion_status status = expansion(ps_ptr, + this_compute_info, cgi, gi, error_info, initial_flag); if (status != expansion_success) - then return status; // *** ERROR RETURN *** + then { + expansion(NULL, + this_compute_info, + cgi, gi, + error_info, false); + return status; // *** ERROR RETURN *** + } + +// compute Theta(h+epsilon/2) +if (active_flag) + then { + ps_ptr->gridfn_copy(gfns::gfn__Theta, gfns::gfn__old_Theta); + ps_ptr->add_to_ghosted_gridfn(epsilon, gfns::gfn__h); + } +const + enum expansion_status status2 = expansion(ps_ptr, + this_compute_info, + cgi, gi, + error_info, initial_flag); +if (status2 != expansion_success) + then return status2; // *** ERROR RETURN *** if (active_flag) then { @@ -575,7 +952,7 @@ if (active_flag) ++isigma) { const int II = ps_ptr->gpn_of_patch_irho_isigma(p, irho,isigma); - const fp old_Theta = p.gridfn(gfns::gfn__save_Theta, + const fp old_Theta = p.gridfn(gfns::gfn__old_Theta, irho,isigma); const fp new_Theta = p.gridfn(gfns::gfn__Theta, irho,isigma); @@ -586,8 +963,14 @@ if (active_flag) } // restore h and Theta - ps_ptr->add_to_ghosted_gridfn(-epsilon, gfns::gfn__h); + ps_ptr->add_to_ghosted_gridfn(-epsilon/2, gfns::gfn__h); ps_ptr->gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta); + ps_ptr->gridfn_copy(gfns::gfn__save_mean_curvature, gfns::gfn__mean_curvature); + if (ps_ptr->N_additional_points()) + then { + const int np = ps_ptr->N_grid_points(); + ps_ptr->gridfn_data(gfns::gfn__Theta)[np] = additional_save_Theta; + } } return expansion_success; // *** NORMAL RETURN *** |