Put lower bound on the conservative density

Examples/Fluid_Kerr/PerfectFluidLevel.cpp

@@ -30,16 +30,16 @@
 #include "InitialFluidData.hpp"
 #include "KerrBH.hpp"
 #include "PerfectFluid.hpp"
+#include "PositiveDensity.hpp"
 #include "SetValue.hpp"
 
 // Things to do at each advance step, after the RK4 is calculated
 void PerfectFluidLevel::specificAdvance()
 {
     // Enforce trace free A_ij and positive chi and alpha
-    BoxLoops::loop(
-        make_compute_pack(TraceARemoval(),
-                          PositiveChiAndAlpha(m_p.min_chi, m_p.min_lapse)),
-        m_state_new, m_state_new, INCLUDE_GHOST_CELLS);
+    BoxLoops::loop(make_compute_pack(TraceARemoval(), PositiveChiAndAlpha(),
+                                     PositiveDensity()),
+                   m_state_new, m_state_new, INCLUDE_GHOST_CELLS);
 
     // Check for nan's
     if (m_p.nan_check)
@@ -87,10 +87,10 @@ void PerfectFluidLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs,
                                         const double a_time)
 {
     // Enforce trace free A_ij and positive chi and alpha
-    BoxLoops::loop(
-        make_compute_pack(TraceARemoval(), PrimitiveRecovery(),
-                          PositiveChiAndAlpha(m_p.min_chi, m_p.min_lapse)),
-        a_soln, a_soln, INCLUDE_GHOST_CELLS);
+    BoxLoops::loop(make_compute_pack(TraceARemoval(), PositiveDensity(),
+                                     PositiveChiAndAlpha(),
+                                     PrimitiveRecovery()),
+                   a_soln, a_soln, INCLUDE_GHOST_CELLS);
 
     // Calculate MatterCCZ4 right hand side with matter_t = ScalarField
     EoS eos(m_p.eos_params);

Examples/Fluid_Kerr/PositiveDensity.hpp

@@ -0,0 +1,46 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+// This compute class enforces the positive chi and alpha condition
+#ifndef POSITIVEDENSITY_HPP_
+#define POSITIVEDENSITY_HPP_
+
+#include "Cell.hpp"
+#include "UserVariables.hpp"
+#include "simd.hpp"
+
+class PositiveDensity
+{
+  private:
+    const double m_min_D;
+    const double m_min_v;
+
+  public:
+    //! Constructor for class
+    PositiveDensity(const double a_min_D = 1e-4, const double a_min_v = 0.)
+        : m_min_D(a_min_D), m_min_v(a_min_v)
+    {
+    }
+
+    template <class data_t> void compute(Cell<data_t> current_cell) const
+    {
+        auto D = current_cell.load_vars(c_D);
+        auto rho = current_cell.load_vars(c_rho);
+        auto tau = current_cell.load_vars(c_tau);
+        Tensor<1, data_t> vi, Sj;
+        vi[0] = current_cell.load_vars(c_vi1);
+        vi[1] = current_cell.load_vars(c_vi2);
+        vi[2] = current_cell.load_vars(c_vi3);
+
+        auto make_zero = simd_compare_lt(D, m_min_D);
+        D = simd_conditional(make_zero, D, m_min_D);
+        FOR(i) vi[i] = simd_conditional(make_zero, vi[i], m_min_v);
+
+        current_cell.store_vars(D, c_D);
+        current_cell.store_vars(vi, GRInterval<c_vi1, c_vi3>());
+    }
+};
+
+#endif /* POSITIVEDENSITY_HPP_ */