aboutsummaryrefslogtreecommitdiff
path: root/src/driver/Newton.cc
diff options
context:
space:
mode:
Diffstat (limited to 'src/driver/Newton.cc')
-rw-r--r--src/driver/Newton.cc694
1 files changed, 642 insertions, 52 deletions
diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc
index 2a54d17..62bab31 100644
--- a/src/driver/Newton.cc
+++ b/src/driver/Newton.cc
@@ -10,15 +10,19 @@
/// Newton_step - take the Newton step, scaling it down if it's too large
//
+#include <iostream>
+
#include <stdio.h>
#include <assert.h>
#include <limits.h>
#include <float.h>
#include <math.h>
+#include <string.h>
#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
#include "config.h"
#include "stdc.h"
@@ -26,6 +30,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"
@@ -48,7 +53,6 @@
// all the code in this file is inside this namespace
namespace AHFinderDirect
{
-using jtutil::error_exit;
//******************************************************************************
@@ -67,6 +71,7 @@ bool broadcast_status(const cGH* GH,
struct iteration_status_buffers& isb);
void broadcast_horizon_data(const cGH* GH,
bool broadcast_flag, bool broadcast_horizon_shape,
+ struct AH_data& AH_data,
struct BH_diagnostics& BH_diagnostics,
patch_system& ps,
struct horizon_buffers& horizon_buffers);
@@ -106,22 +111,22 @@ void Newton(const cGH* GH,
const struct verbose_info& verbose_info,
struct iteration_status_buffers& isb)
{
+cGH const * const cctkGH = GH;
+DECLARE_CCTK_ARGUMENTS;
+DECLARE_CCTK_PARAMETERS;
+
const bool my_active_flag = hs.has_genuine_horizons();
-const int N_horizons = hs.N_horizons();
// print out which horizons we're finding on this processor
-if (verbose_info.print_physics_details)
- then {
- if (hs.has_genuine_horizons())
- then CCTK_VInfo(CCTK_THORNSTRING,
- "proc %d: searching for horizon%s %s/%d",
- my_proc,
- (hs.my_N_horizons() > 1 ? "s" : ""),
- hs.sequence_string(","), N_horizons);
- else CCTK_VInfo(CCTK_THORNSTRING,
- "proc %d: dummy horizon(s) only",
- my_proc);
- }
+if (hs.has_genuine_horizons())
+ then CCTK_VInfo(CCTK_THORNSTRING,
+ "proc %d: searching for horizon%s %s/%d",
+ my_proc,
+ (hs.my_N_horizons() > 1 ? "s" : ""),
+ hs.sequence_string(","), N_horizons);
+ else CCTK_VInfo(CCTK_THORNSTRING,
+ "proc %d: dummy horizon(s) only",
+ my_proc);
//
// each pass through this loop finds a single horizon,
@@ -138,7 +143,10 @@ if (verbose_info.print_physics_details)
"Newton_solve(): processor %d working on horizon %d",
my_proc, hn);
- const bool horizon_is_genuine = hs.is_genuine();
+ // only try to find horizons every find_every time steps
+ const bool horizon_is_genuine =
+ hs.is_genuine() && AH_data_array[hn]->search_flag;
+ // this is only a pessimistic approximation
const bool there_is_another_genuine_horizon = hs.is_next_genuine();
if (verbose_info.print_algorithm_details)
then {
@@ -152,12 +160,262 @@ if (verbose_info.print_physics_details)
struct AH_data* AH_data_ptr
= horizon_is_genuine ? AH_data_array[hn] : NULL;
+
patch_system* const ps_ptr = horizon_is_genuine
? AH_data_ptr->ps_ptr : NULL;
Jacobian* const Jac_ptr = horizon_is_genuine
? AH_data_ptr->Jac_ptr: NULL;
- const fp add_to_expansion = horizon_is_genuine
- ? -AH_data_ptr->surface_expansion : 0.0;
+
+ if (horizon_is_genuine) {
+ // deal with dependent horizons
+ if (AH_data_ptr->depends_on > 0) {
+ assert (AH_data_ptr->depends_on != hn);
+ assert (AH_data_ptr->depends_on < hn);
+ // check that the other horizon is on the same processor!
+ AH_data *AH_other_ptr = AH_data_array[AH_data_ptr->depends_on];
+ assert (AH_other_ptr);
+ AH_data_ptr->compute_info.desired_value
+ = AH_other_ptr->compute_info.desired_value
+ * AH_data_ptr->desired_value_factor
+ + AH_data_ptr->desired_value_offset;
+ const int gnp = ps_ptr->ghosted_N_grid_points();
+ assert (AH_other_ptr->ps_ptr->ghosted_N_grid_points() == gnp);
+ memcpy (ps_ptr->ghosted_gridfn_data(gfns::gfn__h),
+ AH_other_ptr->ps_ptr->ghosted_gridfn_data(gfns::gfn__h),
+ gnp * sizeof(fp));
+ ps_ptr->origin_x (AH_other_ptr->ps_ptr->origin_x());
+ ps_ptr->origin_y (AH_other_ptr->ps_ptr->origin_y());
+ ps_ptr->origin_z (AH_other_ptr->ps_ptr->origin_z());
+ }
+ }
+
+ if (horizon_is_genuine) {
+ if (AH_data_ptr->move_origins
+ && AH_data_ptr->depends_on == 0
+ && AH_data_ptr->found_flag)
+ {
+ // move the origin to the centre
+ patch_system& ps = *ps_ptr;
+ fp cx=0, cy=0, cz=0; // change of origin
+ fp sx=0, sy=0, sz=0; // change of shape
+ if (! predict_origin_movement) {
+ size_t N=0;
+ for (int pn=0; pn<ps.N_patches(); ++pn) {
+ patch& p = ps.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 fp rho = p.rho_of_irho (irho);
+ const fp sigma = p.sigma_of_isigma (isigma);
+ const fp r = p.ghosted_gridfn (gfns::gfn__h, irho, isigma);
+ fp x, y, z;
+ p.xyz_of_r_rho_sigma (r, rho, sigma, x, y, z);
+ cx += x;
+ cy += y;
+ cz += z;
+ ++N;
+ }
+ }
+ }
+ assert (N > 0);
+ cx /= N;
+ cy /= N;
+ cz /= N;
+ sx=cx; sy=cy; sz=cz;
+ } else { // if predict_origin_movement
+ if (ah_centroid_valid[hn-1]) {
+ if (verbose_info.print_algorithm_highlights) {
+ std::cout << "AHF centroid x " << (ah_centroid_x[hn-1]) << " y " << (ah_centroid_y[hn-1]) << " z " << (ah_centroid_z[hn-1]) << " t " << (ah_centroid_t[hn-1]) << std::endl;
+ }
+ if (ah_centroid_valid_p[hn-1]) {
+ // have two previous origins: linear extrapolation
+ if (verbose_info.print_algorithm_highlights) {
+ std::cout << "AHF centroid_p x " << (ah_centroid_x_p[hn-1]) << " y " << (ah_centroid_y_p[hn-1]) << " z " << (ah_centroid_z_p[hn-1]) << " t " << (ah_centroid_t_p[hn-1]) << std::endl;
+ }
+ fp const dt = ah_centroid_t[hn-1] - ah_centroid_t_p[hn-1];
+ fp const dt_n = cctk_time - ah_centroid_t [hn-1];
+ fp const dt_p = cctk_time - ah_centroid_t_p[hn-1];
+ fp const timescale =
+ fabs (dt) + fabs (dt_p) + fabs (dt_n) +
+ 0.001 * fabs (cctk_delta_time);
+ if (10 * fabs (dt ) < timescale ||
+ 10 * fabs (dt_p) < timescale ||
+ 10 * fabs (dt_n) < timescale)
+ {
+ // the times are too similar
+ if (verbose_info.print_algorithm_highlights) {
+ std::cout << "AHF toosim" << std::endl;
+ }
+ cx = ah_centroid_x[hn-1];
+ cy = ah_centroid_y[hn-1];
+ cz = ah_centroid_z[hn-1];
+ } else {
+ fp const fact_n = + dt_p / dt;
+ fp const fact_p = - dt_n / dt;
+ if (fabs (fact_n) > 5 || fabs (fact_p) > 5) {
+ // don't trust a large extrapolation
+ if (verbose_info.print_algorithm_highlights) {
+ std::cout << "AHF notrust" << std::endl;
+ }
+ cx = ah_centroid_x[hn-1];
+ cy = ah_centroid_y[hn-1];
+ cz = ah_centroid_z[hn-1];
+ } else {
+ if (verbose_info.print_algorithm_highlights) {
+ std::cout << "AHF extrap fn " << fact_n << " fp " << fact_p << std::endl;
+ }
+ cx = fact_n * ah_centroid_x[hn-1] + fact_p * ah_centroid_x_p[hn-1];
+ cy = fact_n * ah_centroid_y[hn-1] + fact_p * ah_centroid_y_p[hn-1];
+ cz = fact_n * ah_centroid_z[hn-1] + fact_p * ah_centroid_z_p[hn-1];
+ if (verbose_info.print_algorithm_highlights) {
+ std::cout << "AHF xp " << (ah_centroid_x_p[hn-1]) << " x " << (ah_centroid_x[hn-1]) << " cx " << cx << std::endl;
+ }
+ }
+ }
+ } else {
+ // have one previous origin: constant extrapolation
+ if (verbose_info.print_algorithm_highlights) {
+ std::cout << "AHF const" << std::endl;
+ }
+ cx = ah_centroid_x[hn-1];
+ cy = ah_centroid_y[hn-1];
+ cz = ah_centroid_z[hn-1];
+ }
+ if (verbose_info.print_algorithm_highlights) {
+ std::cout << "AHF predicted position cx " << cx << " cy " << cy << " cz " << cz << std::endl;
+ }
+ cx -= ps.origin_x();
+ cy -= ps.origin_y();
+ cz -= ps.origin_z();
+ }
+ sx = ah_centroid_x[hn-1] - ps.origin_x();
+ sy = ah_centroid_y[hn-1] - ps.origin_y();
+ sz = ah_centroid_z[hn-1] - ps.origin_z();
+ } // if predict_origin_movement
+ switch (ps.type()) {
+ case patch_system::patch_system__full_sphere:
+ break; // do nothing
+ case patch_system::patch_system__plus_z_hemisphere:
+ cz = 0; sz = 0; break;
+ case patch_system::patch_system__plus_xy_quadrant_mirrored:
+ cx = cy = 0; sx = sy = 0; break;
+ case patch_system::patch_system__plus_xy_quadrant_rotating:
+ cx = cy = 0; sx = sy = 0; break;
+ case patch_system::patch_system__plus_xz_quadrant_mirrored:
+ cx = cz = 0; sx = sz = 0; break;
+ case patch_system::patch_system__plus_xz_quadrant_rotating:
+ cx = cz = 0; sx = sz = 0; break;
+ case patch_system::patch_system__plus_xyz_octant_mirrored:
+ cx = cy = cz = 0; sx = sy = sz = 0; break;
+ case patch_system::patch_system__plus_xyz_octant_rotating:
+ cx = cy = cz = 0; sx = sy = sz = 0; break;
+ default: assert(0);
+ }
+ ps.origin_x (ps.origin_x() + cx);
+ ps.origin_y (ps.origin_y() + cy);
+ ps.origin_z (ps.origin_z() + cz);
+ if (reshape_while_moving) {
+ for (int pn=0; pn<ps.N_patches(); ++pn) {
+ patch& p = ps.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 fp rho = p.rho_of_irho (irho);
+ const fp sigma = p.sigma_of_isigma (isigma);
+ fp & r = p.ghosted_gridfn (gfns::gfn__h, irho, isigma);
+ fp x, y, z;
+ p.xyz_of_r_rho_sigma (r, rho, sigma, x, y, z);
+ fp const proj = (sx*x + sy*y + sz*z) / r;
+ r -= proj;
+ }
+ }
+ }
+ ps.synchronize();
+ }
+ }
+
+ // modify the initial guess
+ jtutil::norm<fp> norms;
+ ps_ptr->ghosted_gridfn_norms (gfns::gfn__h, norms);
+ // smoothing:
+ ps_ptr->scale_ghosted_gridfn
+ (1.0 - AH_data_ptr->smoothing_factor, gfns::gfn__h);
+ ps_ptr->add_to_ghosted_gridfn
+ (AH_data_ptr->smoothing_factor * norms.mean(), gfns::gfn__h);
+ // enlarging:
+ ps_ptr->scale_ghosted_gridfn
+ (AH_data_ptr->shiftout_factor, gfns::gfn__h);
+ }
+
+ bool I_am_pretracking = horizon_is_genuine && AH_data_ptr->use_pretracking;
+ bool I_was_pretracking = false;
+ bool pretracking_have_upper_bound = false;
+ bool pretracking_have_lower_bound = false;
+ bool pretracking_was_successful = false;
+ fp const old_pretracking_value = I_am_pretracking ? AH_data_ptr->pretracking_value : 0.0;
+ fp pretracking_upper_bound;
+ fp pretracking_lower_bound;
+ bool pretracking_have_horizon_info;
+ fp pretracking_mean_expansion;
+ for (int pretracking_iter = 0;
+ I_am_pretracking
+ ? pretracking_iter < AH_data_ptr->pretracking_max_iterations
+ : true;
+ ++pretracking_iter)
+ {
+ bool found_this_horizon;
+ if (I_am_pretracking || I_was_pretracking) {
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: iteration %d", pretracking_iter);
+ if (pretracking_have_lower_bound) {
+ if (pretracking_have_upper_bound) {
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: looking for value %g in [%g,%g]",
+ double(AH_data_ptr->pretracking_value),
+ double(pretracking_lower_bound),
+ double(pretracking_upper_bound));
+ } else {
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: looking for value %g in [%g,*]",
+ double(AH_data_ptr->pretracking_value),
+ double(pretracking_lower_bound));
+ }
+ } else {
+ if (pretracking_have_upper_bound) {
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: looking for value %g in [*,%g]",
+ double(AH_data_ptr->pretracking_value),
+ double(pretracking_upper_bound));
+ } else {
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: looking for value %g in [*,*]",
+ double(AH_data_ptr->pretracking_value));
+ }
+ }
+ AH_data_ptr->compute_info.desired_value = AH_data_ptr->pretracking_value;
+ ps_ptr->ghosted_gridfn_copy (gfns::gfn__h, gfns::gfn__save_h);
+ pretracking_have_horizon_info = false;
+ }
+
+ if (! I_am_pretracking) {
+ if (horizon_is_genuine) {
+ if (! AH_data_ptr->initial_guess_info.reset_horizon_after_not_finding) {
+ // save the surface for possible backtracking later
+ ps_ptr->ghosted_gridfn_copy (gfns::gfn__h, gfns::gfn__save_h);
+ }
+ }
+ }
+
+ struct what_to_compute dummy_compute_info;
+ struct what_to_compute & compute_info =
+ horizon_is_genuine
+ ? AH_data_ptr->compute_info
+ : dummy_compute_info;
+
+ if (horizon_is_genuine) {
+ if (ps_ptr->N_additional_points()) {
+ const int gnp = ps_ptr->ghosted_N_grid_points();
+ ps_ptr->ghosted_gridfn_data(gfns::gfn__h)[gnp] = 0;
+ }
+ }
const int max_iterations
= horizon_is_genuine
@@ -166,10 +424,16 @@ if (verbose_info.print_physics_details)
: solver_info.max_Newton_iterations__subsequent)
: INT_MAX;
+ int num_Theta_growth_iterations = 0;
+ fp previous_Theta_norm = 1.0e30;
+ int num_Theta_nonshrink_iterations = 0;
+ fp best_Theta_norm = 1.0e30;
+
//
// each pass through this loop does a single Newton iteration
// on the current horizon (which might be either genuine or dummy)
//
+ bool do_return = false;
for (int iteration = 1 ; ; ++iteration)
{
if (verbose_info.print_algorithm_debug)
@@ -185,20 +449,136 @@ if (verbose_info.print_physics_details)
if (horizon_is_genuine
&& solver_info.debugging_output_at_each_Newton_iteration)
then output_gridfn(*ps_ptr, gfns::gfn__h,
+ "h", GH,
IO_info, IO_info.h_base_file_name,
+ IO_info.h_min_digits,
hn, verbose_info.print_algorithm_highlights,
iteration);
+ // calculate the norms also for a surface a bit further out,
+ // so that we know how they vary in space.
+ // perform this calculation first, so that the real values
+ // do not have to be saved.
+ const fp epsilon = Jacobian_info.perturbation_amplitude;
+ jtutil::norm<fp> shifted_Theta_norms;
+ jtutil::norm<fp> shifted_expansion_Theta_norms;
+ jtutil::norm<fp> shifted_inner_expansion_Theta_norms;
+ jtutil::norm<fp> shifted_product_expansion_Theta_norms;
+ jtutil::norm<fp> shifted_mean_curvature_Theta_norms;
+ fp shifted_area;
+ jtutil::norm<fp> shifted2_Theta_norms;
+ jtutil::norm<fp> shifted2_expansion_Theta_norms;
+ jtutil::norm<fp> shifted2_inner_expansion_Theta_norms;
+ jtutil::norm<fp> shifted2_product_expansion_Theta_norms;
+ jtutil::norm<fp> shifted2_mean_curvature_Theta_norms;
+ fp shifted2_area;
+ if (solver_info.want_expansion_gradients) {
+ if (horizon_is_genuine) {
+ ps_ptr->add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
+ // ps_ptr->scale_ghosted_gridfn(1.0+epsilon, gfns::gfn__h);
+ }
+ const enum expansion_status raw_shifted_expansion_status
+ = expansion(ps_ptr, compute_info,
+ cgi, gi,
+ error_info, (iteration == 1),
+ false, // no, we don't want Jacobian coeffs
+ false,
+ &shifted_Theta_norms,
+ &shifted_expansion_Theta_norms,
+ &shifted_inner_expansion_Theta_norms,
+ &shifted_product_expansion_Theta_norms,
+ &shifted_mean_curvature_Theta_norms);
+ if (horizon_is_genuine) {
+ shifted_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,
+ BH_diagnostics_info.integral_method);
+ ps_ptr->add_to_ghosted_gridfn(-epsilon, gfns::gfn__h);
+ // ps_ptr->scale_ghosted_gridfn(1.0/(1.0+epsilon), gfns::gfn__h);
+
+ ps_ptr->add_to_ghosted_gridfn(-epsilon, gfns::gfn__h);
+ // ps_ptr->scale_ghosted_gridfn(1.0/(1.0+epsilon), gfns::gfn__h);
+ }
+ const enum expansion_status raw_shifted2_expansion_status
+ = expansion(ps_ptr, compute_info,
+ cgi, gi,
+ error_info, (iteration == 1),
+ false, // no, we don't want Jacobian coeffs
+ false,
+ &shifted2_Theta_norms,
+ &shifted2_expansion_Theta_norms,
+ &shifted2_inner_expansion_Theta_norms,
+ &shifted2_product_expansion_Theta_norms,
+ &shifted2_mean_curvature_Theta_norms);
+ if (horizon_is_genuine) {
+ shifted2_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,
+ BH_diagnostics_info.integral_method);
+ ps_ptr->add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
+ // ps_ptr->scale_ghosted_gridfn(1.0+epsilon, gfns::gfn__h);
+ }
+ } // if want_expansion_gradients
+
+ // now calculate the real values.
jtutil::norm<fp> Theta_norms;
+ jtutil::norm<fp> expansion_Theta_norms;
+ jtutil::norm<fp> inner_expansion_Theta_norms;
+ jtutil::norm<fp> product_expansion_Theta_norms;
+ jtutil::norm<fp> mean_curvature_Theta_norms;
const enum expansion_status raw_expansion_status
- = expansion(ps_ptr, add_to_expansion,
+ = expansion(ps_ptr, compute_info,
cgi, gi,
error_info, (iteration == 1),
true, // yes, we want Jacobian coeffs
verbose_info.print_algorithm_details,
- &Theta_norms);
+ &Theta_norms,
+ &expansion_Theta_norms,
+ &inner_expansion_Theta_norms,
+ &product_expansion_Theta_norms,
+ &mean_curvature_Theta_norms);
+ fp area;
+ if (horizon_is_genuine)
+ {
+ 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,
+ BH_diagnostics_info.integral_method);
+ }
const bool Theta_is_ok = (raw_expansion_status == expansion_success);
- const bool norms_are_ok = horizon_is_genuine && Theta_is_ok;
+ const bool norms_are_ok = horizon_is_genuine && Theta_is_ok;
+
+ if (norms_are_ok) {
+ const fp this_norm = Theta_norms.infinity_norm();
+ if (this_norm > previous_Theta_norm) {
+ ++ num_Theta_growth_iterations;
+ } else {
+ num_Theta_growth_iterations = 0;
+ }
+ previous_Theta_norm = this_norm;
+
+ if (this_norm >= best_Theta_norm) {
+ ++ num_Theta_nonshrink_iterations;
+ } else {
+ num_Theta_nonshrink_iterations = 0;
+ best_Theta_norm = this_norm;
+ }
+ }
+
+ if (I_am_pretracking && norms_are_ok) {
+ pretracking_mean_expansion = expansion_Theta_norms.mean();
+ pretracking_have_horizon_info = true;
+ }
+
if (verbose_info.print_algorithm_debug)
then {
CCTK_VInfo(CCTK_THORNSTRING,
@@ -212,27 +592,64 @@ if (verbose_info.print_physics_details)
}
if (horizon_is_genuine && Theta_is_ok
&& solver_info.debugging_output_at_each_Newton_iteration)
- then output_gridfn(*ps_ptr, gfns::gfn__Theta,
+ then {
+ output_gridfn(*ps_ptr, gfns::gfn__Theta,
+ "Theta", GH,
IO_info, IO_info.Theta_base_file_name,
+ IO_info.h_min_digits,
+ hn, verbose_info.print_algorithm_highlights,
+ iteration);
+ output_gridfn(*ps_ptr, gfns::gfn__mean_curvature,
+ "mean_curvature", GH,
+ IO_info, IO_info.mean_curvature_base_file_name,
+ IO_info.h_min_digits,
hn, verbose_info.print_algorithm_highlights,
iteration);
+ }
//
// have we found this horizon?
// if so, compute and output BH diagnostics
//
- const bool found_this_horizon
- = norms_are_ok && (Theta_norms.infinity_norm()
- <= solver_info.Theta_norm_for_convergence);
+ found_this_horizon
+ = (norms_are_ok
+ && (I_was_pretracking
+ ? pretracking_was_successful
+ : Theta_norms.infinity_norm() <= solver_info.Theta_norm_for_convergence));
if (horizon_is_genuine)
then AH_data_ptr->found_flag = found_this_horizon;
- if (found_this_horizon)
+ if (found_this_horizon && ! I_am_pretracking)
then {
struct BH_diagnostics& BH_diagnostics
= AH_data_ptr->BH_diagnostics;
- BH_diagnostics.compute(*ps_ptr, BH_diagnostics_info);
+ const fp mean_expansion = expansion_Theta_norms.mean();
+ const fp mean_inner_expansion = inner_expansion_Theta_norms.mean();
+ const fp mean_product_expansion = product_expansion_Theta_norms.mean();
+ const fp mean_mean_curvature = mean_curvature_Theta_norms.mean();
+ // const fp area_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_area - area ) / epsilon;
+ // const fp mean_expansion_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_expansion_Theta_norms.mean() - expansion_Theta_norms.mean() ) / epsilon;
+ // const fp mean_inner_expansion_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_inner_expansion_Theta_norms.mean() - inner_expansion_Theta_norms.mean()) / epsilon;
+ // const fp mean_product_expansion_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_product_expansion_Theta_norms.mean() - product_expansion_Theta_norms.mean()) / epsilon;
+ // const fp mean_mean_curvature_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_mean_curvature_Theta_norms.mean() - mean_curvature_Theta_norms.mean() ) / epsilon;
+ const fp area_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_area - shifted2_area ) / (2*epsilon);
+ const fp mean_expansion_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_expansion_Theta_norms.mean() - shifted2_expansion_Theta_norms.mean() ) / (2*epsilon);
+ const fp mean_inner_expansion_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_inner_expansion_Theta_norms.mean() - shifted2_inner_expansion_Theta_norms.mean()) / (2*epsilon);
+ const fp mean_product_expansion_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_product_expansion_Theta_norms.mean() - shifted2_product_expansion_Theta_norms.mean()) / (2*epsilon);
+ const fp mean_mean_curvature_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_mean_curvature_Theta_norms.mean() - shifted2_mean_curvature_Theta_norms.mean() ) / (2*epsilon);
+ BH_diagnostics.compute(*ps_ptr,
+ area,
+ mean_expansion,
+ mean_inner_expansion,
+ mean_product_expansion,
+ mean_mean_curvature,
+ area_gradient,
+ mean_expansion_gradient,
+ mean_inner_expansion_gradient,
+ mean_product_expansion_gradient,
+ mean_mean_curvature_gradient,
+ BH_diagnostics_info);
if (IO_info.output_BH_diagnostics)
then {
if (AH_data_ptr->BH_diagnostics_fileptr == NULL)
@@ -251,8 +668,13 @@ if (verbose_info.print_physics_details)
// (if so, we'll give up on this horizon)
//
const bool expansion_is_too_large
- = norms_are_ok && (Theta_norms.infinity_norm()
- > solver_info.max_allowable_Theta);
+ = norms_are_ok
+ && ( Theta_norms.infinity_norm() > solver_info.max_allowable_Theta
+ || ( solver_info.max_allowable_Theta_growth_iterations > 0
+ && num_Theta_growth_iterations > solver_info.max_allowable_Theta_growth_iterations)
+ || ( solver_info.max_allowable_Theta_nonshrink_iterations > 0
+ && num_Theta_nonshrink_iterations > solver_info.max_allowable_Theta_nonshrink_iterations)
+ );
//
@@ -260,10 +682,11 @@ if (verbose_info.print_physics_details)
// then pretend expansion() returned a "surface too large" error status
//
jtutil::norm<fp> h_norms;
- if (horizon_is_genuine)
- then ps_ptr->ghosted_gridfn_norms(gfns::gfn__h, h_norms);
- const fp mean_horizon_radius = horizon_is_genuine ? h_norms.mean()
- : 0.0;
+ if (horizon_is_genuine) {
+ ps_ptr->ghosted_gridfn_norms(gfns::gfn__h, h_norms);
+ }
+ const fp mean_horizon_radius
+ = horizon_is_genuine ? h_norms.mean() : 0.0;
const bool horizon_is_too_large
= (mean_horizon_radius > solver_info
.max_allowable_horizon_radius[hn]);
@@ -279,7 +702,7 @@ if (verbose_info.print_physics_details)
// does *this* horizon need more iterations?
// i.e. has this horizon's Newton iteration not yet converged?
- const bool this_horizon_needs_more_iterations
+ const bool this_horizon_needs_more_iterations
= horizon_is_genuine && Theta_is_ok
&& !found_this_horizon
&& !expansion_is_too_large
@@ -290,7 +713,8 @@ if (verbose_info.print_physics_details)
// on this or a following horizon?
const bool I_need_more_iterations
= this_horizon_needs_more_iterations
- || there_is_another_genuine_horizon;
+ || there_is_another_genuine_horizon
+ || I_am_pretracking;
if (verbose_info.print_algorithm_details)
then {
@@ -305,7 +729,6 @@ if (verbose_info.print_physics_details)
int(I_need_more_iterations));
}
-
//
// broadcast iteration status from each active processor
// to all processors, and inclusive-or the "we need more iterations"
@@ -373,6 +796,7 @@ if (verbose_info.print_physics_details)
broadcast_horizon_data(GH,
my_proc == found_proc,
broadcast_horizon_shape,
+ found_AH_data,
found_AH_data.BH_diagnostics,
*found_AH_data.ps_ptr,
found_AH_data.horizon_buffers);
@@ -387,8 +811,9 @@ if (verbose_info.print_physics_details)
//
// if we found our horizon, maybe output the horizon shape?
//
- if (found_this_horizon)
+ if (found_this_horizon && ! I_am_pretracking)
then {
+ // printf("will output h/Th/mc: %d/%d/%d\n", IO_info.output_h, IO_info.output_Theta, IO_info.output_mean_curvature); //xxxxxxxxxxxx
if (IO_info.output_h)
then {
// if this is the first time we've output h for this
@@ -396,13 +821,24 @@ if (verbose_info.print_physics_details)
if (!AH_data_ptr->h_files_written)
then setup_h_files(*ps_ptr, IO_info, hn);
output_gridfn(*ps_ptr, gfns::gfn__h,
+ "h", GH,
IO_info, IO_info.h_base_file_name,
+ IO_info.h_min_digits,
hn, verbose_info
.print_algorithm_highlights);
}
if (IO_info.output_Theta)
then output_gridfn(*ps_ptr, gfns::gfn__Theta,
+ "Theta", GH,
IO_info, IO_info.Theta_base_file_name,
+ IO_info.h_min_digits,
+ hn, verbose_info
+ .print_algorithm_highlights);
+ if (IO_info.output_mean_curvature)
+ then output_gridfn(*ps_ptr, gfns::gfn__mean_curvature,
+ "mean_curvature", GH,
+ IO_info, IO_info.mean_curvature_base_file_name,
+ IO_info.h_min_digits,
hn, verbose_info
.print_algorithm_highlights);
}
@@ -417,7 +853,8 @@ if (verbose_info.print_physics_details)
if (verbose_info.print_algorithm_details)
then CCTK_VInfo(CCTK_THORNSTRING,
"==> all processors are finished!");
- return; // *** NORMAL RETURN ***
+ do_return = true;
+ break; // *** LOOP EXIT ***
}
if ((N_procs == 1) && !this_horizon_needs_more_iterations)
then {
@@ -441,7 +878,7 @@ if (verbose_info.print_physics_details)
= expansion_Jacobian
(this_horizon_needs_more_iterations ? ps_ptr : NULL,
this_horizon_needs_more_iterations ? Jac_ptr : NULL,
- add_to_expansion,
+ compute_info,
cgi, gi, Jacobian_info,
error_info, (iteration == 1),
verbose_info.print_algorithm_details);
@@ -474,7 +911,7 @@ if (verbose_info.print_physics_details)
verbose_info.print_algorithm_details);
if ((rcond >= 0.0) && (rcond < 100.0*FP_EPSILON))
then {
- CCTK_VWarn(SERIOUS_WARNING, __LINE__, __FILE__, CCTK_THORNSTRING,
+ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"Newton_solve: Jacobian matrix is numerically singular!");
// give up on this horizon
break; // *** LOOP CONTROL ***
@@ -491,7 +928,9 @@ if (verbose_info.print_physics_details)
if (solver_info.debugging_output_at_each_Newton_iteration)
then output_gridfn(*ps_ptr, gfns::gfn__Delta_h,
+ "Delta_h", GH,
IO_info, IO_info.Delta_h_base_file_name,
+ IO_info.h_min_digits,
hn, verbose_info.print_algorithm_details,
iteration);
@@ -507,6 +946,139 @@ if (verbose_info.print_physics_details)
// end of this Newton iteration
}
+ if (! I_am_pretracking) {
+ if (horizon_is_genuine) {
+ if (! AH_data_ptr->initial_guess_info.reset_horizon_after_not_finding) {
+ if (! found_this_horizon) {
+ // the surface failed; backtrack and continue
+ AH_data_ptr->ps_ptr->ghosted_gridfn_copy
+ (gfns::gfn__save_h, gfns::gfn__h);
+ }
+ }
+ }
+ }
+
+ // exit
+ if (do_return) return; // *** NORMAL RETURN ***
+
+ // break out of the pretracking loop if we are not pretracking
+ if (! I_am_pretracking) break;
+
+ if (! found_this_horizon) {
+ // the surface failed; backtrack and continue
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: solving failed; backtracking");
+ AH_data_ptr->ps_ptr->ghosted_gridfn_copy (gfns::gfn__save_h, gfns::gfn__h);
+ if (pretracking_have_lower_bound) {
+ assert (AH_data_ptr->pretracking_value >= pretracking_lower_bound - 1.0e-10 * fabs(pretracking_lower_bound));
+ }
+ pretracking_have_lower_bound = true;
+ pretracking_lower_bound = AH_data_ptr->pretracking_value;
+ if (pretracking_have_upper_bound) {
+ AH_data_ptr->pretracking_delta = 0.5 * (pretracking_upper_bound - pretracking_lower_bound);
+ AH_data_ptr->pretracking_value = pretracking_lower_bound + 0.5 * (pretracking_upper_bound - pretracking_lower_bound);
+ } else {
+ AH_data_ptr->pretracking_delta = fabs(AH_data_ptr->pretracking_delta);
+ AH_data_ptr->pretracking_delta *= 2.0;
+ if (AH_data_ptr->pretracking_delta > AH_data_ptr->pretracking_maximum_delta) {
+ AH_data_ptr->pretracking_delta = AH_data_ptr->pretracking_maximum_delta;
+ }
+ AH_data_ptr->pretracking_value += AH_data_ptr->pretracking_delta;
+ if (AH_data_ptr->pretracking_value > AH_data_ptr->pretracking_maximum_value) {
+ AH_data_ptr->pretracking_value = AH_data_ptr->pretracking_maximum_value;
+ }
+ }
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: new value %g, delta %g",
+ double(AH_data_ptr->pretracking_value),
+ double(AH_data_ptr->pretracking_delta));
+ if (pretracking_lower_bound >= AH_data_ptr->pretracking_maximum_value * 0.9999999999) {
+ // give up
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: upper bound reached, giving up.");
+ I_am_pretracking = false;
+ I_was_pretracking = true;
+ pretracking_was_successful = false;
+ // restore old pretracking goal
+ AH_data_ptr->pretracking_value = old_pretracking_value;
+ } else if (AH_data_ptr->pretracking_delta < AH_data_ptr->pretracking_minimum_delta) {
+ // give up
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: step size too small, giving up.");
+ I_am_pretracking = false;
+ I_was_pretracking = true;
+ pretracking_was_successful = true;
+ }
+ } else {
+ // the surface was okay
+ // get mean expansion
+ assert (pretracking_have_horizon_info);
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: solving succeeded; expansion is now %g",
+ double(pretracking_mean_expansion));
+// if (fabs(AH_data_ptr->pretracking_value) > 2.0 * AH_data_ptr->pretracking_minimum_delta) {
+ if (AH_data_ptr->pretracking_value > AH_data_ptr->pretracking_minimum_value + 1.0e-10 * AH_data_ptr->pretracking_minimum_delta) {
+ // it is not yet a horizon
+ if (pretracking_have_upper_bound) {
+ assert (AH_data_ptr->pretracking_value <= pretracking_upper_bound + 1.0e-10 * fabs(pretracking_upper_bound));
+ }
+ pretracking_have_upper_bound = true;
+ pretracking_upper_bound = AH_data_ptr->pretracking_value;
+ if (pretracking_have_lower_bound) {
+#if 1
+ // TODO
+ // move lower bound further down
+ pretracking_lower_bound -= pretracking_upper_bound - pretracking_lower_bound;
+ if (pretracking_lower_bound < AH_data_ptr->pretracking_minimum_value) pretracking_lower_bound = AH_data_ptr->pretracking_minimum_value;
+#endif
+ AH_data_ptr->pretracking_delta = 0.5 * (pretracking_lower_bound - pretracking_upper_bound);
+ AH_data_ptr->pretracking_value = pretracking_lower_bound + 0.5 * (pretracking_upper_bound - pretracking_lower_bound);
+ } else {
+ AH_data_ptr->pretracking_delta = - fabs(AH_data_ptr->pretracking_delta);
+ AH_data_ptr->pretracking_delta *= 2.0;
+ if (- AH_data_ptr->pretracking_delta > AH_data_ptr->pretracking_maximum_delta) {
+ AH_data_ptr->pretracking_delta = - AH_data_ptr->pretracking_maximum_delta;
+ }
+ AH_data_ptr->pretracking_value += AH_data_ptr->pretracking_delta;
+ if (AH_data_ptr->pretracking_value < AH_data_ptr->pretracking_minimum_value) {
+ AH_data_ptr->pretracking_value = AH_data_ptr->pretracking_minimum_value;
+ }
+ }
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: new value radius %g, delta %g",
+ double(AH_data_ptr->pretracking_value),
+ double(AH_data_ptr->pretracking_delta));
+ if (pretracking_upper_bound <= AH_data_ptr->pretracking_minimum_value * 1.00000000001) {
+ // give up
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: lower bound reached, giving up.");
+ I_am_pretracking = false;
+ I_was_pretracking = true;
+ pretracking_was_successful = false;
+ // restore old pretracking goal
+ AH_data_ptr->pretracking_value = old_pretracking_value;
+ } else if (- AH_data_ptr->pretracking_delta < AH_data_ptr->pretracking_minimum_delta) {
+ // give up
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: step size too small, giving up.");
+ I_am_pretracking = false;
+ I_was_pretracking = true;
+ pretracking_was_successful = true;
+ }
+ } else {
+ // a true horizon was found; we are done
+ CCTK_VInfo(CCTK_THORNSTRING,
+ "Pretracking: done.");
+ I_am_pretracking = false;
+ I_was_pretracking = true;
+ pretracking_was_successful = true;
+ }
+ } // if surface found
+
+ // end of pretracking loop
+ }
+ I_am_pretracking = false;
+
// end of this horizon
}
@@ -703,8 +1275,8 @@ return any_proc_needs_more_iterations;
//
// This function (which must be called on *every* processor) broadcasts
-// the BH diagnostics and optionally also the (ghosted) horizon shape,
-// from a specified processor to all processors.
+// the BH diagnostics and (ghosted) horizon shape from a specified processor
+// to all processors.
//
// The present implementation of this function uses the Cactus reduction
// API. If AHFinderDirect is ported to some other software environment,
@@ -715,15 +1287,11 @@ return any_proc_needs_more_iterations;
// GH --> The Cactus grid hierarchy.
// broadcast_flag = true on the broadcasting processor,
// false on all other processors
-// broadcast_horizon_shape = true to broadcast the (ghosted) horizon shape
-// as well as the BH diagnostics
-// false to only broadcast the BH diagnostics
// BH_diagnostics = On the broadcasting processor, this is the BH diagnostics
// to broadcast; on all other processors, it's set to the
// broadcast BH diagnostics.
-// ps = Used only if broadcast_horizon_shape; in this case...
-// On the broadcasting processor, gfn__h is broadcast;
-// on all other processors, gfn__h is set to the broadcast values.
+// ps = On the broadcasting processor, gfn__h is broadcast;
+// on all other processors, gfn__h is set to the broadcast values.
// horizon_buffers = Internal buffers for use in the broadcast;
// if N_buffer == 0 then we set N_buffer and allocate
// the buffers.
@@ -731,6 +1299,7 @@ return any_proc_needs_more_iterations;
namespace {
void broadcast_horizon_data(const cGH* GH,
bool broadcast_flag, bool broadcast_horizon_shape,
+ struct AH_data& AH_data,
struct BH_diagnostics& BH_diagnostics,
patch_system& ps,
struct horizon_buffers& horizon_buffers)
@@ -751,7 +1320,8 @@ if (horizon_buffers.N_buffer == 0)
// allocate the buffers
horizon_buffers.N_buffer
= BH_diagnostics::N_buffer
- + (broadcast_horizon_shape ? ps.N_grid_points() : 0);
+ + (broadcast_horizon_shape ? ps.N_grid_points() : 0)
+ + 4;
horizon_buffers.send_buffer
= new CCTK_REAL[horizon_buffers.N_buffer];
horizon_buffers.receive_buffer
@@ -762,9 +1332,9 @@ if (broadcast_flag)
then {
// pack the data to be broadcast into the send buffer
BH_diagnostics.copy_to_buffer(horizon_buffers.send_buffer);
+ int posn = BH_diagnostics::N_buffer;
if (broadcast_horizon_shape)
then {
- int posn = BH_diagnostics::N_buffer;
for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
{
const patch& p = ps.ith_patch(pn);
@@ -782,8 +1352,12 @@ if (broadcast_flag)
}
}
}
- assert( posn == horizon_buffers.N_buffer );
}
+ horizon_buffers.send_buffer[posn++] = AH_data.initial_find_flag;
+ horizon_buffers.send_buffer[posn++] = AH_data.really_initial_find_flag;
+ horizon_buffers.send_buffer[posn++] = AH_data.search_flag;
+ horizon_buffers.send_buffer[posn++] = AH_data.found_flag;
+ assert( posn == horizon_buffers.N_buffer );
}
else jtutil::zero_C_array(horizon_buffers.N_buffer,
horizon_buffers.send_buffer);
@@ -814,9 +1388,12 @@ if (!broadcast_flag)
then {
// unpack the data from the receive buffer
BH_diagnostics.copy_from_buffer(horizon_buffers.receive_buffer);
+ ps.origin_x(BH_diagnostics.origin_x);
+ ps.origin_y(BH_diagnostics.origin_y);
+ ps.origin_z(BH_diagnostics.origin_z);
+ int posn = BH_diagnostics::N_buffer;
if (broadcast_horizon_shape)
then {
- int posn = BH_diagnostics::N_buffer;
for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
{
patch& p = ps.ith_patch(pn);
@@ -834,12 +1411,16 @@ if (!broadcast_flag)
}
}
}
- assert( posn == horizon_buffers.N_buffer );
// recover the full ghosted-grid horizon shape
// (we only broadcast the nominal-grid shape)
ps.synchronize();
}
+ AH_data.initial_find_flag = horizon_buffers.receive_buffer[posn++];
+ AH_data.really_initial_find_flag = horizon_buffers.receive_buffer[posn++];
+ AH_data.search_flag = horizon_buffers.receive_buffer[posn++];
+ AH_data.found_flag = horizon_buffers.receive_buffer[posn++];
+ assert( posn == horizon_buffers.N_buffer );
}
}
}
@@ -918,6 +1499,7 @@ const fp scale = (max_Delta_h <= max_allowable_Delta_h)
if (verbose_info.print_algorithm_details)
then {
+
if (scale == 1.0)
then CCTK_VInfo(CCTK_THORNSTRING,
"h += Delta_h (rms-norm=%.1e, infinity-norm=%.1e)",
@@ -948,6 +1530,14 @@ if (verbose_info.print_algorithm_details)
}
}
}
+
+if (ps.N_additional_points())
+ {
+ const int np = ps.N_grid_points();
+ const int gnp = ps.ghosted_N_grid_points();
+ ps.ghosted_gridfn_data(gfns::gfn__h)[gnp]
+ -= scale * ps.gridfn_data(gfns::gfn__Delta_h)[np];
+ }
}
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 04:44:46 +0000