Slowly finding first order violation in Ham on first slice...

Examples/ScalarField/ScalarFieldLevel.cpp

@@ -99,14 +99,14 @@ void ScalarFieldLevel::initialData()
 // Things to do before outputting a checkpoint file
 void ScalarFieldLevel::prePlotLevel()
 {
-    fillAllGhosts();
+    /*fillAllGhosts();
     Potential potential(m_p.potential_params);
     ScalarFieldWithPotential scalar_field(potential);
     BoxLoops::loop(
         MatterConstraints<ScalarFieldWithPotential>(
             scalar_field, m_dx, m_p.G_Newton, c_Ham, Interval(c_Mom, c_Mom), m_p.min_chi, 
 	    c_Ham_abs_terms, Interval(c_Mom_abs_terms, c_Mom_abs_terms)),
-        m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS);
+        m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS);*/
 }
 #endif
 
@@ -184,6 +184,7 @@ void ScalarFieldLevel::specificPostTimeStep()
 	    c_Ham_abs_terms, Interval(c_Mom_abs_terms, c_Mom_abs_terms)),
         m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS);
 
+    MayDay::Error("Constraint calc in sPTS is done.");
     double mass = m_p.potential_params.scalar_mass;
 
     AMRReductions<VariableType::diagnostic> amr_reductions(m_gr_amr);

Source/CCZ4/CCZ4Geometry.hpp

@@ -73,6 +73,7 @@ class CCZ4Geometry
             FOR(k, l) { chris_LLU[i][j][k] += h_UU[k][l] * chris.LLL[i][j][l]; }
         }
 
+        simd<double> ricci_hat_tr(0.);
         FOR(i, j)
         {
             data_t ricci_hat = 0;
@@ -104,6 +105,23 @@ class CCZ4Geometry
 
             out.LL[i][j] =
                 (ricci_chi + vars.chi * ricci_hat + z_terms) / vars.chi;
+
+            simd<double> ricci_hat_tr(0.);
+            if(i==j) { ricci_hat_tr += z_terms; }
+        }
+
+        simd<double> tol(1e-8);
+        simd<double> tr(TensorAlgebra::compute_trace(out.LL, h_UU));
+        if(simd_compare_gt(tr, tol))
+        { 
+            //IntVect m_coords = current_cell.get_int_vect();
+            /*std::cout << "Coords are: \n";
+            for(int s=0; s<3; s++)
+            {
+                std::cout << m_coords[s] << ",";
+            }*/
+            std::cout << "\nRicci LL trace terms at this point: " << tr << "\n";
+            MayDay::Error("RTrace large at the above coords.");
         }
 
         out.scalar = vars.chi * TensorAlgebra::compute_trace(out.LL, h_UU);
@@ -173,6 +191,21 @@ class CCZ4Geometry
             ricci.LL[i][j] += 0.5 * dZ_coeff * z_terms / vars.chi;
         }
         ricci.scalar = vars.chi * TensorAlgebra::compute_trace(ricci.LL, h_UU);
+
+        /*simd<double> tol(1e-8);
+        simd<double> tr(TensorAlgebra::compute_trace(ricci.LL, h_UU));
+        if(simd_compare_gt(tr, tol))
+        { 
+            //IntVect m_coords = current_cell.get_int_vect();
+            /*std::cout << "Coords are: \n";
+            for(int s=0; s<3; s++)
+            {
+                std::cout << m_coords[s] << ",";
+            }
+            std::cout << "\nTrace at this point: " << tr << "\n";
+            MayDay::Error("Trace large at the above coords.");
+        }*/
+
         return ricci;
     }
 

Source/CCZ4/NewConstraints.impl.hpp

@@ -59,10 +59,23 @@ Constraints::Vars<data_t> Constraints::constraint_equations(
         auto A_UU = TensorAlgebra::raise_all(vars.A, h_UU);
         data_t tr_A2 = TensorAlgebra::compute_trace(vars.A, A_UU);
 
-        out.Ham = ricci.scalar +
-                  (GR_SPACEDIM - 1.) * vars.K * vars.K / GR_SPACEDIM - tr_A2;
+        out.Ham = ricci.scalar;// +
+                  //(GR_SPACEDIM - 1.) * vars.K * vars.K / GR_SPACEDIM - tr_A2;
         out.Ham -= 2 * m_cosmological_constant;
 
+        /*simd<double> tol(1e-8);
+        if(simd_compare_gt(ricci.scalar, tol))
+        { 
+            //IntVect m_coords = current_cell.get_int_vect();
+            /*std::cout << "Coords are: \n";
+            for(int s=0; s<3; s++)
+            {
+                std::cout << m_coords[s] << ",";
+            }
+            std::cout << "\nRicci at this point: " << ricci.scalar << "\n";
+            MayDay::Error("Ham constraint large at the above coords.");
+        }*/
+
         out.Ham_abs_terms =
             abs(ricci.scalar) + abs(tr_A2) +
             abs((GR_SPACEDIM - 1.) * vars.K * vars.K / GR_SPACEDIM);

Source/InitialConditions/ScalarFields/RandomField.impl.hpp

@@ -17,7 +17,7 @@
     debug = false;
     kstar = 16.*(2.*M_PI/m_params.L);
     epsilon = 0.5;
-    H0 = -3.0*sqrt((8.0 * M_PI/3.0/m_params.m_pl/m_params.m_pl)
+    H0 = 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(N, 3.);
 
@@ -251,7 +251,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";
 

Source/Matter/NewMatterConstraints.impl.hpp

@@ -31,6 +31,7 @@ void MatterConstraints<matter_t>::compute(Cell<data_t> current_cell) const
     const auto d1 = m_deriv.template diff1<BSSNMatterVars>(current_cell);
     const auto d2 = m_deriv.template diff2<BSSNMatterVars>(current_cell);
 
+
     // Inverse metric and Christoffel symbol
     const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h);
     const auto chris = TensorAlgebra::compute_christoffel(d1.h, h_UU);
@@ -46,6 +47,19 @@ void MatterConstraints<matter_t>::compute(Cell<data_t> current_cell) const
     {
         out.Ham += -16.0 * M_PI * m_G_Newton * emtensor.rho;
         out.Ham_abs_terms += 16.0 * M_PI * m_G_Newton * abs(emtensor.rho);
+
+        /*simd<double> tol(1e-8);
+        if(simd_compare_gt(out.Ham, tol))
+        { 
+            IntVect m_coords = current_cell.get_int_vect();
+            std::cout << "Coords are: \n";
+            for(int s=0; s<3; s++)
+            {
+                std::cout << m_coords[s] << ",";
+            }
+            std::cout << "\nHam at this point: " << out.Ham << "\n";
+            MayDay::Error("Ham constraint large at the above coords.");
+        }*/
     }
 
     // Momentum constraints