exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Merged cummulative developement branch (Hrvoje Jasak). Author: Hrvoje Jasak. Merge: Henrik Rusche

Browse source 
Conflicts:
	src/ODE/sixDOF/finiteRotation/finiteRotation.C
	src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C
	src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H
	src/dynamicMesh/meshMotion/solidBodyMotion/constantVelocity/constantVelocity.C
	src/dynamicMesh/meshMotion/solidBodyMotion/rotatingOscillation/rotatingOscillation.H
	src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.C
	src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C
	src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellMDLimitedGrad/cellMDLimitedGrads.C
	src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C
This commit is contained in:
Henrik Rusche 2016-05-24 12:30:46 +02:00
commit 3d1742c1f0
894 changed files with 29187 additions and 7182 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
ListOfContributors
View file

@ -86,3 +86,4 @@ Contents:
Inno Gatin
Alexey Matveichev
Vuko Vukcevic
Robert Keser

2
applications/solvers/DNS/dnsFoam/globalProperties.H
View file

@ -9,7 +9,7 @@
<< (
0.5*nu*average
(
magSqr(fvc::grad(U) + fvc::grad(U)().T())
magSqr(fvc::grad(U) + T(fvc::grad(U)))
)
).value() << endl;

2
applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
View file

@ -219,7 +219,7 @@ tmp<fvVectorMatrix> PDRkEpsilon::divDevRhoReff(volVectorField& U) const
{
return
(
- fvm::laplacian(muEff(), U) - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
- fvm::laplacian(muEff(), U) - fvc::div(muEff()*dev2(T(fvc::grad(U))))
);
}

4
applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.H
View file

@ -37,7 +37,7 @@ SourceFiles
#ifndef mixedFixedValueSlipFvPatchField_H
#define mixedFixedValueSlipFvPatchField_H
#include "transformFvPatchField.H"
#include "transformFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -45,7 +45,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class mixedFixedValueSlipFvPatch Declaration
Class mixedFixedValueSlipFvPatch Declaration
\*---------------------------------------------------------------------------*/
template<class Type>

2
applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
View file

@ -139,7 +139,7 @@ int main(int argc, char *argv[])
+ aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
+ aSf*p_pos - aSf*p_neg;
volTensorField tauMC("tauMC", mu*dev2(fvc::grad(U)().T()));
volTensorField tauMC("tauMC", mu*dev2(Foam::T(fvc::grad(U))));
// --- Solve density
solve(fvm::ddt(rho) + fvc::div(phi));

2
applications/solvers/compressible/rhopSonicFoam/rhopSonicFoam.C
View file

@ -137,8 +137,6 @@ int main(int argc, char *argv[])
)
+ HbyAblend*fvc::ddtPhiCorr(rrhoUA, rho, rhoU, phi);
p.boundaryField().updateCoeffs();
surfaceScalarField phiGradp =
rrhoUAf*mesh.magSf()*fvc::snGrad(p);

1
applications/solvers/compressible/sonicDyMFoam/eEqn.H
View file

@ -54,4 +54,5 @@
// Recalculate density
rho = thermo.rho();
rho.correctBoundaryConditions();
}

1
applications/solvers/compressible/sonicFoam/eEqn.H
View file

@ -39,4 +39,5 @@
// Recalculate density
rho = thermo.rho();
rho.correctBoundaryConditions();
}

7
applications/solvers/compressible/steadyCompressibleFoam/createFields.H
View file

@ -40,6 +40,13 @@
mesh
);
// Store kinetic energy for total enthalpy equation
volScalarField K
(
"K",
0.5*(magSqr(U))
);
# include "compressibleCreatePhi.H"
Info<< "Creating turbulence model\n" << endl;

15
applications/solvers/compressible/steadyCompressibleFoam/hEqn.H
View file

@ -1,23 +1,16 @@
{
// Solve the enthalpy equation
T.storePrevIter();
// Solve the enthalpy equation in total enthalpy formulation (see K)
// Calculate face velocity from flux
surfaceScalarField faceU
(
"faceU",
phi/fvc::interpolate(rho)
);
K = 0.5*(magSqr(U));
fvScalarMatrix hEqn
(
fvm::ddt(rho, h)
+ fvc::ddt(rho, K)
+ fvm::div(phi, h)
+ fvc::div(phi, K)
- fvm::laplacian(turbulence->alphaEff(), h)
==
// ddt(p) term removed: steady-state. HJ, 27/Apr/2010
fvc::div(faceU, p, "div(U,p)")
- p*fvc::div(faceU)
// Viscous heating: note sign (devRhoReff has a minus in it)
- (turbulence->devRhoReff() && fvc::grad(U))
);

9
applications/solvers/compressible/steadyCompressibleFoam/pEqn.H
View file

@ -21,12 +21,21 @@
p.storePrevIter();
volScalarField divPhid
(
"divPhid",
fvc::div(phid)
);
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::ddt(psis, p)
+ fvm::div(phid, p)
// Convective flux relaxation terms
+ fvm::SuSp(-divPhid, p)
+ divPhid*p
+ fvc::div(phid2)
- fvm::laplacian(rho*rUA, p)
);

1
applications/solvers/compressible/steadyCompressibleFoam/rhoFromP.H
View file

@ -12,4 +12,5 @@
rho = Foam::min(rho, rhoMax);
rho = Foam::max(rho, rhoMin);
rho.relax();
rho.correctBoundaryConditions();
}

9
applications/solvers/compressible/steadyCompressibleMRFFoam/createFields.H
View file

@ -40,6 +40,13 @@
mesh
);
// Store kinetic energy for total enthalpy equation
volScalarField K
(
"K",
0.5*(magSqr(U))
);
# include "compressibleCreatePhi.H"
Info<< "Creating turbulence model\n" << endl;
@ -93,4 +100,4 @@
),
h
);
i -= 0.5*magSqr(Urot);
i == h - 0.5*(magSqr(Urot) - magSqr(Urel));

7
applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H
View file

@ -2,7 +2,12 @@
// Solve the enthalpy equation
T.storePrevIter();
surfaceScalarField faceU = phi/fvc::interpolate(rho);
// Calculate face velocity from flux
surfaceScalarField faceU
(
"faceU",
phi/fvc::interpolate(rho)
);
fvScalarMatrix hEqn
(

29
applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H
View file

@ -1,33 +1,18 @@
{
// Solve the rothalpy equation
T.storePrevIter();
// Create relative velocity
Urel == U;
mrfZones.relativeVelocity(Urel);
// Create rotational velocity (= omega x r)
Urot = U - Urel;
// Calculate face velocity from absolute flux
surfaceScalarField rhof = fvc::interpolate(rho);
surfaceScalarField phiAbs
(
"phiAbs",
phi
);
mrfZones.absoluteFlux(rhof, phiAbs);
surfaceScalarField faceU("faceU", phiAbs/rhof);
Urot == U - Urel;
fvScalarMatrix iEqn
(
fvm::ddt(rho, i)
+ fvm::div(phi, i)
- fvm::laplacian(turbulence->alphaEff(), i)
// u & gradP term (steady-state formulation)
+ fvm::SuSp((fvc::div(faceU, p, "div(U,p)") - fvc::div(faceU)*p)/i, i)
==
// Viscous heating: note sign (devRhoReff has a minus in it)
- (turbulence->devRhoReff() && fvc::grad(Urel))
@ -39,19 +24,9 @@
maxResidual = max(eqnResidual, maxResidual);
// Calculate enthalpy out of rothalpy
h = i + 0.5*magSqr(Urot);
h = i + 0.5*(magSqr(Urot) - magSqr(Urel));
h.correctBoundaryConditions();
// Bound the enthalpy using TMin and TMax
volScalarField Cp = thermo.Cp();
h = Foam::min(h, TMax*Cp);
h = Foam::max(h, TMin*Cp);
h.correctBoundaryConditions();
// Re-initialise rothalpy based on limited enthalpy
i = h - 0.5*magSqr(Urot);
thermo.correct();
psis = thermo.psi()/thermo.Cp()*thermo.Cv();
}

12
applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H
View file

@ -4,9 +4,10 @@
surfaceScalarField psisf = fvc::interpolate(psis);
surfaceScalarField rhof = fvc::interpolate(rho);
// Needs to be outside of loop since p is changing, but psi and rho are not.
// Needs to be outside of loop since p is changing, but psi and rho are not
surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p);
// --- PISO loop
for (int corr = 0; corr < nCorr; corr++)
{
U = rUA*UEqn.H();
@ -25,12 +26,21 @@
p.storePrevIter();
volScalarField divPhid
(
"divPhid",
fvc::div(phid)
);
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::ddt(psis, p)
+ fvm::div(phid, p)
// Convective flux relaxation terms
+ fvm::SuSp(-divPhid, p)
+ divPhid*p
+ fvc::div(phid2)
- fvm::laplacian(rho*rUA, p)
);

1
applications/solvers/compressible/steadyCompressibleMRFFoam/rhoFromP.H
View file

@ -12,4 +12,5 @@
rho = Foam::min(rho, rhoMax);
rho = Foam::max(rho, rhoMin);
rho.relax();
rho.correctBoundaryConditions();
}

7
applications/solvers/compressible/steadyCompressibleSRFFoam/createFields.H
View file

@ -73,6 +73,9 @@
SRF::SRFModel::New(Urel)
);
// Create rotational velocity
volVectorField Urot("Urot", SRF->U());
// Create absolute velocity field
volVectorField Uabs
(
@ -84,7 +87,7 @@
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
Urel + SRF->U()
Urel + Urot
);
// Create rothalpy, in two steps to preserve boundary conditions
@ -100,4 +103,4 @@
),
h
);
i -= 0.5*magSqr(SRF->U());
i == h - 0.5*(magSqr(Urot) - magSqr(Urel));

19
applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H
View file

@ -1,13 +1,5 @@
{
// Solve the rothalpy equation
T.storePrevIter();
// Calculate face velocity from flux
surfaceScalarField faceU
(
"faceU",
phi/fvc::interpolate(rho) + (SRF->faceU() & mesh.Sf())
);
fvScalarMatrix iEqn
(
@ -15,9 +7,6 @@
+ fvm::div(phi, i)
- fvm::laplacian(turbulence->alphaEff(), i)
==
// ddt(p) term removed: steady-state. HJ, 27/Apr/2010
fvc::div(faceU, p, "div(U,p)")
- p*fvc::div(faceU)
// Viscous heating: note sign (devRhoReff has a minus in it)
- (turbulence->devRhoReff() && fvc::grad(Urel))
);
@ -28,15 +17,9 @@
maxResidual = max(eqnResidual, maxResidual);
// Calculate enthalpy out of rothalpy
volVectorField Urot("Urot", SRF->U());
h = i + 0.5*magSqr(Urot);
h = i + 0.5*(magSqr(Urot) - magSqr(Urel));
h.correctBoundaryConditions();
// Re-initialise rothalpy based on limited enthalpy
i = h - 0.5*magSqr(Urot);
// Bounding of enthalpy taken out
thermo.correct();
psis = thermo.psi()/thermo.Cp()*thermo.Cv();
}

12
applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H
View file

