Bugfix/integral symmetry factor (#256)

* Added symmetry correction flag and symmetry factor variable

1. The symmetry correction flag is read in ChomboParameters. It is
defaulted to "off" if the parameter is not provided.
2. If symmetry correction is flag is on, the symmetry factor is
calculated based off of the number of hi or lo boundaries chosen
to be reflective (choice 2).
3. The symmetry factor is a chombo parameter, and must be given to
any AMRReduction function which integrates over the entire box
(currently norm and sum). The symmetry factor is defaulted to 1
if it is not provided to the AMRReduction scheme.
4. The result of sum or norm is multiplied by the symmetry factor
before it is returned, in the case of a single variable.

The symmetry factor is provided correctly to AMRReductions for
the BinaryBH example ONLY in this commit.

* Added symmetry correction flag to all params in all Examples

* Adding params changes, so I can bugfix separately

* Bug fixes to symmetry factor implementation

Fixed a few bugs with how I had first implemented the
symmetry factor:

1. Moved the symmetry factor argument in norm and sum
ahead of the defaulted variables.
2. Changed type of the symmetry factor variable in
sum and norm to int, and cast it to double before applying
to the final result.
3. Changed error in the name of the variable in norm.

* Trying to add symmetry factor from GRAMR but I think I need to start on a new branch...

* Symmetry factor implementation in GRAMR

Created a series of subroutines to directly read in
the symmetry factor associated with using reflective
boundary conditions in the GRAMR class, which can then be
drawn out for use by the sum() and norm() AMRReductions
functions.

1. The symmetry factor is calculated in BoundaryConditions
as a parameter, then assigned to a private BoundaryConditions
variable.
2. A GRAMRLevel subroutine, get_symm(), calls the BoundaryConditions
subroutine to return the private symmetry factor variable.
3. AMRReductions reads in the symmetry factor from get_gramrlevels
in the same style as m_coarsest_dx.
4. The symmetry factor is applied to the norm() and sum() results
immediately before they are returned.

* Removed symmetry argument from uses of AMRReductions

* Fixed clang tidy warnings

* Fixed clang tidy warnings

* Clang formatting changes, after running clang-format locally

* Fixed typo on line 40 of BoundaryConditions.cpp

* Added default to the symmetry_correction flag

* Symmetry factor implementation

Created a series of subroutines to directly read in
the symmetry factor associated with using reflective
boundary conditions in the GRAMR class, which can then be
drawn out for use by the sum() and norm() AMRReductions
functions.

1. If the symmetry_correction flag in params is turn on,
the symmetry factor is calculated in BoundaryConditions
as a parameter, then assigned to a private BoundaryConditions
variable.
2. A GRAMRLevel subroutine, get_symm(), calls the BoundaryConditions
subroutine to return the private symmetry factor variable.
3. AMRReductions reads in the symmetry factor from get_gramrlevels
in the same style as m_coarsest_dx.
4. The symmetry factor is applied to the norm() and sum() results
immediately before they are returned.

* Pre-pull commit

* Removal of old Chombo Parameter function argument from AMRReductions class.

* Fixed Katy's comments on PR

1. Removed symmetry_correction flag, making the symmetry factor
   calculation automatic.
2. Moved application of symmetry factor in AMRReductions::norm function
   into the logic block controlled by the a_normalize_by_volume flag, so
that only when normalising by the volume will the symmetry factor be
applied.
3. Changed get_symm() function names to get_symmetry_factor() (coding
   convention).
4. Removed file path in params_very_cheap.txt in BinarBH example.

* Remove symmetry_correction param from tests

---------

Co-authored-by: dinatraykova <dina.traykova@aei.mpg.de>

Examples/BinaryBH/params_very_cheap.txt

@@ -155,7 +155,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 +183,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
 
 #################################################
 

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_symmetry_factor())
 {
     set_level_data_vect(a_gramr);
     set_ref_ratios_vect(a_gramr);
@@ -107,6 +108,7 @@ Real AMRReductions<var_t>::norm(const Interval &a_vars,
     {
         norm /=
             pow(m_domain_volume, 1.0 / static_cast<double>(a_norm_exponent));
+        norm *= static_cast<double>(m_symmetry_factor);
     }
 
     return norm;
@@ -125,7 +127,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,18 @@ void BoundaryConditions::params_t::read_params(GRParmParse &pp)
 
     if (reflective_boundaries_exist)
     {
+        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 +228,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_symmetry_factor() 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_symmetry_factor(int dir = 0) const
+    {
+        return m_boundaries.get_symmetry_factor();
+    }
 
     /// 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

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

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
+