diff --git a/Examples/Fluid_Kerr_IDfromFile/PositiveDensity.hpp b/Examples/Fluid_Kerr_IDfromFile/PositiveDensity.hpp
index 9fb8c6543..bd1dcd5b4 100644
--- a/Examples/Fluid_Kerr_IDfromFile/PositiveDensity.hpp
+++ b/Examples/Fluid_Kerr_IDfromFile/PositiveDensity.hpp
@@ -20,6 +20,7 @@ class PositiveDensity
         data_t D;
         data_t tau;
         Tensor<1, data_t> Sj;
+        Tensor<1, data_t> vi;
         Tensor<2, data_t> h;
 
         template <typename mapping_function_t>
@@ -32,6 +33,8 @@ class PositiveDensity
             define_enum_mapping(mapping_function, c_tau, tau);
             define_enum_mapping(mapping_function, GRInterval<c_Sj1, c_Sj3>(),
                                 Sj);
+	    define_enum_mapping(mapping_function, GRInterval<c_vi1, c_vi3>(),
+                                vi);
             define_symmetric_enum_mapping(mapping_function,
                                           GRInterval<c_h11, c_h33>(), h);
         }
@@ -43,7 +46,7 @@ class PositiveDensity
 
   public:
     //! Constructor for class
-    PositiveDensity(const double a_min_D = 1e-12, const double a_min_v = 0.)
+    PositiveDensity(const double a_min_D = 1.e-12, const double a_min_v = 0.)
         : m_min_D(a_min_D), m_min_v(a_min_v)
     {
     }
@@ -59,16 +62,25 @@ class PositiveDensity
         // D >= min_D
         vars.D = simd_max(vars.D, m_min_D / pow(chi_regularised, 1.5));
 
-        data_t S2_over_chi = 0.;
-        FOR(i, j) S2_over_chi += vars.h[i][j] * vars.Sj[i] * vars.Sj[j];
+        data_t S2 = 0.;
+        FOR(i, j) S2 += vars.h[i][j] * vars.Sj[i] * vars.Sj[j] / chi_regularised;
 
+	vars.tau = simd_max(vars.tau, 1.e-12/ pow(chi_regularised, 1.5));
+
+	data_t empty;
+        if (!simd_all_false(simd_compare_lt(vars.D, m_min_D
+	    / pow(chi_regularised, 1.5)), empty)) FOR(i) vars.vi[i] = 0.;
         //  S^2 <= (D + tau)^2
         data_t factor =
-            simd_min(1., vars.chi * (vars.D + vars.tau) / sqrt(S2_over_chi));
-        FOR(i) vars.Sj[i] *= factor;
-
+	  simd_min(1., (vars.D + vars.tau) / sqrt(abs(S2)));
+        //FOR(i) vars.Sj[i] *= factor;
+	
+	FOR(i) vars.Sj[i] = simd_min(vars.Sj[i], (vars.D + vars.tau)/sqrt(3.));
+	
         current_cell.store_vars(vars.D, c_D);
+	current_cell.store_vars(vars.tau, c_tau);
         current_cell.store_vars(vars.Sj, GRInterval<c_Sj1, c_Sj3>());
+	current_cell.store_vars(vars.vi, GRInterval<c_vi1, c_vi3>());
     }
 };
 
diff --git a/Examples/Fluid_Kerr_IDfromFile/PrimitiveRecovery.hpp b/Examples/Fluid_Kerr_IDfromFile/PrimitiveRecovery.hpp
index 25fce3b5f..4287cf79d 100644
--- a/Examples/Fluid_Kerr_IDfromFile/PrimitiveRecovery.hpp
+++ b/Examples/Fluid_Kerr_IDfromFile/PrimitiveRecovery.hpp
@@ -46,10 +46,13 @@ class PrimitiveRecovery : public EoS
         data_t S2 = 0.;
         FOR(i, j) S2 += vars.chi * h_UU[i][j] * vars.Sj[i] * vars.Sj[j];
 
+	//data_t D_regularised = simd_max(1e-12, vars.D);
+	//data_t tau_regularised = simd_max(1e-12, vars.tau);
+	
         data_t q = vars.tau / vars.D;
         data_t r = S2 / pow(vars.D, 2.);
         data_t xa = 1.5 * (1. + q);
-        Wa = sqrt(pow(xa, 2.) / (pow(xa, 2.) - r));
+        Wa = sqrt(pow(xa, 2.) / abs(pow(xa, 2.) - r));
 
         vars.rho = pow(vars.chi, 1.5) * vars.D / Wa;
         vars.eps = -1. + xa / Wa * (1. - Wa * Wa) + Wa * (1. + q);
@@ -85,6 +88,11 @@ class PrimitiveRecovery : public EoS
             FOR(j) vars.vi[i] += vars.chi * h_UU[i][j] * vi_D[j];
         }
 
+	data_t chi_regularised = simd_max(1e-6, vars.chi);
+	data_t empty2;
+	if (!simd_all_false(simd_compare_lt(vars.D, 1.e-12
+            / pow(chi_regularised, 1.5)), empty2)) FOR(i) vars.vi[i] = 0.;
+
         current_cell.store_vars(vars);
     }
 };
diff --git a/Examples/Fluid_Kerr_IDfromFile/SimulationParameters.hpp b/Examples/Fluid_Kerr_IDfromFile/SimulationParameters.hpp
index f76955e42..bcffac42f 100644
--- a/Examples/Fluid_Kerr_IDfromFile/SimulationParameters.hpp
+++ b/Examples/Fluid_Kerr_IDfromFile/SimulationParameters.hpp
@@ -36,6 +36,7 @@ class SimulationParameters : public SimulationParametersBase
         pp.load("fluid_width", initial_params.awidth, 0.1);
         pp.load("fluid_delta", initial_params.delta, 0.2);
         pp.load("lambda", lambda, 1.); // eigenvalue for numerical flux
+	pp.load("min_D", min_D, 1.e-6);
 
         // Reading data
         pp.load("spacing", initial_params.spacing);
@@ -98,6 +99,7 @@ class SimulationParameters : public SimulationParametersBase
     // Initial data for matter and potential and BH
     double G_Newton;
     double lambda;
+    double min_D;
     int lines;
     InitialData::params_t initial_params;
     // KerrBH::params_t kerr_params;