@ -4,9 +4,10 @@
surfaceScalarField psisf = fvc::interpolate(psis);
surfaceScalarField rhof = fvc::interpolate(rho);
// Needs to be outside of loop since p is changing, but psi and rho are not.
// Needs to be outside of loop since p is changing, but psi and rho are not
surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p);
// --- PISO loop
for (int corr = 0; corr < nCorr; corr++)
{
Urel = rUrelA*UrelEqn.H();
@ -20,12 +21,21 @@
p.storePrevIter();
volScalarField divPhid
(
"divPhid",
fvc::div(phid)
);
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::ddt(psis, p)
+ fvm::div(phid, p)
// Convective flux relaxation terms
+ fvm::SuSp(-divPhid, p)
+ divPhid*p
+ fvc::div(phid2)
- fvm::laplacian(rho*rUrelA, p)
);

1
applications/solvers/compressible/steadyCompressibleSRFFoam/rhoFromP.H
View file

@ -12,4 +12,5 @@
rho = Foam::min(rho, rhoMax);
rho = Foam::max(rho, rhoMin);
rho.relax();
rho.correctBoundaryConditions();
}

11
applications/solvers/compressible/steadyCompressibleSRFFoam/steadyCompressibleSRFFoam.C
View file

@ -65,15 +65,22 @@ int main(int argc, char *argv[])
# include "initConvergenceCheck.H"
# include "UEqn.H"
# include "iEqn.H"
# include "pEqn.H"
// # include "hEqn.H"
// Solving for rothalpy
# include "iEqn.H"
# include "rhoFromP.H"
// Correct turbulence
turbulence->correct();
Uabs = Urel + SRF->U();
// Update rotational velocity
Urot = SRF->U();
// Update absolute velocity
Uabs = Urel + Urot;
runTime.write();

2
applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C
View file

@ -90,7 +90,7 @@ int main(int argc, char *argv[])
U.correctBoundaryConditions();
p.correctBoundaryConditions();
phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource;
phi = (fvc::interpolate(U) & mesh.Sf()) + tpEqn().flux() + tpresSource;
// Make flux relative in rotating zones
mrfZones.relativeFlux(phi);

97
applications/solvers/coupled/MRFPorousFoam/UEqn.H
View file

@ -1,15 +1,94 @@
// Momentum equation
tmp<fvVectorMatrix> UEqn
{
volScalarField divPhi
(
fvm::ddt(U)
+ fvm::div(phi, U)
"divPhi",
fvc::div(phi)
);
// Momentum equation
fvVectorMatrix UEqn
(
fvm::div(phi, U)
+ turbulence->divDevReff(U)
);
// Add MRF and porous sources
mrfZones.addCoriolis(UEqn());
pZones.addResistance(UEqn());
// Add MRF sources
mrfZones.addCoriolis(UEqn);
UEqn().relax();
// Add porous sources
tmp<volTensorField> tTU;
UpEqn.insertEquation(0, UEqn());
if (addPorosity)
{
tTU = tmp<volTensorField>
(
new volTensorField
(
IOobject
(
"TU",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedTensor("zero", dimless/dimTime, tensor::zero)
)
);
volTensorField& TU = tTU();
pZones.addResistance(UEqn, TU);
trTU = inv(TU + tensor(I)*UEqn.A());
trTU().rename("rAU");
}
else
{
trAU = 1.0/UEqn.A();
trAU().rename("rAU");
}
// Insert the additional components. Note this will destroy the H and A
UEqn += fvm::SuSp(-divPhi, U) + divPhi*U;
UEqn.relax();
// Insert momentum equation
UpEqn.insertEquation(0, UEqn);
if (addPorosity)
{
// Manually over-ride the 3x3 block to handle the off-diagonal
// part of the Ap coefficient
const tensorField& TUIn = tTU().internalField();
CoeffField<vector4>::squareTypeField& DD = UpEqn.diag().asSquare();
const scalarField& V = mesh.V().field();
// Note: this insertion should only happen in porous cell zones
// A rewrite of the porousZones class interface is needed.
// HJ, 14/Mar/2016
forAll (TUIn, cellI)
{
const scalar& cellV = V[cellI];
const tensor& cellTU = TUIn[cellI];
CoeffField<vector4>::squareType& cellDD = DD[cellI];
cellDD(0, 0) += cellV*cellTU.xx();
cellDD(0, 1) += cellV*cellTU.xy();
cellDD(0, 2) += cellV*cellTU.xz();
cellDD(1, 0) += cellV*cellTU.yx();
cellDD(1, 1) += cellV*cellTU.yy();
cellDD(2, 2) += cellV*cellTU.yz();
cellDD(2, 0) += cellV*cellTU.zx();
cellDD(2, 1) += cellV*cellTU.zy();
cellDD(2, 2) += cellV*cellTU.zz();
}
}
}

6
applications/solvers/coupled/MRFPorousFoam/createFields.H
View file

@ -50,3 +50,9 @@ volVector4Field Up
mesh,
dimensionedVector4("zero", dimless, vector4::zero)
);
Info<< "Creating field rAU\n" << endl;
// Note: depending on the porosity treatment, only one of the rAU fields
// will be used. HJ, 1/Mar/2016
tmp<volScalarField> trAU;
tmp<volTensorField> trTU;

7
applications/solvers/coupled/MRFPorousFoam/createPorous.H
View file

@ -1 +1,8 @@
porousZones pZones(mesh);
bool addPorosity = false;
if (!pZones.empty())
{
addPorosity = true;
}

61
applications/solvers/coupled/MRFPorousFoam/pEqn.H
View file

@ -1,25 +1,48 @@
// Pressure parts of the continuity equation
surfaceScalarField rUAf
(
"rUAf",
fvc::interpolate(1.0/UEqn().A())
);
tmp<fvScalarMatrix> tpEqn;
tmp<surfaceScalarField> tpresSource;
UEqn.clear();
if (addPorosity)
{
// Collect pressure source with tensorial 1/Ap
const volTensorField& rTU = trTU();
surfaceScalarField presSource
(
"presSource",
rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf())
);
tpresSource =
(
(mesh.Sf() & fvc::interpolate(rTU))
& fvc::interpolate(fvc::grad(p, "grad(pSource)"))
);
const surfaceScalarField& presSource = tpresSource();
fvScalarMatrix pEqn
(
- fvm::laplacian(rUAf, p)
==
- fvc::div(presSource)
);
// Assemble pressure matrix with tensorial 1/Ap
tpEqn =
(
- fvm::laplacian(rTU, p)
==
- fvc::div(presSource)
);
}
else
{
// Collect pressure source with scalar 1/Ap
const volScalarField& rAU = trAU();
pEqn.setReference(pRefCell, pRefValue);
tpresSource =
(
fvc::interpolate(rAU)*
(mesh.Sf() & fvc::interpolate(fvc::grad(p, "grad(pSource)")))
);
const surfaceScalarField& presSource = tpresSource();
UpEqn.insertEquation(3, pEqn);
// Assemble pressure matrix with tensorial 1/Ap
tpEqn =
(
- fvm::laplacian(rAU, p)
==
- fvc::div(presSource)
);
}
tpEqn().setReference(pRefCell, pRefValue);
UpEqn.insertEquation(3, tpEqn());

3
applications/solvers/coupled/conjugateHeatFoam/Make/options
View file

@ -19,4 +19,5 @@ EXE_LIBS = \
-lradiation \
-lbasicThermophysicalModels \
-lspecie \
-lconjugateHeatTransfer
-lconjugateHeatTransfer \
-llduSolvers

1
applications/solvers/coupled/conjugateHeatFoam/conjugateHeatFoam.C
View file

@ -109,6 +109,7 @@ int main(int argc, char *argv[])
// Update density according to Boussinesq approximation
rhok = 1.0 - beta*(T - TRef);
rhok.correctBoundaryConditions();
runTime.write();

3
applications/solvers/coupled/conjugateHeatSimpleFoam/Make/options
View file

@ -19,4 +19,5 @@ EXE_LIBS = \
-lradiation \
-lbasicThermophysicalModels \
-lspecie \
-lconjugateHeatTransfer
-lconjugateHeatTransfer \
-llduSolvers

9
applications/solvers/coupled/conjugateHeatSimpleFoam/conjugateHeatSimpleFoam.C
View file

@ -81,22 +81,23 @@ int main(int argc, char *argv[])
// Update thermal conductivity in the solid
solidThermo.correct();
ksolid = solidThermo.k();
kSolid = solidThermo.k();
// Coupled patches
# include "attachPatches.H"
kappaEff.correctBoundaryConditions();
ksolid.correctBoundaryConditions();
kSolid.correctBoundaryConditions();
// Interpolate to the faces and add thermal resistance
surfaceScalarField ksolidf = fvc::interpolate(ksolid);
solidThermo.modifyResistance(ksolidf);
surfaceScalarField kSolidf = fvc::interpolate(kSolid);
solidThermo.modifyResistance(kSolidf);
# include "solveEnergy.H"
// Update density according to Boussinesq approximation
rhok = 1.0 - beta*(T - TRef);
rhok.correctBoundaryConditions();
runTime.write();

4
applications/solvers/coupled/conjugateHeatSimpleFoam/createFields.H
View file

@ -100,7 +100,9 @@
(
"rhok",
runTime.timeName(),
mesh
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
1.0 - beta*(T - TRef)
);

2
applications/solvers/coupled/conjugateHeatSimpleFoam/createSolidFields.H
View file

@ -16,7 +16,7 @@
thermalModel solidThermo(Tsolid);
Info<< "Reading solid diffusivity k\n" << endl;
volScalarField ksolid
volScalarField kSolid
(
IOobject
(

7
applications/solvers/coupled/conjugateHeatSimpleFoam/solveEnergy.H
View file

@ -20,21 +20,16 @@
+ 3.0*radiation->Rp()*pow4(T)
);
//Called automatically in 1.6.1-ext
//TFluidEqn->boundaryManipulate(T.boundaryField());
TFluidEqn->relax();
fvScalarMatrix* TSolidEqn = new fvScalarMatrix
(
- fvm::laplacian(ksolidf, Tsolid, "laplacian(k,T)")
- fvm::laplacian(kSolidf, Tsolid, "laplacian(k,T)")
+ fvm::SuSp(-solidThermo.S()/Tsolid, Tsolid)
);
//Called automatically in 1.6.1-ext
//TSolidEqn->boundaryManipulate(Tsolid.boundaryField());
TSolidEqn->relax();
// Add fluid equation
TEqns.set(0, TFluidEqn);

2
applications/solvers/coupled/conjugateHeatTransfer/thermalModel/thermalModel.C
View file

@ -151,7 +151,7 @@ tmp<volScalarField> thermalModel::S() const
(
"zero",
dimEnergy/dimTime/dimVolume,
scalar(0.0)
scalar(0)
)
)
);

36
applications/solvers/coupled/conjugateHeatTransfer/thermalModel/thermalSource/constantThermalSource/constantThermalSource.C
View file

