Resolve Katy's comments

Examples/Cosmo/CosmoDiagnostics.hpp

@@ -51,7 +51,7 @@ class CosmoDiagnostics // public MatterConstraints<matter_t>
     {
     }
 
-    //! The compute member which calculates the constraints at each point in the
+    //! The compute member which calculates the diagnostics at each point in the
     //! box
     template <class data_t> void compute(Cell<data_t> current_cell) const;
 

Examples/Cosmo/CosmoDiagnostics.impl.hpp

@@ -43,14 +43,12 @@ void CosmoDiagnostics<matter_t>::compute(Cell<data_t> current_cell) const
     rho_scaled = emtensor.rho / pow(vars.chi, 3. / 2.);
     S_scaled = emtensor.S / pow(vars.chi, 3. / 2.);
 
-    // Write the constraints into the output FArrayBox
+    // Write the diagnostics into the output vector of cells
     current_cell.store_vars(sqrt_gamma, c_sqrt_gamma);
     current_cell.store_vars(rho, c_rho);
     current_cell.store_vars(rho_scaled, c_rho_scaled);
     current_cell.store_vars(S_scaled, c_S_scaled);
     current_cell.store_vars(K_scaled, c_K_scaled);
-
-    // store_vars(out, current_cell);
 }
 
 #endif /* COSMODIAGNOSTICS_IMPL_HPP_ */

Examples/Cosmo/CosmoLevel.cpp

@@ -69,6 +69,8 @@ void CosmoLevel::initialData()
         m_state_new, m_state_new, INCLUDE_GHOST_CELLS);
 
     fillAllGhosts();
+    // Note that the GammaCaluculator is not necessary since the data is
+    // conformally flat. It is left here for generality.
     BoxLoops::loop(GammaCalculator(m_dx), m_state_new, m_state_new,
                    EXCLUDE_GHOST_CELLS);
 
@@ -78,6 +80,7 @@ void CosmoLevel::initialData()
     InitialK<ScalarFieldWithPotential> my_initial_K(scalar_field, m_dx,
                                                     m_p.G_Newton);
     BoxLoops::loop(my_initial_K, m_state_new, m_state_new, EXCLUDE_GHOST_CELLS);
+
     // Calculate constraints and some diagnostics as we need it in tagging
     // criterion
     BoxLoops::loop(MatterConstraints<ScalarFieldWithPotential>(
@@ -89,14 +92,11 @@ void CosmoLevel::initialData()
         scalar_field, m_dx, m_p.G_Newton);
     BoxLoops::loop(cosmo_diagnostics, m_state_new, m_state_diagnostics,
                    EXCLUDE_GHOST_CELLS);
+
     // Assign initial rho_mean here
-    // from rho = 1/2 m^2 phi^2;
-    // phi0 = A sin(2 pi n x /L);
-    // mean of phi0^2 = A^2 sin^2 = 0.5*A^2;
-    // m = m_mode
-    m_cosmo_amr.set_rho_mean(
-        0.5 * m_p.scalar_field_mode * m_p.scalar_field_mode *
-        (0.5 * m_p.initial_params.amplitude * m_p.initial_params.amplitude));
+    InitialScalarData initial_scalar_data(m_p.initial_params, m_dx, m_p.L,
+                                          m_p.scalar_field_mode);
+    m_cosmo_amr.set_rho_mean(initial_scalar_data.compute_initial_rho_mean());
 }
 
 #ifdef CH_USE_HDF5
@@ -132,22 +132,11 @@ void CosmoLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs,
     ScalarFieldWithPotential scalar_field(potential);
     CosmoMovingPunctureGauge cosmo_moving_puncture_gauge(m_p.ccz4_params);
     cosmo_moving_puncture_gauge.set_K_mean(m_cosmo_amr.get_K_mean());
-    if (m_p.max_spatial_derivative_order == 4)
-    {
-        MatterCCZ4RHS<ScalarFieldWithPotential, CosmoMovingPunctureGauge,
-                      FourthOrderDerivatives>
-            my_ccz4_matter(scalar_field, m_p.ccz4_params, m_dx, m_p.sigma,
-                           m_p.formulation, m_p.G_Newton);
-        BoxLoops::loop(my_ccz4_matter, a_soln, a_rhs, EXCLUDE_GHOST_CELLS);
-    }
-    else if (m_p.max_spatial_derivative_order == 6)
-    {
-        MatterCCZ4RHS<ScalarFieldWithPotential, CosmoMovingPunctureGauge,
-                      SixthOrderDerivatives>
-            my_ccz4_matter(scalar_field, m_p.ccz4_params, m_dx, m_p.sigma,
-                           m_p.formulation, m_p.G_Newton);
-        BoxLoops::loop(my_ccz4_matter, a_soln, a_rhs, EXCLUDE_GHOST_CELLS);
-    }
+    MatterCCZ4RHS<ScalarFieldWithPotential, CosmoMovingPunctureGauge,
+                  FourthOrderDerivatives>
+        my_ccz4_matter(scalar_field, m_p.ccz4_params, m_dx, m_p.sigma,
+                       m_p.formulation, m_p.G_Newton);
+    BoxLoops::loop(my_ccz4_matter, a_soln, a_rhs, EXCLUDE_GHOST_CELLS);
 }
 
 // Things to do at ODE update, after soln + rhs
