Bugfix/integral symmetry factor

Examples/BinaryBH/params.txt

@@ -108,6 +108,11 @@ vars_parity            = 0 0 4 6 0 5 0    #chi and hij
 vars_parity_diagnostic = 0 1 2 3          #Ham and Mom
                          0 7              #Weyl
 
+# Additionally, if reflective boundaries selected, 
+# choose if symmetry factor will be applied to integrals.
+# 0 = false (default), 1 = true
+symmetry_correction = 1
+
 # if sommerfeld boundaries selected, must select
 # asymptotic values
 num_nonzero_asymptotic_vars = 5

Examples/BinaryBH/params_two_punctures.txt

@@ -108,6 +108,11 @@ vars_parity            = 0 0 4 6 0 5 0    #chi and hij
 vars_parity_diagnostic = 0 1 2 3          #Ham and Mom
                          0 7              #Weyl
 
+# Additionally, if reflective boundaries selected, 
+# choose if symmetry factor will be applied to integrals.
+# 0 = false (default), 1 = true
+symmetry_correction = 1
+
 # if sommerfeld boundaries selected, must select
 # asymptotic values
 num_nonzero_asymptotic_vars = 5

Examples/BinaryBH/params_very_cheap.txt

@@ -9,7 +9,7 @@
 verbosity = 0
 
 # location / naming of output files
-# output_path = "" # Main path for all files. Must exist!
+output_path = "/nfs/st01/hpc-gr-epss/eaf49/PR-dump" # Main path for all files. Must exist!
 chk_prefix = BinaryBHChk_
 plot_prefix = BinaryBHPlot_
 # restart_file = BinaryBHChk_000360.3d.hdf5
@@ -108,6 +108,11 @@ vars_parity            = 0 0 4 6 0 5 0    #chi and hij
 vars_parity_diagnostic = 0 1 2 3          #Ham and Mom
                          0 7              #Weyl
 
+# Additionally, if reflective boundaries selected, 
+# choose if symmetry factor will be applied to integrals.
+# 0 = false (default), 1 = true
+symmetry_correction = 1
+
 # if sommerfeld boundaries selected, must select
 # non zero asymptotic values
 num_nonzero_asymptotic_vars = 5
@@ -155,7 +160,7 @@ sigma = 0.3
 track_punctures = 1
 puncture_tracking_level = 1
 
-# calculate_constraint_norms = 0
+calculate_constraint_norms = 1
 
 # min_chi = 1.e-4
 # min_lapse = 1.e-4
@@ -183,14 +188,15 @@ modes = 2 0 # l m for spherical harmonics
 #################################################
 # Apparent Horizon Finder parameters
 
-AH_activate = 1
-AH_num_ranks = 4
-AH_num_points_u = 11
-AH_num_points_v = 10
-AH_solve_interval = 2
-AH_print_interval = 2
+AH_activate = 0
+AH_num_ranks = 65
+AH_num_points_u = 65
+AH_num_points_v = 48
+# AH_solve_interval = 2
+# AH_print_interval = 2
 # AH_merger_search_factor = 1.
 # AH_merger_pre_factor = 1.
+AH_set_origins_to_punctures = 1
 
 #################################################
 

Examples/KerrBH/params.txt

@@ -96,6 +96,11 @@ vars_parity            = 0 0 4 6 0 5 0    #chi and hij
                          0 1 2 3 1 2 3    #lapse shift and B
 vars_parity_diagnostic = 0 1 2 3          #Ham and Mom
 
+# Additionally, if reflective boundaries selected, 
+# choose if symmetry factor will be applied to integrals.
+# 0 = false (default), 1 = true
+symmetry_correction = 1
+
 # if sommerfeld boundaries selected, must select
 # non zero asymptotic values
 num_nonzero_asymptotic_vars = 5

Examples/KerrBH/params_cheap.txt

