Pre-pull request update

Examples/FRWScalarField/FRWScalarLevel.cpp

@@ -16,8 +16,8 @@
 #include "FixedGridsTaggingCriterion.hpp"
 
 // Problem specific includes
-#include "InitialScalarData.hpp"
 #include "FRW.hpp"
+#include "InitialScalarData.hpp"
 #include "MatterEvolution.hpp"
 #include "ScalarField.hpp"
 #include "ScalarPotential.hpp"
@@ -32,18 +32,17 @@ void FRWScalarLevel::initialData()
     // First set everything to zero ... we don't want undefined values in
     // constraints etc, then initial conditions for fields
     SetValue set_zero(0.0);
-    FRW frw_bg(m_p.bg_params, m_dx, m_time); 
+    FRW frw_bg(m_p.bg_params, m_dx, m_time);
     InitialScalarData initial_sf(m_p.initial_params, m_dx);
     auto compute_pack = make_compute_pack(set_zero, frw_bg);
     BoxLoops::loop(compute_pack, m_state_diagnostics, m_state_diagnostics,
                    SKIP_GHOST_CELLS);
     BoxLoops::loop(initial_sf, m_state_new, m_state_new, FILL_GHOST_CELLS);
-
 }
 
 // Things to do in RHS update, at each RK4 step
 void FRWScalarLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs,
-                                        const double a_time)
+                                     const double a_time)
 {
     // Calculate MatterCCZ4 right hand side with matter_t = ScalarField
     // We don't want undefined values floating around in the constraints so
@@ -54,7 +53,6 @@ void FRWScalarLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs,
     MatterEvolution<ScalarFieldWithPotential, FRW> my_matter(
         scalar_field, frw_bg, m_p.sigma, m_dx, m_p.center);
     BoxLoops::loop(my_matter, a_soln, a_rhs, SKIP_GHOST_CELLS);
-
 }
 
 void FRWScalarLevel::specificPostTimeStep()
@@ -65,25 +63,23 @@ void FRWScalarLevel::specificPostTimeStep()
                        disable_simd());
 
     // At any level, but after the min_level timestep
-    int min_level = 0; 
+    int min_level = 0;
     bool calculate_diagnostics = at_level_timestep_multiple(min_level);
     if (calculate_diagnostics)
     {
         fillAllGhosts();
         ScalarPotential potential(m_p.initial_params);
         ScalarFieldWithPotential scalar_field(potential);
         FRW frw_bg(m_p.bg_params, m_dx, m_time);
-        BoxLoops::loop(frw_bg, m_state_new,
-                       m_state_diagnostics, SKIP_GHOST_CELLS);
-
+        BoxLoops::loop(frw_bg, m_state_new, m_state_diagnostics,
+                       SKIP_GHOST_CELLS);
     }
 
     // write out the integral after each timestep on minimum level
-
 }
 
 void FRWScalarLevel::computeTaggingCriterion(FArrayBox &tagging_criterion,
-                                                const FArrayBox &current_state)
+                                             const FArrayBox &current_state)
 {
     BoxLoops::loop(FixedGridsTaggingCriterion(m_dx, m_level, m_p.L, m_p.center),
                    current_state, tagging_criterion);

Examples/FRWScalarField/Main_FRWScalar.cpp

@@ -28,7 +28,7 @@ int runGRChombo(int argc, char *argv[])
     // and an associated LevelFactory)
     GRAMR gr_amr;
     DefaultLevelFactory<FRWScalarLevel> scalar_field_level_fact(gr_amr,
-                                                                   sim_params);
+                                                                sim_params);
     setupAMRObject(gr_amr, scalar_field_level_fact);
 
     // call this after amr object setup so grids known

Examples/FRWScalarField/SimulationParameters.hpp

@@ -11,8 +11,8 @@
 #include "GRParmParse.hpp"
 
 // Problem specific includes:
-#include "InitialScalarData.hpp"
 #include "FRW.hpp"
+#include "InitialScalarData.hpp"
 
 class SimulationParameters : public FixedBGSimulationParametersBase
 {
@@ -35,8 +35,8 @@ class SimulationParameters : public FixedBGSimulationParametersBase
         // Initial and Kerr data
         pp.load("frw_rho0", bg_params.rho0, 0.0);
         pp.load("frw_omega", bg_params.omega, 0.0);
+        pp.load("frw_alpha", bg_params.alpha, 0.0);
         pp.load("center", bg_params.center, center);
-
     }
 
     void check_params()
@@ -45,7 +45,6 @@ class SimulationParameters : public FixedBGSimulationParametersBase
                        initial_params.mass < 0.2 / coarsest_dx / dt_multiplier,
                        "oscillations of scalar field do not appear to be "
                        "resolved on coarsest level");
-
     }
 
     // Problem specific parameters

Examples/FRWScalarField/params.txt

@@ -42,7 +42,8 @@ scalar_amplitude = 0.1
 
 
 frw_rho0 = 1.0 # Initial energy density
-frw_omega = 0.0 # Equation of state parameter
+frw_omega = 0.0 # Background equation of state 
+frw_alpha = 0.0  # = 0 (cosmic time), 1 (conformal time)
 
 
 

Source/Background/FRW.hpp

@@ -17,123 +17,109 @@
 
 //! Class which computes an FRW background
 
