...
 
Commits (1)
  • Kutalmis Bercin's avatar
    ENH: Refactor dnsFoam · a3506b91
    Kutalmis Bercin authored
        - If applied: This commit will restructure dnsFoam in line with other
                      solvers, e.g. pisoFoam, without changing its external
                      behaviour.
    
        - Why: Prior to this change, dnsFoam structure reflected v.2.x and older
               solver style.
    
        - How: This change collects naked dnsFoam code parts under general file
               structure, e.g. via UEqn.H.
    
        - Verification: No change in any output of dnsFoam tutorial case.
    a3506b91
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));
#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));
......@@ -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
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();
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 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);