diff --git a/Examples/Fluid_Kerr/InitialFluidData.hpp b/Examples/Fluid_Kerr/InitialFluidData.hpp
index 9bcc612d4..24cdcb20e 100644
--- a/Examples/Fluid_Kerr/InitialFluidData.hpp
+++ b/Examples/Fluid_Kerr/InitialFluidData.hpp
@@ -17,7 +17,6 @@
 #include "simd.hpp"
 
 //! Class which sets the initial fluid matter config
-// aka Kevin data
 class InitialFluidData
 {
   public:
@@ -33,61 +32,20 @@ class InitialFluidData
     };
 
     //! The constructor
-    InitialFluidData(params_t a_params, double a_dx)
-        : m_dx(a_dx), m_params(a_params)
-    {
-    }
+    InitialFluidData(params_t a_params, double a_dx);
+    //              : m_dx(a_dx), m_params(a_params)
+    //{
+    //}
 
     //! Function to compute the value of all the initial vars on the grid
-    template <class data_t> void compute(Cell<data_t> current_cell) const
-    {
-        // load vars
-        FluidCCZ4RHS<PerfectFluid<>>::Vars<data_t> vars;
-        VarsTools::assign(vars, 0.);
-        // where am i?
-        Coordinates<data_t> coords(current_cell, m_dx, m_params.center);
-        data_t rr = coords.get_radius();
-        data_t rr2 = rr * rr;
-
-        data_t x = coords.x;
-        double y = coords.y;
-        double z = coords.z;
-
-        // Tensor<1, data_t> vi, Sj;
-        data_t chi_regularised = simd_max(vars.chi, 1e-6);
-
-        // vi[0] = 0.;
-        // vi[1] = 0.;
-        // vi[2] = 0.;
-
-        // calculate the field value
-        vars.rho = m_params.rho0 * (exp(-pow(rr / m_params.awidth, 2.0))) +
-                   m_params.delta;
-        data_t v2 = 0.;
-        FOR(i, j)
-        v2 += vars.h[i][j] * vars.vi[i] * vars.vi[j] / chi_regularised;
-        vars.eps = 0.;
-        data_t P = vars.rho * (1. + vars.eps) / 3.;
-        data_t WW = 1. / (1. - v2);
-        data_t hh = 1. + vars.eps + P / vars.rho;
-
-        vars.D = vars.rho * sqrt(WW);
-        vars.tau = vars.rho * hh * WW - P - vars.D;
-        FOR(i)
-        {
-            vars.Sj[i] = 0.;
-            FOR(j)
-            vars.Sj[i] += vars.rho * hh * WW * vars.h[i][j] * vars.vi[j] /
-                          chi_regularised;
-        }
-
-        // store the vars
-        current_cell.store_vars(vars);
-    }
+    //! The constructor
+    template <class data_t> void compute(Cell<data_t> current_cell) const;
 
   protected:
     double m_dx;
     const params_t m_params;
 };
 
+#include "InitialFluidData.impl.hpp"
+
 #endif /* INITIALFLUIDDATA_HPP_ */
diff --git a/Examples/Fluid_Kerr/InitialFluidData.impl.hpp b/Examples/Fluid_Kerr/InitialFluidData.impl.hpp
new file mode 100644
index 000000000..d87e2743a
--- /dev/null
+++ b/Examples/Fluid_Kerr/InitialFluidData.impl.hpp
@@ -0,0 +1,67 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#if !defined(INITIALFLUIDDATA_HPP_)
+#error "This file should only be included through PrimitiveRecovery.hpp"
+#endif
+
+#ifndef INITIALFLUIDDATA_IMPL_HPP_
+#define INITIALFLUIDDATA_IMPL_HPP_
+
+inline InitialFluidData::InitialFluidData(params_t a_params, double a_dx)
+    : m_dx(a_dx), m_params(a_params)
+{
+}
+
+template <class data_t>
+void InitialFluidData::compute(Cell<data_t> current_cell) const
+{
+    // Set up EoS
+    // data_t P_over_rho = 0.;
+    // my_eos.compute_eos(P_over_rho, vars);
+    // load vars
+    FluidCCZ4RHS<PerfectFluid<>>::Vars<data_t> vars;
+    VarsTools::assign(vars, 0.);
+    // where am i?
+    Coordinates<data_t> coords(current_cell, m_dx, m_params.center);
+    data_t rr = coords.get_radius();
+    data_t rr2 = rr * rr;
+
+    data_t x = coords.x;
+    double y = coords.y;
+    double z = coords.z;
+
+    // Tensor<1, data_t> vi, Sj;
+    data_t chi_regularised = simd_max(vars.chi, 1e-6);
+
+    // vi[0] = 0.;
+    // vi[1] = 0.;
+    // vi[2] = 0.;
+
+    // calculate the field value
+    vars.rho =
+        m_params.rho0 * (exp(-pow(rr / m_params.awidth, 2.0))) + m_params.delta;
+    data_t v2 = 0.;
+    FOR(i, j) v2 += vars.h[i][j] * vars.vi[i] * vars.vi[j] / chi_regularised;
+    vars.eps = 0.;
+    data_t P_over_rho = vars.rho * (1. + vars.eps) / 3.;
+    data_t WW = 1. / (1. - v2);
+    data_t hh = 1. + vars.eps + P_over_rho;
+
+    vars.D = vars.rho * sqrt(WW);
+    vars.tau = vars.rho * (hh * WW - P_over_rho) - vars.D;
+    FOR(i)
+    {
+        vars.Sj[i] = 0.;
+        FOR(j)
+        vars.Sj[i] +=
+            vars.rho * hh * WW * vars.h[i][j] * vars.vi[j] / chi_regularised;
+    }
+
+    // store the vars
+    current_cell.store_vars(vars);
+}
+
+#endif /* INITIALFLUIDDATA_IMPL_HPP_ */
diff --git a/Examples/Fluid_Kerr/PerfectFluid.hpp b/Examples/Fluid_Kerr/PerfectFluid.hpp
index 9e6481e2d..139682702 100644
--- a/Examples/Fluid_Kerr/PerfectFluid.hpp
+++ b/Examples/Fluid_Kerr/PerfectFluid.hpp
@@ -25,13 +25,10 @@
      This class is an example of a matter_t object which calculates the
      matter type specific elements for the RHS update and the evaluation
      of the constraints. This includes the Energy Momentum Tensor, and
-     the matter evolution terms. In this case, a scalar field,
-     the matter elements are phi and (minus) its conjugate momentum, Pi.
-     It is templated over a potential function potential_t which the
+     the matter evolution terms.
+     It is templated over an equation of state function eos_t, which the
      user must specify in a class, although a default is provided which
-     sets dVdphi and V_of_phi to zero.
-     It assumes minimal coupling of the field to gravity.
-     \sa MatterCCZ4(), ConstraintsMatter()
+     sets P_over_rho to (1 + eps) / 3.
 */
 template <class eos_t = DefaultEoS> class PerfectFluid
 {