@ -31,7 +31,12 @@ License
namespace Foam
{
defineTypeNameAndDebug(constantThermalSource, 0);
addToRunTimeSelectionTable(thermalSource, constantThermalSource, dictionary);
addToRunTimeSelectionTable
(
thermalSource,
constantThermalSource,
dictionary
);
}
@ -78,9 +83,34 @@ void Foam::constantThermalSource::addSource(volScalarField& source) const
const labelList& cells = mesh().cellZones()[zoneID];
forAll(cells, cellI)
// Sum up the volumes
scalar sumVolume = 0;
const scalarField& V = mesh().V().field();
forAll (cells, cellI)
{
source[cells[cellI]] = S_.value();
sumVolume += V[cells[cellI]];
}
reduce(sumVolume, sumOp<scalar>());
if (sumVolume < SMALL)
{
FatalErrorIn
(
"constantThermalSource::addSourcex()\n"
) << "Zone " << zones_[zoneI]
<< " specified in source " << name()
<< " has zero volume: " << sumVolume
<< abort(FatalError);
}
const scalar SoverV = S_.value()/sumVolume;
forAll (cells, cellI)
{
source[cells[cellI]] = SoverV;
}
}
}

9
applications/solvers/coupled/conjugateHeatTransfer/thermalModel/thermalSource/constantThermalSource/constantThermalSource.H
View file

@ -25,7 +25,9 @@ Class
constantThermalSource
Description
Constant thermal properties
Constant thermal source, given in Watts for the component
consisting of face zones. Therefore, S from dictionary is divided by
total volume of the component
Author
Henrik Rusche, Wikki GmbH. All rights reserved.
@ -46,7 +48,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class constantThermalSource Declaration
Class constantThermalSource Declaration
\*---------------------------------------------------------------------------*/
class constantThermalSource
@ -58,9 +60,10 @@ class constantThermalSource
//- Strengh of the source
dimensionedScalar S_;
//- list of cell zones
//- List of cell zones
const wordList zones_;
// Private Member Functions
//- Disallow default bitwise copy construct

30
applications/solvers/coupled/pUCoupledFoam/UEqn.H
View file

@ -1,11 +1,23 @@
// Momentum equation
tmp<fvVectorMatrix> UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
);
{
volScalarField divPhi
(
"divPhi",
fvc::div(phi)
);
UEqn().relax();
// Momentum equation
fvVectorMatrix UEqn
(
fvm::div(phi, U)
+ turbulence->divDevReff(U)
);
UpEqn.insertEquation(0, UEqn());
rAU = 1.0/UEqn.A();
// Insert the additional components. Note this will destroy the H and A
UEqn += fvm::SuSp(-divPhi, U) + divPhi*U;
UEqn.relax();
UpEqn.insertEquation(0, UEqn);
}

15
applications/solvers/coupled/pUCoupledFoam/createFields.H
View file

@ -50,3 +50,18 @@ volVector4Field Up
mesh,
dimensionedVector4("zero", dimless, vector4::zero)
);
Info<< "Creating field rAU\n" << endl;
volScalarField rAU
(
IOobject
(
"rAU",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
runTime.deltaT()
);

13
applications/solvers/coupled/pUCoupledFoam/pEqn.H
View file

@ -1,21 +1,14 @@
// Pressure parts of the continuity equation
surfaceScalarField rUAf
(
"rUAf",
fvc::interpolate(1.0/UEqn().A())
);
UEqn.clear();
surfaceScalarField presSource
(
"presSource",
rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf())
fvc::interpolate(rAU)*
(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf())
);
fvScalarMatrix pEqn
(
- fvm::laplacian(rUAf, p)
- fvm::laplacian(rAU, p)
==
- fvc::div(presSource)
);

5
applications/solvers/immersedBoundary/porousSimpleIbFoam/pEqn.H
View file

@ -1,5 +1,4 @@
p.boundaryField().updateCoeffs();
{
volScalarField AU = UEqn().A();
U = UEqn().H()/AU;
@ -48,4 +47,4 @@
// Momentum corrector
U -= fvc::grad(p)/AU;
U.correctBoundaryConditions();
}

5
applications/solvers/immersedBoundary/simpleIbFoam/pEqn.H
View file

@ -1,5 +1,4 @@
p.boundaryField().updateCoeffs();
{
volScalarField AU = UEqn().A();
U = UEqn().H()/AU;
@ -48,4 +47,4 @@
// Momentum corrector
U -= fvc::grad(p)/AU;
U.correctBoundaryConditions();
}

1
applications/solvers/incompressible/MRFSimpleFoam/MRFSimpleFoam.C
View file

@ -74,7 +74,6 @@ int main(int argc, char *argv[])
solve(UEqn() == -fvc::grad(p));
p.boundaryField().updateCoeffs();
volScalarField rAU = 1.0/UEqn().A();
U = rAU*UEqn().H();
UEqn.clear();

3
applications/solvers/incompressible/consistentIcoFoam/Make/files Normal file
View file

@ -0,0 +1,3 @@
consistentIcoFoam.C
EXE = $(FOAM_APPBIN)/consistentIcoFoam

6
applications/solvers/incompressible/consistentIcoFoam/Make/options Normal file
View file

@ -0,0 +1,6 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-llduSolvers

125
applications/solvers/incompressible/consistentIcoFoam/consistentIcoFoam.C Normal file
View file