@@ -97,6 +97,11 @@ vars_parity            = 0 0 4 6 0 5 0    #chi and hij
                          0 1 2 3 1 2 3    #lapse shift and B
 vars_parity_diagnostic = 0 1 2 3          #Ham and Mom
 
+# Additionally, if reflective boundaries selected, 
+# choose if symmetry factor will be applied to integrals.
+# 0 = false (default), 1 = true
+symmetry_correction = 1
+
 # if sommerfeld boundaries selected, must select
 # non zero asymptotic values
 num_nonzero_asymptotic_vars = 5

Examples/ScalarField/params.txt

@@ -114,6 +114,11 @@ vars_parity            = 0 0 4 6 0 5 0    #chi and hij
                          0 0              #phi and Pi
 vars_parity_diagnostic = 0 1 2 3          #Ham and Mom
 
+# Additionally, if reflective boundaries selected, 
+# choose if symmetry factor will be applied to integrals.
+# 0 = false (default), 1 = true
+symmetry_correction = 1
+
 # if sommerfeld boundaries selected, must select
 # non zero asymptotic values
 num_nonzero_asymptotic_vars = 5

Source/BoxUtils/AMRReductions.hpp

@@ -33,6 +33,7 @@ template <VariableType var_t> class AMRReductions
     double m_domain_volume;
     Vector<LevelData<FArrayBox> *> m_level_data_ptrs;
     Vector<int> m_ref_ratios;
+    const int m_symmetry_factor;
 
     //! constructs a Vector of LevelData<FArrayBox> pointers and stores them
     void set_level_data_vect(const GRAMR &a_gramr);

Source/BoxUtils/AMRReductions.impl.hpp

@@ -14,7 +14,8 @@ template <VariableType var_t>
 AMRReductions<var_t>::AMRReductions(const GRAMR &a_gramr,
                                     const int a_base_level)
     : m_base_level(a_base_level),
-      m_coarsest_dx(a_gramr.get_gramrlevels()[0]->get_dx())
+      m_coarsest_dx(a_gramr.get_gramrlevels()[0]->get_dx()),
+      m_symmetry_factor(a_gramr.get_gramrlevels()[0]->get_symm())
 {
     set_level_data_vect(a_gramr);
     set_ref_ratios_vect(a_gramr);
@@ -109,7 +110,7 @@ Real AMRReductions<var_t>::norm(const Interval &a_vars,
             pow(m_domain_volume, 1.0 / static_cast<double>(a_norm_exponent));
     }
 
-    return norm;
+    return norm * static_cast<double>(m_symmetry_factor);
 }
 
 template <VariableType var_t>
@@ -125,7 +126,8 @@ Real AMRReductions<var_t>::sum(const Interval &a_vars) const
     CH_assert(a_vars.begin() >= 0 && a_vars.end() < m_num_vars);
     CH_TIME("AMRReductions::sum");
     return computeSum(m_level_data_ptrs, m_ref_ratios, m_coarsest_dx, a_vars,
-                      m_base_level);
+                      m_base_level) *
+           static_cast<double>(m_symmetry_factor);
 }
 
 template <VariableType var_t>

Source/GRChomboCore/BoundaryConditions.cpp

@@ -34,6 +34,7 @@ BoundaryConditions::params_t::params_t()
     vars_parity_diagnostic.fill(BoundaryConditions::EXTRAPOLATING_BC);
     vars_asymptotic_values.fill(0.0);
     extrapolation_order = 1;
