reactingTwoPhaseEulerFoam: Add fvOption handling to the continuity error

correction in MovingPhaseModel<BasePhaseModel>::correct()

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/correctContErrs.H

-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);

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);

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));

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);
 }
 
 

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>

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_;
 }

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)