Feature/adaptive moving puncture

InstallNotes/MakeDefsLocalExamples/Cosma-Durham.Make.defs.local

@@ -1,18 +1,19 @@
 DIM           = 3
 CXX           = icpc -std=c++14 -qopenmp -mkl=sequential -traditional -xCORE-AVX512
 FC            = ifort -qopenmp -mkl=sequential
-OPT           = TRUE
+OPT           = HIGH
 DEBUG         = FALSE
 MPI           = TRUE
-OMP           = TRUE
-OPT           = HIGH
+OPENMPCC           = TRUE
 USE_64         = TRUE
 USE_HDF        = TRUE
 MPICXX        = mpiicpc -std=c++14 -qopenmp -mkl=sequential -xCORE-AVX512
 HDFINCFLAGS    = -I${HDF5_HOME}/include
 HDFLIBFLAGS    = -L${HDF5_HOME}/lib -lhdf5 -lz
 HDFMPIINCFLAGS = -I${HDF5_HOME}/include
 HDFMPILIBFLAGS = -L${HDF5_HOME}/lib -lhdf5 -lz
-cxxdbgflags      = -g -Wl,--eh-frame-hdr
-cxxoptflags      = -O3 -qoverride-limits
-fdbgflags        = -g -Wl,--eh-frame-hdr
+cxxdbgflags      = -g
+cxxoptflags      = -O3
+fdbgflags        = -g
+
+

InstallNotes/MakeDefsLocalExamples/Cosma-Durham.modules.sh

@@ -1,4 +1,7 @@
-module load intel_comp/2019 intel_mpi/2019 parallel_hdf5/1.10.3
-module load pythonconda3/4.5.4
-module unload python
+module load intel_comp/2019 
+module load intel_mpi/2019 
+module load parallel_hdf5/1.10.3
+
+
+
 

Source/CCZ4/CCZ4RHS.hpp

@@ -18,6 +18,7 @@
 #include "UserVariables.hpp" //This files needs NUM_VARS - total number of components
 
 #include <array>