-
 class FRW
 {
   public:
     //! Struct for the params of the  BH
     struct params_t
     {
-        double rho0; // initial energy density
+        double rho0;  // initial energy density
         double omega; // equation of state parameter = 0 (MD), 1/3 (RD)
-        std::array<double, CH_SPACEDIM> center; 
-
+        double alpha; // can set to conformal time
+        std::array<double, CH_SPACEDIM> center;
     };
 
     template <class data_t> using Vars = ADMFixedBGVars::Vars<data_t>;
 
-
     const params_t m_params;
 
     const double m_time;
     const double m_dx;
-
-
-   FRW(params_t a_params, double a_dx, double a_time) : m_params(a_params), m_dx(a_dx), m_time(a_time){}
 
+    FRW(params_t a_params, double a_dx, double a_time)
+        : m_params(a_params), m_dx(a_dx), m_time(a_time)
+    {
+    }
 
     template <class data_t> void compute(Cell<data_t> current_cell) const
     {
 
-
-
         const Coordinates<data_t> coords(current_cell, m_dx, m_params.center);
         Vars<data_t> metric_vars;
 
-
         compute_metric_background(metric_vars, coords);
-
-
-
 
-       //Reconstruct FRW quantities to check as diagnostics
-       data_t scalefac = sqrt(metric_vars.gamma[0][0]);
-       data_t hubbleparam = -metric_vars.K_tensor[0][0];
+        // Reconstruct FRW quantities to check as diagnostics
+        data_t scalefac = sqrt(metric_vars.gamma[0][0]);
+        data_t hubbleparam = -metric_vars.K_tensor[0][0];
 
-       current_cell.store_vars(scalefac, c_scalefac);
-       current_cell.store_vars(hubbleparam, c_hubbleparam);
-
-
-
+        current_cell.store_vars(scalefac, c_scalefac);
+        current_cell.store_vars(hubbleparam, c_hubbleparam);
     }
 
-
-    //Analytical background
-     template <class data_t, template <typename> class vars_t>
+    // Analytical background
+    template <class data_t, template <typename> class vars_t>
     void compute_metric_background(vars_t<data_t> &vars,
                                    const Coordinates<data_t> &coords) const
     {
 
-        data_t scalefac; 
-        data_t hubparam; 
-
-        // Zero components - to do this in level?
-        vars.K = 0.0;
-
-        FOR(i) vars.shift[i] = 0.0;
-        FOR(i,j) vars.gamma[i][j] = 0.0;
-        FOR(i,j) vars.K_tensor[i][j] = 0.0;
-
-
-        FOR(i,j, k) vars.d1_gamma[i][j][k] = 0.0;
-        FOR(i) vars.d1_lapse[i] = 0.0;
-        FOR(i,j) vars.d1_shift[i][j] = 0.0;
+        data_t scalefac;
+        data_t hubparam;
 
-        // ds^2 = -dt^2 + a^2 dx_i dx^i assumin flat - TODO: curvature, conformal time
+        // ds^2 = -dt^2 + a^2 dx_i dx^i
 
-        vars.lapse = 1.0; // Conformal time - a^2
-
-        if (m_params.rho0 == 0.0){
+        if (m_params.rho0 == 0.0)
+        {
 
             scalefac = 1.0;
 
             FOR(i) vars.gamma[i][i] = scalefac; // Minkowski
 
+            vars.lapse = 1.0;
         }
-
-
-        else{
-
-        data_t K0 = -sqrt(24.0*M_PI*m_params.rho0); //Ham constraint
-
-        data_t t0 = -1.0/((1.0+m_params.omega)*K0); //Defining t0 from time evolution of scale factor
-
-
-
-        scalefac = pow((m_time/t0)+1, 2/(3*(1+m_params.omega)));  //Scale factor evolution
-
-        hubparam = 2.0/(3.0*(1+m_params.omega)*(m_time+t0)); //Hubble parameter evolution
-
-
-
-
-        vars.K = -3.0*hubparam; //K = -3H (not used but included for completeness)
-
-        FOR(i) vars.gamma[i][i] = pow(scalefac, 2.0); //gamma_ij = a^2
-
-        FOR(i) vars.K_tensor[i][i] = -hubparam; //Kij = -H gamma_ij
-
-        }
-
 
-    }
+        else
+        {
+
+            data_t K0 = -sqrt(24.0 * M_PI * m_params.rho0); // Ham constraint
+
+            data_t t0 =
+                -1.0 / ((1.0 + m_params.omega) *
+                        K0); // Defining t0 from time evolution of scale factor
+
+            if (m_params.omega == -1)
+            { // Cosmo constant
 
+                hubparam = -K0 / 3.0;
 
-
-
+                scalefac = exp(hubparam * m_time);
+            }
 
+            else
+            {
 
+                scalefac = pow(
+                    (m_time / t0) + 1,
+                    2 / (3 * (1 + m_params.omega))); // Scale factor evolution
 
+                hubparam = 2.0 / (3.0 * (1 + m_params.omega) *
+                                  (m_time + t0)); // Hubble parameter evolution
+            }
+
+            vars.lapse = pow(
+                scalefac,
+                m_params
+                    .alpha); // Cosmic time alpha = 0, conformal time alpha = 1
+
+            vars.K =
+                -3.0 *
+                hubparam; // K = -3H (not used but included for completeness)
+
+            FOR(i) vars.gamma[i][i] = pow(scalefac, 2.0); // gamma_ij = a^2
+
+            FOR(i) vars.K_tensor[i][i] = -hubparam; // Kij = -H gamma_ij
+        }
+    }
 };
 
 #endif /* FRW_HPP_ */