Add proca constraint excision within the horizon

Examples/BoostedBHProca/BoostedBHProcaLevel.cpp

@@ -22,7 +22,7 @@
 // Problem specific includes
 #include "EnergyConservation.hpp"
 #include "ExcisionDiagnostics.hpp"
-#include "ExcisionEvolution.hpp"
+#include "ExcisionProcaEvolution.hpp"
 #include "FluxExtraction.hpp"
 #include "InitialProcaData.hpp"
 #include "LinearMomConservation.hpp"
@@ -40,26 +40,37 @@ void BoostedBHProcaLevel::initialData()
     // This is just for the diagnostics
     SetValue set_zero(0.0);
     BoostedBH boosted_bh(m_p.bg_params, m_dx); // just calculates chi
-    auto compute_pack = make_compute_pack(set_zero, boosted_bh); //pack up the classes into a single BoxLoops call
-    BoxLoops::loop(compute_pack, m_state_diagnostics, m_state_diagnostics,
-                   SKIP_GHOST_CELLS); //Loop over box cells, skipping the ghost cells
+    auto compute_pack = make_compute_pack(
+        set_zero,
+        boosted_bh); // pack up the classes into a single BoxLoops call
+    BoxLoops::loop(
+        compute_pack, m_state_diagnostics, m_state_diagnostics,
+        SKIP_GHOST_CELLS); // Loop over box cells, skipping the ghost cells
 
     // Now set the actual evolution variables
-    InitialProcaData initial_pf(m_p.proca_amplitude, m_p.center,
-                                m_p.bg_params, m_dx, m_p.proca_initial_data_profile); //Initial proca field data
-    BoxLoops::loop(initial_pf, m_state_new, m_state_new, FILL_GHOST_CELLS); //Loop over box cells, filling the ghost cells
+    InitialProcaData initial_pf(
+        m_p.proca_amplitude, m_p.center, m_p.bg_params, m_dx,
+        m_p.proca_initial_data_profile); // Initial proca field data
+    BoxLoops::loop(
+        initial_pf, m_state_new, m_state_new,
+        FILL_GHOST_CELLS); // Loop over box cells, filling the ghost cells
 
     // now the gauss constraint
-    fillAllGhosts(); //Is this necessary here?
-    ProcaConstraint enforce_constraint(m_p.center, m_p.bg_params,
-                                       m_p.proca_mass, m_dx); //This sets the scalar part of the Proca field, using the gauss constraint
-    BoxLoops::loop(enforce_constraint, m_state_new, m_state_new,
-                   EXCLUDE_GHOST_CELLS); //Loop over box cells, excluding the ghost cells
+    fillAllGhosts(); // Is this necessary here?
+    ProcaConstraint enforce_constraint(
+        m_p.center, m_p.bg_params, m_p.proca_mass,
+        m_dx); // This sets the scalar part of the Proca field, using the gauss
+               // constraint
+    BoxLoops::loop(
+        enforce_constraint, m_state_new, m_state_new,
+        EXCLUDE_GHOST_CELLS); // Loop over box cells, excluding the ghost cells
 
     // excise evolution vars within horizon, turn off simd vectorisation
-    BoxLoops::loop(
-        ExcisionEvolution<ProcaField, BoostedBH>(m_dx, m_p.center, boosted_bh),
-        m_state_new, m_state_new, SKIP_GHOST_CELLS, disable_simd()); //Now excise the evolution variables inside the horizon
+    BoxLoops::loop(ExcisionProcaEvolution<ProcaField, BoostedBH>(
+                       m_dx, m_p.center, boosted_bh),
+                   m_state_new, m_state_new, SKIP_GHOST_CELLS,
+                   disable_simd()); // Now excise the evolution variables inside
+                                    // the horizon
 }
 
 void BoostedBHProcaLevel::specificPostTimeStep()
