Midway through debug, just to share with Amelia

Examples/ScalarField/params.txt

@@ -128,7 +128,7 @@ stop_time = 40
 max_steps = 5000
 
 # Spatial derivative order (only affects CCZ4 RHS)
-max_spatial_derivative_order = 4 # can be 4 or 6
+max_spatial_derivative_order = 6 # can be 4 or 6
 
 nan_check = 1
 

Source/CCZ4/CCZ4Geometry.hpp

@@ -82,15 +82,29 @@ class CCZ4Geometry
                 // We call this ricci_hat rather than ricci_tilde as we have
                 // replaced what should be \tilde{Gamma} with \hat{Gamma} in
                 // order to avoid adding terms that cancel later on
-                ricci_hat += 0.5 * (vars.h[k][i] * d1.Gamma[k][j] +
-                                    vars.h[k][j] * d1.Gamma[k][i]);
+                //ricci_hat += (vars.h[k][i] * d1.Gamma[k][j] +
+                //                    vars.h[k][j] * d1.Gamma[k][i]);
+                //if(i==j) { ricci_hat_tr += vars.h[k][i]; }// * d1.Gamma[k][j]; }
+                /*if(i==0 && j==0) 
+                { 
+                    //std::cout << "The i,j flag works.\n";
+                    simd<double> tol(1e-8);
+                    test = (ricci_hat); 
+                    if(simd_compare_gt(test, tol))
+                    {
+                        std::cout << "R hat term at this point: " << test << "\n";
+                        MayDay::Error("LL part is large at this point.");
+                    }
+                }*/
+
                 ricci_hat += 0.5 * vars.Gamma[k] * d1.h[i][j][k];
+
                 FOR(l)
                 {
-                    ricci_hat += -0.5 * h_UU[k][l] * d2.h[i][j][k][l] +
-                                 (chris.ULL[k][l][i] * chris_LLU[j][k][l] +
-                                  chris.ULL[k][l][j] * chris_LLU[i][k][l] +
-                                  chris.ULL[k][i][l] * chris_LLU[k][j][l]);
+                    ricci_hat += -0.5 * h_UU[k][l] * d2.h[i][j][k][l];// +
+                                 //(chris.ULL[k][l][i] * chris_LLU[j][k][l] +
+                                  //chris.ULL[k][l][j] * chris_LLU[i][k][l] +
+                                  //chris.ULL[k][i][l] * chris_LLU[k][j][l]);
                 }
             }
 
@@ -106,22 +120,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> ricci_hat_tr(0.);
+            if(i==j) { ricci_hat_tr += ricci_hat; }
         }
 
         simd<double> tol(1e-8);
-        simd<double> tr(TensorAlgebra::compute_trace(out.LL, h_UU));
-        if(simd_compare_gt(tr, tol))
+        //simd<double> LL(out.LL[0][0]);
+        //simd<double> tr(TensorAlgebra::compute_trace(out.LL, h_UU));
+        if(simd_compare_gt(ricci_hat_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.");
+            std::cout << "\nRhat tr term at this point: " << ricci_hat_tr << "\n";
+            MayDay::Error("R hat tr large at the above coords.");
         }
 
         out.scalar = vars.chi * TensorAlgebra::compute_trace(out.LL, h_UU);

Source/CCZ4/GammaCalculator.hpp

@@ -32,7 +32,7 @@ class GammaCalculator
     };
 
   protected:
-    const FourthOrderDerivatives
+    const SixthOrderDerivatives
         m_deriv; //!< An object for calculating derivatives of the variables
 
   public:
@@ -48,6 +48,23 @@ class GammaCalculator
         using namespace TensorAlgebra;
         const auto h_UU = compute_inverse_sym(vars.h);
         const auto chris = compute_christoffel(d1.h, h_UU);
+        MayDay::Error("Chris calc ended.");
+
+        /*simd<double> tol(1e-8);
+        if(simd_compare_gt(chris.contracted[0], tol))
+        {
+            std::cout << "Christoffel comp 0 is large here: " << chris.contracted[0] << "\n";
+            MayDay::Error("Chris failed.");
+        }*/
+
+        /*simd<double> tol(1e-8);
+        simd<double> chrs(0.);
+        for(int s=0; s<3; s++) { chrs += chris.contracted[s]; }
+        if(simd_compare_gt(chrs, tol))
+        {
+            std::cout << "Chris has this magnitude in GammaCalc: " << chrs << "\n";
+            MayDay::Error("Check chris error.");
+        }*/
 
         // assign values of Gamma^k = h_UU^ij * \tilde{Gamma}^k_ij in the output
         // FArrayBox

