Pulled from feature branch for constraint investigation

Examples/ScalarField/DiagnosticVariables.hpp

@@ -13,6 +13,7 @@ enum
     c_sfd,
     c_sf2,
     c_ch2,
+    c_Ham_abs_pbp,
 
     c_a,
     c_H,
@@ -34,6 +35,8 @@ static const std::array<std::string, NUM_DIAGNOSTIC_VARS> variable_names = {
     "Phi",
     "Pi",
     "PhiSq",
+    "ChiSq",
+    "AbsHamPBP",
 
     "ScaleFactor",
     "HubbleFactor",

Examples/ScalarField/Main_ScalarField.cpp

@@ -15,6 +15,7 @@
 #include "GRParmParse.hpp"
 #include "SetupFunctions.hpp"
 #include "SimulationParameters.hpp"
+#include "MultiLevelTask.hpp"
 
 // Problem specific includes:
 #include "ScalarFieldLevel.hpp"
@@ -53,6 +54,19 @@ int runGRChombo(int argc, char *argv[])
             level->specificPostTimeStep();
     };
 
+    // (!!) only un-comment one of the code blocks below
+    // at a time.
+
+    // call specificPostTimeStep at every time step
+    //MultiLevelTaskPtr<> call_task(task);
+    //call_task.execute(gr_amr);
+
+    // call specificPostTimeStep at some_interval
+    int some_interval = 1;
+    bool reverse_levels = true;
+    MultiLevelTaskPtr<> call_task(task, reverse_levels, some_interval);
+    gr_amr.schedule(call_task);
+
     // Engage! Run the evolution
     gr_amr.run(sim_params.stop_time, sim_params.max_steps);
     gr_amr.conclude();

Examples/ScalarField/MeansVars.hpp

@@ -28,8 +28,7 @@ class MeansVars
             data_t Pi;
             data_t chi;
             data_t K;
-
-            data_t h11;
+            data_t Ham;
 
             template <typename mapping_function_t>
             void enum_mapping(mapping_function_t mapping_function);

Examples/ScalarField/MeansVars.impl.hpp

@@ -30,12 +30,14 @@ inline
      data_t phisq = vars.phi*vars.phi;
      data_t chisq = vars.chi*vars.chi;
      data_t kin = vars.Pi*vars.Pi;
+     data_t hamabspbp = abs(vars.Ham);
 
     //store class (Vars) variables as diagnostic variables on the grid
      current_cell.store_vars(vars.phi, c_sf);
      current_cell.store_vars(vars.Pi, c_sfd);
      current_cell.store_vars(vars.chi, c_a);
      current_cell.store_vars(vars.K, c_H);
+     current_cell.store_vars(hamabspbp, c_Ham_abs_pbp);
      current_cell.store_vars(phisq, c_sf2);
      current_cell.store_vars(chisq, c_ch2);
      current_cell.store_vars(kin, c_kin);
@@ -50,7 +52,7 @@ inline
      define_enum_mapping(mapping_function, c_Pi, Pi);
      define_enum_mapping(mapping_function, c_chi, chi);
      define_enum_mapping(mapping_function, c_K, K);
-     define_enum_mapping(mapping_function, c_h11, h11);
+     define_enum_mapping(mapping_function, c_Ham, Ham);
  }
 
  #endif /* MEANSVARS_IMPL_HPP_ */

Examples/ScalarField/ScalarFieldLevel.cpp

@@ -74,7 +74,6 @@ void ScalarFieldLevel::initialData()
 
     pfield.clear_data();
     pout() << "Calculating position ICs ended.\n";
-   // MayDay::Error("Calculating position ICs ended.");
 
     std::vector<std::vector<double> > h(std::pow(N, 3.), std::vector<double>(6, 0.)); // input array memory allocation
     std::vector<std::vector<double> > hdot(std::pow(N, 3.), std::vector<double>(6, 0.));
@@ -90,9 +89,7 @@ void ScalarFieldLevel::initialData()
         make_compute_pack(InitialScalarData(m_p.initial_params)),
     m_state_new, m_state_new, INCLUDE_GHOST_CELLS,disable_simd());
 
-    cout << "IC set-up ended.\n";
-
-    //MayDay::Error("Check for crash.");
+    pout() << "IC set-up ended.\n";
 
     fillAllGhosts();
     BoxLoops::loop(GammaCalculator(m_dx), m_state_new, m_state_new,
@@ -197,6 +194,7 @@ void ScalarFieldLevel::specificPostTimeStep()
     double H = -amr_reductions.sum(c_H)/vol/3.;
 
     double hambar = amr_reductions.sum(c_Ham)/vol;
+    double hamabspbpSum = amr_reductions.sum(c_Ham_abs_pbp)/vol;
     double mombar = amr_reductions.sum(c_Mom)/vol;
     double habsbar = amr_reductions.sum(c_Ham_abs_terms)/vol;
     double mabsbar = amr_reductions.sum(c_Mom_abs_terms)/vol;
@@ -222,7 +220,8 @@ void ScalarFieldLevel::specificPostTimeStep()
     if(first_step) 
     {
         means_file.write_header_line({"Scalar field mean","Scalar field variance","Pi mean","Scale factor","Conformal factor variance","Hubble factor",
-            "Kinetic ED","Potential ED","First SRP","Second SRP","Avg Ham constr","Avg Mom constr","Avg Ham abs term","Avg Mom abs term","Avg lapse"});
+            "Kinetic ED","Potential ED","First SRP","Second SRP","Avg Ham constr","Avg |Ham| constr (point by point)","Avg Mom constr",
+            "Avg Ham abs term","Avg Mom abs term","Avg lapse"});
     }