@@ -105,8 +116,8 @@ void BoostedBHProcaLevel::specificPostTimeStep()
             double rhoLinMom_sum = amr_reductions.sum(c_rhoLinMom);
             double sourceLinMom_sum = amr_reductions.sum(c_sourceLinMom);
 
-            SmallDataIO integral_file(m_p.data_path + m_p.integrals_filename, m_dt,
-                                      m_time, m_restart_time,
+            SmallDataIO integral_file(m_p.data_path + m_p.integrals_filename,
+                                      m_dt, m_time, m_restart_time,
                                       SmallDataIO::APPEND, first_step);
             // remove any duplicate data if this is post restart
             integral_file.remove_duplicate_time_data();
@@ -117,9 +128,8 @@ void BoostedBHProcaLevel::specificPostTimeStep()
             // write data
             if (first_step)
             {
-                integral_file.write_header_line({"Energy",
-                                                 "Lin. Mom.",
-                                                 "Lin. Mom. Source"});
+                integral_file.write_header_line(
+                    {"Energy", "Lin. Mom.", "Lin. Mom. Source"});
             }
             integral_file.write_time_data_line(data_for_writing);
 
@@ -150,9 +160,9 @@ void BoostedBHProcaLevel::specificEvalRHS(GRLevelData &a_soln,
     BoxLoops::loop(my_evolution, a_soln, a_rhs, SKIP_GHOST_CELLS);
 
     // Do excision within horizon, don't use vectorisation
-    BoxLoops::loop(
-        ExcisionEvolution<ProcaField, BoostedBH>(m_dx, m_p.center, boosted_bh),
-        a_soln, a_rhs, SKIP_GHOST_CELLS, disable_simd());
+    BoxLoops::loop(ExcisionProcaEvolution<ProcaField, BoostedBH>(
+                       m_dx, m_p.center, boosted_bh),
+                   a_soln, a_rhs, SKIP_GHOST_CELLS, disable_simd());
 }
 
 // Note that for the fixed grids this only happens on the initial timestep

Examples/BoostedBHProca/ExcisionProcaEvolution.hpp

@@ -0,0 +1,57 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef EXCISIONPROCAEVOLUTION_HPP_
+#define EXCISIONPROCAEVOLUTION_HPP_
+
+#include "Cell.hpp"
+#include "Coordinates.hpp"
+#include "GRInterval.hpp"
+#include "Tensor.hpp"
+#include "UserVariables.hpp" //This files needs NUM_VARS - total number of components
+#include "VarsTools.hpp"
+
+//! Does excision for fixed BG BH solutions
+template <class matter_t, class background_t> class ExcisionProcaEvolution
+{
+    // Use matter_t class
+    using Vars = typename matter_t::template Vars<double>;
+
+  protected:
+    const double m_dx;                              //!< The grid spacing
+    const std::array<double, CH_SPACEDIM> m_center; //!< The BH center
+    const FourthOrderDerivatives m_deriv;
+    const background_t m_background;
+
+  public:
+    ExcisionProcaEvolution(const double a_dx,
+                           const std::array<double, CH_SPACEDIM> a_center,
+                           background_t a_background)
+        : m_dx(a_dx), m_deriv(m_dx), m_center(a_center),
+          m_background(a_background)
+    {
+    }
+
+    void compute(const Cell<double> current_cell) const
+    {
+        const Coordinates<double> coords(current_cell, m_dx, m_center);
+        bool is_excised = m_background.check_if_excised(coords);
+        if (is_excised)
+        {
+            // the matter rhs vars within the excision zone
+            // recalculate them - for now set to decay to zero
+            Vars vars;
+            VarsTools::assign(vars, 0.0);
+            // assign values of rhs or vars in output box
+            current_cell.store_vars(vars);
+        } // else do nothing
+        if (coords.get_radius() < 1.5)
+        {
+            current_cell.store_vars(0.0, c_Zvec);
+        }
+    }
+};
+
+#endif /* EXCISIONPROCAEVOLUTION_HPP_ */