Source/CCZ4/NewConstraints.impl.hpp

@@ -59,20 +59,20 @@ 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))
+        if(simd_compare_gt(out.Ham, tol))
         { 
-            //IntVect m_coords = current_cell.get_int_vect();
-            /*std::cout << "Coords are: \n";
+            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";
+            std::cout << "\nRicci at this point: " << out.Ham << "\n";
             MayDay::Error("Ham constraint large at the above coords.");
         }*/
 

Source/InitialConditions/ScalarFields/RandomField.impl.hpp

@@ -70,6 +70,36 @@ void RandomField::compute(Cell<data_t> current_cell) const
         MayDay::Error("RandomField: Cell index greater than resolution^3 at coarsest level.");
     }
 
+    if(hx[0][r] + hx[4][r] + hx[8][r] > 1e-12) 
+    { 
+        std::cout << "Trace of hij is large here: \n";
+        std::cout << "(" << i << "," << j << "," << k << ")\n";
+        std::cout << hx[0][r] + hx[4][r] + hx[8][r] << "\n";
+        MayDay::Error();
+    }
+
+    double det = (1. + m_params.A * hx[0][r]) * (1. + m_params.A * hx[4][r]) * (1. + m_params.A * hx[8][r]) +
+                 2. * pow(m_params.A, 3.) * hx[1][r] * hx[2][r] * hx[5][r] -
+                 (1. + m_params.A * hx[0][r]) * pow(m_params.A, 2.) * hx[5][r] * hx[5][r] -
+                 (1. + m_params.A * hx[4][r]) * pow(m_params.A, 2.) * hx[2][r] * hx[2][r] -
+                 (1. + m_params.A * hx[8][r]) * pow(m_params.A, 2.) * hx[1][r] * hx[1][r];
+
+    if((det - 1.) > 1e-12)
+    {
+        std::cout << "Determinant of hij is large here: \n";
+        std::cout << "(" << i << "," << j << "," << k << ")\n";
+        std::cout << det - 1. << "\n";
+        MayDay::Error();
+    }
+
+    std::ofstream det_print("./det-RF.dat", ios::app);
+
+    det_print << i << "," << j << "," << k << ",";
+    det_print << setprecision(12) << "," << det << "\n";// ": " << (1. + m_params.A * hx[0][r]) << "," << (1. + m_params.A * hx[0][r]) << "," << (1. + m_params.A * hx[0][r]) << ": ";
+    //det_print << setprecision(12) << m_params.A * hx[5][r] << "," << m_params.A * hx[2][r] << "," << m_params.A * hx[1][r] << "\n";
+
+    det_print.close();
+
     if(m_spec_type == "position")
     {
         //store tensor metric variables, g_ij = delta_ij + 1/2 h_ij
@@ -99,6 +129,21 @@ void RandomField::compute(Cell<data_t> current_cell) const
     }
 
     else { MayDay::Error("RandomField: Spec type entered is not a viable option."); }
+
+    /*double det = c_h11 * c_h22 * c_h33 +
+                 2. * c_h12 * c_h13 * c_h23 -
+                 c_h11 * c_h23 * c_h23 -
+                 c_h22 * c_h13 * c_h13 -
+                 c_h33 * c_h12 * c_h12;
+
+    //if((det) > 1e-12)
+    //{
+        std::cout << "Determinant of hij is large here: \n";
+        std::cout << "(" << i << "," << j << "," << k << ")\n";
+        std::cout << c_h11 << "," << c_h22 << "," << c_h33 << "\n";
+        std::cout << det << "\n";
+        MayDay::Error();
+    //}*/
 }
 
 void RandomField::clear_data()

Source/Matter/NewMatterConstraints.impl.hpp

@@ -48,7 +48,7 @@ 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);
+        simd<double> tol(1e-8);
         if(simd_compare_gt(out.Ham, tol))
         { 
             IntVect m_coords = current_cell.get_int_vect();
@@ -59,7 +59,7 @@ void MatterConstraints<matter_t>::compute(Cell<data_t> current_cell) const
             }
             std::cout << "\nHam at this point: " << out.Ham << "\n";
             MayDay::Error("Ham constraint large at the above coords.");
-        }*/
+        }
     }
 
     // Momentum constraints

Source/utils/TensorAlgebra.hpp

@@ -9,8 +9,12 @@
 #include "AlwaysInline.hpp"
 #include "DimensionDefinitions.hpp"
 #include "Tensor.hpp"
