diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/correctContErrs.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/correctContErrs.H
deleted file mode 100644
index 6aeb5d024bfcfc48f5747caa6a7703012a14bba1..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/correctContErrs.H
+++ /dev/null
@@ -1,7 +0,0 @@
-phase1.continuityError() =
-    fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
-  - (fvOptions(alpha1, rho1)&rho1);
-
-phase2.continuityError() =
-    fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
-  - (fvOptions(alpha2, rho2)&rho2);
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
index 338b0d88b67ffc8a2352d12ae98ea1601f52d684..b9029e4190f4aff6b8c2fdd197ca547b1135f7d4 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -89,7 +89,7 @@ tmp<surfaceScalarField> phiF2;
 while (pimple.correct())
 {
     // Update continuity errors due to temperature changes
-    #include "correctContErrs.H"
+    fluid.correct();
 
     // Correct fixed-flux BCs to be consistent with the velocity BCs
     MRF.correctBoundaryFlux(U1, phi1);
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index 8736d0d5a3fba16bc0605cb7b58837878c10d38c..8bd69f606c39812e518010192f690b92755a8acf 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -93,7 +93,7 @@ tmp<surfaceScalarField> Ff2;
 while (pimple.correct())
 {
     // Update continuity errors due to temperature changes
-    #include "correctContErrs.H"
+    fluid.correct();
 
     surfaceScalarField rhof1(fvc::interpolate(rho1));
     surfaceScalarField rhof2(fvc::interpolate(rho2));
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
index 8adb86740163a9aaf50b3f736760376bbc42de37..3d2e2750806c6a3e51737fc6bdf24b8ed96edda8 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
@@ -217,8 +217,11 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correct()
 
     this->fluid().MRF().correctBoundaryVelocity(U_);
 
+    volScalarField& rho = this->thermo().rho();
+
     continuityError_ =
-        fvc::ddt(*this, this->thermo().rho()) + fvc::div(alphaRhoPhi_);
+        fvc::ddt(*this, rho) + fvc::div(alphaRhoPhi_)
+      - (this->fluid().fvOptions()(*this, rho)&rho);
 }
 
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
index e12f52df926c5fb80672d6399cd3e37f604af8d7..691a663ffbe4d604c1d438780153e40b09237755 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
@@ -179,7 +179,7 @@ protected:
         IOMRFZoneList MRF_;
 
         //- Optional FV-options
-        fv::IOoptionList fvOptions_;
+        mutable fv::IOoptionList fvOptions_;
 
         //- Blending methods
         blendingMethodTable blendingMethods_;
@@ -371,7 +371,7 @@ public:
             inline const IOMRFZoneList& MRF() const;
 
             //- Optional FV-options
-            inline fv::IOoptionList& fvOptions();
+            inline fv::IOoptionList& fvOptions() const;
 
             //- Access a sub model between a phase pair
             template <class modelType>
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
index 0613ba47acbcb1309c5fd20b48380bf3d61f2b01..3d7e837533ac4372487e7a5ec09b4e4558557d7b 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
@@ -61,7 +61,7 @@ inline const Foam::IOMRFZoneList& Foam::phaseSystem::MRF() const
 }
 
 
-inline Foam::fv::IOoptionList& Foam::phaseSystem::fvOptions()
+inline Foam::fv::IOoptionList& Foam::phaseSystem::fvOptions() const
 {
     return fvOptions_;
 }
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
index 087434ba5b141036c64cc42dce7858f01ee0c45d..f74e3e93f081252a48dbf30da6697af8120980c7 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
@@ -94,8 +94,6 @@ int main(int argc, char *argv[])
             fluid.solve();
             fluid.correct();
 
-            #include "correctContErrs.H"
-
             #include "YEqns.H"
 
             if (faceMomentum)