aboutsummaryrefslogtreecommitdiff
path: root/src/driver/setup.cc
diff options
context:
space:
mode:
Diffstat (limited to 'src/driver/setup.cc')
-rw-r--r--src/driver/setup.cc392
1 files changed, 225 insertions, 167 deletions
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index d010d32..5217097 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -5,24 +5,23 @@
// <<<access to persistent data>>>
//
// AHFinderDirect_setup - top-level driver to setup persistent data structures
-// AHFinderDirect_update - top-level driver to update from steerable parameters
-//
-/// set_mask_pars - set internal data structures from mask parameters
///
/// decode_method - decode the method parameter
/// decode_verbose_level - decode the verbose_level parameter
-/// decode_horizon_file_format - decode the horizon_file_format parameter
///
/// allocate_horizons_to_processor - choose which horizons this proc will find
///
/// choose_patch_system_type - choose patch system type
///
+#include <assert.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <string.h>
+#include <vector>
+
#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"
@@ -36,6 +35,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"
@@ -58,7 +58,6 @@
// all the code in this file is inside this namespace
namespace AHFinderDirect
{
-using jtutil::error_exit;
//******************************************************************************
@@ -67,17 +66,14 @@ using jtutil::error_exit;
//
namespace {
-void set_mask_pars(CCTK_ARGUMENTS);
-
enum method
decode_method(const char method_string[]);
enum verbose_level
decode_verbose_level(const char verbose_level_string[]);
-enum horizon_file_format
- decode_horizon_file_format(const char horizon_file_format_string[]);
int allocate_horizons_to_processor(int N_procs, int my_proc,
int N_horizons, bool multiproc_flag,
+ const CCTK_INT depends_on[],
horizon_sequence& my_hs,
const struct verbose_info& verbose_info);
@@ -113,9 +109,36 @@ CCTK_VInfo(CCTK_THORNSTRING,
//
+// check parameters
+//
+int need_zones = 0;
+for (int n=1; n<N_horizons; ++n) {
+ need_zones = jtutil::max(need_zones, N_zones_per_right_angle[n]);
+}
+if (need_zones > max_N_zones_per_right_angle) {
+ CCTK_VWarn (FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "AHFinderDirect_setup(): "
+ "The parameter max_N_zones_per_right_angle must be at least the maximum of all N_zones_per_right_angle[] parameters. "
+ "Set max_N_zones_per_right_angle to %d or higher to continue.",
+ need_zones); /*NOTREACHED*/
+ }
+for (int n=1; n<N_horizons; ++n) {
+ if (depends_on[n] != 0) {
+ assert (depends_on[n] >= 1 && depends_on[n] < n);
+ if (N_zones_per_right_angle[n] != N_zones_per_right_angle[depends_on[n]]) {
+ CCTK_VWarn (FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "AHFinderDirect_setup(): "
+ "The parameter N_zones_per_right_angle must be the same for a horizon as for the horizon on which it depends. "
+ "Horizon %d depends on horizon %d, but they have different resolutions.",
+ n, depends_on[n]); /*NOTREACHED*/
+ }
+ }
+}
+
+
+//
// basic setup
//
-state.find_every = find_every;
state.method = decode_method(method);
state.error_info.warn_level__point_outside__initial
@@ -144,9 +167,7 @@ verbose_info.print_algorithm_details
verbose_info.print_algorithm_debug
= (state.verbose_info.verbose_level >= verbose_level__algorithm_debug);
-state.metric_type = /* ADMBase:: */ metric_type;
-
-state.timer_handle = (print_timing_stats != 0) ? CCTK_TimerCreateI() : -1;
+state.timer_handle = (print_timing_stats != 0) ? CCTK_TimerCreate("finding apparent horizons") : -1;
state.N_procs = CCTK_nProcs(cctkGH);
state.my_proc = CCTK_MyProc(cctkGH);
@@ -173,6 +194,7 @@ if (cgi.coord_system_handle < 0)
coordinate_system_name); /*NOTREACHED*/
cgi.use_Cactus_conformal_metric = false; // dummy value, may change later
+cgi.mask_varindex = Cactus_gridfn_varindex("AHFinderDirect::ahmask");
cgi.g_dd_11_varindex = Cactus_gridfn_varindex("ADMBase::gxx");
cgi.g_dd_12_varindex = Cactus_gridfn_varindex("ADMBase::gxy");
cgi.g_dd_13_varindex = Cactus_gridfn_varindex("ADMBase::gxz");
@@ -222,9 +244,6 @@ gi.check_that_geometry_is_finite = (check_that_geometry_is_finite != 0);
//
// Jacobian info
//
-state.test_all_Jacobian_compute_methods
- = (test_all_Jacobian_compute_methods != 0);
-
struct Jacobian_info& Jac_info = state.Jac_info;
Jac_info.Jacobian_compute_method
= decode_Jacobian_compute_method(Jacobian_compute_method);
@@ -243,6 +262,8 @@ solver_info.linear_solver_pars.ILUCG_pars.error_tolerance
= ILUCG__error_tolerance;
solver_info.linear_solver_pars.ILUCG_pars.limit_CG_iterations
= (ILUCG__limit_CG_iterations != 0);
+solver_info.linear_solver_pars.UMFPACK_pars.N_II_iterations
+ = UMFPACK__N_II_iterations;
solver_info.max_Newton_iterations__initial
= max_Newton_iterations__initial;
solver_info.max_Newton_iterations__subsequent
@@ -250,9 +271,13 @@ solver_info.max_Newton_iterations__subsequent
solver_info.max_allowable_Delta_h_over_h = max_allowable_Delta_h_over_h;
solver_info.Theta_norm_for_convergence = Theta_norm_for_convergence;
solver_info.max_allowable_Theta = max_allowable_Theta;
+solver_info.max_allowable_Theta_growth_iterations
+ = max_allowable_Theta_growth_iterations;
+solver_info.max_allowable_Theta_nonshrink_iterations
+ = max_allowable_Theta_nonshrink_iterations;
// ... horizon numbers run from 1 to N_horizons inclusive
// so the array size is N_horizons+1
-solver_info.max_allowable_horizon_radius = new double[state.N_horizons+1];
+solver_info.max_allowable_horizon_radius = new fp[state.N_horizons+1];
{
for (int hn = 0 ; hn <= N_horizons ; ++hn)
{
@@ -260,18 +285,22 @@ solver_info.max_allowable_horizon_radius = new double[state.N_horizons+1];
= max_allowable_horizon_radius[hn];
}
}
+solver_info.want_expansion_gradients = want_expansion_gradients;
//
// I/O info
//
struct IO_info& IO_info = state.IO_info;
-IO_info.horizon_file_format = decode_horizon_file_format(horizon_file_format);
+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
@@ -290,7 +319,9 @@ IO_info.h_directory
: 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;
@@ -304,8 +335,6 @@ IO_info.time = 0.0;
state.BH_diagnostics_info.integral_method
= patch::decode_integration_method(integral_method);
-state.always_broadcast_horizon_shape = (always_broadcast_horizon_shape != 0);
-
//
// mask parameters
@@ -326,33 +355,43 @@ mask_info.set_mask_for_this_horizon = new bool[N_horizons+1];
}
}
if (mask_info.set_mask_for_any_horizon)
- then set_mask_pars(CCTK_PASS_CTOC);
-
-
-//
-// announce parameters
-//
-state.announce_centroid_flag = (which_horizon_to_announce_centroid != 0);
-if (state.announce_centroid_flag)
then {
- state.which_horizon_to_announce_centroid
- = which_horizon_to_announce_centroid;
- if ( (state.which_horizon_to_announce_centroid < 1)
- || (state.which_horizon_to_announce_centroid > state.N_horizons) )
- then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-" AHFinderDirect_setup():\n"
-" invalid which_horizon_to_announce_centroid = %d!\n"
-" (valid range is [1,N_horizons=%d])!\n"
- ,
- state.which_horizon_to_announce_centroid,
- state.N_horizons); /*NOTREACHED*/
+ mask_info.radius_multiplier = mask_radius_multiplier;
+ mask_info.radius_offset = mask_radius_offset;
+ mask_info.buffer_thickness = mask_buffer_thickness;
+ mask_info.mask_is_noshrink = mask_is_noshrink;
+ mask_info.min_horizon_radius_points_for_mask
+ = min_horizon_radius_points_for_mask;
+ mask_info.set_old_style_mask = (set_old_style_mask != 0);
+ mask_info.set_new_style_mask = (set_new_style_mask != 0);
+ if (mask_info.set_old_style_mask)
+ then {
+ struct mask_info::old_style_mask_info& osmi
+ = mask_info.old_style_mask_info;
+ osmi.gridfn_name = old_style_mask_gridfn_name;
+ osmi.gridfn_varindex = Cactus_gridfn_varindex(osmi.gridfn_name);
+ osmi.gridfn_dataptr = NULL; // dummy value; fixup later
+ osmi.inside_value = old_style_mask_inside_value;
+ osmi.buffer_value = old_style_mask_buffer_value;
+ osmi.outside_value = old_style_mask_outside_value;
+ }
+ if (mask_info.set_new_style_mask)
+ then {
+ struct mask_info::new_style_mask_info& nsmi
+ = mask_info.new_style_mask_info;
+ nsmi.gridfn_name = new_style_mask_gridfn_name;
+ nsmi.gridfn_varindex = Cactus_gridfn_varindex(nsmi.gridfn_name);
+ nsmi.gridfn_dataptr = NULL; // dummy value; fixup later
+ nsmi.bitfield_name = new_style_mask_bitfield_name;
+ nsmi.bitfield_bitmask = 0; // dummy value; fixup later
+ nsmi.inside_value = new_style_mask_inside_value;
+ nsmi.buffer_value = new_style_mask_buffer_value;
+ nsmi.outside_value = new_style_mask_outside_value;
+ nsmi.inside_bitvalue = 0; // dummy value; fixup later
+ nsmi.buffer_bitvalue = 0; // dummy value; fixup later
+ nsmi.outside_bitvalue = 0; // dummy value; fixup later
+ }
}
- else state.which_horizon_to_announce_centroid = 0; // dummy value; unused
-
-// initial value, will get each horizon's store_info_in_SS_vars
-// flag inclusive-ored into it below
-state.store_info_in_SS_vars_for_any_horizon = false;
//
@@ -371,6 +410,7 @@ const bool multiproc_flag = (state.method == method__find_horizons);
state.N_active_procs
= allocate_horizons_to_processor(state.N_procs, state.my_proc,
state.N_horizons, multiproc_flag,
+ depends_on,
hs,
verbose_info);
@@ -474,6 +514,8 @@ if (strlen(surface_interpolator_name) > 0)
true, verbose_info.print_algorithm_details);
patch_system& ps = *AH_data.ps_ptr;
if (genuine_flag)
+ then ps.set_gridfn_to_constant(0.0, gfns::gfn__zero);
+ if (genuine_flag)
then ps.set_gridfn_to_constant(1.0, gfns::gfn__one);
AH_data.Jac_ptr = genuine_flag
@@ -482,9 +524,67 @@ if (strlen(surface_interpolator_name) > 0)
verbose_info.print_algorithm_details)
: NULL;
- AH_data.surface_expansion = surface_expansion[hn];
-
- AH_data.initial_find_flag = genuine_flag;
+ AH_data.compute_info.surface_definition =
+ STRING_EQUAL(surface_definition[hn], "expansion")
+ ? definition_expansion
+ : STRING_EQUAL(surface_definition[hn], "inner expansion")
+ ? definition_inner_expansion
+ : STRING_EQUAL(surface_definition[hn], "mean curvature")
+ ? definition_mean_curvature
+ : STRING_EQUAL(surface_definition[hn], "expansion product")
+ ? definition_expansion_product
+ : (CCTK_WARN (0, "internal error"), definition_error);
+ AH_data.compute_info.surface_modification =
+ STRING_EQUAL(surface_modification[hn], "none")
+ ? modification_none
+ : STRING_EQUAL(surface_modification[hn], "radius")
+ ? modification_radius
+ : STRING_EQUAL(surface_modification[hn], "radius^2")
+ ? modification_radius2
+#if 0
+ : STRING_EQUAL(surface_modification[hn], "mean radius")
+ ? modification_mean_radius
+ : STRING_EQUAL(surface_modification[hn], "areal radius")
+ ? modification_areal_radius
+#endif
+ : (CCTK_WARN (0, "internal error"), modification_error);
+ AH_data.compute_info.surface_selection =
+ STRING_EQUAL(surface_selection[hn], "definition")
+ ? selection_definition
+ : STRING_EQUAL(surface_selection[hn], "mean coordinate radius")
+ ? selection_mean_coordinate_radius
+ : STRING_EQUAL(surface_selection[hn], "areal radius")
+ ? selection_areal_radius
+ : STRING_EQUAL(surface_selection[hn], "expansion times mean coordinate radius")
+ ? selection_expansion_mean_coordinate_radius
+ : STRING_EQUAL(surface_selection[hn], "expansion times areal radius")
+ ? selection_expansion_areal_radius
+ : (CCTK_WARN (0, "internal error"), selection_error);
+ AH_data.compute_info.desired_value = desired_value[hn];
+
+ AH_data.move_origins = move_origins;
+
+ AH_data.use_pretracking = use_pretracking[hn];
+ AH_data.pretracking_max_iterations = pretracking_max_iterations[hn];
+
+ AH_data.pretracking_value = pretracking_value[hn];
+ AH_data.pretracking_minimum_value = pretracking_minimum_value[hn];
+ AH_data.pretracking_maximum_value = pretracking_maximum_value[hn];
+ AH_data.pretracking_delta = pretracking_delta[hn];
+ AH_data.pretracking_minimum_delta = pretracking_minimum_delta[hn];
+ AH_data.pretracking_maximum_delta = pretracking_maximum_delta[hn];
+
+ AH_data.depends_on = depends_on[hn];
+ AH_data.desired_value_factor = desired_value_factor[hn];
+ AH_data.desired_value_offset = desired_value_offset[hn];
+
+ AH_data.shiftout_factor = shiftout_factor[hn];
+ AH_data.smoothing_factor = smoothing_factor[hn];
+
+ // AH_data.initial_find_flag = genuine_flag;
+ // AH_data.really_initial_find_flag = AH_data.initial_find_flag;
+ AH_data.initial_find_flag = true;
+ AH_data.really_initial_find_flag = AH_data.initial_find_flag;
if (genuine_flag)
then {
@@ -493,6 +593,8 @@ if (strlen(surface_interpolator_name) > 0)
" setting initial guess parameters etc");
AH_data.initial_guess_info.method
= decode_initial_guess_method(initial_guess_method[hn]);
+ AH_data.initial_guess_info.reset_horizon_after_not_finding
+ = reset_horizon_after_not_finding[hn];
// ... read from named file
AH_data.initial_guess_info.read_from_named_file_info.file_name
= initial_guess__read_from_named_file__file_name[hn];
@@ -542,106 +644,47 @@ if (strlen(surface_interpolator_name) > 0)
= initial_guess__coord_ellipsoid__z_radius[hn];
}
+ AH_data.search_flag = false;
AH_data.found_flag = false;
AH_data.h_files_written = false;
AH_data.BH_diagnostics_fileptr = NULL;
-
- AH_data.store_info_in_SS_vars = (which_surface_to_store_info[hn] >= 0);
- if (AH_data.store_info_in_SS_vars)
- then {
- AH_data.SS_surface_number = which_surface_to_store_info[hn];
- if (AH_data.SS_surface_number
- >= /* SphericalSurface:: */ nsurfaces)
- then CCTK_VWarn(FATAL_ERROR,
- __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-" AHFinderDirect_setup():\n"
-" which_surface_to_store_info[%d] = %d\n"
-" is outside the legal range [0,SphericalSurface::nsurfaces=%d)!\n"
- ,
- hn, AH_data.SS_surface_number,
- int(/* SphericalSurface:: */ nsurfaces));
- /*NOTREACHED*/
- }
- else AH_data.SS_surface_number = 0; // dummy value
-
- state.store_info_in_SS_vars_for_any_horizon
- |= AH_data.store_info_in_SS_vars;
}
}
-}
-//******************************************************************************
-
-//
-// This function is called by the Cactus scheduler to update some of our
-// persistent data structures (stored in struct state ) from parameters
-// which may have been steered. It is also called by AHFinderDirect_setup()
-// (above) as part of the initial setup of these data structures.
-//
-extern "C"
- void AHFinderDirect_update(CCTK_ARGUMENTS)
-{
-set_mask_pars(CCTK_PASS_CTOC);
+ // Initialise the centroids. These values may later be overwritten
+ // when they are recovered from a checkpoint. However, if new
+ // horizons are enabled during recovery, they are then correctly
+ // initialised.
+ for (int n = 0; n < N_horizons; ++ n) {
+ // Horizon centroids are not valid
+ ah_centroid_x[n] = 0.0;
+ ah_centroid_y[n] = 0.0;
+ ah_centroid_z[n] = 0.0;
+ ah_centroid_t[n] = 0.0;
+ ah_centroid_valid[n] = 0;
+ ah_centroid_iteration[n] = -1;
+
+ ah_centroid_x_p[n] = 0.0;
+ ah_centroid_y_p[n] = 0.0;
+ ah_centroid_z_p[n] = 0.0;
+ ah_centroid_t_p[n] = 0.0;
+ ah_centroid_valid_p[n] = 0;
+ ah_centroid_iteration_p[n] = -1;
+
+ struct AH_data& AH_data = *state.AH_data_array[n+1];
+ ah_initial_find_flag[n] = AH_data.initial_find_flag;
+ ah_really_initial_find_flag[n] = AH_data.really_initial_find_flag;
+ ah_search_flag[n] = AH_data.search_flag;
+ ah_found_flag[n] = AH_data.found_flag;
+ if (verbose_info.print_algorithm_highlights) {
+ printf ("AHF setup %d initial_find_flag=%d\n", n+1, (int) AH_data.initial_find_flag);
+ printf ("AHF setup %d really_initial_find_flag=%d\n", n+1, (int) AH_data.really_initial_find_flag);
+ printf ("AHF setup %d search_flag=%d\n", n+1, (int) AH_data.search_flag);
+ printf ("AHF setup %d found_flag=%d\n", n+1, (int) AH_data.found_flag);
+ }
+ }
}
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
-// This function sets our internal data structures from the the mask
-// parameters.
-//
-namespace {
-void set_mask_pars(CCTK_ARGUMENTS)
-{
-DECLARE_CCTK_ARGUMENTS
-DECLARE_CCTK_PARAMETERS
-
-struct mask_info& mask_info = state.mask_info;
-mask_info.set_old_style_mask = (set_old_style_mask != 0);
-mask_info.set_new_style_mask = (set_new_style_mask != 0);
-if (mask_info.set_mask_for_any_horizon)
- then {
- mask_info.radius_multiplier = mask_radius_multiplier;
- mask_info.radius_offset = mask_radius_offset;
- mask_info.buffer_thickness = mask_buffer_thickness;
- mask_info.mask_is_noshrink = mask_is_noshrink;
- mask_info.min_horizon_radius_points_for_mask
- = min_horizon_radius_points_for_mask;
- mask_info.set_old_style_mask = (set_old_style_mask != 0);
- mask_info.set_new_style_mask = (set_new_style_mask != 0);
- if (mask_info.set_old_style_mask)
- then {
- struct mask_info::old_style_mask_info& osmi
- = mask_info.old_style_mask_info;
- osmi.gridfn_name = old_style_mask_gridfn_name;
- osmi.gridfn_varindex = Cactus_gridfn_varindex(osmi.gridfn_name);
- osmi.gridfn_dataptr = NULL; // dummy value; fixup later
- osmi.inside_value = old_style_mask_inside_value;
- osmi.buffer_value = old_style_mask_buffer_value;
- osmi.outside_value = old_style_mask_outside_value;
- }
- if (mask_info.set_new_style_mask)
- then {
- struct mask_info::new_style_mask_info& nsmi
- = mask_info.new_style_mask_info;
- nsmi.gridfn_name = new_style_mask_gridfn_name;
- nsmi.gridfn_varindex = Cactus_gridfn_varindex(nsmi.gridfn_name);
- nsmi.gridfn_dataptr = NULL; // dummy value; fixup later
- nsmi.bitfield_name = new_style_mask_bitfield_name;
- nsmi.bitfield_bitmask = 0; // dummy value; fixup later
- nsmi.inside_value = new_style_mask_inside_value;
- nsmi.buffer_value = new_style_mask_buffer_value;
- nsmi.outside_value = new_style_mask_outside_value;
- nsmi.inside_bitvalue = 0; // dummy value; fixup later
- nsmi.buffer_bitvalue = 0; // dummy value; fixup later
- nsmi.outside_bitvalue = 0; // dummy value; fixup later
- }
- }
-}
- }
//******************************************************************************
//******************************************************************************
@@ -677,9 +720,7 @@ namespace {
enum verbose_level
decode_verbose_level(const char verbose_level_string[])
{
-if (STRING_EQUAL(verbose_level_string, "no output"))
- then return verbose_level__no_output;
-else if (STRING_EQUAL(verbose_level_string, "physics highlights"))
+if (STRING_EQUAL(verbose_level_string, "physics highlights"))
then return verbose_level__physics_highlights;
else if (STRING_EQUAL(verbose_level_string, "physics details"))
then return verbose_level__physics_details;
@@ -696,26 +737,6 @@ else CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
}
//******************************************************************************
-
-//
-// This function decodes the horizon_file_format parameter (string) into an
-// internal enum for future use.
-//
-namespace {
-enum horizon_file_format
- decode_horizon_file_format(const char horizon_file_format_string[])
-{
-if (STRING_EQUAL(horizon_file_format_string, "ASCII (gnuplot)"))
- then return horizon_file_format__ASCII_gnuplot;
-else if (STRING_EQUAL(horizon_file_format_string, "HDF5"))
- then return horizon_file_format__HDF5;
-else CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"decode_horizon_file_format(): unknown horizon_file_format_string=\"%s\"!",
- horizon_file_format_string); /*NOTREACHED*/
-}
- }
-
-//******************************************************************************
//******************************************************************************
//******************************************************************************
@@ -743,6 +764,7 @@ else CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
namespace {
int allocate_horizons_to_processor(int N_procs, int my_proc,
int N_horizons, bool multiproc_flag,
+ const CCTK_INT depends_on[],
horizon_sequence& my_hs,
const struct verbose_info& verbose_info)
{
@@ -753,15 +775,38 @@ const int N_active_procs = multiproc_flag ? jtutil::min(N_procs, N_horizons)
// Implementation note:
// We allocate the horizons to active processors in round-robin order.
//
+std::vector<int> proc_of_horizon (N_horizons+1);
+ for (int hn = 1 ; hn <= N_horizons ; ++hn)
+ {
+ proc_of_horizon.at(hn) = -1;
+ }
+
int proc = 0;
for (int hn = 1 ; hn <= N_horizons ; ++hn)
{
+ if (depends_on[hn] < 0 || depends_on[hn] > N_horizons) {
+ CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "horizon %d depends on a horizon with the illegal index %d",
+ hn, int(depends_on[hn]));
+ } else if (depends_on[hn] == hn) {
+ CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "horizon %d depends on itself",
+ hn);
+ } else if (depends_on[hn] > hn) {
+ CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "horizon %d depends on a horizon with a larger index %d",
+ hn, int(depends_on[hn]));
+ }
+ const int this_horizons_proc
+ = depends_on[hn] == 0 ? proc : proc_of_horizon.at(depends_on[hn]);
+ assert (this_horizons_proc >= 0 && this_horizons_proc < N_procs);
+ proc_of_horizon.at(hn) = this_horizons_proc;
if (verbose_info.print_algorithm_highlights)
then CCTK_VInfo(CCTK_THORNSTRING,
" allocating horizon %d to processor #%d",
- hn, proc);
- if (proc == my_proc)
- then my_hs.append_hn(hn);
+ hn, this_horizons_proc);
+ if (this_horizons_proc == my_proc)
+ then my_hs.append_hn(hn);
if (++proc >= N_active_procs)
then proc = 0;
}
@@ -875,6 +920,19 @@ else if (STRING_EQUAL(grid_domain, "octant"))
" ==> using \"+xz quadrant (mirrored)\" patch system");
return patch_system::patch_system__plus_xz_quadrant_mirrored;
}
+ else if ((origin_y == 0.0) && (origin_z == 0.0))
+ then {
+ CCTK_VInfo(CCTK_THORNSTRING,
+ " Cactus has mirrored octant mode");
+ CCTK_VInfo(CCTK_THORNSTRING,
+ " but patch system origin_x=%g != 0",
+ double(origin_z));
+ CCTK_VInfo(CCTK_THORNSTRING,
+ " ==> using \"+yz quadrant (mirrored)\" patch system");
+ //return patch_system::patch_system__plus_yz_quadrant_mirrored;
+ CCTK_WARN (0, "This patch system type is not implemented");
+ return patch_system::patch_system__full_sphere;
+ }
else {
CCTK_VInfo(CCTK_THORNSTRING,
" Cactus has mirrored octant mode");
@@ -888,7 +946,7 @@ else if (STRING_EQUAL(grid_domain, "octant"))
}
else {
- CCTK_VWarn(SERIOUS_WARNING, __LINE__, __FILE__, CCTK_THORNSTRING,
+ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
" choose_patch_system_type()\n"
" grid::domain = \"%s\" not supported!\n"