@ -0,0 +1,125 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 3.2
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Application
consistentIcoFoam
Description
Transient solver for incompressible, laminar flow of Newtonian fluids.
Consistent formulation without time-step and relaxation dependence by Jasak
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readPISOControls.H"
# include "CourantNo.H"
// Convection-diffusion matrix
fvVectorMatrix HUEqn
(
fvm::div(phi, U)
- fvm::laplacian(nu, U)
);
// ddt matrix
fvVectorMatrix ddtUEqn(fvm::ddt(U));
solve(ddtUEqn + HUEqn == -fvc::grad(p));
// Prepare clean Ap without time derivative contribution
// HJ, 26/Oct/2015
volScalarField aU = HUEqn.A();
// --- PISO loop
for (int corr = 0; corr < nCorr; corr++)
{
U = HUEqn.H()/aU;
phi = (fvc::interpolate(U) & mesh.Sf());
adjustPhi(phi, U, p);
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(1/aU, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
# include "continuityErrs.H"
// Note: cannot call H(U) here because the velocity is not complete
// HJ, 22/Jan/2016
U = 1.0/(aU + ddtUEqn.A())*
(
U*aU - fvc::grad(p) + ddtUEqn.H()
);
U.correctBoundaryConditions();
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

55
applications/solvers/incompressible/consistentIcoFoam/createFields.H Normal file
View file

@ -0,0 +1,55 @@
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar nu
(
transportProperties.lookup("nu")
);
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);

3
applications/solvers/incompressible/consistentPimpleDyMFoam/Make/files Normal file
View file

@ -0,0 +1,3 @@
consistentPimpleDyMFoam.C
EXE = $(FOAM_APPBIN)/consistentPimpleDyMFoam

21
applications/solvers/incompressible/consistentPimpleDyMFoam/Make/options Normal file
View file

@ -0,0 +1,21 @@
EXE_INC = \
-I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-ldynamicFvMesh \
-ltopoChangerFvMesh \
-ldynamicMesh \
-lmeshTools \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume \
-llduSolvers \
-L$(MESQUITE_LIB_DIR) -lmesquite

33
applications/solvers/incompressible/consistentPimpleDyMFoam/UEqn.H Normal file
View file

@ -0,0 +1,33 @@
// Convection-diffusion matrix
fvVectorMatrix HUEqn
(
fvm::div(phi, U)
+ turbulence->divDevReff(U)
);
// ddt matrix
fvVectorMatrix ddtUEqn(fvm::ddt(U));
// Get under-relaxation factor and under-relax the diagonal
scalar UUrf = mesh.solutionDict().relaxationFactor(U.name());
if (oCorr == nOuterCorr - 1)
{
if (mesh.solutionDict().relax("UFinal"))
{
UUrf = mesh.solutionDict().relaxationFactor("UFinal");
}
else
{
UUrf = 1;
}
}
// Solve momentum predictor
solve
(
ddtUEqn
+ relax(HUEqn, UUrf)
==
-fvc::grad(p)
);

128
applications/solvers/incompressible/consistentPimpleDyMFoam/consistentPimpleDyMFoam.C Normal file
View file

@ -0,0 +1,128 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 3.2
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Application
pimpleDyMFoam.C
Description
Transient solver for incompressible, flow of Newtonian fluids
with dynamic mesh using the PIMPLE (merged PISO-SIMPLE) algorithm.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulenceModel.H"
#include "dynamicFvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createDynamicFvMesh.H"
# include "readPIMPLEControls.H"
# include "initContinuityErrs.H"
# include "createFields.H"
# include "readTimeControls.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
# include "readControls.H"
# include "CourantNo.H"
# include "setDeltaT.H"
// Make the fluxes absolute
fvc::makeAbsolute(phi, U);
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
bool meshChanged = mesh.update();
# include "volContinuity.H"
if (checkMeshCourantNo)
{
# include "meshCourantNo.H"
}
// Mesh motion update
if (correctPhi && meshChanged)
{
// Fluxes will be corrected to absolute velocity
// HJ, 6/Feb/2009
# include "correctPhi.H"
}
if (meshChanged)
{
# include "CourantNo.H"
}
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
// --- PIMPLE loop
label oCorr = 0;
do
{
if (nOuterCorr != 1)
{
p.storePrevIter();
}
# include "UEqn.H"
// --- PISO loop
for (int corr = 0; corr < nCorr; corr++)
{
# include "pEqn.H"
}
turbulence->correct();
} while (++oCorr < nOuterCorr);
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

58
applications/solvers/incompressible/consistentPimpleDyMFoam/correctPhi.H Normal file
View file

@ -0,0 +1,58 @@
{
# include "continuityErrs.H"
wordList pcorrTypes
(
p.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
for (label i = 0; i < p.boundaryField().size(); i++)
{
if (p.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}
volScalarField pcorr
(
IOobject
(
"pcorr",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("pcorr", p.dimensions(), 0),
pcorrTypes
);
// Initialise flux with interpolated velocity
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, pcorr);
for(int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(1/aU, pcorr) == fvc::div(phi)
);
pcorrEqn.setReference(pRefCell, pRefValue);
pcorrEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phi -= pcorrEqn.flux();
}
// Fluxes are corrected to absolute velocity and further corrected
// later. HJ, 6/Feb/2009
}
# include "continuityErrs.H"
}

58
applications/solvers/incompressible/consistentPimpleDyMFoam/createFields.H Normal file
View file

@ -0,0 +1,58 @@
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
Info<< "Reading field aU if present\n" << endl;
volScalarField aU
(
IOobject
(
"aU",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
1/runTime.deltaT(),
zeroGradientFvPatchScalarField::typeName
);

69
applications/solvers/incompressible/consistentPimpleDyMFoam/pEqn.H Normal file
View file

@ -0,0 +1,69 @@
{
// Prepare clean Ap without time derivative contribution and
// without contribution from under-relaxation
// HJ, 26/Oct/2015
aU = HUEqn.A();
// Store velocity under-relaxation point before using U for
// the flux precursor
U.storePrevIter();
U = HUEqn.H()/aU;
phi = (fvc::interpolate(U) & mesh.Sf());
// Global flux balance
adjustPhi(phi, U, p);
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(1/aU, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
if
(
// oCorr == nOuterCorr - 1
corr == nCorr - 1
&& nonOrth == nNonOrthCorr
)
{
pEqn.solve
(
mesh.solutionDict().solver(p.name() + "Final")
);
}
else
{
pEqn.solve(mesh.solutionDict().solver(p.name()));
}
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
// Explicitly relax pressure for momentum corrector
if (oCorr != nOuterCorr - 1)
{
p.relax();
}
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
# include "movingMeshContinuityErrs.H"
U = UUrf*
(
1.0/(aU + ddtUEqn.A())*
(
U*aU - fvc::grad(p) + ddtUEqn.H()
)
)
+ (1 - UUrf)*U.prevIter();
U.correctBoundaryConditions();
}

14
applications/solvers/incompressible/consistentPimpleDyMFoam/readControls.H Normal file
View file

@ -0,0 +1,14 @@
# include "readTimeControls.H"
# include "readPIMPLEControls.H"
bool correctPhi = false;
if (pimple.found("correctPhi"))
{
correctPhi = Switch(pimple.lookup("correctPhi"));
}
bool checkMeshCourantNo = false;
if (pimple.found("checkMeshCourantNo"))
{
checkMeshCourantNo = Switch(pimple.lookup("checkMeshCourantNo"));
}

3
applications/solvers/incompressible/consistentSimpleFoam/Make/files Normal file
View file

@ -0,0 +1,3 @@
consistentSimpleFoam.C
EXE = $(FOAM_APPBIN)/consistentSimpleFoam

12
applications/solvers/incompressible/consistentSimpleFoam/Make/options Normal file
View file

@ -0,0 +1,12 @@
EXE_INC = -g \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lincompressibleRASModels \
-lincompressibleTransportModels \
-lfiniteVolume \
-llduSolvers

20
applications/solvers/incompressible/consistentSimpleFoam/UEqn.H Normal file
View file

@ -0,0 +1,20 @@
// Solve the momentum equation
tmp<fvVectorMatrix> HUEqn
(
fvm::div(phi, U)
+ turbulence->divDevReff(U)
);
// Get under-relaxation factor and under-relax the diagonal
const scalar UUrf = mesh.solutionDict().relaxationFactor(U.name());
// Momentum solution
eqnResidual = solve
(
relax(HUEqn(), UUrf)
==
-fvc::grad(p)
).initialResidual();
maxResidual = max(eqnResidual, maxResidual);

86
applications/solvers/incompressible/consistentSimpleFoam/consistentSimpleFoam.C Normal file
View file

@ -0,0 +1,86 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 3.2
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Application
consistentSimpleFoam
Description
Steady-state solver for incompressible, turbulent flow
Consistent formulation without time-step and relaxation dependence by Jasak
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "RASModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSIMPLEControls.H"
# include "initConvergenceCheck.H"
p.storePrevIter();
// Pressure-velocity SIMPLE corrector
{
# include "UEqn.H"
# include "pEqn.H"
}
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
# include "convergenceCheck.H"
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

9
applications/solvers/incompressible/consistentSimpleFoam/convergenceCheck.H Normal file
View file

@ -0,0 +1,9 @@
// check convergence
if (maxResidual < convergenceCriterion)
{
Info<< "reached convergence criterion: " << convergenceCriterion << endl;
runTime.writeAndEnd();
Info<< "latestTime = " << runTime.timeName() << endl;
}

42
applications/solvers/incompressible/consistentSimpleFoam/createFields.H Normal file
View file

@ -0,0 +1,42 @@
Info << "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info << "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::RASModel> turbulence
(
incompressible::RASModel::New(U, phi, laminarTransport)
);

7
applications/solvers/incompressible/consistentSimpleFoam/initConvergenceCheck.H Normal file
View file

@ -0,0 +1,7 @@
// initialize values for convergence checks
scalar eqnResidual = 1, maxResidual = 0;
scalar convergenceCriterion = 0;
simple.readIfPresent("convergence", convergenceCriterion);

58
applications/solvers/incompressible/consistentSimpleFoam/pEqn.H Normal file
View file

@ -0,0 +1,58 @@
{
// Prepare clean 1/Ap without contribution from under-relaxation
// HJ, 26/Oct/2015
volScalarField rUA
(
"(1|A(U))",
1/HUEqn().A()
);
// Store velocity under-relaxation point before using U for
// the flux precursor
U.storePrevIter();
U = rUA*HUEqn().H();
HUEqn.clear();
phi = fvc::interpolate(U) & mesh.Sf();
// Global flux balance
adjustPhi(phi, U, p);
// Non-orthogonal pressure corrector loop
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
// Retain the residual from the first iteration
if (nonOrth == 0)
{
eqnResidual = pEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
pEqn.solve();
}
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
# include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
// Momentum corrector
// Note: since under-relaxation does not change aU, H/a in U can be
// re-used. HJ, 22/Jan/2016
U = UUrf*(U - rUA*fvc::grad(p)) + (1 - UUrf)*U.prevIter();
U.correctBoundaryConditions();
}

6
applications/solvers/incompressible/simpleFoam/pEqn.H
View file

@ -1,9 +1,9 @@
p.boundaryField().updateCoeffs();
{
volScalarField AU = UEqn().A();
U = UEqn().H()/AU;
UEqn.clear();
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p);
// Non-orthogonal pressure corrector loop
@ -41,4 +41,4 @@
// Momentum corrector
U -= fvc::grad(p)/AU;
U.correctBoundaryConditions();
}

1
applications/solvers/incompressible/simpleSRFFoam/simpleSRFFoam.C
View file

@ -75,7 +75,6 @@ int main(int argc, char *argv[])
solve(UrelEqn() == -fvc::grad(p));
p.boundaryField().updateCoeffs();
volScalarField AUrel = UrelEqn().A();
Urel = UrelEqn().H()/AUrel;
UrelEqn.clear();

2
applications/solvers/multiphase/cavitatingFoam/pEqn.H
View file

@ -18,8 +18,6 @@
phiv = (fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phiv);
p.boundaryField().updateCoeffs();
surfaceScalarField phiGradp = rUAf*mesh.magSf()*fvc::snGrad(p);
phiv -= phiGradp/rhof;

2
applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H
View file

@ -11,7 +11,7 @@
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muEff, U)
- (fvc::grad(U) & fvc::grad(muEff))
//- fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T())))
- fvc::div(muEff*(mesh.Sf() & fvc::interpolate(T(fvc::grad(U)))))
);
UEqn.relax();

3
applications/solvers/multiphase/interFoam/alphaEqn.H
View file

@ -33,3 +33,6 @@
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
}
solve(fvm::ddt(rho) + fvc::div(rhoPhi));

2
applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
View file

@ -6,5 +6,5 @@
+ fvm::div(mixture.rhoPhi(), U)
- fvm::laplacian(muf, U)
- (fvc::grad(U) & fvc::grad(muf))
//- fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T())))
//- fvc::div(muf*(mesh.Sf() & fvc::interpolate(T(fvc::grad(U)))))
);

1
applications/solvers/viscoelastic/viscoelasticFluidFoam/viscoelasticFluidFoam.C
View file

@ -79,7 +79,6 @@ int main(int argc, char *argv[])
solve(UEqn() == -fvc::grad(p));
p.boundaryField().updateCoeffs();
volScalarField rUA = 1.0/UEqn().A();
U = rUA*UEqn().H();
UEqn.clear();

16
applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
View file

@ -250,6 +250,22 @@ Foam::label Foam::checkTopology
rs
);
ctr.write();
// Count number of cells in all regions
labelList nCellsInRegions(rs.nRegions(), 0);
forAll (rs, rsI)
{
nCellsInRegions[rs[rsI]]++;
}
Info<< "Nuumber of cells per region: " << nl;
forAll (nCellsInRegions, regionI)
{
Info<< tab << regionI << tab << nCellsInRegions[regionI] << nl;
}
Info<< endl;
}
}

3
applications/utilities/mesh/manipulation/splitMesh/splitMesh.C
View file

@ -263,6 +263,9 @@ int main(int argc, char *argv[])
splitter.changeMesh();
// Remove zones - not needed
mesh.removeZones();
Info<< "Writing mesh to " << runTime.timeName() << endl;
if (!mesh.write())
{

21
applications/utilities/mesh/manipulation/transformPoints/transformPoints.C
View file

@ -165,6 +165,14 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
if (args.options().empty())
{
FatalErrorIn(args.executable())
<< "No options supplied, please use one or more of "
"-translate, -rotate, -scale, or -cylToCart options."
<< exit(FatalError);
}
word regionName = polyMesh::defaultRegion;
fileName meshDir;
@ -191,14 +199,7 @@ int main(int argc, char *argv[])
)
);
if (args.options().empty())
{
FatalErrorIn(args.executable())
<< "No options supplied, please use one or more of "
"-translate, -rotate, -scale, or -cylToCart options."
<< exit(FatalError);
}
// Translation options
if (args.optionFound("translate"))
{
@ -209,6 +210,8 @@ int main(int argc, char *argv[])
points += transVector;
}
// Rotation options
if (args.optionFound("rotate"))
{
Pair<vector> n1n2(args.optionLookup("rotate")());
@ -298,7 +301,7 @@ int main(int argc, char *argv[])
}
}
// Scale options
if (args.optionFound("scale"))
{

4
applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C
View file

@ -122,8 +122,8 @@ int main(int argc, char *argv[])
Info<< "Patch " << patchi
<< " named " << currPatch.name()
<< " y+ : min: " << min(Yp) << " max: " << max(Yp)
<< " average: " << average(Yp) << nl << endl;
<< " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
<< " average: " << gAverage(Yp) << nl << endl;
}
}

12
applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C
View file

@ -91,8 +91,8 @@ void calcIncompressibleYPlus
Info<< "Patch " << patchi
<< " named " << nutPw.patch().name()
<< " y+ : min: " << min(Yp) << " max: " << max(Yp)
<< " average: " << average(Yp) << nl << endl;
<< " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
<< " average: " << gAverage(Yp) << nl << endl;
}
}
@ -171,8 +171,8 @@ void calcCompressibleYPlus
Info<< "Patch " << patchi
<< " named " << mutPw.patch().name()
<< " y+ : min: " << min(Yp) << " max: " << max(Yp)
<< " average: " << average(Yp) << nl << endl;
<< " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
<< " average: " << gAverage(Yp) << nl << endl;
}
}
@ -243,8 +243,8 @@ void calcTwoPhaseYPlus
Info<< "Patch " << patchi
<< " named " << nutPw.patch().name()
<< " y+ : min: " << min(Yp) << " max: " << max(Yp)
<< " average: " << average(Yp) << nl << endl;
<< " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
<< " average: " << gAverage(Yp) << nl << endl;
}
}

27
applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C
View file

