aboutsummaryrefslogtreecommitdiff
path: root/src/driver/find_horizons.cc
diff options
context:
space:
mode:
Diffstat (limited to 'src/driver/find_horizons.cc')
-rw-r--r--src/driver/find_horizons.cc258
1 files changed, 168 insertions, 90 deletions
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index bf3f25a..45cc780 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -6,8 +6,6 @@
// AHFinderDirect_find_horizons - top-level driver to find apparent horizons
///
/// find_horizon - find a horizon
-/// print_summary_of_which_horizons_found
-///
/// do_evaluate_expansions
/// do_test_expansion_Jacobian
///
@@ -15,10 +13,12 @@
#include <stdio.h>
#include <assert.h>
#include <math.h>
+#include <string.h>
-#include "util_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 +26,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 +49,6 @@
// all the code in this file is inside this namespace
namespace AHFinderDirect
{
-using jtutil::error_exit;
//******************************************************************************
@@ -63,8 +63,6 @@ extern struct state state;
// ***** prototypes for functions local to this file
//
namespace {
-void print_summary_of_which_horizons_found
- (int N_horizons, const struct AH_data* const* AH_data_array);
void do_evaluate_expansions(int my_proc, int N_horizons,
horizon_sequence& hs,
struct AH_data* const AH_data_array[],
@@ -89,6 +87,36 @@ void do_test_expansion_Jacobians(int my_proc, int N_horizons,
//******************************************************************************
//
+// This function is called by the Cactus scheduler to import
+// the excision mask.
+//
+extern "C"
+ void AHFinderDirect_import_mask(CCTK_ARGUMENTS)
+{
+DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
+
+for (int k=0; k<cctk_lsh[2]; ++k)
+for (int j=0; j<cctk_lsh[1]; ++j)
+for (int i=0; i<cctk_lsh[0]; ++i)
+{
+ const int ind = CCTK_GFINDEX3D(cctkGH,i,j,k);
+ // zero means: point can be used,
+ // non-zero means: point must be avoided
+ ahmask[ind] = 0;
+ if (use_mask)
+ // grid points with mask values of 1.0 can be used,
+ // values of 0.0 and 0.5 must be avoided.
+ // the excision boundary cannot be used because
+ // (a) it is inaccurate
+ // (b) it does not respect the symmetries e.g. in Kerr.
+ then ahmask[ind] = fabs(emask[ind] - 1.0) > 0.01;
+}
+}
+
+//******************************************************************************
+
+//
// This function is called by the Cactus scheduler to find the apparent
// horizon or horizons in the current slice.
//
@@ -96,22 +124,44 @@ extern "C"
void AHFinderDirect_find_horizons(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
-// only try to find horizons every find_every time steps
-if (! state.find_now(cctk_iteration))
- then return; // *** NO-OP RETURN ***
+// determine whether a horizon should be found at this iteration
+bool find_any = false;
+for (int hn = 1; hn <= state.my_hs->N_horizons(); ++ hn)
+{
+ // only try to find horizons every find_every time steps
+ const int my_find_after = find_after_individual[hn];
+ const int my_dont_find_after = dont_find_after_individual[hn];
+ const fp my_find_after_time = find_after_individual_time[hn];
+ const fp my_dont_find_after_time = dont_find_after_individual_time[hn];
+ const int my_find_every = (find_every_individual[hn] >= 0
+ ? find_every_individual[hn]
+ : find_every);
+ const bool find_this = cctk_iteration >= my_find_after
+ && (my_dont_find_after < 0
+ ? true
+ : cctk_iteration <= my_dont_find_after)
+ && cctk_time >= my_find_after_time
+ && (my_dont_find_after_time <= my_find_after_time
+ ? true
+ : cctk_time <= my_dont_find_after_time)
+ && my_find_every > 0
+ && cctk_iteration % my_find_every == 0
+ && ! disable_horizon[hn];
+ struct AH_data& AH_data = *state.AH_data_array[hn];
+ AH_data.search_flag = find_this;
+ find_any = find_any || find_this;
+}
+if (! find_any) return;
if (state.timer_handle >= 0)
then CCTK_TimerResetI(state.timer_handle);
const int my_proc = state.my_proc;
-const int N_horizons = state.N_horizons;
horizon_sequence& hs = *state.my_hs;
const bool active_flag = hs.has_genuine_horizons();
-const bool broadcast_horizon_shape
- = state.always_broadcast_horizon_shape
- || state.mask_info.set_mask_for_any_horizon
- || state.store_info_in_SS_vars_for_any_horizon;
+const bool broadcast_horizon_shape = true;
struct cactus_grid_info& cgi = state.cgi;
const struct geometry_info& gi = state.gi;
@@ -122,13 +172,50 @@ const struct verbose_info& verbose_info = state.verbose_info;
// what are the semantics of the Cactus gxx variables? (these may
// change from one call to another, so we have to re-check each time)
-if (CCTK_Equals(state.metric_type, "physical"))
+if (CCTK_Equals(metric_type, "physical"))
then cgi.use_Cactus_conformal_metric = false;
-else if (CCTK_Equals(state.metric_type, "static conformal"))
+else if (CCTK_Equals(metric_type, "static conformal"))
then cgi.use_Cactus_conformal_metric = (*conformal_state > 0);
else CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"AHFinderDirect_find_horizons(): unknown ADMBase::metric_type=\"%s\"!",
- state.metric_type); /*NOTREACHED*/
+"AHFinderDirect_find_horizons(): unknown metric_type=\"%s\"!",
+ metric_type); /*NOTREACHED*/
+
+// update parameters
+IO_info.output_ASCII_files = (output_ASCII_files != 0);
+IO_info.output_HDF5_files = (output_HDF5_files != 0);
+IO_info.output_initial_guess = (output_initial_guess != 0);
+IO_info.output_h_every = output_h_every;
+IO_info.output_Theta_every = output_Theta_every;
+IO_info.output_mean_curvature_every = output_mean_curvature_every;
+IO_info.output_h = false; // dummy value
+IO_info.output_Theta = false; // dummy value
+IO_info.output_mean_curvature = false; // dummy value
+
+IO_info.output_BH_diagnostics = (output_BH_diagnostics != 0);
+IO_info.BH_diagnostics_directory
+ = (strlen(BH_diagnostics_directory) == 0)
+ ? /* IO:: */ out_dir
+ : BH_diagnostics_directory;
+IO_info.BH_diagnostics_base_file_name = BH_diagnostics_base_file_name;
+IO_info.BH_diagnostics_file_name_extension = BH_diagnostics_file_name_extension;
+
+IO_info.output_ghost_zones_for_h = (output_ghost_zones_for_h != 0);
+IO_info.ASCII_gnuplot_file_name_extension = ASCII_gnuplot_file_name_extension;
+IO_info.HDF5_file_name_extension = HDF5_file_name_extension;
+IO_info.h_directory
+ = (strlen(h_directory) == 0)
+ ? /* IO:: */ out_dir
+ : h_directory;
+IO_info.h_base_file_name = h_base_file_name;
+IO_info.Theta_base_file_name = Theta_base_file_name;
+IO_info.mean_curvature_base_file_name = mean_curvature_base_file_name;
+IO_info.Delta_h_base_file_name = Delta_h_base_file_name;
+IO_info.h_min_digits = h_min_digits;
+IO_info.Jacobian_base_file_name = Jacobian_base_file_name;
+IO_info.output_OpenDX_control_files = (output_OpenDX_control_files != 0);
+IO_info.OpenDX_control_file_name_extension = OpenDX_control_file_name_extension;
+IO_info.time_iteration = 0;
+IO_info.time = 0.0;
// get the Cactus time step and decide if we want to output h and/or Theta now
IO_info.time_iteration = cctk_iteration;
@@ -139,6 +226,9 @@ IO_info.output_h
IO_info.output_Theta
= (IO_info.output_Theta_every > 0)
&& ((IO_info.time_iteration % IO_info.output_Theta_every) == 0);
+IO_info.output_mean_curvature
+ = (IO_info.output_mean_curvature_every > 0)
+ && ((IO_info.time_iteration % IO_info.output_mean_curvature_every) == 0);
// set initial guess for any (genuine) horizons that need it,
// i.e. for any (genuine) horizons where we didn't find the horizon previously
@@ -146,21 +236,39 @@ IO_info.output_Theta
{
assert( state.AH_data_array[hn] != NULL );
struct AH_data& AH_data = *state.AH_data_array[hn];
+ if (verbose_info.print_algorithm_highlights) {
+ printf ("AHF find_horizons[%d] initial_find_flag=%d\n", hn, (int) AH_data.initial_find_flag);
+ printf ("AHF find_horizons[%d] really_initial_find_flag=%d\n", hn, (int) AH_data.really_initial_find_flag);
+ printf ("AHF find_horizons[%d] search_flag=%d\n", hn, (int) AH_data.search_flag);
+ printf ("AHF find_horizons[%d] found_flag=%d\n", hn, (int) AH_data.found_flag);
+ }
if (AH_data.found_flag)
- then AH_data.initial_find_flag = false;
+ then {
+ AH_data.initial_find_flag = false;
+ AH_data.really_initial_find_flag = false;
+ }
else {
- patch_system& ps = *AH_data.ps_ptr;
- setup_initial_guess(ps,
- AH_data.initial_guess_info,
- IO_info,
- hn, N_horizons, verbose_info);
- if (active_flag && IO_info.output_initial_guess)
- then output_gridfn(ps, gfns::gfn__h,
- IO_info, IO_info.h_base_file_name,
- hn, verbose_info
- .print_algorithm_highlights);
- AH_data.initial_find_flag = true;
- }
+ if (AH_data.really_initial_find_flag
+ || AH_data.initial_guess_info.reset_horizon_after_not_finding)
+ then {
+ patch_system& ps = *AH_data.ps_ptr;
+ if (verbose_info.print_algorithm_highlights) {
+ printf ("AHF find_horizons[%d] setup_initial_guess\n", hn);
+ }
+ setup_initial_guess(ps,
+ AH_data.initial_guess_info,
+ IO_info,
+ hn, N_horizons, verbose_info);
+ if (active_flag && IO_info.output_initial_guess)
+ then output_gridfn(ps, gfns::gfn__h,
+ "h", cgi.GH,
+ IO_info, IO_info.h_base_file_name,
+ IO_info.h_min_digits,
+ hn, verbose_info
+ .print_algorithm_highlights);
+ AH_data.initial_find_flag = true;
+ }
+ }
}
//
@@ -180,7 +288,7 @@ case method__test_expansion_Jacobians:
do_test_expansion_Jacobians(my_proc, N_horizons,
state.AH_data_array,
cgi, gi, Jac_info,
- state.test_all_Jacobian_compute_methods,
+ (test_all_Jacobian_compute_methods != 0),
IO_info, error_info, verbose_info,
state.timer_handle);
break;
@@ -198,10 +306,6 @@ case method__find_horizons:
state.isb);
if (state.timer_handle >= 0)
then CCTK_TimerStopI(state.timer_handle);
- if (verbose_info.print_physics_highlights
- && !verbose_info.print_physics_details)
- then print_summary_of_which_horizons_found(N_horizons,
- state.AH_data_array);
break;
}
@@ -223,53 +327,6 @@ if (state.timer_handle >= 0)
}
//******************************************************************************
-
-//
-// This function prints (using CCTK_VInfo()) a 1-line summary of
-// which AHs were / were not found.
-//
-namespace {
-void print_summary_of_which_horizons_found
- (int N_horizons, const struct AH_data* const* AH_data_array)
-{
-const int hn_buffer_size = 10; // buffer for single "%d" for hn
-char hn_buffer[hn_buffer_size];
-const int msg_buffer_size = 200; // buffer for the entire line
-char msg_buffer[msg_buffer_size];
-
-msg_buffer[0] = '\0';
-Util_Strlcpy(msg_buffer, "found horizon(s) ", msg_buffer_size);;
-
-bool already_found_flag = false; // running inclusive-or of found_flag
- for (int hn = 1 ; hn <= N_horizons ; ++hn)
- {
- assert(AH_data_array[hn] != NULL);
- bool found_flag = AH_data_array[hn]->found_flag;
- if (found_flag)
- then {
- if (already_found_flag)
- then Util_Strlcat(msg_buffer, ",", msg_buffer_size);
- already_found_flag = true;
- snprintf(hn_buffer, hn_buffer_size, "%d", hn);
- Util_Strlcat(msg_buffer, hn_buffer, msg_buffer_size);
- }
- }
-
-if (already_found_flag)
- then {
- // we found one or more horizons
- Util_Strlcat(msg_buffer, " of ", msg_buffer_size);
-
- snprintf(hn_buffer, hn_buffer_size, "%d", N_horizons);
- Util_Strlcat(msg_buffer, hn_buffer, msg_buffer_size);
- }
- else Util_Strlcpy(msg_buffer, "no horizons found", msg_buffer_size);
-
-CCTK_VInfo(CCTK_THORNSTRING, msg_buffer);
-}
- }
-
-//******************************************************************************
//******************************************************************************
//******************************************************************************
@@ -317,7 +374,8 @@ if (active_flag)
if (timer_handle >= 0)
then CCTK_TimerStartI(timer_handle);
jtutil::norm<fp> Theta_norms;
- const bool Theta_ok = expansion(&ps, -AH_data.surface_expansion,
+ const bool Theta_ok = expansion(&ps,
+ AH_data.compute_info,
cgi, gi,
error_info, true,// initial eval
false, // no Jacobian coeffs
@@ -328,7 +386,9 @@ if (active_flag)
if (IO_info.output_h)
then output_gridfn(ps, gfns::gfn__h,
+ "h", cgi.GH,
IO_info, IO_info.h_base_file_name,
+ IO_info.h_min_digits,
hn, verbose_info.print_algorithm_details);
if (Theta_ok)
@@ -338,17 +398,28 @@ if (active_flag)
Theta_norms.rms_norm(), Theta_norms.infinity_norm());
if (IO_info.output_Theta)
then output_gridfn(ps, gfns::gfn__Theta,
+ "Theta", cgi.GH,
IO_info, IO_info
.Theta_base_file_name,
+ IO_info.h_min_digits,
+ hn, verbose_info
+ .print_algorithm_details);
+ if (IO_info.output_mean_curvature)
+ then output_gridfn(ps, gfns::gfn__mean_curvature,
+ "mean_curvature", cgi.GH,
+ IO_info, IO_info
+ .mean_curvature_base_file_name,
+ IO_info.h_min_digits,
hn, verbose_info
.print_algorithm_details);
}
}
}
else {
+ struct what_to_compute new_compute_info;
for (int i = 0 ; i < N_horizons ; ++i)
{
- expansion(NULL, 0.0, // dummy computation
+ expansion(NULL, new_compute_info,
cgi, gi,
error_info, true); // initial evaluation
}
@@ -402,17 +473,22 @@ const int hn = 1;
struct AH_data* const AH_data_ptr = active_flag ? AH_data_array[hn] : NULL;
patch_system* const ps_ptr = active_flag ? AH_data_ptr->ps_ptr : NULL;
-const fp add_to_expansion = active_flag ? -AH_data_ptr->surface_expansion : 0.0;
+struct what_to_compute dummy_compute_info;
+struct what_to_compute & compute_info =
+ active_flag
+ ? AH_data_ptr->compute_info
+ : dummy_compute_info;
//
// numerical-perturbation Jacobian
//
Jacobian* Jac_NP_ptr = active_flag ? AH_data_ptr->Jac_ptr : NULL;
-expansion(ps_ptr, add_to_expansion,
+expansion(ps_ptr, compute_info,
cgi, gi,
error_info, true); // initial evaluation
Jac_info.Jacobian_compute_method = Jacobian__numerical_perturbation;
-expansion_Jacobian(ps_ptr, Jac_NP_ptr, add_to_expansion,
+expansion_Jacobian(ps_ptr, Jac_NP_ptr,
+ compute_info,
cgi, gi, Jac_info,
error_info, true, // initial evaluation
print_msg_flag);
@@ -426,12 +502,13 @@ if (test_all_Jacobian_compute_methods)
*ps_ptr,
verbose_info.print_algorithm_details)
: NULL;
- expansion(ps_ptr, add_to_expansion,
+ expansion(ps_ptr, compute_info,
cgi, gi,
error_info, true, // initial evaluation
true); // compute SD Jacobian coeffs
Jac_info.Jacobian_compute_method = Jacobian__symbolic_diff_with_FD_dr;
- expansion_Jacobian(ps_ptr, Jac_SD_FDdr_ptr, add_to_expansion,
+ expansion_Jacobian(ps_ptr, Jac_SD_FDdr_ptr,
+ compute_info,
cgi, gi, Jac_info,
error_info, true, // initial evaluation
print_msg_flag);
@@ -441,6 +518,7 @@ if (active_flag)
then output_Jacobians(*ps_ptr,
Jac_NP_ptr, Jac_SD_FDdr_ptr,
IO_info, IO_info.Jacobian_base_file_name,
+ IO_info.h_min_digits,
hn, print_msg_flag);
}
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 03:45:44 +0000