applications/solvers/DNS/dnsFoam/globalProperties.H → applications/solvers/DNS/dnsFoam/UEqn.H

+force.primitiveFieldRef() = ReImSum
+(
+    fft::reverseTransform
+    (
+        K/(mag(K) + 1.0e-6) ^ forceGen.newField(), K.nn()
+    )*recRootN
+);
+
 {
     Info<< "k(" << runTime.timeName() << ") = "
         << 3.0/2.0*average(magSqr(U)).value() << endl;
@@ -8,3 +16,14 @@
     Info<< "U.f(" << runTime.timeName() << ") = "
         << 181.0*average(U & force).value() << endl;
 }
+
+fvVectorMatrix UEqn
+(
+    fvm::ddt(U)
+  + fvm::div(phi, U)
+  - fvm::laplacian(nu, U)
+ ==
+    force
+);
+
+solve(UEqn == -fvc::grad(p));

applications/solvers/DNS/dnsFoam/createFields.H

-#include "readTransportProperties.H"
+Info<< "Reading transportProperties\n" << endl;
+
+IOdictionary transportProperties
+(
+    IOobject
+    (
+        "transportProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+dimensionedScalar nu
+(
+    "nu",
+    dimViscosity,
+    transportProperties
+);
 
 Info<< "Reading field p\n" << endl;
 volScalarField p
@@ -32,4 +51,31 @@ volVectorField U
 
 mesh.setFluxRequired(p.name());
 
-#include "readTurbulenceProperties.H"
+Info<< "Reading turbulenceProperties\n" << endl;
+
+IOdictionary turbulenceProperties
+(
+    IOobject
+    (
+        "turbulenceProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+volVectorField force
+(
+    U/dimensionedScalar("dt", dimTime, runTime.deltaTValue())
+);
+
+Kmesh K(mesh);
+UOprocess forceGen(K, runTime.deltaTValue(), turbulenceProperties);
+
+label ntot = 1;
+forAll(K.nn(), idim)
+{
+    ntot *= K.nn()[idim];
+}
+const scalar recRootN = 1.0/Foam::sqrt(scalar(ntot));

applications/solvers/DNS/dnsFoam/dnsFoam.C

@@ -32,6 +32,9 @@ Group
 Description
     Direct numerical simulation solver for boxes of isotropic turbulence.
 
+Note
+    The solver requires the optional FFTW library.
+
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
@@ -62,70 +65,18 @@ int main(int argc, char *argv[])
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    label ntot = 1;
-    forAll(K.nn(), idim)
-    {
-        ntot *= K.nn()[idim];
-    }
-    const scalar recRootN = 1.0/Foam::sqrt(scalar(ntot));
-
     Info<< nl << "Starting time loop" << endl;
 
     while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        force.primitiveFieldRef() = ReImSum
-        (
-            fft::reverseTransform
-            (
-                K/(mag(K) + 1.0e-6) ^ forceGen.newField(), K.nn()
-            )*recRootN
-        );
-
-        #include "globalProperties.H"
-
-        fvVectorMatrix UEqn
-        (
-            fvm::ddt(U)
-          + fvm::div(phi, U)
-          - fvm::laplacian(nu, U)
-         ==
-            force
-        );
-
-        solve(UEqn == -fvc::grad(p));
-
+        #include "UEqn.H"
 
         // --- PISO loop
         while (piso.correct())
         {
-            volScalarField rAU(1.0/UEqn.A());
-            surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-            surfaceScalarField phiHbyA
-            (
-                "phiHbyA",
-                fvc::flux(HbyA)
-              + rAUf*fvc::ddtCorr(U, phi)
-            );
-
-            // Update the pressure BCs to ensure flux consistency
-            constrainPressure(p, U, phiHbyA, rAUf);
-
-            fvScalarMatrix pEqn
-            (
-                fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
-            );
-
-            pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
-
-            phi = phiHbyA - pEqn.flux();
-
-            #include "continuityErrs.H"
-
-            U = HbyA - rAU*fvc::grad(p);
-            U.correctBoundaryConditions();
+            #include "pEqn.H"
         }
 
         runTime.write();
@@ -149,4 +100,4 @@ int main(int argc, char *argv[])
 }
 
 
-// ************************************************************************* //
+// ************************************************************************* //
\ No newline at end of file

applications/solvers/DNS/dnsFoam/pEqn.H

+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    fvc::flux(HbyA)
+  + rAUf*fvc::ddtCorr(U, phi)
+);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, U, phiHbyA, rAUf);
+
+fvScalarMatrix pEqn
+(
+    fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
+);
+
+pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+
+phi = phiHbyA - pEqn.flux();
+
+#include "continuityErrs.H"
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();

applications/solvers/DNS/dnsFoam/readTransportProperties.H

-    Info<< "Reading transportProperties\n" << endl;
-
-    IOdictionary transportProperties
-    (
-        IOobject
-        (
-            "transportProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar nu
-    (
-        "nu",
-        dimViscosity,
-        transportProperties
-    );

applications/solvers/DNS/dnsFoam/readTurbulenceProperties.H

-    Info<< "Reading turbulenceProperties\n" << endl;
-
-    IOdictionary turbulenceProperties
-    (
-        IOobject
-        (
-            "turbulenceProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    );
-
-    volVectorField force
-    (
-        U/dimensionedScalar("dt", dimTime, runTime.deltaTValue())
-    );
-
-    Kmesh K(mesh);
-    UOprocess forceGen(K, runTime.deltaTValue(), turbulenceProperties);
-