-    means_file.write_time_data_line({phibar, phivar, pibar, a, chivar, H, kinb, potb, epsilon, delta, hambar, mombar, habsbar, mabsbar, lapse});
+    means_file.write_time_data_line({phibar, phivar, pibar, a, chivar, H, kinb, potb, epsilon, delta, hambar, hamabspbpSum, mombar, habsbar, mabsbar, lapse});
 }

Source/InitialConditions/ScalarFields/RandomField.impl.hpp

@@ -15,8 +15,8 @@
 {
     kstar = 32.*(2.*M_PI/m_params.L);
     epsilon = 0.05;
-    H0 = 0.204692;//-3.0*sqrt((8.0 * M_PI/3.0/m_params.m_pl/m_params.m_pl)
-            //*(0.5*m_params.velocity*m_params.velocity + 0.5*pow(m_params.m * m_params.amplitude, 2.0)));
+    H0 = -3.0*sqrt((8.0 * M_PI/3.0/m_params.m_pl/m_params.m_pl)
+            *(0.5*m_params.velocity*m_params.velocity + 0.5*pow(m_params.m * m_params.amplitude, 2.0)));
     norm = pow(m_params.N, 3.);
 
     calc_spectrum();
@@ -197,8 +197,8 @@ void RandomField::calc_spectrum()
 
         // Start of with random numbers filling the entire array
         // Real parts of h+, hx and hij
-        hplus[k + (N/2+1)*(j + N*i)][0] = find_rayleigh_factor(kmag, m_spec_type, sigma_dist(engine), 0);// * cos(theta_dist(engine));
-        hcross[k + (N/2+1)*(j + N*i)][0] = find_rayleigh_factor(kmag, m_spec_type, sigma_dist(engine), 0);// * cos(theta_dist(engine));
+        hplus[k + (N/2+1)*(j + N*i)][0] = find_rayleigh_factor(kmag, m_spec_type, sigma_dist(engine), 0) * cos(theta_dist(engine));
+        hcross[k + (N/2+1)*(j + N*i)][0] = find_rayleigh_factor(kmag, m_spec_type, sigma_dist(engine), 0) * cos(theta_dist(engine));
 
         calc_transferse_vectors(i, j, k, mhat, nhat);
 
@@ -218,8 +218,8 @@ void RandomField::calc_spectrum()
         // Else, fill the imaginary part of each field appropriately
         else
         {
-            hplus[k + (N/2+1)*(j + N*i)][1] = find_rayleigh_factor(kmag, m_spec_type, sigma_dist(engine), 1);// * sin(theta_dist(engine));
-            hcross[k + (N/2+1)*(j + N*i)][1] = find_rayleigh_factor(kmag, m_spec_type, sigma_dist(engine), 1);// * sin(theta_dist(engine));
+            hplus[k + (N/2+1)*(j + N*i)][1] = find_rayleigh_factor(kmag, m_spec_type, sigma_dist(engine), 1) * sin(theta_dist(engine));
+            hcross[k + (N/2+1)*(j + N*i)][1] = find_rayleigh_factor(kmag, m_spec_type, sigma_dist(engine), 1) * sin(theta_dist(engine));
             for (int l=0; l<3; l++) for (int p=0; p<3; p++)
             {
                 hk[lut[l][p]][k + (N/2+1)*(j + N*i)][1] = ((mhat[l]*mhat[p] - nhat[l]*nhat[p]) * hplus[k + (N/2+1)*(j + N*i)][1]
@@ -246,7 +246,7 @@ void RandomField::calc_spectrum()
     }
 
     hkprint.close();
-    MayDay::Error("Printed file for comparison with stand-alone IC generator.");
+    //MayDay::Error("Printed file for comparison with stand-alone IC generator.");
 
     pout() << "Moving to configuration space.\n";
 
@@ -310,7 +310,7 @@ double RandomField::find_rayleigh_factor(double km, std::string spec_type, doubl
     }
 
     // Apply the tanh window function and the uniform draw
-    windowed_value *= 0.5 * (1.0 - tanh(epsilon * (km - kstar)));// * sqrt(-2. * log(uniform_draw));
+    windowed_value *= 0.5 * (1.0 - tanh(epsilon * (km - kstar))) * sqrt(-2. * log(uniform_draw));
     return windowed_value;
 }