From 8995550787a14dc4163476b129315dc7a6e6efd9 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Tue, 1 Dec 2015 10:04:56 +0000
Subject: [PATCH] PDRFoam: Added support for fvOptions

---
 applications/solvers/combustion/PDRFoam/EaEqn.H       | 11 ++++++++++-
 applications/solvers/combustion/PDRFoam/EauEqn.H      | 10 +++++++++-
 applications/solvers/combustion/PDRFoam/PDRFoam.C     |  3 +++
 .../solvers/combustion/PDRFoam/PDRFoamAutoRefine.C    |  1 +
 applications/solvers/combustion/PDRFoam/UEqn.H        |  4 ++++
 applications/solvers/combustion/PDRFoam/bEqn.H        |  8 ++++++++
 applications/solvers/combustion/PDRFoam/pEqn.H        |  5 +++++
 7 files changed, 40 insertions(+), 2 deletions(-)

diff --git a/applications/solvers/combustion/PDRFoam/EaEqn.H b/applications/solvers/combustion/PDRFoam/EaEqn.H
index 1baaa7180f5..8b844a81591 100644
--- a/applications/solvers/combustion/PDRFoam/EaEqn.H
+++ b/applications/solvers/combustion/PDRFoam/EaEqn.H
@@ -1,7 +1,7 @@
 {
     volScalarField& hea = thermo.he();
 
-    solve
+    fvScalarMatrix EaEqn
     (
         betav*fvm::ddt(rho, hea) + mvConvection->fvmDiv(phi, hea)
       + betav*fvc::ddt(rho, K) + fvc::div(phi, K)
@@ -16,7 +16,16 @@
           : -betav*dpdt
         )
       - fvm::laplacian(Db, hea)
+      + betav*fvOptions(rho, hea)
     );
 
+    EaEqn.relax();
+
+    fvOptions.constrain(EaEqn);
+
+    EaEqn.solve();
+
+    fvOptions.correct(hea);
+
     thermo.correct();
 }
diff --git a/applications/solvers/combustion/PDRFoam/EauEqn.H b/applications/solvers/combustion/PDRFoam/EauEqn.H
index f031d907323..91bb49b9136 100644
--- a/applications/solvers/combustion/PDRFoam/EauEqn.H
+++ b/applications/solvers/combustion/PDRFoam/EauEqn.H
@@ -2,7 +2,7 @@ if (ign.ignited())
 {
     volScalarField& heau = thermo.heu();
 
-    solve
+    fvScalarMatrix heauEqn
     (
         betav*fvm::ddt(rho, heau) + mvConvection->fvmDiv(phi, heau)
       + (betav*fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
@@ -23,5 +23,13 @@ if (ign.ignited())
         // A possible solution would be to solve for ftu as well as ft.
         //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau)
         //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau)
+     ==
+        betav*fvOptions(rho, heau)
     );
+
+    fvOptions.constrain(heauEqn);
+
+    heauEqn.solve();
+
+    fvOptions.correct(heau);
 }
diff --git a/applications/solvers/combustion/PDRFoam/PDRFoam.C b/applications/solvers/combustion/PDRFoam/PDRFoam.C
index fa8518ba495..65f52e2ee4c 100644
--- a/applications/solvers/combustion/PDRFoam/PDRFoam.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoam.C
@@ -77,6 +77,7 @@ Description
 #include "Switch.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "fvIOoptionList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -93,11 +94,13 @@ int main(int argc, char *argv[])
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createMRF.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "createTimeControls.H"
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    turbulence->validate();
     scalar StCoNum = 0.0;
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C b/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
index 057286f1e94..4e561052fd0 100644
--- a/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
@@ -87,6 +87,7 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    turbulence->validate();
     scalar StCoNum = 0.0;
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/combustion/PDRFoam/UEqn.H b/applications/solvers/combustion/PDRFoam/UEqn.H
index 1d949d17b7c..983e020cab7 100644
--- a/applications/solvers/combustion/PDRFoam/UEqn.H
+++ b/applications/solvers/combustion/PDRFoam/UEqn.H
@@ -7,13 +7,17 @@
       + turbulence->divDevRhoReff(U)
      ==
         betav*rho*g
+      + betav*fvOptions(rho, U)
     );
 
+    fvOptions.constrain(UEqn);
+
     volSymmTensorField invA(inv(I*UEqn.A() + drag->Dcu()));
 
     if (pimple.momentumPredictor())
     {
         U = invA & (UEqn.H() - betav*fvc::grad(p));
         U.correctBoundaryConditions();
+        fvOptions.correct(U);
         K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/combustion/PDRFoam/bEqn.H b/applications/solvers/combustion/PDRFoam/bEqn.H
index abec1938c32..8ad11618713 100644
--- a/applications/solvers/combustion/PDRFoam/bEqn.H
+++ b/applications/solvers/combustion/PDRFoam/bEqn.H
@@ -73,6 +73,8 @@ if (ign.ignited())
       + fvm::div(phiSt, b)
       - fvm::Sp(fvc::div(phiSt), b)
       - fvm::laplacian(Db, b)
+     ==
+        betav*fvOptions(rho, b)
     );
 
 
@@ -82,8 +84,14 @@ if (ign.ignited())
 
     // Solve for b
     // ~~~~~~~~~~~
+    bEqn.relax();
+
+    fvOptions.constrain(bEqn);
+
     bEqn.solve();
 
+    fvOptions.correct(b);
+
     Info<< "min(b) = " << min(b).value() << endl;
 
     if (composition.contains("ft"))
diff --git a/applications/solvers/combustion/PDRFoam/pEqn.H b/applications/solvers/combustion/PDRFoam/pEqn.H
index 714a3e84e8b..499885875bd 100644
--- a/applications/solvers/combustion/PDRFoam/pEqn.H
+++ b/applications/solvers/combustion/PDRFoam/pEqn.H
@@ -25,6 +25,8 @@ if (pimple.transonic())
             betav*fvm::ddt(psi, p)
           + fvm::div(phid, p)
           - fvm::laplacian(rho*invA, p)
+         ==
+            betav*fvOptions(psi, p, rho.name())
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -53,6 +55,8 @@ else
             betav*fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(rho*invA, p)
+         ==
+            betav*fvOptions(psi, p, rho.name())
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -69,6 +73,7 @@ else
 
 U = HbyA - (invA & (betav*fvc::grad(p)));
 U.correctBoundaryConditions();
+fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())
-- 
GitLab