+#include <type_traits>
 
 /// Base parameter struct for CCZ4
 /** This struct collects the gauge independent CCZ4 parameters i.e. the damping
@@ -78,15 +79,33 @@ class CCZ4RHS
     const deriv_t m_deriv;
 
   public:
-    /// Constructor
+    /// Constructor that allows manual construction of the gauge class
     CCZ4RHS(
         params_t a_params,            //!< The CCZ4 parameters
+        const gauge_t &a_gauge,       //!< The gauge class 
         double a_dx,                  //!< The grid spacing
         double a_sigma,               //!< Kreiss-Oliger dissipation coefficient
         int a_formulation = USE_CCZ4, //!< Switches between CCZ4, BSSN,...
         double a_cosmological_constant = 0 //!< Value of the cosmological const.
     );
 
+    /// Simple constructor when gauge class can be constructed with just a
+    /// params struct. The dummy template parameter and final argument are just
+    /// to remove this constructor from overload resolution when the gauge class
+    /// does not support this type of construction.
+    // MR: There are much nicer ways to do this using C++20
+    template <typename dummy_t = gauge_t>
+    CCZ4RHS(params_t a_params, //!< The CCZ4 parameters
+            double a_dx,       //!< The grid spacing
+            double a_sigma,    //!< Kreiss-Oliger dissipation coefficient
+            int a_formulation = USE_CCZ4, //!< Switches between CCZ4, BSSN,...
+            double a_cosmological_constant =
+                0, //!< Value of the cosmological const.
+            typename std::enable_if_t<
+                std::is_constructible<dummy_t, params_t>::value> * =
+                0 //!< dummy argument for SFINAE
+    );
+
     /// Compute function
     /** This function orchestrates the calculation of the rhs for one specific
      * grid cell. This function is called by the BoxLoops::loop for each grid
@@ -112,8 +131,10 @@ class CCZ4RHS
         const diff2_vars_t<Tensor<2, data_t>>
             &d2, //!< The second derivative the variables
         const vars_t<data_t>
-            &advec //!< The advection derivatives of the variables
-    ) const;
+            &advec, //!< The advection derivatives of the variables
+        const Cell<data_t> &current_cell //!< Cell class that can be manipulated in
+                                         //!< in the RHS and gauge class.
+        ) const;
 };
 
 #include "CCZ4RHS.impl.hpp"

Source/CCZ4/CCZ4RHS.impl.hpp

@@ -16,9 +16,11 @@
 
 template <class gauge_t, class deriv_t>
 inline CCZ4RHS<gauge_t, deriv_t>::CCZ4RHS(
-    CCZ4_params_t<typename gauge_t::params_t> a_params, double a_dx,
-    double a_sigma, int a_formulation, double a_cosmological_constant)
-    : m_params(a_params), m_gauge(a_params), m_sigma(a_sigma),
+    CCZ4_params_t<typename gauge_t::params_t> a_params,
+    const gauge_t &a_gauge,
+    double a_dx, double a_sigma, int a_formulation,
+    double a_cosmological_constant)
+    : m_params(a_params), m_gauge(a_gauge), m_sigma(a_sigma),
       m_formulation(a_formulation),
       m_cosmological_constant(a_cosmological_constant), m_deriv(a_dx)
 {
@@ -36,6 +38,18 @@ inline CCZ4RHS<gauge_t, deriv_t>::CCZ4RHS(
         MayDay::Error("The requested formulation is not supported");
 }
 
+template <class gauge_t, class deriv_t>
+template <typename dummy_t>
+inline CCZ4RHS<gauge_t, deriv_t>::CCZ4RHS(
+    CCZ4_params_t<typename gauge_t::params_t> a_params, double a_dx,
+    double a_sigma, int a_formulation, double a_cosmological_constant,
+    typename std::enable_if_t<std::is_constructible<dummy_t, params_t>::value>
+        *)
+    : CCZ4RHS(a_params, gauge_t(a_params), a_dx, a_sigma, a_formulation,
+              a_cosmological_constant)
+{
+}
+
 template <class gauge_t, class deriv_t>
 template <class data_t>
 void CCZ4RHS<gauge_t, deriv_t>::compute(Cell<data_t> current_cell) const
@@ -47,7 +61,7 @@ void CCZ4RHS<gauge_t, deriv_t>::compute(Cell<data_t> current_cell) const
         m_deriv.template advection<Vars>(current_cell, vars.shift);
 
     Vars<data_t> rhs;
-    rhs_equation(rhs, vars, d1, d2, advec);
+    rhs_equation(rhs, vars, d1, d2, advec, current_cell);
 
     m_deriv.add_dissipation(rhs, current_cell, m_sigma);
 
@@ -61,7 +75,8 @@ void CCZ4RHS<gauge_t, deriv_t>::rhs_equation(
     vars_t<data_t> &rhs, const vars_t<data_t> &vars,
     const vars_t<Tensor<1, data_t>> &d1,
     const diff2_vars_t<Tensor<2, data_t>> &d2,
-    const vars_t<data_t> &advec) const
+    const vars_t<data_t> &advec,
+    const Cell<data_t> &current_cell) const
 {
     using namespace TensorAlgebra;
 
@@ -223,7 +238,7 @@ void CCZ4RHS<gauge_t, deriv_t>::rhs_equation(
 
     FOR1(i) { rhs.Gamma[i] = advec.Gamma[i] + Gammadot[i]; }
 
-    m_gauge.rhs_gauge(rhs, vars, d1, d2, advec);
+    m_gauge.rhs_gauge(rhs, vars, d1, d2, advec, current_cell);
 }
 
 #endif /* CCZ4RHS_IMPL_HPP_ */

Source/CCZ4/IntegratedMovingPunctureGauge.hpp