@ -30,6 +30,8 @@ Description
- pitch (rotation about y) followed by
- yaw (rotation about z)
or -rotateAlongVector (vector and angle in degrees)
The yawPitchRoll does yaw followed by pitch followed by roll.
\*---------------------------------------------------------------------------*/
@ -42,11 +44,11 @@ Description
#include "transformField.H"
#include "Pair.H"
#include "quaternion.H"
#include "RodriguesRotation.H"
using namespace Foam;
using namespace Foam::mathematicalConstant;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
@ -59,6 +61,7 @@ int main(int argc, char *argv[])
argList::validArgs.append("output surface file");
argList::validOptions.insert("translate", "vector");
argList::validOptions.insert("rotate", "(vector vector)");
argList::validOptions.insert("rotateAlongVector", "(vector angleInDegree)");
argList::validOptions.insert("scale", "vector");
argList::validOptions.insert("rollPitchYaw", "(roll pitch yaw)");
argList::validOptions.insert("yawPitchRoll", "(yaw pitch roll)");
@ -72,7 +75,6 @@ int main(int argc, char *argv[])
Info<< "Writing surf to " << outFileName << " ..." << endl;
if (args.options().empty())
{
FatalErrorIn(args.executable())
@ -85,6 +87,8 @@ int main(int argc, char *argv[])
pointField points(surf1.points());
// Translation options
if (args.optionFound("translate"))
{
vector transVector(args.optionLookup("translate")());
@ -94,6 +98,8 @@ int main(int argc, char *argv[])
points += transVector;
}
// Rotation options
if (args.optionFound("rotate"))
{
Pair<vector> n1n2(args.optionLookup("rotate")());
@ -133,7 +139,6 @@ int main(int argc, char *argv[])
<< " pitch " << v.y() << nl
<< " roll " << v.z() << endl;
// Convert to radians
v *= pi/180.0;
@ -148,6 +153,22 @@ int main(int argc, char *argv[])
Info<< "Rotating points by quaternion " << R << endl;
points = transform(R, points);
}
else if (args.optionFound("rotateAlongVector"))
{
vector rotationAxis;
scalar rotationAngle;
args.optionLookup("rotateAlongVector")()
>> rotationAxis
>> rotationAngle;
tensor T = RodriguesRotation(rotationAxis, rotationAngle);
Info << "Rotating points by " << T << endl;
points = transform(T, points);
}
// Scale options
if (args.optionFound("scale"))
{

3
src/ODE/sixDOF/finiteRotation/finiteRotation.C
View file

@ -84,7 +84,8 @@ Foam::vector Foam::finiteRotation::eulerAngles(const tensor& rotT)
// Use mag to avoid negative value due to round-off
// HJ, 24/Feb/2016
const scalar c2 = sqrt(Foam::max(0, rotT.xx() + rotT.xy()));
// Bugfix: sqr. SS, 18/Apr/2016
const scalar c2 = sqrt(sqr(rotT.xx()) + sqr(rotT.xy()));
// Calculate pitch angle
pitchAngle = atan2(-rotT.xz(), c2);

174
src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C
View file

@ -207,6 +207,69 @@ void Foam::sixDOFqODE::constrainTranslation(vector& vec) const
}
void Foam::sixDOFqODE::aitkensRelaxation
(
const scalar min,
const scalar max
)
{
// Calculate translational relax factor
const scalar saveOldRelFacT = oldRelaxFactorT_;
oldRelaxFactorT_ = relaxFactorT_;
if(magSqr(A_[0] - An_[1] - A_[1] - An_[2]) > SMALL)
{
relaxFactorT_ = saveOldRelFacT + (saveOldRelFacT - 1)*
(
(A_[1] - An_[2])
& (A_[0] - An_[1] - A_[1] - An_[2])
)/
magSqr(A_[0] - An_[1] - A_[1] - An_[2]);
}
else
{
relaxFactorT_ = min;
}
// Bound the relaxation factor for stability
if(relaxFactorT_ > max)
{
relaxFactorT_ = max;
}
else if(relaxFactorT_ < min)
{
relaxFactorT_ = min;
}
// Calculate rotational relax factor
const scalar saveOldRelFacR = oldRelaxFactorR_;
oldRelaxFactorR_ = relaxFactorR_;
if(magSqr(OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2]) > SMALL)
{
relaxFactorR_ =
saveOldRelFacR + (saveOldRelFacR - 1)*
((OmegaDot_[1] - OmegaDotn_[2]) &
(OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2]))/
magSqr(OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2]);
}
else
{
relaxFactorR_ = min;
}
// Bound the relaxation factor for stability
if(relaxFactorR_ > max)
{
relaxFactorR_ = max;
}
else if(relaxFactorR_ < min)
{
relaxFactorR_ = min;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
@ -220,6 +283,10 @@ Foam::sixDOFqODE::sixDOFqODE(const IOobject& io)
Xequilibrium_(lookup("equilibriumPosition")),
linSpringCoeffs_(lookup("linearSpring")),
linDampingCoeffs_(lookup("linearDamping")),
relaxFactorT_(1.0),
relaxFactorR_(1.0),
oldRelaxFactorT_(1.0),
oldRelaxFactorR_(1.0),
Xrel_(lookup("Xrel")),
U_(lookup("U")),
@ -233,6 +300,10 @@ Foam::sixDOFqODE::sixDOFqODE(const IOobject& io)
omegaAverage_(omega_),
omegaAverageAbsolute_(omega_),
A_(3, vector::zero),
OmegaDot_(3, vector::zero),
An_(3, vector::zero),
OmegaDotn_(3, vector::zero),
force_(lookup("force")),
moment_(lookup("moment")),
forceRelative_(lookup("forceRelative")),
@ -286,6 +357,10 @@ Foam::sixDOFqODE::sixDOFqODE
Xequilibrium_(sd.Xequilibrium_.name(), sd.Xequilibrium_),
linSpringCoeffs_(sd.linSpringCoeffs_.name(), sd.linSpringCoeffs_),
linDampingCoeffs_(sd.linDampingCoeffs_.name(), sd.linDampingCoeffs_),
relaxFactorT_(sd.relaxFactorT_),
relaxFactorR_(sd.relaxFactorR_),
oldRelaxFactorT_(sd.oldRelaxFactorT_),
oldRelaxFactorR_(sd.oldRelaxFactorR_),
Xrel_(sd.Xrel_.name(), sd.Xrel_),
U_(sd.U_.name(), sd.U_),
@ -299,6 +374,10 @@ Foam::sixDOFqODE::sixDOFqODE
sd.omegaAverageAbsolute_
),
A_(sd.A_),
OmegaDot_(sd.OmegaDot_),
An_(sd.An_),
OmegaDotn_(sd.OmegaDotn_),
force_(sd.force_.name(), sd.force_),
moment_(sd.moment_.name(), sd.moment_),
forceRelative_(sd.forceRelative_.name(), sd.forceRelative_),
@ -327,6 +406,8 @@ Foam::sixDOFqODE::~sixDOFqODE()
void Foam::sixDOFqODE::setState(const sixDOFqODE& sd)
{
// Set state does not copy AList_, AOld_, relaxFactor_ and relaxFactorOld_ to
// allow Aitkens relaxation (IG 5/May/2016)
mass_ = sd.mass_;
momentOfInertia_ = sd.momentOfInertia_;
@ -388,23 +469,6 @@ void Foam::sixDOFqODE::derivatives
dydx[4] = accel.y();
dydx[5] = accel.z();
// // Add translational constraints by setting RHS of given components to zero
// if (fixedSurge_)
// {
// dydx[0] = 0; // Surge velocity
// dydx[3] = 0; // Surge acceleration
// }
// if (fixedSway_)
// {
// dydx[1] = 0; // Sway velocity
// dydx[4] = 0; // Sway acceleration
// }
// if (fixedHeave_)
// {
// dydx[2] = 0; // Heave velocity
// dydx[5] = 0; // Heave acceleration
// }
// Set the derivatives for rotation
dimensionedVector curOmega
(
@ -413,25 +477,6 @@ void Foam::sixDOFqODE::derivatives
vector(y[6], y[7], y[8])
);
// dimensionedVector curGlobalOmega = curRotation.invR() & curOmega;
//
// // Add rotational constraints by setting RHS of given components to zero
// if (fixedRoll_)
// {
// curGlobalOmega.value().x() = 0;
// }
// if (fixedPitch_)
// {
// curGlobalOmega.value().y() = 0;
// }
// if (fixedYaw_)
// {
// curGlobalOmega.value().z() = 0;
// }
//
//
// curOmega = curRotation.R() & curGlobalOmega;
const vector omegaDot = OmegaDot(curRotation, curOmega).value();
dydx[6] = omegaDot.x();
@ -543,6 +588,63 @@ void Foam::sixDOFqODE::update(const scalar delta)
}
void Foam::sixDOFqODE::relaxAcceleration
(
const scalar minRelFactor,
const scalar maxRelFactor
)
{
if (minRelFactor - maxRelFactor < SMALL)
{
// Fixed relaxation
relaxFactorT_ = minRelFactor;
relaxFactorR_ = minRelFactor;
}
else
{
// Use Aitkens relaxation
// Update Aitkens relaxation factor
aitkensRelaxation(minRelFactor, maxRelFactor);
// Update non relaxed accelerations
An_[1] = An_[2];
An_[2] = (force()/mass_).value();
OmegaDotn_[1] = OmegaDotn_[2];
OmegaDotn_[2] = (inv(momentOfInertia_) & moment()).value();
}
const dimensionedVector Aold
(
"",
dimensionSet(0, 1, -2, 0, 0, 0, 0),
A_[2]
);
const dimensionedVector OmegaDotold
(
"",
dimensionSet(0, 0, -2, 0, 0, 0, 0),
OmegaDot_[2]
);
force() =
Aold*mass_ + relaxFactorT_*(force() - Aold*mass_);
moment() =
(momentOfInertia_ & OmegaDotold)
+ relaxFactorR_*(moment() - (momentOfInertia_ & OmegaDotold));
// Update relaxed old accelerations
A_[0] = A_[1];
A_[1] = A_[2];
A_[2] = (force()/mass_).value();
OmegaDot_[0] = OmegaDot_[1];
OmegaDot_[1] = OmegaDot_[2];
OmegaDot_[2] = (inv(momentOfInertia_) & moment()).value();
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
bool Foam::sixDOFqODE::writeData(Ostream& os) const

32
src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H
View file

@ -84,6 +84,18 @@ class sixDOFqODE
//- Linear damping coeffs
dimensionedDiagTensor linDampingCoeffs_;
//- Translational relaxation factor
scalar relaxFactorT_;
//- Rotational relaxation factor
scalar relaxFactorR_;
//- Old translational relaxation factor
scalar oldRelaxFactorT_;
//- Old rotational relaxation factor
scalar oldRelaxFactorR_;
// Initial body state variables
@ -114,6 +126,15 @@ class sixDOFqODE
// External forces
//- Accelerations from previous iterations
// A_[2] is the old value, A_[1] old old, A_[0] old old old
List<vector> A_;
List<vector> OmegaDot_;
//- Previos iteration non relaxed accelerations
List<vector> An_;
List<vector> OmegaDotn_;
//- Force driving the motion in inertial coord. sys.
dimensionedVector force_;
@ -205,6 +226,9 @@ class sixDOFqODE
// system
void constrainTranslation(vector& vec) const;
//- Update Aitkens relaxation parameters
void aitkensRelaxation(const scalar min, const scalar max);
public:
@ -344,13 +368,19 @@ public:
//- Return access to moment in relative coord. sys.
inline dimensionedVector& momentRelative();
//- Return total force in inertial coord. sys.
inline dimensionedVector forceTotal() const;
//- Return total moment in inertial coord. sys.
inline dimensionedVector momentTotal() const;
//- Relax the force(acceleration) using Aitkens or fixed relaxation
void relaxAcceleration
(
const scalar minRelFactor,
const scalar maxRelFactor
);
// Rotations

4
src/POD/scalarTransportPOD/scalarTransportPOD.C
View file

@ -234,13 +234,11 @@ void Foam::scalarTransportPOD::calcDerivativeCoeffs() const
{
const volScalarField& snapI = b.orthoField(i);
volVectorField gradSnapI = fvc::grad(snapI);
for (label j = 0; j < b.baseSize(); j++)
{
const volScalarField& snapJ = b.orthoField(j);
volVectorField gradSnapJ = fvc::grad(snapJ);
const tmp<volVectorField> gradSnapJ = fvc::grad(snapJ);
derivative[i][j] =
DT.value()*POD::projection(fvc::div(gradSnapJ), snapI)

24
src/coupledMatrix/coupledFvMatrices/coupledFvMatrix/coupledFvMatrix.C
View file

@ -58,7 +58,7 @@ Foam::word Foam::coupledFvMatrix<Type>::coupledPsiName() const
if (rowI < matrices.size() - 1)
{
cpn += "+";
cpn += "_";
}
}
@ -198,7 +198,14 @@ Foam::coupledFvMatrix<Type>::solve(const dictionary& solverControls)
{
fvMatrix<Type>& curMatrix =
static_cast<fvMatrix<Type>& >(matrices[rowI]);
curMatrix.psi().internalField().replace(cmpt, psiCmpt[rowI]);
GeometricField<Type, fvPatchField, volMesh>& psiRef =
const_cast<GeometricField<Type, fvPatchField, volMesh>&>
(
curMatrix.psi()
);
psiRef.internalField().replace(cmpt, psiCmpt[rowI]);
curMatrix.diag() = saveDiag[rowI];
}
}
@ -206,9 +213,16 @@ Foam::coupledFvMatrix<Type>::solve(const dictionary& solverControls)
// Correct boundary conditions
forAll (matrices, rowI)
{
fvMatrix<Type>& curMatrix =
static_cast<fvMatrix<Type>& >(matrices[rowI]);
curMatrix.psi().correctBoundaryConditions();
const fvMatrix<Type>& curMatrix =
static_cast<const fvMatrix<Type>&>(matrices[rowI]);
GeometricField<Type, fvPatchField, volMesh>& psiRef =
const_cast<GeometricField<Type, fvPatchField, volMesh>&>
(
curMatrix.psi()
);
psiRef.correctBoundaryConditions();
}
return solverPerfVec;

34
src/coupledMatrix/coupledFvMatrices/coupledFvScalarMatrix/coupledFvScalarMatrix.C
View file

@ -69,10 +69,13 @@ coupledSolverPerformance coupledFvMatrix<scalar>::solve
static_cast<fvScalarMatrix&>(matrices[rowI]);
saveDiag.set(rowI, new scalarField(curMatrix.diag()));
// HR 17/Feb/2013
// Need to be able to compare references to support hacks such as in jumpCyclic
// psi.set(rowI, new scalarField(curMatrix.psi()));
psi.set(rowI, &curMatrix.psi());
// Need to be able to compare references to support hacks
// such as in jumpCyclic. HR 17/Feb/2013
scalarField& psiRef =
const_cast<scalarField&>(curMatrix.psi().internalField());
psi.set(rowI, &psiRef);
source.set(rowI, new scalarField(curMatrix.source()));
curMatrix.addBoundarySource(source[rowI], 0);
@ -113,27 +116,20 @@ coupledSolverPerformance coupledFvMatrix<scalar>::solve
solverPerf.print();
// HR 17/Feb/2013
// Not needed since reference is used
// Update solution
//forAll (matrices, rowI)
//{
// fvScalarMatrix& curMatrix =
// static_cast<fvScalarMatrix&>(matrices[rowI]);
//
// curMatrix.psi().internalField() = psi[rowI];
//}
// Correct boundary conditions
forAll (matrices, rowI)
{
fvScalarMatrix& curMatrix =
static_cast<fvScalarMatrix&>(matrices[rowI]);
const fvScalarMatrix& curMatrix =
static_cast<const fvScalarMatrix&>(matrices[rowI]);
curMatrix.psi().correctBoundaryConditions();
volScalarField& psiRef =
const_cast<volScalarField&>(curMatrix.psi());
psiRef.correctBoundaryConditions();
}
//HR 17.2.2013: Clear references to internal field without deleting the objects
// Clear references to internal field without deleting the objects
// HR 17.2.2013
forAll (matrices, rowI)
{
psi.set(rowI, NULL).ptr();

232
src/dbns/limiter/BarthJespersenLimiter.H
View file

@ -1,232 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 3.2
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
BarthJespersenLimiter
Description
Barth-Jespersen limiter
Author
Aleksandar Jemcov, Hrvoje Jasak
\*---------------------------------------------------------------------------*/
#ifndef BarthJespersenLimiter_H
#define BarthJespersenLimiter_H
#include "vector.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class BarthJespersenLimiter Declaration
\*---------------------------------------------------------------------------*/
class BarthJespersenLimiter
{
public:
// Constructors
//- Construct null
BarthJespersenLimiter()
{}
// Desctructor - default
// Member functions
//- Return limiter value
inline scalar limiter
(
const scalar& cellVolume,
const scalar& deltaOneMax,
const scalar& deltaOneMin,
const scalar& deltaTwo
)
{
if (mag(deltaTwo) < SMALL)
{
return 1;
}
else
{
scalar deltaTwoStab = stabilise(deltaTwo, VSMALL);
return min
(
max
(
max(deltaOneMax/deltaTwoStab, 0),
max(deltaOneMin/deltaTwoStab, 0)
),
1
);
}
}
//- Return scalar limiter
inline scalar limiter
(
const scalar& cellVolume,
const scalar& faceFlux,
const scalar& deltaOneMax,
const scalar& deltaOneMin,
const vector& gradPhiP,
const vector& gradPhiN,
const vector& d
)
{
return limiter
(
cellVolume,
deltaOneMax,
deltaOneMin,
(d & gradPhiP)
);
}
//- Return vector limiter
inline vector limiter
(
const scalar& cellVolume,
const scalar& faceFlux,
const vector& deltaOneMax,
const vector& deltaOneMin,
const tensor& gradPhiP,
const tensor& gradPhiN,
const vector& d
)
{
vector deltaTwo = d & gradPhiP;
return vector
(
limiter
(
cellVolume,
deltaOneMax[vector::X],
deltaOneMin[vector::X],
deltaTwo[vector::X]
),
limiter
(
cellVolume,
deltaOneMax[vector::Y],
deltaOneMin[vector::Y],
deltaTwo[vector::Y]
),
limiter
(
cellVolume,
deltaOneMax[vector::Z],
deltaOneMin[vector::Z],
deltaTwo[vector::Z]
)
);
}
//- Return limiter for a symmTensor
inline symmTensor limiter
(
const scalar& cellVolume,
const scalar& faceFlux,
const symmTensor& deltaOneMax,
const symmTensor& deltaOneMin,
const vector& gradPhiPXX,
const vector& gradPhiNXX,
const vector& gradPhiPXY,
const vector& gradPhiNXY,
const vector& gradPhiPXZ,
const vector& gradPhiNXZ,
const vector& gradPhiPYY,
const vector& gradPhiNYY,
const vector& gradPhiPYZ,
const vector& gradPhiNYZ,
const vector& gradPhiPZZ,
const vector& gradPhiNZZ,
const vector& d
)
{
return symmTensor
(
limiter
(
cellVolume,
deltaOneMax[0],
deltaOneMin[0],
d & gradPhiPXX
),
limiter
(
cellVolume,
deltaOneMax[1],
deltaOneMin[1],
d & gradPhiPXY
),
limiter
(
cellVolume,
deltaOneMax[2],
deltaOneMin[2],
d & gradPhiPXZ
),
limiter
(
cellVolume,
deltaOneMax[3],
deltaOneMin[3],
d & gradPhiPYY
),
limiter
(
cellVolume,
deltaOneMax[4],
deltaOneMin[4],
d & gradPhiPYZ
),
limiter
(
cellVolume,
deltaOneMax[5],
deltaOneMin[5],
d & gradPhiPZZ
)
);
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
#endif
// ************************************************************************* //

44
src/dbns/limiter/MDLimiter.H
View file

@ -191,10 +191,6 @@ public:
const GradFieldType& gradPhiIn = gradPhi.internalField();
// Compute limiter values, internal faces
// Dummy flux
scalar upwindFlux = 1;
forAll (owner, faceI)
{
const label& own = owner[faceI];
@ -203,34 +199,27 @@ public:
vector deltaRLeft = faceCentre[faceI] - cellCentre[own];
vector deltaRRight = faceCentre[faceI] - cellCentre[nei];
Type phiOwnerLimiter = limitFunction.limiter
// Find minimal limiter value in each cell
// Owner side
limitFunction.limiter
(
phiLimiterIn[own],
cellVolume[own],
upwindFlux,
phiMaxIn[own] - phi_[own],
phiMinIn[own] - phi_[own],
gradPhiIn[own],
gradPhiIn[nei],
deltaRLeft
(deltaRLeft & gradPhiIn[own])
);
Type phiNeighbourLimiter = limitFunction.limiter
// Neighbour side
limitFunction.limiter
(
phiLimiterIn[nei],
cellVolume[nei],
upwindFlux,
phiMaxIn[nei] - phi_[nei],
phiMinIn[nei] - phi_[nei],
gradPhiIn[nei],
gradPhiIn[own],
deltaRRight
(deltaRRight & gradPhiIn[nei])
);
// find minimal limiter value in each cell
phiLimiterIn[own] =
min(phiLimiterIn[own], phiOwnerLimiter);
phiLimiterIn[nei] =
min(phiLimiterIn[nei], phiNeighbourLimiter);
}
// Coupled boundaries
@ -255,24 +244,19 @@ public:
const GradFieldType gradPhiRight =
gradPhi.boundaryField()[patchI].patchNeighbourField();
// Find minimal limiter value in each cell
forAll (fc, faceI)
{
const label& curFC = fc[faceI];
Type phiFCLimiter = limitFunction.limiter
limitFunction.limiter
(
phiLimiterIn[curFC],
cellVolume[curFC],
upwindFlux,
phiMaxIn[curFC] - phi_[curFC],
phiMinIn[curFC] - phi_[curFC],
gradPhiLeft[faceI],
gradPhiRight[faceI],
deltaR[faceI]
(deltaR[faceI] & gradPhiLeft[faceI])
);
// Find minimal limiter value in each cell
phiLimiterIn[curFC] =
min(phiLimiterIn[curFC], phiFCLimiter);
}
}
}

273
src/dbns/limiter/VenkatakrishnanLimiter.H
View file

@ -1,273 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 3.2
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
VenkatakrishnanLimiter
Description
Venkatakrishnan differentiable limiter
Author
Aleksandar Jemcov
Updated by Hrvoje Jasak
\*---------------------------------------------------------------------------*/
#ifndef VenkatakrishnanLimiter_H
#define VenkatakrishnanLimiter_H
#include "vector.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class VenkatakrishnanLimiter Declaration
\*---------------------------------------------------------------------------*/
class VenkatakrishnanLimiter
{
// Private data
//- Limiter value
scalar k_;
public:
// Constructor
//- Construct null
VenkatakrishnanLimiter()
:
k_(5)
{}
//- Return limiter
inline scalar limiter
(
const scalar& cellVolume,
const scalar& deltaOneMax,
const scalar& deltaOneMin,
const scalar& deltaTwo
)
{
scalar epsilonSquare = pow3(k_)*cellVolume;
if (deltaTwo > 0)
{
return max
(
min
(
(
(sqr(deltaOneMax) + epsilonSquare)*deltaTwo
+ 2*sqr(deltaTwo)*deltaOneMax
)/
stabilise
(
deltaTwo*
(
sqr(deltaOneMax)
+ 2.0*sqr(deltaTwo)
+ deltaOneMax*deltaTwo
+ epsilonSquare
),
SMALL
),
1
),
0
);
}
else if (deltaTwo < 0)
{
return max
(
min
(
(
(sqr(deltaOneMin) + epsilonSquare)*deltaTwo
+ 2*sqr(deltaTwo)*deltaOneMin
)/
stabilise
(
deltaTwo*
(
sqr(deltaOneMin)
+ 2.0*sqr(deltaTwo) + deltaOneMin*deltaTwo
+ epsilonSquare
),
SMALL
),
1
),
0
);
}
else
{
return 1;
}
}
//- Return limiter for a scalar
inline scalar limiter
(
const scalar& cellVolume,
const scalar& faceFlux,
const scalar& deltaOneMax,
const scalar& deltaOneMin,
const vector& gradPhiP,
const vector& gradPhiN,
const vector& d
)
{
return limiter
(
cellVolume,
deltaOneMax,
deltaOneMin,
d & gradPhiP
);
}
//- Return limiter for a vector
inline vector limiter
(
const scalar& cellVolume,
const scalar& faceFlux,
const vector& deltaOneMax,
const vector& deltaOneMin,
const tensor& gradPhiP,
const tensor& gradPhiN,
const vector& d
)
{
vector deltaTwo = d & gradPhiP;
return vector
(
limiter
(
cellVolume,
deltaOneMax[vector::X],
deltaOneMin[vector::X],
deltaTwo[vector::X]
),
limiter
(
cellVolume,
deltaOneMax[vector::Y],
deltaOneMin[vector::Y],
deltaTwo[vector::Y]
),
limiter
(
cellVolume,
deltaOneMax[vector::Z],
deltaOneMin[vector::Z],
deltaTwo[vector::Z]
)
);
}
//- Return limiter for a symmTensor
inline symmTensor limiter
(
const scalar& cellVolume,
const scalar& faceFlux,
const symmTensor& deltaOneMax,
const symmTensor& deltaOneMin,
const vector& gradPhiPXX,
const vector& gradPhiNXX,
const vector& gradPhiPXY,
const vector& gradPhiNXY,
const vector& gradPhiPXZ,
const vector& gradPhiNXZ,
const vector& gradPhiPYY,
const vector& gradPhiNYY,
const vector& gradPhiPYZ,
const vector& gradPhiNYZ,
const vector& gradPhiPZZ,
const vector& gradPhiNZZ,
const vector& d
)
{
return symmTensor
(
limiter
(
cellVolume,
deltaOneMax[0],
deltaOneMin[0],
d & gradPhiPXX
),
limiter
(
cellVolume,
deltaOneMax[1],
deltaOneMin[1],
d & gradPhiPXY
),
limiter
(
cellVolume,
deltaOneMax[2],
deltaOneMin[2],
d & gradPhiPXZ
),
limiter
(
cellVolume,
deltaOneMax[3],
deltaOneMin[3],
d & gradPhiPYY
),
limiter
(
cellVolume,
deltaOneMax[4],
deltaOneMin[4],
d & gradPhiPYZ
),
limiter
(
cellVolume,
deltaOneMax[5],
deltaOneMin[5],
d & gradPhiPZZ
)
);
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
#endif
// ************************************************************************* //

72
src/dbns/limiter/firstOrderLimiter.H
View file

@ -43,7 +43,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class firstOrderLimiter Declaration
Class firstOrderLimiter Declaration
\*---------------------------------------------------------------------------*/
class firstOrderLimiter
@ -56,71 +56,37 @@ public:
firstOrderLimiter()
{}
//- Return limiter
inline scalar limiter
// Destructor - default
// Member functions
//- Set scalar limiter value
inline void limiter
(
scalar& lim,
const scalar& cellVolume,
const scalar& deltaOneMax,
const scalar& deltaOneMin,
const scalar& deltaTwo
)
{
return 0;
lim = 0;
}
//- Return limiter for a scalar
inline scalar limiter
//- Set Type limiter
template<class Type>
inline void limiter
(
Type& lim,
const scalar& cellVolume,
const scalar& faceFlux,
const scalar& deltaOneMax,
const scalar& deltaOneMin,
const vector& gradPhiP,
const vector& gradPhiN,
const vector& d
const Type& deltaOneMax,
const Type& deltaOneMin,
const Type& extrapolate
)
{
return 0;
}
//- Return limiter for a vector
inline vector limiter
(
const scalar& cellVolume,
const scalar& faceFlux,
const vector& deltaOneMax,
const vector& deltaOneMin,
const tensor& gradPhiP,
const tensor& gradPhiN,
const vector& d
)
{
return vector::zero;
}
//- Return limiter for a symmTensor
inline symmTensor limiter
(
const scalar& cellVolume,
const scalar& faceFlux,
const symmTensor& deltaOneMax,
const symmTensor& deltaOneMin,
const vector& gradPhiPXX,
const vector& gradPhiNXX,
const vector& gradPhiPXY,
const vector& gradPhiNXY,
const vector& gradPhiPXZ,
const vector& gradPhiNXZ,
const vector& gradPhiPYY,
const vector& gradPhiNYY,
const vector& gradPhiPYZ,
const vector& gradPhiNYZ,
const vector& gradPhiPZZ,
const vector& gradPhiNZZ,
const vector& d
)
{
return symmTensor::zero;
lim = pTraits<Type>::zero;
}
};

46
src/dbns/numericFlux/fineNumericFlux.C
View file

@ -44,10 +44,7 @@ Foam::fineNumericFlux<Flux, Limiter>::fineNumericFlux
thermo_(thermo),
rhoFlux_(fieldLevel.rhoFlux()),
rhoUFlux_(fieldLevel.rhoUFlux()),
rhoEFlux_(fieldLevel.rhoEFlux()),
gradP_(fvc::grad(p_)),
gradU_(fvc::grad(U_)),
gradT_(fvc::grad(T_))
rhoEFlux_(fieldLevel.rhoEFlux())
{}
@ -67,31 +64,36 @@ void Foam::fineNumericFlux<Flux, Limiter>::computeInteriorFlux()
const scalarField& Cv = fieldLevel_.CvVar();
const scalarField& R = fieldLevel_.R();
gradP_ = fvc::grad(p_);
gradP_.correctBoundaryConditions();
// Get gradients
// Coupled patch update on gradients moved into gradScheme.C
// HJ, 22/Apr/2016;
gradU_ = fvc::grad(U_);
gradU_.correctBoundaryConditions();
// Changed return type for gradient cacheing. HJ, 22/Apr/2016
const tmp<volVectorField> tgradP = fvc::grad(p_);
const volVectorField& gradP = tgradP();
gradT_ = fvc::grad(T_);
gradT_.correctBoundaryConditions();
const tmp<volTensorField> tgradU = fvc::grad(U_);
const volTensorField& gradU = tgradU();
const tmp<volVectorField> tgradT = fvc::grad(T_);
const volVectorField& gradT = tgradT();
MDLimiter<scalar, Limiter> scalarPLimiter
(
this->p_,
this->gradP_
gradP
);
MDLimiter<vector, Limiter> vectorULimiter
(
this->U_,
this->gradU_
gradU
);
MDLimiter<scalar, Limiter> scalarTLimiter
(
this->T_,
this->gradT_
gradT
);
// Get limiters
@ -114,12 +116,12 @@ void Foam::fineNumericFlux<Flux, Limiter>::computeInteriorFlux()
rhoFlux_[faceI],
rhoUFlux_[faceI],
rhoEFlux_[faceI],
p_[own] + pLimiter[own]*(deltaRLeft & gradP_[own]),
p_[nei] + pLimiter[nei]*(deltaRRight & gradP_[nei]),
U_[own] + cmptMultiply(ULimiter[own], (deltaRLeft & gradU_[own])),
U_[nei] + cmptMultiply(ULimiter[nei], (deltaRRight & gradU_[nei])),
T_[own] + TLimiter[own]*(deltaRLeft & gradT_[own]),
T_[nei] + TLimiter[nei]*(deltaRRight & gradT_[nei]),
p_[own] + pLimiter[own]*(deltaRLeft & gradP[own]),
p_[nei] + pLimiter[nei]*(deltaRRight & gradP[nei]),
U_[own] + cmptMultiply(ULimiter[own], (deltaRLeft & gradU[own])),
U_[nei] + cmptMultiply(ULimiter[nei], (deltaRRight & gradU[nei])),
T_[own] + TLimiter[own]*(deltaRLeft & gradT[own]),
T_[nei] + TLimiter[nei]*(deltaRRight & gradT[nei]),
R[own],
R[nei],
Cv[own],
@ -147,9 +149,9 @@ void Foam::fineNumericFlux<Flux, Limiter>::computeInteriorFlux()
const scalarField& pR = fieldLevel_.patchR(patchi);
// Gradients
const fvPatchVectorField& pGradP = gradP_.boundaryField()[patchi];
const fvPatchTensorField& pGradU = gradU_.boundaryField()[patchi];
const fvPatchVectorField& pGradT = gradT_.boundaryField()[patchi];
const fvPatchVectorField& pGradP = gradP.boundaryField()[patchi];
const fvPatchTensorField& pGradU = gradU.boundaryField()[patchi];
const fvPatchVectorField& pGradT = gradT.boundaryField()[patchi];
// Limiters
const fvPatchScalarField& pPatchLim = pLimiter.boundaryField()[patchi];

32
src/dbns/numericFlux/fineNumericFlux.H
View file

@ -98,16 +98,7 @@ class fineNumericFlux
surfaceScalarField& rhoEFlux_;
// Gradients
//- Static pressure gradient
volVectorField gradP_;
//- Velocity gradient
volTensorField gradU_;
//- Static temperature gradient
volVectorField gradT_;
// Gradients: no longer stored. HJ, 22/Apr/2016
// Private Member Functions
@ -226,27 +217,6 @@ public:
}
// Return Gradients
//- Return pressure gradient
const volVectorField& gradP() const
{
return gradP_;
}
//- Return pressure gradient
const volTensorField& gradU() const
{
return gradU_;
}
//- Return Temperature gradient
const volVectorField& gradT() const
{
return gradT_;
}
// Update fluxes based on current state
//- Compute flux

46
src/dbns/numericFlux/numericFlux.C
View file

@ -78,10 +78,7 @@ Foam::numericFlux<Flux, Limiter>::numericFlux
IOobject::NO_WRITE
),
rhoFlux_*linearInterpolate(thermo.Cv()*T_ + 0.5*magSqr(U_))
),
gradP_(fvc::grad(p_)),
gradU_(fvc::grad(U_)),
gradT_(fvc::grad(T_))
)
{}
@ -105,31 +102,36 @@ void Foam::numericFlux<Flux, Limiter>::computeFlux()
const volScalarField Cv = thermo_.Cv();
const volScalarField R = thermo_.Cp() - Cv;
gradP_ = fvc::grad(p_);
gradP_.correctBoundaryConditions();
// Get gradients
// Coupled patch update on gradients moved into gradScheme.C
// HJ, 22/Apr/2016;
gradU_ = fvc::grad(U_);
gradU_.correctBoundaryConditions();
// Changed return type for gradient cacheing. HJ, 22/Apr/2016
const tmp<volVectorField> tgradP = fvc::grad(p_);
const volVectorField& gradP = tgradP();
gradT_ = fvc::grad(T_);
gradT_.correctBoundaryConditions();
const tmp<volTensorField> tgradU = fvc::grad(U_);
const volTensorField& gradU = tgradU();
const tmp<volVectorField> tgradT = fvc::grad(T_);
const volVectorField& gradT = tgradT();
MDLimiter<scalar, Limiter> scalarPLimiter
(
this->p_,
this->gradP_
gradP
);
MDLimiter<vector, Limiter> vectorULimiter
(
this->U_,
this->gradU_
gradU
);
MDLimiter<scalar, Limiter> scalarTLimiter
(
this->T_,
this->gradT_
gradT
);
// Get limiters
@ -152,12 +154,12 @@ void Foam::numericFlux<Flux, Limiter>::computeFlux()
rhoFlux_[faceI],
rhoUFlux_[faceI],
rhoEFlux_[faceI],
p_[own] + pLimiter[own]*(deltaRLeft & gradP_[own]),
p_[nei] + pLimiter[nei]*(deltaRRight & gradP_[nei]),
U_[own] + cmptMultiply(ULimiter[own], (deltaRLeft & gradU_[own])),
U_[nei] + cmptMultiply(ULimiter[nei], (deltaRRight & gradU_[nei])),
T_[own] + TLimiter[own]*(deltaRLeft & gradT_[own]),
T_[nei] + TLimiter[nei]*(deltaRRight & gradT_[nei]),
p_[own] + pLimiter[own]*(deltaRLeft & gradP[own]),
p_[nei] + pLimiter[nei]*(deltaRRight & gradP[nei]),
U_[own] + cmptMultiply(ULimiter[own], (deltaRLeft & gradU[own])),
U_[nei] + cmptMultiply(ULimiter[nei], (deltaRRight & gradU[nei])),
T_[own] + TLimiter[own]*(deltaRLeft & gradT[own]),
T_[nei] + TLimiter[nei]*(deltaRRight & gradT[nei]),
R[own],
R[nei],
Cv[own],
@ -186,9 +188,9 @@ void Foam::numericFlux<Flux, Limiter>::computeFlux()
const scalarField& pR = R.boundaryField()[patchi];
// Gradients
const fvPatchVectorField& pGradP = gradP_.boundaryField()[patchi];
const fvPatchTensorField& pGradU = gradU_.boundaryField()[patchi];
const fvPatchVectorField& pGradT = gradT_.boundaryField()[patchi];
const fvPatchVectorField& pGradP = gradP.boundaryField()[patchi];
const fvPatchTensorField& pGradU = gradU.boundaryField()[patchi];
const fvPatchVectorField& pGradT = gradT.boundaryField()[patchi];
// Limiters
const fvPatchScalarField& pPatchLim = pLimiter.boundaryField()[patchi];