@@ -179,8 +168,9 @@ void CosmoLevel::computeTaggingCriterion(
     FArrayBox &tagging_criterion, const FArrayBox &current_state,
     const FArrayBox &current_state_diagnostics)
 {
-    BoxLoops::loop(CosmoHamTaggingCriterion(m_dx, m_p.center_tag, m_p.rad,
-                                            m_cosmo_amr.get_rho_mean()),
+    double rho_mean = m_cosmo_amr.get_rho_mean();
+    BoxLoops::loop(CosmoHamTaggingCriterion(m_dx, m_p.tagging_center,
+                                            m_p.tagging_radius, rho_mean),
                    current_state_diagnostics, tagging_criterion);
 }
 void CosmoLevel::specificPostTimeStep()
@@ -211,23 +201,23 @@ void CosmoLevel::specificPostTimeStep()
         if (m_level == min_level)
         {
             // AMRReductions for diagnostic variables
-            AMRReductions<VariableType::diagnostic> amr_reductions_diag(
+            AMRReductions<VariableType::diagnostic> amr_reductions_diagnostic(
                 m_cosmo_amr);
-            double phys_vol = amr_reductions_diag.sum(c_sqrt_gamma);
-            double L2_Ham = amr_reductions_diag.norm(c_Ham);
-            double L2_Mom = amr_reductions_diag.norm(c_Mom);
-            double K_total = amr_reductions_diag.sum(c_K_scaled);
-            m_cosmo_amr.set_rho_mean(amr_reductions_diag.sum(c_rho_scaled) /
-                                     phys_vol);
-            m_cosmo_amr.set_S_mean(amr_reductions_diag.sum(c_S_scaled) /
+            double phys_vol = amr_reductions_diagnostic.sum(c_sqrt_gamma);
+            double L2_Ham = amr_reductions_diagnostic.norm(c_Ham);
+            double L2_Mom = amr_reductions_diagnostic.norm(c_Mom);
+            double K_total = amr_reductions_diagnostic.sum(c_K_scaled);
+            m_cosmo_amr.set_rho_mean(
+                amr_reductions_diagnostic.sum(c_rho_scaled) / phys_vol);
+            m_cosmo_amr.set_S_mean(amr_reductions_diagnostic.sum(c_S_scaled) /
                                    phys_vol);
             m_cosmo_amr.set_K_mean(K_total / phys_vol);
 
             // AMRReductions for evolution variables
-            AMRReductions<VariableType::evolution> amr_reductions_evo(
+            AMRReductions<VariableType::evolution> amr_reductions_evolution(
                 m_cosmo_amr);
 
-            double chi_mean = amr_reductions_evo.sum(c_chi) / phys_vol;
+            double chi_mean = amr_reductions_evolution.sum(c_chi) / phys_vol;
 
             // Write output file
             SmallDataIO constraints_file(m_p.data_path + "data_out", m_dt,
@@ -256,11 +246,11 @@ void CosmoLevel::specificPostTimeStep()
             interpolator.refresh();
 
             // set up the query and execute it
-            std::array<double, CH_SPACEDIM> extr_point = {
+            std::array<double, CH_SPACEDIM> extraction_points = {
                 0., m_p.L / 2, m_p.L / 2}; // specified point {x \in [0,L],y \in
                                            // [0,L], z \in [0,L]}
             CustomExtraction extraction(c_rho, m_p.lineout_num_points, m_p.L,
-                                        extr_point, m_dt, m_time);
+                                        extraction_points, m_dt, m_time);
             extraction.execute_query(&interpolator, m_p.data_path + "lineout");
         }
     }

Examples/Cosmo/CosmoLevel.hpp

@@ -3,25 +3,24 @@
  * Please refer to LICENSE in GRChombo's root directory.
  */
 
-#ifndef INFLATIONLEVEL_HPP_
-#define INFLATIONLEVEL_HPP_
+#ifndef COSMOLEVEL_HPP_
+#define COSMOLEVEL_HPP_
 
-// #include "BHAMR.hpp"
 #include "CosmoAMR.hpp"
 #include "DefaultLevelFactory.hpp"
 #include "GRAMRLevel.hpp"
 // Problem specific includes
 #include "Potential.hpp"
 #include "ScalarField.hpp"
 
-//!  A class for the evolution of a scalar field, minimally coupled to gravity
+//!  A class for the evolution of a scalar field in cosmological spacetime
 /*!
-     The class takes some initial data for a scalar field (variables phi and Pi)
-     and evolves it using the CCZ4 equations. It is possible to specify an
-   initial period of relaxation for the conformal factor chi, for non analytic
-   initial conditions (for example, a general field configuration at a moment of
-   time symmetry assuming conformal flatness). \sa MatterCCZ4(),
-   ConstraintsMatter(), ScalarField(), RelaxationChi()
+     The class takes initial data for a scalar field (variables phi and Pi)
+     and evolves it using the CCZ4 equations. An example of specifying initial
+   scalar field data is provided as an example. The given initial data example
+   is obtained by analytically solving the constraint using a sinusoidal scalar
+   field profile with a quadratic potential. See initial scalar field details in
+   InitialData(), InitialScalarData.hpp and Potential.hpp
 */
 class CosmoLevel : public GRAMRLevel
 {

Examples/Cosmo/DiagnosticVariables.hpp

@@ -9,15 +9,16 @@
 // assign an enum to each variable
 enum
 {
-    c_Ham,
-    c_Mom,
-    c_Ham_abs_sum, // c_Ham_abs_sum
-    c_Mom_abs_sum, // c_Mom_abs_sum
-    c_rho,
-    c_sqrt_gamma, // sqrt(gamma) = pow(chi,-3/2) volume factor of spatial metric
-    c_rho_scaled,
-    c_S_scaled,
-    c_K_scaled,
+    c_Ham,         // Hamiltonian constraint
+    c_Mom,         // Momentum constraints
+    c_Ham_abs_sum, // Sum of absolute value of each term in Ham
+    c_Mom_abs_sum, // Sum of absolute value of each term in Mom
+    c_rho,         // Energy density
+    c_sqrt_gamma,  // \sqrt(\gamma) = pow(chi,-3/2) volume factor of spatial
+                   // metric
+    c_rho_scaled,  // \rho * volume factor
+    c_S_scaled,    // S * volume factor
+    c_K_scaled,    // K * volume factor
 
     NUM_DIAGNOSTIC_VARS
 };

Examples/Cosmo/InitialScalarData.hpp

@@ -26,12 +26,11 @@ class InitialScalarData
     {
         double amplitude; //!< Amplitude of bump in initial SF bubble
         std::array<double, CH_SPACEDIM>
-            center;   //!< Centre of perturbation in initial SF bubble
-        double width; //!< Width of bump in initial SF bubble
+            center; //!< Centre of perturbation in initial SF bubble
     };
 
     //! The constructor
-    InitialScalarData(params_t a_params, double a_dx, double a_L, double a_mode)
+    InitialScalarData(params_t a_params, double a_dx, double a_L, int a_mode)
         : m_dx(a_dx), m_L(a_L), m_mode(a_mode), m_params(a_params)
     {
     }
@@ -44,7 +43,8 @@ class InitialScalarData
 
         // From rho = 1/2 (dphi/dx)^2 + V(phi) with V(phi) = 1/2 m^2 phi^2
         // ,we choose phi = A sin(2 n pi x/L) and m = 2 n pi/L such that initial
-        // rho = constant Calculate the field value
+        // rho = constant
+        // Calculate the field value
         data_t phi =
             m_params.amplitude * sin(2 * m_mode * M_PI * coords.x / m_L);
 
@@ -59,11 +59,23 @@ class InitialScalarData
         current_cell.store_vars(1.0, c_chi);
     }
 
+    // Function to analytically compute initial rho mean
+    double compute_initial_rho_mean()
+    {
+        // Compute initial rho_mean here:
+        // from rho = 1/2 m^2 phi^2
+        // phi0 = A sin(2 pi n x /L)
+        // mean of phi0^2 = A^2 sin^2 = 0.5*A^2
+        double rho_mean = 0.5 * m_mode * m_mode *
+                          (0.5 * m_params.amplitude * m_params.amplitude);
+        return rho_mean;
+    }
+
   protected:
     double m_dx;
     const params_t m_params; //!< The matter initial condition params
     double m_L;              // box length
-    double m_mode;           // SF mode
+    int m_mode;              // SF mode
 };
 
 #endif /* INITIALSCALARDATA_HPP_ */

Examples/Cosmo/Potential.hpp

@@ -18,12 +18,12 @@ class Potential
 
   private:
     params_t m_params;
-    double m_mode;
+    int m_mode;
     double m_L;
 
   public:
     //! The constructor
-    Potential(params_t a_params, double a_L, double a_mode)
+    Potential(params_t a_params, double a_L, int a_mode)
         : m_params(a_params), m_L(a_L), m_mode(a_mode)
     {
     }

Examples/Cosmo/SimulationParameters.hpp

@@ -32,16 +32,15 @@ class SimulationParameters : public SimulationParametersBase
             center; // already read in SimulationParametersBase
         pp.load("G_Newton", G_Newton, 1.0);
         pp.load("scalar_amplitude", initial_params.amplitude, 0.1);
-        pp.load("scalar_width", initial_params.width, 1.0);
         pp.load("scalar_mass", potential_params.scalar_mass, 0.1);
-        pp.load("scalar_field_mode", scalar_field_mode, 1.0);
+        pp.load("scalar_field_mode", scalar_field_mode, 1);
 
         // Lineout params
         pp.load("lineout_num_points", lineout_num_points, 10);
 
         // Tagging params
-        pp.load("center_tag", center_tag, center);
-        pp.load("rad", rad, L);
+        pp.load("tagging_center", tagging_center, center);
+        pp.load("tagging_radius", tagging_radius, L);
 
 #ifdef USE_AHFINDER
         double AH_guess =
@@ -57,16 +56,12 @@ class SimulationParameters : public SimulationParametersBase
                            0.2 / coarsest_dx / dt_multiplier,
                        "oscillations of scalar field do not appear to be "
                        "resolved on coarsest level");
-        warn_parameter("scalar_width", initial_params.width,
-                       initial_params.width < 0.5 * L,
-                       "is greater than half the domain size");
     }
 
     // Initial data for matter and potential and BH
-    double G_Newton;
-    double scalar_field_mode, rad;
-    int lineout_num_points;
-    std::array<double, CH_SPACEDIM> center_tag;
+    double G_Newton, tagging_radius;
+    int lineout_num_points, scalar_field_mode;
+    std::array<double, CH_SPACEDIM> tagging_center;
     InitialScalarData::params_t initial_params;
     Potential::params_t potential_params;
     KerrBH::params_t kerr_params;

Examples/Cosmo/params.txt

@@ -44,14 +44,18 @@ print_progress_only_to_rank_0 = 1
 G_Newton = 1.
 
 # Scalar field initial data
+# This given scalar field params makes Hubble length 10 times wavelength of scalar field 
 scalar_amplitude = 0.0777635
-scalar_width = 1.0 # not use
 scalar_mass = 1.0
-scalar_field_mode = 1.0
+scalar_field_mode = 1 # must be an integer
 
-# lineout params
+# Lineout params
 lineout_num_points = 10
 
+# Tagging params
+# tagging_center = 0 0 0
+# tagging_radius = 7
+
 #################################################
 # Grid parameters
 
@@ -66,23 +70,21 @@ N_full = 128
 L_full = 1
 
 # Maximum number of times you can regrid above coarsest level
-max_level = 0 # There are (max_level+1) grids, so min is zero
+max_level = 1 # There are (max_level+1) grids, so min is zero
 
 # Frequency of regridding at each level and thresholds on the tagging
 # Need one for each level except the top one, ie max_level items
 # Generally you do not need to regrid frequently on every level
 # in this example turn off regridding on all levels
 # Level Regridding: 0   1   2   3   4   5
 regrid_interval   = 0   0   0   0   0   0
-# regrid_threshold  = 0.5
-
-regrid_interval = 0 0 0 0 0 0
+regrid_threshold  = 1e6
 
 # Max and min box sizes
-max_box_size = 8
-min_box_size = 8
+max_box_size = 16
+min_box_size = 16
 
-tag_buffer_size = 0 # this example uses a fixed grid
+# tag_buffer_size = 0
 
 # grid_buffer_size = 8
 # fill_ratio = 0.7
@@ -134,7 +136,8 @@ stop_time = 100.
 # max_steps = 4
 
 # Spatial derivative order (only affects CCZ4 RHS)
-max_spatial_derivative_order = 4 # can be 4 or 6
+# In cosmology sixth order stencils are rarely use so the example only have fourth order
+max_spatial_derivative_order = 4 # can be 6 but need to add the stencils in CosmoLevel.cpp (See other examples)
 
 nan_check = 1