+    symmetry_factor = 1;
 }
 
 void BoundaryConditions::params_t::set_is_periodic(
@@ -141,6 +142,23 @@ void BoundaryConditions::params_t::read_params(GRParmParse &pp)
 
     if (reflective_boundaries_exist)
     {
+        int symmetry_correction = 0;
+        pp.load("symmetry_correction", symmetry_correction, 0);
+        if (symmetry_correction != 0)
+        {
+            for (int d = 0; d < CH_SPACEDIM; d++)
+            {
+                if (lo_boundary[d] == 2)
+                {
+                    symmetry_factor *= 2;
+                }
+                else if (hi_boundary[d] == 2)
+                {
+                    symmetry_factor *= 2;
+                }
+            }
+        }
+
         pp.load("vars_parity", vars_parity);
         if (pp.contains("vars_parity_diagnostic"))
             pp.load("vars_parity_diagnostic", vars_parity_diagnostic);
@@ -215,6 +233,7 @@ void BoundaryConditions::define(double a_dx,
 {
     m_dx = a_dx;
     m_params = a_params;
+    m_symmetry_factor = a_params.symmetry_factor;
     m_domain = a_domain;
     m_domain_box = a_domain.domainBox();
     m_num_ghosts = a_num_ghosts;

Source/GRChomboCore/BoundaryConditions.hpp

@@ -73,6 +73,7 @@ class BoundaryConditions
         bool sommerfeld_boundaries_exist;
         bool extrapolating_boundaries_exist;
         bool mixed_boundaries_exist;
+        int symmetry_factor;
 
         std::array<int, NUM_VARS> vars_parity;
         std::array<int, NUM_DIAGNOSTIC_VARS>
@@ -97,11 +98,14 @@ class BoundaryConditions
     ProblemDomain m_domain; // the problem domain (excludes boundary cells)
     Box m_domain_box;       // The box representing the domain
     bool is_defined; // whether the BoundaryConditions class members are defined
+    int m_symmetry_factor;
 
   public:
     /// Default constructor - need to call define afterwards
     BoundaryConditions() { is_defined = false; }
 
+    int get_symm() const { return m_symmetry_factor; }
+
     /// define function sets members and is_defined set to true
     void define(double a_dx, std::array<double, CH_SPACEDIM> a_center,
                 const params_t &a_params, ProblemDomain a_domain,

Source/GRChomboCore/GRAMRLevel.hpp

@@ -177,6 +177,10 @@ class GRAMRLevel : public AMRLevel, public InterpSource<>
 
     // direction irrelevant, but relevant for InterpSource
     ALWAYS_INLINE double get_dx(int dir = 0) const { return m_dx; };
+    ALWAYS_INLINE int get_symm(int dir = 0) const
+    {
+        return m_boundaries.get_symm();
+    }
 
     /// Returns true if m_time is the same as the time at the end of the current
     /// timestep on level a_level and false otherwise

Tests/AMRInterpolatorTest/AMRInterpolatorTest.inputs

@@ -26,3 +26,5 @@ lo_boundary = 2 0 2
 # 1,2,3 = odd x, y, z
 # 4,5,6 = odd xy, yz, xz
 vars_parity = 0 1 			# A and B
+
+symmetry_correction = 0

Tests/ApparentHorizonFinderTest2D/ApparentHorizonTest2D.inputs

@@ -18,6 +18,8 @@ lo_boundary = 1 2
 vars_parity            = 2  # V
 vars_asymptotic_values = 0. # V
 
+symmetry_correction = 0
+
 #Max and min box sizes
 max_grid_size = 32
 block_factor = 4

Tests/ApparentHorizonFinderTest3D/ApparentHorizonTest3D.inputs

@@ -23,6 +23,8 @@ vars_asymptotic_values = 1. 1. 0. 0. 1. 0. 1. #chi and hij
                          1. 0. 0. 0. 0. 0. 0. #lapse shift and B
                          0. 0. 0. 0. 		  # Ham and Mom
 
+symmetry_correction = 0
+
 #Max and min box sizes
 max_grid_size = 32
 block_factor = 16

run_clang_format

@@ -1,2 +1,3 @@
 for f in $(find . -name "*.hpp"); do clang-format -style=file -i $f; done
 for f in $(find . -name "*.cpp"); do clang-format -style=file -i $f; done
+