Bugfix/234 (#262)

* first attempt * Compiles and passes checks * IMP in KerrBH --------- Co-authored-by: K Clough <katyclough@yahoo.com>
GRTLCollaboration · Nov 11, 2024 · 3e6efab · 3e6efab
1 parent d73b71d
commit 3e6efab
Show file tree

Hide file tree

Showing 4 changed files with 52 additions and 7 deletions.
diff --git a/Examples/KerrBH/KerrBHLevel.cpp b/Examples/KerrBH/KerrBHLevel.cpp
@@ -8,6 +8,7 @@
 #include "CCZ4RHS.hpp"
 #include "ChiTaggingCriterion.hpp"
 #include "ComputePack.hpp"
+#include "IntegratedMovingPunctureGauge.hpp"
 #include "KerrBHLevel.hpp"
 #include "NanCheck.hpp"
 #include "NewConstraints.hpp"
@@ -50,6 +51,10 @@ void KerrBHLevel::initialData()
     BoxLoops::loop(GammaCalculator(m_dx), m_state_new, m_state_new,
                    EXCLUDE_GHOST_CELLS);
 
+    fillAllGhosts();
+    BoxLoops::loop(IntegratedMovingPunctureGauge(m_p.ccz4_params), m_state_new,
+                   m_state_new, EXCLUDE_GHOST_CELLS);
+
 #ifdef USE_AHFINDER
     // Diagnostics needed for AHFinder
     BoxLoops::loop(Constraints(m_dx, c_Ham, Interval(c_Mom1, c_Mom3)),
@@ -82,15 +87,17 @@ void KerrBHLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs,
     // Calculate CCZ4 right hand side
     if (m_p.max_spatial_derivative_order == 4)
     {
-        BoxLoops::loop(CCZ4RHS<MovingPunctureGauge, FourthOrderDerivatives>(
-                           m_p.ccz4_params, m_dx, m_p.sigma, m_p.formulation),
-                       a_soln, a_rhs, EXCLUDE_GHOST_CELLS);
+        BoxLoops::loop(
+            CCZ4RHS<IntegratedMovingPunctureGauge, FourthOrderDerivatives>(
+                m_p.ccz4_params, m_dx, m_p.sigma, m_p.formulation),
+            a_soln, a_rhs, EXCLUDE_GHOST_CELLS);
     }
     else if (m_p.max_spatial_derivative_order == 6)
     {
-        BoxLoops::loop(CCZ4RHS<MovingPunctureGauge, SixthOrderDerivatives>(
-                           m_p.ccz4_params, m_dx, m_p.sigma, m_p.formulation),
-                       a_soln, a_rhs, EXCLUDE_GHOST_CELLS);
+        BoxLoops::loop(
+            CCZ4RHS<IntegratedMovingPunctureGauge, SixthOrderDerivatives>(
+                m_p.ccz4_params, m_dx, m_p.sigma, m_p.formulation),
+            a_soln, a_rhs, EXCLUDE_GHOST_CELLS);
     }
 }
 

diff --git a/Source/CCZ4/IntegratedMovingPunctureGauge.hpp b/Source/CCZ4/IntegratedMovingPunctureGauge.hpp
@@ -91,6 +91,18 @@ class IntegratedMovingPunctureGauge
             rhs.B[i] = 0.; // static, stays the same to save initial condition
         }
     }
+
+    // Add the matter terms to the rhs of the gauge equations
+    // For the IMP gauge, this function should do nothing
+    template <class data_t, template <typename> class vars_t>
+    inline void rhs_gauge_add_matter_terms(vars_t<data_t> &matter_rhs,
+                                           const vars_t<data_t> &matter_vars,
+                                           Tensor<2, data_t, 3> h_UU,
+                                           const emtensor_t<data_t> emtensor,
+                                           const double G_Newton) const
+    {
+        // no changes as B should remain constant
+    }
 };
 
 #endif /* INTEGRATEDMOVINGPUNCTUREGAUGE_HPP_ */
diff --git a/Source/CCZ4/MovingPunctureGauge.hpp b/Source/CCZ4/MovingPunctureGauge.hpp
@@ -64,6 +64,29 @@ class MovingPunctureGauge
                        rhs.Gamma[i] - m_params.eta * vars.B[i];
         }
     }
+
+    // Add the matter terms to the rhs of the gauge equations
+    // For the MP gauge, this only changes the rhs of B through the rhs of Gamma
+    template <class data_t, template <typename> class vars_t>
+    inline void rhs_gauge_add_matter_terms(vars_t<data_t> &matter_rhs,
+                                           const vars_t<data_t> &matter_vars,
+                                           Tensor<2, data_t, 3> h_UU,
+                                           const emtensor_t<data_t> emtensor,
+                                           const double G_Newton) const
+    {
+        FOR(i)
+        {
+            data_t matter_term_Gamma = 0.0;
+            FOR(j)
+            {
+                matter_term_Gamma += -16.0 * M_PI * G_Newton *
+                                     matter_vars.lapse * h_UU[i][j] *
+                                     emtensor.Si[j];
+            }
+
+            matter_rhs.B[i] += matter_term_Gamma;
+        }
+    }
 };
 
 #endif /* MOVINGPUNCTUREGAUGE_HPP_ */
diff --git a/Source/Matter/MatterCCZ4RHS.impl.hpp b/Source/Matter/MatterCCZ4RHS.impl.hpp
@@ -101,8 +101,11 @@ void MatterCCZ4RHS<matter_t, gauge_t, deriv_t>::add_emtensor_rhs(
         }
 
         matter_rhs.Gamma[i] += matter_term_Gamma;
-        matter_rhs.B[i] += matter_term_Gamma;
     }
+
+    // Add matter contribution to RHS of gauge evolution
+    this->m_gauge.rhs_gauge_add_matter_terms(matter_rhs, matter_vars, h_UU,
+                                             emtensor, m_G_Newton);
 }
 
 #endif /* MATTERCCZ4RHS_IMPL_HPP_ */