@@ -77,7 +77,8 @@ class IntegratedMovingPunctureGauge
     inline void rhs_gauge(vars_t<data_t> &rhs, const vars_t<data_t> &vars,
                           const vars_t<Tensor<1, data_t>> &d1,
                           const diff2_vars_t<Tensor<2, data_t>> &d2,
-                          const vars_t<data_t> &advec) const
+                          const vars_t<data_t> &advec,
+                          const Cell<data_t> &current_cell) const
     {
         rhs.lapse = m_params.lapse_advec_coeff * advec.lapse -
                     m_params.lapse_coeff *

Source/CCZ4/MovingPunctureGauge.hpp

@@ -6,6 +6,8 @@
 #ifndef MOVINGPUNCTUREGAUGE_HPP_
 #define MOVINGPUNCTUREGAUGE_HPP_
 
+
+#include "Cell.hpp"
 #include "DimensionDefinitions.hpp"
 #include "Tensor.hpp"
 
@@ -49,7 +51,8 @@ class MovingPunctureGauge
     inline void rhs_gauge(vars_t<data_t> &rhs, const vars_t<data_t> &vars,
                           const vars_t<Tensor<1, data_t>> &d1,
                           const diff2_vars_t<Tensor<2, data_t>> &d2,
-                          const vars_t<data_t> &advec) const
+                          const vars_t<data_t> &advec,
+                          const Cell<data_t> &current_cell) const
     {
         rhs.lapse = m_params.lapse_advec_coeff * advec.lapse -
                     m_params.lapse_coeff *

Source/GRChomboCore/ChomboParameters.hpp

@@ -55,8 +55,10 @@ class ChomboParameters
         // L's, N's and center
         read_grid_params(pp);
 
+#ifdef CH_USE_HDF5
         pp.load("ignore_checkpoint_name_mismatch",
                 ignore_checkpoint_name_mismatch, false);
+#endif
 
         pp.load("max_level", max_level, 0);
         // the reference ratio is hard coded to 2
@@ -93,12 +95,13 @@ class ChomboParameters
         pp.load("plot_interval", plot_interval, 0);
         pp.load("stop_time", stop_time, 1.0);
         pp.load("max_steps", max_steps, 1000000);
+#ifdef CH_USE_HDF5
         pp.load("write_plot_ghosts", write_plot_ghosts, false);
 
         // load vars to write to plot files
         UserVariables::load_vars_to_vector(pp, "plot_vars", "num_plot_vars",
                                            plot_vars, num_plot_vars);
-
+#endif
         // alias the weird chombo names to something more descriptive
         // for these box params, and default to some reasonable values
         if (pp.contains("max_grid_size"))
@@ -129,15 +132,15 @@ class ChomboParameters
     {
         // In this function, cannot use default value - it may print a 'default
         // message' to pout and a 'setPoutBaseName' must happen before
-
         restart_from_checkpoint = pp.contains("restart_file");
+#ifdef CH_USE_HDF5
         if (restart_from_checkpoint)
         {
             pp.load("restart_file", restart_file);
         }
-
         pp.load("chk_prefix", checkpoint_prefix);
         pp.load("plot_prefix", plot_prefix);
+#endif
 
 #ifdef CH_MPI
         // Again, cannot use default value
@@ -378,13 +381,15 @@ class ChomboParameters
 
         // check the restart_file exists and can be read if restarting from a
         // checkpoint
+#ifdef CH_USE_HDF5
         if (restart_from_checkpoint)
         {
             bool restart_file_exists =
-                (access(restart_file.c_str(), R_OK) == 0);
+                (access((hdf5_path + restart_file).c_str(), R_OK) == 0);
             check_parameter("restart_file", restart_file, restart_file_exists,
                             "file cannot be opened for reading");
         }
+#endif
 
         check_parameter("dt_multiplier", dt_multiplier, dt_multiplier > 0.0,
                         "must be > 0.0");
@@ -419,9 +424,11 @@ class ChomboParameters
         // (MR); while this would technically work (any plot files would just
         // overwrite a checkpoint file), I think a user would only ever do
         // this unintentinally
+#ifdef CH_USE_HDF5
         check_parameter("plot_prefix", plot_prefix,
                         plot_interval <= 0 || plot_prefix != checkpoint_prefix,
                         "should be different to checkpoint_prefix");
+#endif
 
         check_parameter("output_path", output_path,
                         FilesystemTools::directory_exists(output_path),
@@ -474,25 +481,30 @@ class ChomboParameters
     Vector<int> regrid_interval; // steps between regrid at each level
     int max_steps;
     bool restart_from_checkpoint; // whether or not to restart or start afresh
-    std::string restart_file;     // The path to the restart_file
-    bool ignore_checkpoint_name_mismatch;   // ignore mismatch of variable names
-                                            // between restart file and program
+#ifdef CH_USE_HDF5
+    std::string restart_file;             // The path to the restart_file
+    bool ignore_checkpoint_name_mismatch; // ignore mismatch of variable names
+                                          // between restart file and program
+#endif
     double dt_multiplier, stop_time;        // The Courant factor and stop time
     int checkpoint_interval, plot_interval; // Steps between outputs
     int max_grid_size, block_factor;        // max and min box sizes
     double fill_ratio; // determines how fussy the regridding is about tags
-    std::string checkpoint_prefix, plot_prefix, pout_prefix; // naming of files
+#ifdef CH_USE_HDF5
+    std::string checkpoint_prefix, plot_prefix; // naming of files
+#endif
     std::string output_path; // base path to use for all files
 #ifdef CH_MPI
-    std::string pout_path; // base path for pout files
+    std::string pout_prefix; // pout file prefix
+    std::string pout_path;   // base path for pout files
 #endif
 #ifdef CH_USE_HDF5
     std::string hdf5_path; // base path for pout files
-#endif
     bool write_plot_ghosts;
     int num_plot_vars;
     std::vector<std::pair<int, VariableType>>
         plot_vars; // vars to write to plot file
+#endif
 
     std::array<double, CH_SPACEDIM> origin,
         dx; // location of coarsest origin and dx