Merge pull request #106 from maorz1998/master

add strang splitting in dfLowMachFoam
deepmodeling · Oct 28, 2022 · 7368ed7 · 7368ed7
2 parents 10dc117 + 0ab5c7a
commit 7368ed7
Show file tree

Hide file tree

Showing 4 changed files with 96 additions and 9 deletions.
diff --git a/applications/solvers/dfLowMachFoam/YEqn.H b/applications/solvers/dfLowMachFoam/YEqn.H
@@ -20,13 +20,16 @@ forAll(Y, i)
 const surfaceScalarField phiUc = linearInterpolate(sumYDiffError) & mesh.Sf();
 
 {
-    start = std::clock();
-    combustion->correct();
-    label flag_mpi_init;
-    MPI_Initialized(&flag_mpi_init);
-    if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
-    end = std::clock();
-    time_monitor_chem += double(end - start) / double(CLOCKS_PER_SEC);
+    if (!splitting)
+    {
+        start = std::clock();
+        combustion->correct();
+        label flag_mpi_init;
+        MPI_Initialized(&flag_mpi_init);
+        if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
+        end = std::clock();
+        time_monitor_chem += double(end - start) / double(CLOCKS_PER_SEC);
+    }
 
     volScalarField Yt(0.0*Y[0]);
 
@@ -44,10 +47,13 @@ const surfaceScalarField phiUc = linearInterpolate(sumYDiffError) & mesh.Sf();
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(DEff(), Yi)
               + mvConvection->fvmDiv(phiUc, Yi)
              ==
-                combustion->R(Yi)
+                (
+                    splitting
+                ?   fvm::laplacian(DEff(), Yi)
+                :  (fvm::laplacian(DEff(), Yi) + combustion->R(Yi))
+                )
             );
 
             YiEqn.relax();

diff --git a/applications/solvers/dfLowMachFoam/YEqn_RR.H b/applications/solvers/dfLowMachFoam/YEqn_RR.H
@@ -0,0 +1,61 @@
+if (!(timeIndex % 2))
+{
+    volScalarField Yt(0.0*Y[0]);
+
+    scalar dtSave = runTime.deltaT().value();
+    runTime.setDeltaT(dtSave * 2);
+
+    start = std::clock();
+    combustion->correct();
+
+    label flag_mpi_init;
+    MPI_Initialized(&flag_mpi_init);
+    if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
+    end = std::clock();
+    time_monitor_chem += double(end - start) / double(CLOCKS_PER_SEC);
+
+    forAll(Y, i)
+    {
+        volScalarField& Yi = Y[i];
+
+        if (i != inertIndex)
+        {
+            volScalarField& Yi = Y[i];
+            fvScalarMatrix YiEqn
+            (
+                fvm::ddt(rho, Yi)
+                ==
+                combustion->R(Yi)
+            );
+
+            YiEqn.relax();
+
+            YiEqn.solve("Yi");
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+    }
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+
+    forAll (Y, i)
+    {
+        volScalarField& tYi = Y[i].oldTime();
+
+        forAll(tYi, celli)
+        {
+            tYi[celli] = Y[i][celli];
+        }
+        volScalarField::Boundary& Bf = tYi.boundaryFieldRef(); 
+        forAll(Bf, patchi)
+        {
+            forAll(Bf[patchi], facei)   
+            {
+                Bf[patchi][facei] = Y[i].boundaryField()[patchi][facei];
+            }
+        }
+    }
+
+    runTime.setDeltaT(dtSave);
+}
diff --git a/applications/solvers/dfLowMachFoam/createFields.H b/applications/solvers/dfLowMachFoam/createFields.H
@@ -132,3 +132,16 @@ volVectorField sumYDiffError
     mesh,
     dimensionedVector("sumYDiffError", dimDynamicViscosity/dimLength, Zero)
 );
+
+IOdictionary CanteraTorchProperties
+(
+    IOobject
+    (
+        "CanteraTorchProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+const Switch splitting = CanteraTorchProperties.lookupOrDefault("splittingStrategy", false);
diff --git a/applications/solvers/dfLowMachFoam/dfLowMachFoam.C b/applications/solvers/dfLowMachFoam/dfLowMachFoam.C
@@ -81,6 +81,7 @@ int main(int argc, char *argv[])
     double time_monitor_Y=0;
     double time_monitor_E=0;
     double time_monitor_corrThermo=0;
+    label timeIndex = 0;
     clock_t start, end;
 
     turbulence->validate();
@@ -97,6 +98,8 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
+        timeIndex ++;
+
         #include "readDyMControls.H"
 
         if (LTS)
@@ -116,6 +119,10 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            if (splitting) 
+            {
+                #include "YEqn_RR.H"
+            }
             if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
             {
                 // Store momentum to set rhoUf for introduced faces.