32
src/dbns/numericFlux/numericFlux.H
View file

@ -90,16 +90,7 @@ class numericFlux
surfaceScalarField rhoEFlux_;
// Gradients
//- Static pressure gradient
volVectorField gradP_;
//- Velocity gradient
volTensorField gradU_;
//- Static temperature gradient
volVectorField gradT_;
// Gradients: no longer stored. HJ, 22/Apr/2016
// Private Member Functions
@ -218,27 +209,6 @@ public:
}
// Return Gradients
//- Return pressure gradient
const volVectorField& gradP() const
{
return gradP_;
}
//- Return pressure gradient
const volTensorField& gradU() const
{
return gradU_;
}
//- Return Temperature gradient
const volVectorField& gradT() const
{
return gradT_;
}
// Update fluxes based on current state
//- Compute flux

6
src/dbns/timeStepping/EulerLocalDdtScheme/EulerLocalDdtScheme.C
View file

@ -296,7 +296,7 @@ template<class Type>
tmp<fvMatrix<Type> >
EulerLocalDdtScheme<Type>::fvmDdt
(
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
const objectRegistry& registry = this->mesh();
@ -342,7 +342,7 @@ tmp<fvMatrix<Type> >
EulerLocalDdtScheme<Type>::fvmDdt
(
const dimensionedScalar& rho,
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
const objectRegistry& registry = this->mesh();
@ -386,7 +386,7 @@ tmp<fvMatrix<Type> >
EulerLocalDdtScheme<Type>::fvmDdt
(
const volScalarField& rho,
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
const objectRegistry& registry = this->mesh();

6
src/dbns/timeStepping/EulerLocalDdtScheme/EulerLocalDdtScheme.H
View file

@ -133,19 +133,19 @@ public:
tmp<fvMatrix<Type> > fvmDdt
(
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
tmp<fvMatrix<Type> > fvmDdt
(
const dimensionedScalar&,
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
tmp<fvMatrix<Type> > fvmDdt
(
const volScalarField&,
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;

6
src/dbns/timeStepping/backwardDualDdtScheme/backwardDualDdtScheme.C
View file

@ -399,7 +399,7 @@ template<class Type>
tmp<fvMatrix<Type> >
backwardDualDdtScheme<Type>::fvmDdt
(
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
const objectRegistry& registry = this->mesh();
@ -479,7 +479,7 @@ tmp<fvMatrix<Type> >
backwardDualDdtScheme<Type>::fvmDdt
(
const dimensionedScalar& rho,
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
const objectRegistry& registry = this->mesh();
@ -557,7 +557,7 @@ tmp<fvMatrix<Type> >
backwardDualDdtScheme<Type>::fvmDdt
(
const volScalarField& rho,
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
const objectRegistry& registry = this->mesh();

6
src/dbns/timeStepping/backwardDualDdtScheme/backwardDualDdtScheme.H
View file

@ -165,19 +165,19 @@ public:
tmp<fvMatrix<Type> > fvmDdt
(
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
tmp<fvMatrix<Type> > fvmDdt
(
const dimensionedScalar&,
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
tmp<fvMatrix<Type> > fvmDdt
(
const volScalarField&,
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;

3
src/dynamicMesh/dynamicMesh/polyMeshModifiers/repatchCoverage/repatchCoverage.C
View file

@ -292,7 +292,8 @@ bool Foam::repatchCoverage::changeTopology() const
slave,
tensorField::zero, // forwardT
tensorField::zero, // reverseT
vectorField::zero // separation
vectorField::zero, // separation
false // Patches are not global
);
// Check uncovered master and slave faces

2
src/dynamicMesh/dynamicTopoFvMesh/fieldMapping/topoMapperTemplates.C
View file

@ -86,6 +86,8 @@ void topoMapper::storeGradients
// Register field under a name that's unique
word registerName("remapGradient(" + field.name() + ')');
// Note: potential issue with cached gradients. HJ, 22/Apr/2016
// Make a new entry
if (mesh_.schemesDict().subDict("gradSchemes").found(gradName))
{

9
src/dynamicMesh/dynamicTopoFvMesh/lengthScaleEstimator/lengthScaleEstimator.C
View file

@ -32,7 +32,7 @@ Author
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
\*---------------------------------------------------------------------------*/
#include "fvc.H"
#include "volFields.H"
@ -1039,9 +1039,6 @@ void lengthScaleEstimator::calculateLengthScale
// Lookup various field types, and evaluate the gradient
bool invalidObject = true;
// Evaluate using gradient scheme
word gradName("grad(" + field_ + ')');
// Register field under a name that's unique
word registerName("lengthScaleGradient(" + field_ + ')');
@ -1065,7 +1062,7 @@ void lengthScaleEstimator::calculateLengthScale
IOobject::NO_WRITE,
false
),
mag(fvc::grad(field, gradName))
mag(fvc::grad(field))
)
);
@ -1092,7 +1089,7 @@ void lengthScaleEstimator::calculateLengthScale
IOobject::NO_WRITE,
false
),
mag(fvc::grad(field, gradName))
mag(fvc::grad(field))
)
);

3
src/dynamicMesh/meshMotion/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
View file

@ -136,7 +136,8 @@ void Foam::displacementSBRStressFvMotionSolver::solve()
surfaceScalarField Df = diffusivityPtr_->operator()();
volTensorField gradCd = fvc::grad(cellDisplacement_);
const tmp<volTensorField> tgradCd = fvc::grad(cellDisplacement_);
const volTensorField& gradCd = tgradCd();
Foam::solve
(

10
src/dynamicMesh/meshMotion/solidBodyMotion/constantVelocity/constantVelocity.C
View file

@ -57,11 +57,11 @@ Foam::solidBodyMotionFunctions::constantVelocity::calcTransformation
const vector translation = transVelocity_*(t - startMotionTime_);
const vector rotation = rotVelocity_*(t - startMotionTime_);
const quaternion R(rotation.x(), rotation.y(), rotation.z());
const septernion TR
(
septernion(origin_ + translation)*R*septernion(-origin_)
);
const quaternion R(rotation.x(), rotation.y(), rotation.z());
const septernion TR
(
septernion(origin_ + translation)*R*septernion(-origin_)
);
return TR;
}

17
src/dynamicMesh/meshMotion/solidBodyMotion/rotatingOscillation/rotatingOscillation.H
View file

@ -47,7 +47,7 @@ namespace solidBodyMotionFunctions
{
/*---------------------------------------------------------------------------*\
Class rotatingOscillation Declaration
Class rotatingOscillation Declaration
\*---------------------------------------------------------------------------*/
class rotatingOscillation
@ -97,19 +97,6 @@ public:
const Time& runTime
);
/* //- Construct and return a clone
virtual autoPtr<solidBodyMotionFunction> clone() const
{
return autoPtr<solidBodyMotionFunction>
(
new rotatingOscillation
(
SBMFCoeffs_,
time_
)
);
}
*/
//- Destructor
virtual ~rotatingOscillation();
@ -119,7 +106,7 @@ public:
//- Return the solid-body motion transformation septernion
virtual septernion transformation() const;
//- Return the solid-body motion velocity
virtual septernion velocity() const;

Some files were not shown because too many files have changed in this diff Show more