+#include "Cell.hpp"
 #include <array>
 
+/// Computes the (i,j) component of the Kronecker delta
+constexpr int delta(int i, int j) { return (i == j); }
+
 template <class data_t> struct chris_t
 {
     Tensor<3, data_t> ULL;        //!< standard christoffel symbols
@@ -25,7 +29,7 @@ template <class data_t>
 ALWAYS_INLINE data_t compute_determinant_sym(const Tensor<2, data_t, 3> &matrix)
 {
     data_t det = matrix[0][0] * matrix[1][1] * matrix[2][2] +
-                 2 * matrix[0][1] * matrix[0][2] * matrix[1][2] -
+                 2. * matrix[0][1] * matrix[0][2] * matrix[1][2] -
                  matrix[0][0] * matrix[1][2] * matrix[1][2] -
                  matrix[1][1] * matrix[0][2] * matrix[0][2] -
                  matrix[2][2] * matrix[0][1] * matrix[0][1];
@@ -53,7 +57,15 @@ template <class data_t>
 Tensor<2, data_t> compute_inverse_sym(const Tensor<2, data_t, 3> &matrix)
 {
     data_t deth = compute_determinant_sym(matrix);
-    data_t deth_inverse = 1. / deth;
+
+    /*std::ofstream det_print("./det-GammaC.dat", ios::app);
+
+    det_print << setprecision(12) << deth << ": " << matrix[0][0] << "," << matrix[1][1] << "," << matrix[2][2] << ": ";
+    det_print << matrix[1][2] << "," << matrix[0][2] << "," << matrix[0][1] << "\n";
+
+    det_print.close();*/
+
+    data_t deth_inverse = ((data_t) 1.) / deth;
     Tensor<2, data_t> h_UU;
     h_UU[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[1][2]) *
                  deth_inverse;
@@ -71,6 +83,32 @@ Tensor<2, data_t> compute_inverse_sym(const Tensor<2, data_t, 3> &matrix)
     h_UU[2][0] = h_UU[0][2];
     h_UU[2][1] = h_UU[1][2];
 
+    simd<double> tr_UU(0.);
+    simd<double> tr_matrix(0.);
+    simd<double> one(1.);
+    simd<double> three(3.);
+    simd<double> tol(1e-8);
+    simd<double> tol2(1e-12);
+    Tensor<2, data_t> dotprod;
+
+    //std::cout << h_UU[0][0] << "\n";// << "\n" << matrix[1][1] * matrix[2][2] << "\n";
+    //std::cout << matrix[1][1] * matrix[2][2] - one << "\n" << (matrix[0][0] - one) * (matrix[0][0] - one) << "\n";
+    //std::cout << matrix[1][1] * matrix[2][2] + matrix[0][0] * matrix[2][2] + matrix[0][0] * matrix[1][1] - three << "\n";
+    //std::cout << matrix[1][2] << "\n" << matrix[0][2] << "\n" << matrix[0][1] << "\n";
+    //std::cout << matrix[1][2] * matrix[1][2] << "\n" << matrix[0][2] * matrix[0][2] << "\n" << matrix[0][1] * matrix[0][1] << "\n";
+    //std::cout << deth_inverse - one << "\n";
+
+    /*for(int i=0; i<3; i++) for(int j=0; j<3; j++)
+    {
+        dotprod[i][j] = matrix[i][j] - h_UU[i][j];
+        //dotprod[i][j] -= delta(i, j);
+        std::cout << dotprod[i][j] << "\n";
+    }*/
+
+    for(int s=0; s<3; s++) { tr_UU += h_UU[s][s]; tr_matrix += matrix[s][s]; }// - one - (matrix[s][s] - one); }
+    if(simd_compare_gt(tr_UU - three, tol) || simd_compare_gt(tr_matrix - three, tol2)) { std::cout << tr_UU-three << "\n" << tr_matrix-three << "\n"; MayDay::Error("Trace cond. failed here (in UU fn)."); }
+    //MayDay::Error("Check tensor dot product.");
+
     return h_UU;
 }
 
@@ -231,9 +269,6 @@ ALWAYS_INLINE Tensor<2, data_t> lower_all(const Tensor<2, data_t> &tensor_UU,
     return raise_all(tensor_UU, metric);
 }
 
-/// Computes the (i,j) component of the Kronecker delta
-constexpr int delta(int i, int j) { return (i == j); }
-
 /// Computes the levi-civita symbol (3D, NB, symbol, not the Tensor)
 inline Tensor<3, double> epsilon()
 {