zhyang-dev/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2012-01-03 12:00:33 +00:00
parent d941a69831 e1c876f59a
commit f2f769e299
268 changed files with 8881 additions and 2833 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
View File

@ -109,6 +109,9 @@ int main(int argc, char *argv[])
surfaceScalarField phiv_pos(U_pos & mesh.Sf());
surfaceScalarField phiv_neg(U_neg & mesh.Sf());
fvc::makeRelative(phiv_pos, U);
fvc::makeRelative(phiv_neg, U);
volScalarField c(sqrt(thermo.Cp()/thermo.Cv()*rPsi));
surfaceScalarField cSf_pos
(
@ -161,17 +164,9 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
mesh.movePoints(motionPtr->newPoints());
phiv_pos = U_pos & mesh.Sf();
phiv_neg = U_neg & mesh.Sf();
fvc::makeRelative(phiv_pos, U);
fvc::makeRelative(phiv_neg, U);
phiv_neg -= mesh.phi();
phiv_pos *= a_pos;
phiv_neg *= a_neg;
aphiv_pos = phiv_pos - aSf;
aphiv_neg = phiv_neg + aSf;
surfaceScalarField phi("phi", aphiv_pos*rho_pos + aphiv_neg*rho_neg);
phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
Info<< phi.boundaryField()[0] << endl;
surfaceVectorField phiUp
(
@ -183,6 +178,7 @@ int main(int argc, char *argv[])
(
aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
+ aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
+ mesh.phi()*(a_pos*p_pos + a_neg*p_neg)
+ aSf*p_pos - aSf*p_neg
);

12
applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H Normal file
View File

@ -0,0 +1,12 @@
{
solve
(
fvm::ddt(rho, e)
+ fvm::div(phi, e)
- fvm::laplacian(turbulence->alphaEff(), e)
==
- p*fvc::div(phi/fvc::interpolate(rho) + mesh.phi())
);
thermo.correct();
}

41
applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
View File

@ -76,34 +76,23 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
volScalarField bp(pow(beta, -2.65));
volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{
if (Re[celli] > 1000.0)
{
Cds[celli] = 0.44;
}
}
volScalarField Cds
(
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ pos(Re - 1000)*0.44
);
// Wen and Yu (1966)
tmp<volScalarField> tKWenYu(0.75*Cds*phase2_.rho()*Ur*bp/d);
volScalarField& KWenYu = tKWenYu();
// Ergun
forAll (beta, cellj)
{
if (beta[cellj] <= 0.8)
{
KWenYu[cellj] =
150.0*alpha_[cellj]*phase2_.nu().value()*phase2_.rho().value()
/sqr(beta[cellj]*d[cellj])
+ 1.75*phase2_.rho().value()*Ur[cellj]
/(beta[cellj]*d[cellj]);
}
}
return tKWenYu;
return
(
pos(beta - 0.8)
*(0.75*Cds*phase2_.rho()*Ur*bp/d)
+ neg(beta - 0.8)
*(
150.0*alpha_*phase2_.nu()*phase2_.rho()/(sqr(beta*d))
+ 1.75*phase2_.rho()*Ur/(beta*d)
)
);
}

14
applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
View File

@ -75,15 +75,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowSchillerNaumann::K
volScalarField bp(pow(beta, -2.65));
volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{
if (Re[celli] > 1000.0)
{
Cds[celli] = 0.44;
}
}
volScalarField Cds
(
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ pos(Re - 1000)*0.44
);
return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
}

14
applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
View File

@ -72,15 +72,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SchillerNaumann::K
) const
{
volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{
if (Re[celli] > 1000.0)
{
Cds[celli] = 0.44;
}
}
volScalarField Cds
(
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ pos(Re - 1000)*0.44
);
return 0.75*Cds*phase2_.rho()*Ur/phase1_.d();
}

14
applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
View File

@ -73,15 +73,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SyamlalOBrien::K
{
volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6)));
volScalarField A(pow(beta, 4.14));
volScalarField B(0.8*pow(beta, 1.28));
forAll (beta, celli)
{
if (beta[celli] > 0.85)
{
B[celli] = pow(beta[celli], 2.65);
}
}
volScalarField B
(
neg(beta - 0.85)*(0.8*pow(beta, 1.28))
+ pos(beta - 0.85)*(pow(beta, 2.65))
);
volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));

14
applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
View File

@ -75,15 +75,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::WenYu::K
volScalarField bp(pow(beta, -2.65));
volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{
if (Re[celli] > 1000.0)
{
Cds[celli] = 0.44;
}
}
volScalarField Cds
(
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ pos(Re - 1000)*0.44
);
return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
}

2
applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
View File

@ -149,6 +149,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
}
}
muff.correctBoundaryConditions();
return tmuf;
}

31
applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H
View File

@ -21,16 +21,27 @@ forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+ fvm::div(phase.phiAlpha(), U)
- fvm::Sp(fvc::ddt(alpha) + fvc::div(phase.phiAlpha()), U)
)
- fvm::laplacian(alpha*nuEff, U)
- fvc::div
(
alpha*(nuEff*dev(T(fvc::grad(U))) /*- ((2.0/3.0)*I)*k*/),
"div(Rc)"
)
==
- fvm::Sp(fluid.dragCoeff(phase, dragCoeffs())/phase.rho(), U)
//- (alpha*phase.rho())*fluid.lift(phase)
+ (alpha/phase.rho())*fluid.Svm(phase)
- fvm::laplacian(alpha*nuEff, U)
- fvc::div
(
alpha*(nuEff*dev(T(fvc::grad(U))) /*- ((2.0/3.0)*I)*k*/),
"div(Rc)"
)
==
- fvm::Sp(fluid.dragCoeff(phase, dragCoeffs())/phase.rho(), U)
//- (alpha*phase.rho())*fluid.lift(phase)
+ (alpha/phase.rho())*fluid.Svm(phase)
- fvm::Sp
(
slamDampCoeff
*max
(
mag(U.dimensionedInternalField()) - maxSlamVelocity,
dimensionedScalar("U0", dimVelocity, 0)
)
/pow(mesh.V(), 1.0/3.0),
U
)
)
);
mrfZones.addCoriolis(alpha, UEqns[phasei]);

12
applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
View File

@ -53,6 +53,18 @@
phi += fvc::interpolate(alpha)*phase.phi();
}
scalar slamDampCoeff
(
fluid.lookupOrDefault<scalar>("slamDampCoeff", 1)
);
dimensionedScalar maxSlamVelocity
(
"maxSlamVelocity",
dimVelocity,
fluid.lookupOrDefault<scalar>("maxSlamVelocity", GREAT)
);
// dimensionedScalar pMin
// (
// "pMin",

80
applications/solvers/multiphase/multiphaseEulerFoam/interfacialCoeffs.H
View File

@ -1,80 +0,0 @@
volScalarField dragCoeff
(
IOobject
(
"dragCoeff",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("dragCoeff", dimensionSet(1, -3, -1, 0, 0), 0)
);
volVectorField liftForce
(
IOobject
(
"liftForce",
runTime.timeName(),
mesh
),
mesh,
dimensionedVector("liftForce", dimensionSet(1, -2, -2, 0, 0), vector::zero)
);
volScalarField heatTransferCoeff
(
IOobject
(
"heatTransferCoeff",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("heatTransferCoeff", dimensionSet(1, -1, -3, -1, 0), 0)
);
{
volVectorField Ur = U1 - U2;
volScalarField magUr = mag(Ur);
if (dispersedPhase == "1")
{
dragCoeff = drag1->K(magUr);
heatTransferCoeff = heatTransfer1->K(magUr);
}
else if (dispersedPhase == "2")
{
dragCoeff = drag2->K(magUr);
heatTransferCoeff = heatTransfer2->K(magUr);
}
else if (dispersedPhase == "both")
{
dragCoeff =
(
alpha2*drag1->K(magUr)
+ alpha1*drag2->K(magUr)
);
heatTransferCoeff =
(
alpha2*heatTransfer1->K(magUr)
+ alpha1*heatTransfer2->K(magUr)
);
}
else
{
FatalErrorIn(args.executable())
<< "dispersedPhase: " << dispersedPhase << " is incorrect"
<< exit(FatalError);
}
volScalarField alphaCoeff
(
(alpha1 + minInterfaceAlpha)*(alpha2 + minInterfaceAlpha)
);
dragCoeff *= alphaCoeff;
heatTransferCoeff *= alphaCoeff;
liftForce = Cl*(alpha1*rho1 + alpha2*rho2)*(Ur ^ fvc::curl(U));
}

41
applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
View File

@ -75,34 +75,23 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
volScalarField bp(pow(beta, -2.65));
volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{
if (Re[celli] > 1000.0)
{
Cds[celli] = 0.44;
}
}
volScalarField Cds
(
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ pos(Re - 1000)*0.44
);
// Wen and Yu (1966)
tmp<volScalarField> tKWenYu = 0.75*Cds*phase2_.rho()*Ur*bp/d;
volScalarField& KWenYu = tKWenYu();
// Ergun
forAll (beta, cellj)
{
if (beta[cellj] <= 0.8)
{
KWenYu[cellj] =
150.0*phase1_[cellj]*phase2_.nu().value()*phase2_.rho().value()
/sqr(beta[cellj]*d[cellj])
+ 1.75*phase2_.rho().value()*Ur[cellj]
/(beta[cellj]*d[cellj]);
}
}
return tKWenYu;
return
(
pos(beta - 0.8)
*(0.75*Cds*phase2_.rho()*Ur*bp/d)
+ neg(beta - 0.8)
*(
150.0*phase1_*phase2_.nu()*phase2_.rho()/(sqr(beta*d))
+ 1.75*phase2_.rho()*Ur/(beta*d)
)
);
}

14
applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
View File

@ -74,15 +74,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowSchillerNaumann::K
volScalarField bp(pow(beta, -2.65));
volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{
if (Re[celli] > 1000.0)
{
Cds[celli] = 0.44;
}
}
volScalarField Cds
(
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ pos(Re - 1000)*0.44
);
return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
}

14
applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
View File

@ -71,15 +71,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SchillerNaumann::K
) const
{
volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{
if (Re[celli] > 1000.0)
{
Cds[celli] = 0.44;
}
}
volScalarField Cds
(
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ pos(Re - 1000)*0.44
);
return 0.75*Cds*phase2_.rho()*Ur/phase1_.d();
}

14
applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
View File

@ -72,15 +72,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SyamlalOBrien::K
{
volScalarField beta(max(phase2_, scalar(1.0e-6)));
volScalarField A(pow(beta, 4.14));
volScalarField B(0.8*pow(beta, 1.28));
forAll (beta, celli)
{
if (beta[celli] > 0.85)
{
B[celli] = pow(beta[celli], 2.65);
}
}
volScalarField B
(
neg(beta - 0.85)*(0.8*pow(beta, 1.28))
+ pos(beta - 0.85)*(pow(beta, 2.65))
);
volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));

14
applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
View File

@ -74,15 +74,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::WenYu::K
volScalarField bp(pow(beta, -2.65));
volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{
if (Re[celli] > 1000.0)
{
Cds[celli] = 0.44;
}
}
volScalarField Cds
(
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ pos(Re - 1000)*0.44
);
return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
}

2
applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
View File

@ -149,6 +149,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
}
}
muff.correctBoundaryConditions();
return tmuf;
}

1
applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
View File

@ -82,7 +82,6 @@ int main(int argc, char *argv[])
rho = fluid.rho();
#include "zonePhaseVolumes.H"
//#include "interfacialCoeffs.H"
//#include "TEqns.H"
#include "UEqns.H"

5
applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
View File

@ -629,8 +629,9 @@ Foam::multiphaseSystem::dragCoeffs() const
(
max
(
fvc::average(dm.phase1())*fvc::average(dm.phase2()),
//dm.phase1()*dm.phase2(),
//fvc::average(dm.phase1()*dm.phase2()),
//fvc::average(dm.phase1())*fvc::average(dm.phase2()),
dm.phase1()*dm.phase2(),
dm.residualPhaseFraction()
)
*dm.K

23
applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
View File

@ -37,8 +37,10 @@
U0s.set(phasei, new volVectorField(phase.U()));
phi0s.set(phasei, new surfaceScalarField(phase.phi()));
mrfZones.absoluteFlux(phi0s[phasei]);
phasei++;
}
surfaceScalarField phi0
@ -67,6 +69,8 @@
phase.U() = rAUs[phasei]*UEqns[phasei].H();
phase.phi().oldTime();
mrfZones.absoluteFlux(phase.phi().oldTime());
mrfZones.absoluteFlux(phase.phi());
phase.phi() =
(
@ -74,7 +78,8 @@
+ fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi())
);
mrfZones.relativeFlux(phase.phi());
surfaceScalarField pphi0("pphi0", phase.phi());
mrfZones.relativeFlux(phase.phi().oldTime());
surfaceScalarField pphi0("pphi0", 1.0*phase.phi());
pphi0 += rAlphaAUfs[phasei]*(g & mesh.Sf());
multiphaseSystem::dragModelTable::const_iterator dmIter =
@ -196,7 +201,7 @@
p.relax();
mSfGradp = pEqnIncomp.flux()/Dp;
U = dimensionedVector("U", dimVelocity, vector::zero);
U = dimensionedVector("U", dimVelocity, vector::zero);
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
@ -254,9 +259,17 @@
phasej++;
}
phase.U() +=
(1.0/phase.rho())
*rAUs[phasei]*(*dcIter())*U0s[phasej];
// phase.U() +=
// (1.0/phase.rho())*rAUs[phasei]*(*dcIter())
// *U0s[phasej];
phase.U() += fvc::reconstruct
(
fvc::interpolate
(
(1.0/phase.rho())*rAUs[phasei]*(*dcIter())
)*phi0s[phasej]
);
}
}

8
applications/solvers/multiphase/multiphaseEulerFoam/phaseModel/phaseModel/phaseModel.C
View File

@ -212,10 +212,10 @@ bool Foam::phaseModel::read(const dictionary& phaseDict)
//if (nuModel_->read(phaseDict_))
{
phaseDict_.lookup("nu") >> nu_;
phaseDict_.lookup("kappa") >> kappa_;
phaseDict_.lookup("Cp") >> Cp_;
phaseDict_.lookup("rho") >> rho_;
phaseDict_.lookup("nu") >> nu_.value();
phaseDict_.lookup("kappa") >> kappa_.value();
phaseDict_.lookup("Cp") >> Cp_.value();
phaseDict_.lookup("rho") >> rho_.value();
return true;
}

158
applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C
View File

@ -32,6 +32,7 @@ Description
#include "volFields.H"
#include "PatchEdgeFaceWave.H"
#include "patchEdgeFaceInfo.H"
#include "patchDist.H"
using namespace Foam;
@ -49,80 +50,115 @@ int main(int argc, char *argv[])
// Get name of patch
const word patchName = args[1];
const polyPatch& patch = patches[patchName];
// Data on all edges and faces
List<patchEdgeFaceInfo> allEdgeInfo(patch.nEdges());
List<patchEdgeFaceInfo> allFaceInfo(patch.size());
// Initial seed
DynamicList<label> initialEdges;
DynamicList<patchEdgeFaceInfo> initialEdgesInfo;
// Just set an edge on the master
if (Pstream::master())
// 1. Walk from a single edge
{
label edgeI = 0;
Info<< "Starting walk on edge " << edgeI << endl;
// Data on all edges and faces
List<patchEdgeFaceInfo> allEdgeInfo(patch.nEdges());
List<patchEdgeFaceInfo> allFaceInfo(patch.size());
initialEdges.append(edgeI);
const edge& e = patch.edges()[edgeI];
initialEdgesInfo.append
(
patchEdgeFaceInfo
// Initial seed
DynamicList<label> initialEdges;
DynamicList<patchEdgeFaceInfo> initialEdgesInfo;
// Just set an edge on the master
if (Pstream::master())
{
label edgeI = 0;
Info<< "Starting walk on edge " << edgeI << endl;
initialEdges.append(edgeI);
const edge& e = patch.edges()[edgeI];
initialEdgesInfo.append
(
e.centre(patch.localPoints()),
0.0
)
);
}
patchEdgeFaceInfo
(
e.centre(patch.localPoints()),
0.0
)
);
}
// Walk
PatchEdgeFaceWave
<
primitivePatch,
patchEdgeFaceInfo
> calc
(
mesh,
patch,
initialEdges,
initialEdgesInfo,
allEdgeInfo,
allFaceInfo,
returnReduce(patch.nEdges(), sumOp<label>())
);
// Extract as patchField
volScalarField vsf
(
IOobject
// Walk
PatchEdgeFaceWave
<
primitivePatch,
patchEdgeFaceInfo
> calc
(
"patchDist",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("patchDist", dimLength, 0.0)
);
scalarField pf(vsf.boundaryField()[patch.index()].size());
forAll(pf, faceI)
{
pf[faceI] = Foam::sqrt(allFaceInfo[faceI].distSqr());
patch,
initialEdges,
initialEdgesInfo,
allEdgeInfo,
allFaceInfo,
returnReduce(patch.nEdges(), sumOp<label>())
);
// Extract as patchField
volScalarField vsf
(
IOobject
(
"patchDist",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("patchDist", dimLength, 0.0)
);
scalarField pf(vsf.boundaryField()[patch.index()].size());
forAll(pf, faceI)
{
pf[faceI] = Foam::sqrt(allFaceInfo[faceI].distSqr());
}
vsf.boundaryField()[patch.index()] = pf;
Info<< "Writing patchDist volScalarField to " << runTime.value()
<< endl;
vsf.write();
}
vsf.boundaryField()[patch.index()] = pf;
Info<< "Writing patchDist volScalarField to " << runTime.value()
<< endl;
vsf.write();
// 2. Use a wrapper to walk from all boundary edges on selected patches
{
labelHashSet otherPatchIDs(identity(mesh.boundaryMesh().size()));
otherPatchIDs.erase(patch.index());
Info<< "Walking on patch " << patch.index()
<< " from edges shared with patches " << otherPatchIDs
<< endl;
patchDist pwd(patch, otherPatchIDs);
// Extract as patchField
volScalarField vsf
(
IOobject
(
"otherPatchDist",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("otherPatchDist", dimLength, 0.0)
);
vsf.boundaryField()[patch.index()] = pwd;
Info<< "Writing otherPatchDist volScalarField to " << runTime.value()
<< endl;
vsf.write();
}
Info<< "\nEnd\n" << endl;
return 0;

4
applications/test/tensor/Test-tensor.C
View File

@ -50,10 +50,14 @@ int main()
<< (t1 & t7 & t1.T()) << " " << transform(t1, t7) << endl;
symmTensor st1(1, 2, 3, 4, 5, 6);
symmTensor st2(7, 8, 9, 10, 11, 12);
Info<< "Check symmetric transformation "
<< transform(t1, st1) << endl;
Info<< "Check for dot product of symmetric tensors "
<< (st1 & st2) << endl;
vector v1(1, 2, 3);
Info<< sqr(v1) << endl;

2
applications/utilities/mesh/conversion/gambitToFoam/gambitToFoam.L
View File

@ -433,7 +433,7 @@ mtype {space}"MTYPE:"{space}
}
else
{
curGroupID = readLabel(groupStream);;
curGroupID = readLabel(groupStream);
}
BEGIN(cellStreams);

2
applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.H
View File

@ -237,8 +237,6 @@ public:
Ostream& os,
const cv2DControls& s
);
};

1
applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControlsI.H
View File

@ -23,7 +23,6 @@ License
\*---------------------------------------------------------------------------*/
inline Foam::scalar Foam::cv2DControls::minCellSize() const
{
return minCellSize_;

2
applications/utilities/mesh/generation/cv2DMesh/indexedFaceI.H
View File

@ -106,3 +106,5 @@ inline int CGAL::indexedFace<Gt, Fb>::faceIndex() const
return index_;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

3
applications/utilities/mesh/generation/cv2DMesh/indexedVertexI.H
View File

@ -228,3 +228,6 @@ bool CGAL::outsideTriangle
|| (v1.farPoint() || v1.ppSlave())
|| (v2.farPoint() || v2.ppSlave());
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

2
applications/utilities/mesh/manipulation/renumberMesh/Make/options
View File

@ -11,4 +11,4 @@ EXE_LIBS = \
-lfiniteVolume \
-lgenericPatchFields \
-lrenumberMethods \
-ldecompositionMethods
-ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lmetisDecomp -lscotchDecomp

319
applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
View File

@ -29,7 +29,7 @@ Description
renumbering all fields from all the time directories.
By default uses bandCompression (CuthillMcKee) but will
read system/renumberMeshDict if present and use the method from there.
read system/renumberMeshDict if -dict option is present
\*---------------------------------------------------------------------------*/
@ -45,6 +45,7 @@ Description
#include "renumberMethod.H"
#include "zeroGradientFvPatchFields.H"
#include "CuthillMcKeeRenumber.H"
#include "fvMeshSubset.H"
using namespace Foam;
@ -455,30 +456,86 @@ autoPtr<mapPolyMesh> reorderMesh
}
// Return new to old cell numbering
labelList regionRenumber
(
const renumberMethod& method,
const fvMesh& mesh,
const labelList& cellToRegion
)
{
Info<< "Determining cell order:" << endl;
labelList cellOrder(cellToRegion.size());
label nRegions = max(cellToRegion)+1;
labelListList regionToCells(invertOneToMany(nRegions, cellToRegion));
label cellI = 0;
forAll(regionToCells, regionI)
{
Info<< " region " << regionI << " starts at " << cellI << endl;
// Make sure no parallel comms
bool oldParRun = UPstream::parRun();
UPstream::parRun() = false;
// Per region do a reordering.
fvMeshSubset subsetter(mesh);
subsetter.setLargeCellSubset(cellToRegion, regionI);
const fvMesh& subMesh = subsetter.subMesh();
labelList subReverseCellOrder = method.renumber
(
subMesh,
subMesh.cellCentres()
);
labelList subCellOrder
(
invert
(
subMesh.nCells(),
subReverseCellOrder
)
);
// Restore state
UPstream::parRun() = oldParRun;
const labelList& cellMap = subsetter.cellMap();
forAll(subCellOrder, i)
{
cellOrder[cellI++] = cellMap[subCellOrder[i]];
}
}
Info<< endl;
return cellOrder;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addOption
argList::addNote
(
"blockSize",
"block size",
"order cells into blocks (using decomposition) before ordering"
);
argList::addBoolOption
(
"orderPoints",
"order points into internal and boundary points"
);
argList::addBoolOption
(
"writeMaps",
"write cellMap, faceMap, pointMap in polyMesh/"
"Renumber mesh to minimise bandwidth"
);
# include "addRegionOption.H"
# include "addOverwriteOption.H"
# include "addTimeOptions.H"
argList::addBoolOption
(
"dict",
"renumber according to system/renumberMeshDict"
);
# include "setRootCase.H"
# include "createTime.H"
@ -495,39 +552,7 @@ int main(int argc, char *argv[])
# include "createNamedMesh.H"
const word oldInstance = mesh.pointsInstance();
label blockSize = 0;
args.optionReadIfPresent("blockSize", blockSize, 0);
if (blockSize > 0)
{
Info<< "Ordering cells into regions of size " << blockSize
<< " (using decomposition);"
<< " ordering faces into region-internal and region-external." << nl
<< endl;
if (blockSize < 0 || blockSize >= mesh.nCells())
{
FatalErrorIn(args.executable())
<< "Block size " << blockSize << " should be positive integer"
<< " and less than the number of cells in the mesh."
<< exit(FatalError);
}
}
const bool orderPoints = args.optionFound("orderPoints");
if (orderPoints)
{
Info<< "Ordering points into internal and boundary points." << nl
<< endl;
}
const bool writeMaps = args.optionFound("writeMaps");
if (writeMaps)
{
Info<< "Writing renumber maps (new to old) to polyMesh." << nl
<< endl;
}
const bool readDict = args.optionFound("dict");
const bool overwrite = args.optionFound("overwrite");
@ -538,29 +563,85 @@ int main(int argc, char *argv[])
<< returnReduce(band, maxOp<label>()) << nl << endl;
// Construct renumberMethod
IOobject io
(
"renumberMeshDict",
runTime.system(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
);
bool sortCoupledFaceCells = false;
bool writeMaps = false;
bool orderPoints = false;
label blockSize = 0;
// Construct renumberMethod
autoPtr<IOdictionary> renumberDictPtr;
autoPtr<renumberMethod> renumberPtr;
if (io.headerOk())
if (readDict)
{
Info<< "Detected local " << runTime.system()/io.name() << "." << nl
<< "Using this to select renumberMethod." << nl << endl;
renumberPtr = renumberMethod::New(IOdictionary(io));
Info<< "Renumber according to renumberMeshDict." << nl << endl;
renumberDictPtr.reset
(
new IOdictionary
(
IOobject
(
"renumberMeshDict",
runTime.system(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
)
);
const IOdictionary renumberDict = renumberDictPtr();
renumberPtr = renumberMethod::New(renumberDict);
sortCoupledFaceCells = renumberDict.lookupOrDefault
(
"sortCoupledFaceCells",
false
);
if (sortCoupledFaceCells)
{
Info<< "Sorting cells on coupled boundaries to be last." << nl
<< endl;
}
blockSize = renumberDict.lookupOrDefault("blockSize", 0);
if (blockSize > 0)
{
Info<< "Ordering cells into regions of size " << blockSize
<< " (using decomposition);"
<< " ordering faces into region-internal and region-external."
<< nl << endl;
if (blockSize < 0 || blockSize >= mesh.nCells())
{
FatalErrorIn(args.executable())
<< "Block size " << blockSize
<< " should be positive integer"
<< " and less than the number of cells in the mesh."
<< exit(FatalError);
}
}
orderPoints = renumberDict.lookupOrDefault("orderPoints", false);
if (orderPoints)
{
Info<< "Ordering points into internal and boundary points." << nl
<< endl;
}
renumberDict.lookup("writeMaps") >> writeMaps;
if (writeMaps)
{
Info<< "Writing renumber maps (new to old) to polyMesh." << nl
<< endl;
}
}
else
{
Info<< "No local " << runTime.system()/io.name()
<< " dictionary found. Using default renumberMethod." << nl
<< endl;
Info<< "Using default renumberMethod." << nl << endl;
dictionary renumberDict;
renumberPtr.reset(new CuthillMcKeeRenumber(renumberDict));
}
@ -671,22 +752,15 @@ int main(int argc, char *argv[])
// fields is done correctly!
label nBlocks = mesh.nCells() / blockSize;
Info<< "nBlocks = " << nBlocks << endl;
Info<< "nBlocks = " << nBlocks << endl;
// Read decomposePar dictionary
IOdictionary decomposeDict
(
IOobject
(
"decomposeParDict",
runTime.system(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
// Read decompositionMethod dictionary
dictionary decomposeDict(renumberDictPtr().subDict("blockCoeffs"));
decomposeDict.set("numberOfSubdomains", nBlocks);
bool oldParRun = UPstream::parRun();
UPstream::parRun() = false;
autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
(
decomposeDict
@ -701,6 +775,9 @@ int main(int argc, char *argv[])
)
);
// Restore state
UPstream::parRun() = oldParRun;
// For debugging: write out region
createScalarField
(
@ -714,23 +791,7 @@ int main(int argc, char *argv[])
<< nl << endl;
// Find point per region
pointField regionPoints(nBlocks, vector::zero);
forAll(cellToRegion, cellI)
{
regionPoints[cellToRegion[cellI]] = mesh.cellCentres()[cellI];
}
// Use block based renumbering.
// Detemines old to new cell ordering
labelList reverseCellOrder = renumberPtr().renumber
(
mesh,
cellToRegion,
regionPoints
);
cellOrder = invert(mesh.nCells(), reverseCellOrder);
cellOrder = regionRenumber(renumberPtr(), mesh, cellToRegion);
// Determine new to old face order with new cell numbering
faceOrder = getRegionFaceOrder
@ -751,6 +812,76 @@ int main(int argc, char *argv[])
cellOrder = invert(mesh.nCells(), reverseCellOrder);
if (sortCoupledFaceCells)
{
// Change order so all coupled patch faceCells are at the end.
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
// Collect all boundary cells on coupled patches
label nBndCells = 0;
forAll(pbm, patchI)
{
if (pbm[patchI].coupled())
{
nBndCells += pbm[patchI].size();
}
}
labelList bndCellMap(nBndCells);
labelList bndCells(bndCellMap);
nBndCells = 0;
forAll(pbm, patchI)
{
if (pbm[patchI].coupled())
{
const labelUList& faceCells = pbm[patchI].faceCells();
forAll(faceCells, i)
{
label cellI = faceCells[i];
bndCells[nBndCells] = cellI;
bndCellMap[nBndCells++] = reverseCellOrder[cellI];
}
}
}
bndCells.setSize(nBndCells);
bndCellMap.setSize(nBndCells);
// Sort
labelList order;
sortedOrder(bndCellMap, order);
// Redo newReverseCellOrder
labelList newReverseCellOrder(mesh.nCells(), -1);
label sortedI = mesh.nCells();
forAllReverse(order, i)
{
label origCellI = bndCells[order[i]];
newReverseCellOrder[origCellI] = --sortedI;
}
Info<< "Ordered all " << nBndCells << " cells with a coupled face"
<< " to the end of the cell list, starting at " << sortedI
<< endl;
// Compact
sortedI = 0;
forAll(cellOrder, newCellI)
{
label origCellI = cellOrder[newCellI];
if (newReverseCellOrder[origCellI] == -1)
{
newReverseCellOrder[origCellI] = sortedI++;
}
}
// Update sorted back to original (unsorted) map
cellOrder = invert(mesh.nCells(), newReverseCellOrder);
}
// Determine new to old face order with new cell numbering
faceOrder = getFaceOrder
(

44
applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
View File

@ -15,16 +15,37 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Write maps from renumbered back to original mesh
writeMaps true;
// Optional entry: sort cells on coupled boundaries to last for use with
// e.g. nonBlockingGaussSeidel.
sortCoupledFaceCells false;
// Optional entry: renumber on a block-by-block basis. It uses a
// blockCoeffs dictionary to construct a decompositionMethod to do
// a block subdivision) and then applies the renumberMethod to each
// block in turn. This can be used in large cases to keep the blocks
// fitting in cache with all the the cache misses bunched at the end.
// This number is the approximate size of the blocks - this gets converted
// to a number of blocks that is the input to the decomposition method.
//blockSize 1000;
// Optional entry: sort points into internal and boundary points
//orderPoints false;
method CuthillMcKee;
//method manual;
//method random;
//method spring;
CuthillMcKeeCoeffs
{
// Reverse CuthillMcKee (RCM) or plain
reverse true;
}
//CuthillMcKeeCoeffs
//{
// // Reverse CuthillMcKee (RCM) or plain
// reverse true;
//}
manualCoeffs
@ -48,4 +69,17 @@ springCoeffs
}
blockCoeffs
{
method scotch;
//method hierarchical;
//hierarchicalCoeffs
//{
// n (1 2 1);
// delta 0.001;
// order xyz;
//}
}
// ************************************************************************* //

3
applications/utilities/mesh/manipulation/topoSet/topoSetDict
View File

@ -211,7 +211,8 @@ FoamFile
// source patchToFace;
// sourceInfo
// {
// name ".*Wall"; // Name of patch, regular expressions allowed
// name ".*Wall"; // Name of patch or patch group,
// // (regular expressions allowed)
// }
//
// // All faces of faceZone

10
applications/utilities/parallelProcessing/decomposePar/decomposePar.C
View File

@ -528,7 +528,7 @@ int main(int argc, char *argv[])
<< "Cell number should be between 0 and "
<< mesh.nCells()-1 << nl
<< "On this mesh the particle should be in cell "
<< mesh.findCell(iter().position())
<< mesh.findCell(iter().position())
<< exit(FatalError);
}
@ -789,6 +789,14 @@ int main(int argc, char *argv[])
// Point fields
if
(
pointScalarFields.size()
|| pointVectorFields.size()
|| pointSphericalTensorFields.size()
|| pointSymmTensorFields.size()
|| pointTensorFields.size()
)
{
labelIOList pointProcAddressing
(

7
applications/utilities/parallelProcessing/redistributePar/redistributePar.C
View File

@ -795,9 +795,9 @@ int main(int argc, char *argv[])
);
fvMesh& mesh = meshPtr();
//Pout<< "Read mesh:" << endl;
//printMeshData(Pout, mesh);
//Pout<< endl;
// Print some statistics
Info<< "Before distribution:" << endl;
printMeshData(mesh);
@ -1022,7 +1022,6 @@ int main(int argc, char *argv[])
//map().distributeFaceData(faceCc);
// Print a bit
// Print some statistics
Info<< "After distribution:" << endl;
printMeshData(mesh);

4
applications/utilities/postProcessing/graphics/PV3Readers/Allwmake
View File

@ -1,6 +1,6 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
#set -x
if [ -d "$ParaView_DIR" -a -r "$ParaView_DIR" ]
then
@ -15,6 +15,8 @@ then
wmake libso vtkPV3Readers
PV3blockMeshReader/Allwmake
PV3FoamReader/Allwmake
else
echo "ERROR: ParaView not found in $ParaView_DIR"
fi
# ----------------------------------------------------------------- end-of-file

2
applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_get_number_of_files_in_dataset.H
View File

@ -9,5 +9,5 @@ int USERD_get_number_of_files_in_dataset(void)
// use 1 insted of 0 which gives an un-necessary warning.
Num_dataset_files = 1;
return Num_dataset_files;;
return Num_dataset_files;
}

1
applications/utilities/postProcessing/sampling/sample/sampleDict
View File

@ -20,6 +20,7 @@ FoamFile
// gnuplot
// raw
// vtk
// ensight
// csv
setFormat raw;

7
applications/utilities/postProcessing/wall/wallShearStress/Make/options
View File

@ -1,11 +1,14 @@
EXE_INC = \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lincompressibleRASModels \
-lincompressibleTransportModels \
-lincompressibleRASModels \
-lbasicThermophysicalModels \
-lspecie \
-lcompressibleRASModels \
-lfiniteVolume \
-lgenericPatchFields

22
applications/utilities/postProcessing/wall/wallShearStress/createFields.H
View File

@ -1,22 +0,0 @@
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::RASModel> RASModel
(
incompressible::RASModel::New(U, phi, laminarTransport)
);

153
applications/utilities/postProcessing/wall/wallShearStress/wallShearStress.C
View File

@ -25,23 +25,130 @@ Application
wallShearStress
Description
Calculates and writes the wall shear stress, for the specified times.
Calculates and reports wall shear stress for all patches, for the
specified times when using RAS turbulence models.
Default behaviour assumes operating in incompressible mode.
Use the -compressible option for compressible RAS cases.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
#include "RASModel.H"
#include "incompressible/RAS/RASModel/RASModel.H"
#include "basicPsiThermo.H"
#include "compressible/RAS/RASModel/RASModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void calcIncompressible
(
const fvMesh& mesh,
const Time& runTime,
const volVectorField& U,
volVectorField& wallShearStress
)
{
#include "createPhi.H"
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::RASModel> model
(
incompressible::RASModel::New(U, phi, laminarTransport)
);
const volSymmTensorField Reff(model->devReff());
forAll(wallShearStress.boundaryField(), patchI)
{
wallShearStress.boundaryField()[patchI] =
(
-mesh.Sf().boundaryField()[patchI]
/mesh.magSf().boundaryField()[patchI]
) & Reff.boundaryField()[patchI];
}
}
void calcCompressible
(
const fvMesh& mesh,
const Time& runTime,
const volVectorField& U,
volVectorField& wallShearStress
)
{
IOobject rhoHeader
(
"rho",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
);
if (!rhoHeader.headerOk())
{
Info<< " no rho field" << endl;
return;
}
Info<< "Reading field rho\n" << endl;
volScalarField rho(rhoHeader, mesh);
#include "compressibleCreatePhi.H"
autoPtr<basicPsiThermo> pThermo
(
basicPsiThermo::New(mesh)
);
basicPsiThermo& thermo = pThermo();
autoPtr<compressible::RASModel> model
(
compressible::RASModel::New
(
rho,
U,
phi,
thermo
)
);
const volSymmTensorField Reff(model->devRhoReff());
forAll(wallShearStress.boundaryField(), patchI)
{
wallShearStress.boundaryField()[patchI] =
(
-mesh.Sf().boundaryField()[patchI]
/mesh.magSf().boundaryField()[patchI]
) & Reff.boundaryField()[patchI];
}
}
int main(int argc, char *argv[])
{
timeSelector::addOptions();
#include "addRegionOption.H"
argList::addBoolOption
(
"compressible",
"calculate compressible wall shear stress"
);
#include "setRootCase.H"
#include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
#include "createMesh.H"
#include "createNamedMesh.H"
const bool compressible = args.optionFound("compressible");
forAll(timeDirs, timeI)
{
@ -49,10 +156,6 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << endl;
mesh.readUpdate();
#include "createFields.H"
volSymmTensorField Reff(RASModel->devReff());
volVectorField wallShearStress
(
IOobject
@ -67,19 +170,41 @@ int main(int argc, char *argv[])
dimensionedVector
(
"wallShearStress",
Reff.dimensions(),
sqr(dimLength)/sqr(dimTime),
vector::zero
)
);
forAll(wallShearStress.boundaryField(), patchi)
IOobject UHeader
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
);
if (UHeader.headerOk())
{
wallShearStress.boundaryField()[patchi] =
(
-mesh.Sf().boundaryField()[patchi]
/mesh.magSf().boundaryField()[patchi]
) & Reff.boundaryField()[patchi];
Info<< "Reading field U\n" << endl;
volVectorField U(UHeader, mesh);
if (compressible)
{
calcCompressible(mesh, runTime, U, wallShearStress);
}
else
{
calcIncompressible(mesh, runTime, U, wallShearStress);
}
}
else
{
Info<< " no U field" << endl;
}
Info<< "Writing wall shear stress to field " << wallShearStress.name()
<< nl << endl;
wallShearStress.write();
}

20
applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
View File

@ -35,12 +35,13 @@ License
namespace Foam
{
template<class Type>
template<class Type, class CombineOp>
void MapConsistentVolFields
(
const IOobjectList& objects,
const meshToMesh& meshToMeshInterp,
const meshToMesh::order& mapOrder
const meshToMesh::order& mapOrder,
const CombineOp& cop
)
{
const fvMesh& meshSource = meshToMeshInterp.fromMesh();
@ -84,7 +85,13 @@ void MapConsistentVolFields
);
// Interpolate field
meshToMeshInterp.interpolate(fieldTarget, fieldSource, mapOrder);
meshToMeshInterp.interpolate//<Type, eqOp<Type> >
(
fieldTarget,
fieldSource,
mapOrder,
cop
);
// Write field
fieldTarget.write();
@ -97,7 +104,12 @@ void MapConsistentVolFields
GeometricField<Type, fvPatchField, volMesh> fieldTarget
(
fieldTargetIOobject,
meshToMeshInterp.interpolate(fieldSource, mapOrder)
meshToMeshInterp.interpolate//<Type, eqOp<Type> >
(
fieldSource,
mapOrder,
cop
)
);
// Write field

263
applications/utilities/preProcessing/mapFields/MapMeshes.H Normal file
View File

@ -0,0 +1,263 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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.
OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#ifndef MapMeshes_H
#define MapMeshes_H
#include "MapVolFields.H"
#include "MapConsistentVolFields.H"
#include "mapLagrangian.H"
#include "UnMapped.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<template<class> class CombineOp>
void MapConsistentMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh::order& mapOrder
)
{
// Create the interpolation scheme
meshToMesh meshToMeshInterp(meshSource, meshTarget);
Info<< nl
<< "Consistently creating and mapping fields for time "
<< meshSource.time().timeName() << nl << endl;
{
// Search for list of objects for this time
IOobjectList objects(meshSource, meshSource.time().timeName());
// Map volFields
// ~~~~~~~~~~~~~
MapConsistentVolFields<scalar>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<scalar>()
);
MapConsistentVolFields<vector>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<vector>()
);
MapConsistentVolFields<sphericalTensor>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<sphericalTensor>()
);
MapConsistentVolFields<symmTensor>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<symmTensor>()
);
MapConsistentVolFields<tensor>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<tensor>()
);
}
{
// Search for list of target objects for this time
IOobjectList objects(meshTarget, meshTarget.time().timeName());
// Mark surfaceFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<surfaceScalarField>(objects);
UnMapped<surfaceVectorField>(objects);
UnMapped<surfaceSphericalTensorField>(objects);
UnMapped<surfaceSymmTensorField>(objects);
UnMapped<surfaceTensorField>(objects);
// Mark pointFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<pointScalarField>(objects);
UnMapped<pointVectorField>(objects);
UnMapped<pointSphericalTensorField>(objects);
UnMapped<pointSymmTensorField>(objects);
UnMapped<pointTensorField>(objects);
}
mapLagrangian(meshToMeshInterp);
}
template<template<class> class CombineOp>
void MapSubMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const HashTable<word>& patchMap,
const wordList& cuttingPatches,
const meshToMesh::order& mapOrder
)
{
// Create the interpolation scheme
meshToMesh meshToMeshInterp
(
meshSource,
meshTarget,
patchMap,
cuttingPatches
);
Info<< nl
<< "Mapping fields for time " << meshSource.time().timeName()
<< nl << endl;
{
// Search for list of source objects for this time
IOobjectList objects(meshSource, meshSource.time().timeName());
// Map volFields
// ~~~~~~~~~~~~~
MapVolFields<scalar>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<scalar>()
);
MapVolFields<vector>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<vector>()
);
MapVolFields<sphericalTensor>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<sphericalTensor>()
);
MapVolFields<symmTensor>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<symmTensor>()
);
MapVolFields<tensor>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<tensor>()
);
}
{
// Search for list of target objects for this time
IOobjectList objects(meshTarget, meshTarget.time().timeName());
// Mark surfaceFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<surfaceScalarField>(objects);
UnMapped<surfaceVectorField>(objects);
UnMapped<surfaceSphericalTensorField>(objects);
UnMapped<surfaceSymmTensorField>(objects);
UnMapped<surfaceTensorField>(objects);
// Mark pointFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<pointScalarField>(objects);
UnMapped<pointVectorField>(objects);
UnMapped<pointSphericalTensorField>(objects);
UnMapped<pointSymmTensorField>(objects);
UnMapped<pointTensorField>(objects);
}
mapLagrangian(meshToMeshInterp);
}
template<template<class> class CombineOp>
void MapConsistentSubMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh::order& mapOrder
)
{
HashTable<word> patchMap;
HashTable<label> cuttingPatchTable;
forAll(meshTarget.boundary(), patchi)
{
if (!isA<processorFvPatch>(meshTarget.boundary()[patchi]))
{
patchMap.insert
(
meshTarget.boundary()[patchi].name(),
meshTarget.boundary()[patchi].name()
);
}
else
{
cuttingPatchTable.insert
(
meshTarget.boundaryMesh()[patchi].name(),
-1
);
}
}
MapSubMesh<CombineOp>
(
meshSource,
meshTarget,
patchMap,
cuttingPatchTable.toc(),
mapOrder
);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

13
applications/utilities/preProcessing/mapFields/MapVolFields.H
View File

@ -35,12 +35,13 @@ License
namespace Foam
{
template<class Type>
template<class Type, class CombineOp>
void MapVolFields
(
const IOobjectList& objects,
const meshToMesh& meshToMeshInterp,
const meshToMesh::order& mapOrder
const meshToMesh::order& mapOrder,
const CombineOp& cop
)
{
const fvMesh& meshSource = meshToMeshInterp.fromMesh();
@ -84,7 +85,13 @@ void MapVolFields
);
// Interpolate field
meshToMeshInterp.interpolate(fieldTarget, fieldSource, mapOrder);
meshToMeshInterp.interpolate//<Type, eqOp<Type> >
(
fieldTarget,
fieldSource,
mapOrder,
cop
);
// Write field
fieldTarget.write();

217
applications/utilities/preProcessing/mapFields/mapFields.C
View File

@ -34,11 +34,8 @@ Description
#include "fvCFD.H"
#include "meshToMesh.H"
#include "MapVolFields.H"
#include "MapConsistentVolFields.H"
#include "UnMapped.H"
#include "processorFvPatch.H"
#include "mapLagrangian.H"
#include "MapMeshes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -46,56 +43,28 @@ void mapConsistentMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh::order& mapOrder
const meshToMesh::order& mapOrder,
const bool subtract
)
{
// Create the interpolation scheme
meshToMesh meshToMeshInterp(meshSource, meshTarget);
Info<< nl
<< "Consistently creating and mapping fields for time "
<< meshSource.time().timeName() << nl << endl;
if (subtract)
{
// Search for list of objects for this time
IOobjectList objects(meshSource, meshSource.time().timeName());
// Map volFields
// ~~~~~~~~~~~~~
MapConsistentVolFields<scalar>(objects, meshToMeshInterp, mapOrder);
MapConsistentVolFields<vector>(objects, meshToMeshInterp, mapOrder);
MapConsistentVolFields<sphericalTensor>
MapConsistentMesh<minusEqOp>
(
objects,
meshToMeshInterp,
meshSource,
meshTarget,
mapOrder
);
MapConsistentVolFields<symmTensor>(objects, meshToMeshInterp, mapOrder);
MapConsistentVolFields<tensor>(objects, meshToMeshInterp, mapOrder);
}
else
{
// Search for list of target objects for this time
IOobjectList objects(meshTarget, meshTarget.time().timeName());
// Mark surfaceFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<surfaceScalarField>(objects);
UnMapped<surfaceVectorField>(objects);
UnMapped<surfaceSphericalTensorField>(objects);
UnMapped<surfaceSymmTensorField>(objects);
UnMapped<surfaceTensorField>(objects);
// Mark pointFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<pointScalarField>(objects);
UnMapped<pointVectorField>(objects);
UnMapped<pointSphericalTensorField>(objects);
UnMapped<pointSymmTensorField>(objects);
UnMapped<pointTensorField>(objects);
MapConsistentMesh<eqOp>
(
meshSource,
meshTarget,
mapOrder
);
}
mapLagrangian(meshToMeshInterp);
}
@ -105,57 +74,32 @@ void mapSubMesh
const fvMesh& meshTarget,
const HashTable<word>& patchMap,
const wordList& cuttingPatches,
const meshToMesh::order& mapOrder
const meshToMesh::order& mapOrder,
const bool subtract
)
{
// Create the interpolation scheme
meshToMesh meshToMeshInterp
(
meshSource,
meshTarget,
patchMap,
cuttingPatches
);
Info<< nl
<< "Mapping fields for time " << meshSource.time().timeName()
<< nl << endl;
if (subtract)
{
// Search for list of source objects for this time
IOobjectList objects(meshSource, meshSource.time().timeName());
// Map volFields
// ~~~~~~~~~~~~~
MapVolFields<scalar>(objects, meshToMeshInterp, mapOrder);
MapVolFields<vector>(objects, meshToMeshInterp, mapOrder);
MapVolFields<sphericalTensor>(objects, meshToMeshInterp, mapOrder);
MapVolFields<symmTensor>(objects, meshToMeshInterp, mapOrder);
MapVolFields<tensor>(objects, meshToMeshInterp, mapOrder);
MapSubMesh<minusEqOp>
(
meshSource,
meshTarget,
patchMap,
cuttingPatches,
mapOrder
);
}
else
{
// Search for list of target objects for this time
IOobjectList objects(meshTarget, meshTarget.time().timeName());
// Mark surfaceFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<surfaceScalarField>(objects);
UnMapped<surfaceVectorField>(objects);
UnMapped<surfaceSphericalTensorField>(objects);
UnMapped<surfaceSymmTensorField>(objects);
UnMapped<surfaceTensorField>(objects);
// Mark pointFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<pointScalarField>(objects);
UnMapped<pointVectorField>(objects);
UnMapped<pointSphericalTensorField>(objects);
UnMapped<pointSymmTensorField>(objects);
UnMapped<pointTensorField>(objects);
MapSubMesh<eqOp>
(
meshSource,
meshTarget,
patchMap,
cuttingPatches,
mapOrder
);
}
mapLagrangian(meshToMeshInterp);
}
@ -163,40 +107,28 @@ void mapConsistentSubMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh::order& mapOrder
const meshToMesh::order& mapOrder,
const bool subtract
)
{
HashTable<word> patchMap;
HashTable<label> cuttingPatchTable;
forAll(meshTarget.boundary(), patchi)
if (subtract)
{
if (!isA<processorFvPatch>(meshTarget.boundary()[patchi]))
{
patchMap.insert
(
meshTarget.boundary()[patchi].name(),
meshTarget.boundary()[patchi].name()
);
}
else
{
cuttingPatchTable.insert
(
meshTarget.boundaryMesh()[patchi].name(),
-1
);
}
MapConsistentSubMesh<minusEqOp>
(
meshSource,
meshTarget,
mapOrder
);
}
else
{
MapConsistentSubMesh<eqOp>
(
meshSource,
meshTarget,
mapOrder
);
}
mapSubMesh
(
meshSource,
meshTarget,
patchMap,
cuttingPatchTable.toc(),
mapOrder
);
}
@ -288,6 +220,11 @@ int main(int argc, char *argv[])
"word",
"specify the mapping method"
);
argList::addBoolOption
(
"subtract",
"subtract mapped source from target"
);
argList args(argc, argv);
@ -350,6 +287,13 @@ int main(int argc, char *argv[])
Info<< "Mapping method: " << mapMethod << endl;
}
const bool subtract = args.optionFound("subtract");
if (subtract)
{
Info<< "Subtracting mapped source field from target" << endl;
}
#include "createTimes.H"
HashTable<word> patchMap;
@ -431,7 +375,13 @@ int main(int argc, char *argv[])
if (consistent)
{
mapConsistentSubMesh(meshSource, meshTarget, mapOrder);
mapConsistentSubMesh
(
meshSource,
meshTarget,
mapOrder,
subtract
);
}
else
{
@ -441,7 +391,8 @@ int main(int argc, char *argv[])
meshTarget,
patchMap,
cuttingPatches,
mapOrder
mapOrder,
subtract
);
}
}
@ -503,7 +454,13 @@ int main(int argc, char *argv[])
if (consistent)
{
mapConsistentSubMesh(meshSource, meshTarget, mapOrder);
mapConsistentSubMesh
(
meshSource,
meshTarget,
mapOrder,
subtract
);
}
else
{
@ -513,7 +470,8 @@ int main(int argc, char *argv[])
meshTarget,
patchMap,
addProcessorPatches(meshTarget, cuttingPatches),
mapOrder
mapOrder,
subtract
);
}
}
@ -629,7 +587,8 @@ int main(int argc, char *argv[])
(
meshSource,
meshTarget,
mapOrder
mapOrder,
subtract
);
}
else
@ -640,7 +599,8 @@ int main(int argc, char *argv[])
meshTarget,
patchMap,
addProcessorPatches(meshTarget, cuttingPatches),
mapOrder
mapOrder,
subtract
);
}
}
@ -679,7 +639,7 @@ int main(int argc, char *argv[])
if (consistent)
{
mapConsistentMesh(meshSource, meshTarget, mapOrder);
mapConsistentMesh(meshSource, meshTarget, mapOrder, subtract);
}
else
{
@ -689,7 +649,8 @@ int main(int argc, char *argv[])
meshTarget,
patchMap,
cuttingPatches,
mapOrder
mapOrder,
subtract
);
}
}

2
src/Allwmake
View File

@ -46,6 +46,8 @@ wmake $makeType genericPatchFields
# Build the proper scotchDecomp, metisDecomp etc.
parallel/Allwmake $*
renumber/Allwmake $*
wmake $makeType conversion
wmake $makeType sampling

1
src/OpenFOAM/Make/files
View File

@ -250,6 +250,7 @@ $(lduMatrix)/solvers/ICCG/ICCG.C
$(lduMatrix)/solvers/BICCG/BICCG.C
$(lduMatrix)/smoothers/GaussSeidel/GaussSeidelSmoother.C
$(lduMatrix)/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.C
$(lduMatrix)/smoothers/DIC/DICSmoother.C
$(lduMatrix)/smoothers/DICGaussSeidel/DICGaussSeidelSmoother.C
$(lduMatrix)/smoothers/DILU/DILUSmoother.C

3
src/OpenFOAM/fields/FieldFields/symmTensorFieldField/symmTensorFieldField.C
View File

@ -55,6 +55,9 @@ UNARY_FUNCTION(symmTensor, symmTensor, inv)
UNARY_OPERATOR(vector, symmTensor, *, hdual)
BINARY_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
BINARY_TYPE_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

3
src/OpenFOAM/fields/FieldFields/symmTensorFieldField/symmTensorFieldField.H
View File

@ -65,6 +65,9 @@ UNARY_FUNCTION(symmTensor, symmTensor, inv)
UNARY_OPERATOR(vector, symmTensor, *, hdual)
BINARY_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
BINARY_TYPE_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

3
src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
View File

@ -164,6 +164,9 @@ tmp<Field<symmTensor> > transformFieldMask<symmTensor>
UNARY_OPERATOR(vector, symmTensor, *, hdual)
BINARY_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
BINARY_TYPE_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

4
src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.H
View File

@ -39,6 +39,7 @@ SourceFiles
#include "vectorField.H"
#include "sphericalTensor.H"
#include "symmTensor.H"
#include "tensor.H"
#define TEMPLATE
#include "FieldFunctionsM.H"
@ -69,6 +70,9 @@ UNARY_FUNCTION(symmTensor, symmTensor, inv)
UNARY_OPERATOR(vector, symmTensor, *, hdual)
BINARY_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
BINARY_TYPE_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

3
src/OpenFOAM/fields/GeometricFields/GeometricSymmTensorField/GeometricSymmTensorField.C
View File

@ -53,6 +53,9 @@ UNARY_FUNCTION(symmTensor, symmTensor, inv, inv)
UNARY_OPERATOR(vector, symmTensor, *, hdual, transform)
BINARY_OPERATOR(tensor, symmTensor, symmTensor, &, '&', dot)
BINARY_TYPE_OPERATOR(tensor, symmTensor, symmTensor, &, '&', dot)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

4
src/OpenFOAM/fields/GeometricFields/GeometricSymmTensorField/GeometricSymmTensorField.H
View File

@ -65,6 +65,10 @@ UNARY_FUNCTION(symmTensor, symmTensor, inv, inv)
UNARY_OPERATOR(vector, symmTensor, *, hdual, transform)
BINARY_OPERATOR(tensor, symmTensor, symmTensor, &, '&', dot)
BINARY_TYPE_OPERATOR(tensor, symmTensor, symmTensor, &, '&', dot)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

28
src/OpenFOAM/matrices/lduMatrix/smoothers/GaussSeidel/GaussSeidelSmoother.C
View File

@ -107,16 +107,29 @@ void Foam::GaussSeidelSmoother::smooth
// To compensate for this, it is necessary to turn the
// sign of the contribution.
FieldField<Field, scalar> mBouCoeffs(interfaceBouCoeffs_.size());
//FieldField<Field, scalar> mBouCoeffs(interfaceBouCoeffs_.size());
//
//forAll(mBouCoeffs, patchi)
//{
// if (interfaces_.set(patchi))
// {
// mBouCoeffs.set(patchi, -interfaceBouCoeffs_[patchi]);
// }
//}
FieldField<Field, scalar>& mBouCoeffs =
const_cast<FieldField<Field, scalar>&>
(
interfaceBouCoeffs_
);
forAll(mBouCoeffs, patchi)
{
if (interfaces_.set(patchi))
{
mBouCoeffs.set(patchi, -interfaceBouCoeffs_[patchi]);
mBouCoeffs[patchi].negate();
}
}
for (label sweep=0; sweep<nSweeps; sweep++)
{
bPrime = source;
@ -170,6 +183,15 @@ void Foam::GaussSeidelSmoother::smooth
psiPtr[cellI] = curPsi;
}
}
// Restore interfaceBouCoeffs_
forAll(mBouCoeffs, patchi)
{
if (interfaces_.set(patchi))
{
mBouCoeffs[patchi].negate();
}
}
}

266
src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.C Normal file
View File

@ -0,0 +1,266 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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.
OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "nonBlockingGaussSeidelSmoother.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(nonBlockingGaussSeidelSmoother, 0);
lduMatrix::smoother::
addsymMatrixConstructorToTable<nonBlockingGaussSeidelSmoother>
addnonBlockingGaussSeidelSmootherSymMatrixConstructorToTable_;
lduMatrix::smoother::
addasymMatrixConstructorToTable<nonBlockingGaussSeidelSmoother>
addnonBlockingGaussSeidelSmootherAsymMatrixConstructorToTable_;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::nonBlockingGaussSeidelSmoother::nonBlockingGaussSeidelSmoother
(
const word& fieldName,
const lduMatrix& matrix,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces
)
:
lduMatrix::smoother
(
fieldName,
matrix,
interfaceBouCoeffs,
interfaceIntCoeffs,
interfaces
)
{
// Check that all interface addressing is sorted to be after the
// non-interface addressing.
register const label nCells = matrix.diag().size();
blockStart_ = nCells;
labelList startCellI(interfaceBouCoeffs.size(), -1);
forAll(interfaces, patchi)
{
if (interfaces.set(patchi))
{
const labelUList& faceCells = matrix_.lduAddr().patchAddr(patchi);
blockStart_ = min(blockStart_, min(faceCells));
}
}
if (debug)
{
Pout<< "nonBlockingGaussSeidelSmoother :"
<< " Starting block on cell " << blockStart_
<< " out of " << nCells << endl;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::nonBlockingGaussSeidelSmoother::smooth
(
const word& fieldName_,
scalarField& psi,
const lduMatrix& matrix_,
const label blockStart,
const scalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs_,
const lduInterfaceFieldPtrsList& interfaces_,
const direction cmpt,
const label nSweeps
)
{
register scalar* __restrict__ psiPtr = psi.begin();
register const label nCells = psi.size();
scalarField bPrime(nCells);
register scalar* __restrict__ bPrimePtr = bPrime.begin();
register const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
register const scalar* const __restrict__ upperPtr =
matrix_.upper().begin();
register const scalar* const __restrict__ lowerPtr =
matrix_.lower().begin();
register const label* const __restrict__ uPtr =
matrix_.lduAddr().upperAddr().begin();
register const label* const __restrict__ ownStartPtr =
matrix_.lduAddr().ownerStartAddr().begin();
// Parallel boundary initialisation. The parallel boundary is treated
// as an effective jacobi interface in the boundary.
// Note: there is a change of sign in the coupled
// interface update. The reason for this is that the
// internal coefficients are all located at the l.h.s. of
// the matrix whereas the "implicit" coefficients on the
// coupled boundaries are all created as if the
// coefficient contribution is of a source-kind (i.e. they
// have a sign as if they are on the r.h.s. of the matrix.
// To compensate for this, it is necessary to turn the
// sign of the contribution.
FieldField<Field, scalar>& mBouCoeffs =
const_cast<FieldField<Field, scalar>&>
(
interfaceBouCoeffs_
);
forAll(mBouCoeffs, patchi)
{
if (interfaces_.set(patchi))
{
mBouCoeffs[patchi].negate();
}
}
for (label sweep=0; sweep<nSweeps; sweep++)
{
bPrime = source;
matrix_.initMatrixInterfaces
(
mBouCoeffs,
interfaces_,
psi,
bPrime,
cmpt
);
register scalar curPsi;
register label fStart;
register label fEnd = ownStartPtr[0];
for (register label cellI=0; cellI<blockStart; cellI++)
{
// Start and end of this row
fStart = fEnd;
fEnd = ownStartPtr[cellI + 1];
// Get the accumulated neighbour side
curPsi = bPrimePtr[cellI];
// Accumulate the owner product side
for (register label curFace=fStart; curFace<fEnd; curFace++)
{
curPsi -= upperPtr[curFace]*psiPtr[uPtr[curFace]];
}
// Finish current psi
curPsi /= diagPtr[cellI];
// Distribute the neighbour side using current psi
for (register label curFace=fStart; curFace<fEnd; curFace++)
{
bPrimePtr[uPtr[curFace]] -= lowerPtr[curFace]*curPsi;
}
psiPtr[cellI] = curPsi;
}
matrix_.updateMatrixInterfaces
(
mBouCoeffs,
interfaces_,
psi,
bPrime,
cmpt
);
// Update rest of the cells
for (label cellI=blockStart; cellI < nCells; cellI++)
{
// Start and end of this row
fStart = fEnd;
fEnd = ownStartPtr[cellI + 1];
// Get the accumulated neighbour side
curPsi = bPrimePtr[cellI];
// Accumulate the owner product side
for (register label curFace=fStart; curFace<fEnd; curFace++)
{
curPsi -= upperPtr[curFace]*psiPtr[uPtr[curFace]];
}
// Finish current psi
curPsi /= diagPtr[cellI];
// Distribute the neighbour side using current psi
for (register label curFace=fStart; curFace<fEnd; curFace++)
{
bPrimePtr[uPtr[curFace]] -= lowerPtr[curFace]*curPsi;
}
psiPtr[cellI] = curPsi;
}
}
// Restore interfaceBouCoeffs_
forAll(mBouCoeffs, patchi)
{
if (interfaces_.set(patchi))
{
mBouCoeffs[patchi].negate();
}
}
}
void Foam::nonBlockingGaussSeidelSmoother::smooth
(
scalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
) const
{
smooth
(
fieldName_,
psi,
matrix_,
blockStart_,
source,
interfaceBouCoeffs_,
interfaces_,
cmpt,
nSweeps
);
}
// ************************************************************************* //

120
src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.H Normal file
View File

@ -0,0 +1,120 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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.
OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::nonBlockingGaussSeidelSmoother
Description
Variant of gaussSeidelSmoother that expects processor boundary
cells to be sorted last and so can block later. Only when the
cells are actually visited does it need the results to be present.
It is expected that there is little benefit to be gained from doing
this on a patch by patch basis since the number of processor interfaces
is quite small and the overhead of checking whether a processor interface
is finished might be quite high (call into mpi). Also this would
require a dynamic memory allocation to store the state of the outstanding
requests.
SourceFiles
nonBlockingGaussSeidelSmoother.C
\*---------------------------------------------------------------------------*/
#ifndef nonBlockingGaussSeidelSmoother_H
#define nonBlockingGaussSeidelSmoother_H
#include "lduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class nonBlockingGaussSeidelSmoother Declaration
\*---------------------------------------------------------------------------*/
class nonBlockingGaussSeidelSmoother
:
public lduMatrix::smoother
{
// Private data
//- Starting cell when to block
label blockStart_;
public:
//- Runtime type information
TypeName("nonBlockingGaussSeidel");
// Constructors
//- Construct from components
nonBlockingGaussSeidelSmoother
(
const word& fieldName,
const lduMatrix& matrix,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces
);
// Member Functions
//- Smooth for the given number of sweeps
static void smooth
(
const word& fieldName,
scalarField& psi,
const lduMatrix& matrix,
const label blockStart,
const scalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const direction cmpt,
const label nSweeps
);
//- Smooth the solution for a given number of sweeps
virtual void smooth
(
scalarField& psi,
const scalarField& Source,
const direction cmpt,
const label nSweeps
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

6
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
View File

@ -144,11 +144,13 @@ public:
inline tetPointRef oldTet(const polyMesh& mesh) const;
//- Return the geometry corresponding to the tri on the
// mesh face for this tet from the supplied mesh
// mesh face for this tet from the supplied mesh. Normal of
// the tri points out of the cell.
inline triPointRef faceTri(const polyMesh& mesh) const;
//- Return the point indices corresponding to the tri on the mesh
// face for this tet from the supplied mesh
// face for this tet from the supplied mesh. Normal of
// the tri points out of the cell.
inline triFace faceTriIs(const polyMesh& mesh) const;
//- Return the geometry corresponding to the tri on the

4
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
View File

@ -266,7 +266,7 @@ void Foam::cyclicPolyPatch::calcTransforms
// Calculate using the given rotation axis and centre. Do not
// use calculated normals.
vector n0 = findFaceMaxRadius(half0Ctrs);
vector n1 = findFaceMaxRadius(half1Ctrs);
vector n1 = -findFaceMaxRadius(half1Ctrs);
n0 /= mag(n0) + VSMALL;
n1 /= mag(n1) + VSMALL;
@ -424,7 +424,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
case ROTATIONAL:
{
vector n0 = findFaceMaxRadius(half0Ctrs);
vector n1 = findFaceMaxRadius(half1Ctrs);
vector n1 = -findFaceMaxRadius(half1Ctrs);
n0 /= mag(n0) + VSMALL;
n1 /= mag(n1) + VSMALL;

6
src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
View File

@ -277,9 +277,9 @@ Foam::scalar Foam::triangle<Point, PointRef>::barycentric
bary.setSize(3);
bary[0] = (d11*d20 - d01*d21)/denom;
bary[1] = (d00*d21 - d01*d20)/denom;
bary[2] = 1.0 - bary[0] - bary[1];
bary[1] = (d11*d20 - d01*d21)/denom;
bary[2] = (d00*d21 - d01*d20)/denom;
bary[0] = 1.0 - bary[1] - bary[2];
return denom;
}

8
src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H
View File

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "Vector.H"
#include "Tensor.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -184,18 +185,21 @@ inline Vector<Cmpt> operator*(const SymmTensor<Cmpt>& st)
//- Inner-product between two symmetric tensors
template <class Cmpt>
inline SymmTensor<Cmpt>
inline Tensor<Cmpt>
operator&(const SymmTensor<Cmpt>& st1, const SymmTensor<Cmpt>& st2)
{
return SymmTensor<Cmpt>
return Tensor<Cmpt>
(
st1.xx()*st2.xx() + st1.xy()*st2.xy() + st1.xz()*st2.xz(),
st1.xx()*st2.xy() + st1.xy()*st2.yy() + st1.xz()*st2.yz(),
st1.xx()*st2.xz() + st1.xy()*st2.yz() + st1.xz()*st2.zz(),
st1.xy()*st2.xx() + st1.yy()*st2.xy() + st1.yz()*st2.xz(),
st1.xy()*st2.xy() + st1.yy()*st2.yy() + st1.yz()*st2.yz(),
st1.xy()*st2.xz() + st1.yy()*st2.yz() + st1.yz()*st2.zz(),
st1.xz()*st2.xx() + st1.yz()*st2.xy() + st1.zz()*st2.xz(),
st1.xz()*st2.xy() + st1.yz()*st2.yy() + st1.zz()*st2.yz(),
st1.xz()*st2.xz() + st1.yz()*st2.yz() + st1.zz()*st2.zz()
);
}

4
src/OpenFOAM/primitives/Tensor/Tensor.H
View File

@ -40,13 +40,15 @@ SourceFiles
#include "Vector.H"
#include "SphericalTensor.H"
#include "SymmTensor.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<class Cmpt>
class SymmTensor;
/*---------------------------------------------------------------------------*\
Class Tensor Declaration
\*---------------------------------------------------------------------------*/

2
src/OpenFOAM/primitives/Tensor/TensorI.H
View File

@ -23,6 +23,8 @@ License
\*---------------------------------------------------------------------------*/
#include "SymmTensor.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam

86
src/OpenFOAM/primitives/functions/DataEntry/CompatibilityConstant/CompatibilityConstant.C Normal file
View File

@ -0,0 +1,86 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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.
OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "CompatibilityConstant.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
Foam::CompatibilityConstant<Type>::CompatibilityConstant
(
const word& entryName, const dictionary& dict
)
:
DataEntry<Type>(entryName),
value_(pTraits<Type>::zero)
{
dict.lookup(entryName) >> value_;
}
template<class Type>
Foam::CompatibilityConstant<Type>::CompatibilityConstant
(
const CompatibilityConstant<Type>& cnst
)
:
DataEntry<Type>(cnst),
value_(cnst.value_)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
Foam::CompatibilityConstant<Type>::~CompatibilityConstant()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Type Foam::CompatibilityConstant<Type>::value(const scalar x) const
{
return value_;
}
template<class Type>
Type Foam::CompatibilityConstant<Type>::integrate
(
const scalar x1,
const scalar x2
) const
{
return (x2 - x1)*value_;
}
// * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
#include "CompatibilityConstantIO.C"
// ************************************************************************* //

143
src/OpenFOAM/primitives/functions/DataEntry/CompatibilityConstant/CompatibilityConstant.H Normal file
View File

@ -0,0 +1,143 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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.
OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::CompatibilityConstant
Description
Templated basic entry that holds a constant value for backwards
compatibility (when DataEntry type is not present)
Usage - for entry \<entryName\> having the value <value>:
\verbatim
<entryName> <value>
\endverbatim
SourceFiles
CompatibilityConstant.C
\*---------------------------------------------------------------------------*/
#ifndef CompatibilityConstant_H
#define CompatibilityConstant_H
#include "DataEntry.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<class Type>
class CompatibilityConstant;
template<class Type>
Ostream& operator<<(Ostream&, const CompatibilityConstant<Type>&);
/*---------------------------------------------------------------------------*\
Class CompatibilityConstant Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class CompatibilityConstant
:
public DataEntry<Type>
{
// Private data
//- Constant value
Type value_;
// Private Member Functions
//- Disallow default bitwise assignment
void operator=(const CompatibilityConstant<Type>&);
public:
// Runtime type information
TypeName("CompatibilityConstant");
// Constructors
//- Construct from entry name and Istream
CompatibilityConstant(const word& entryName, const dictionary& dict);
//- Copy constructor
CompatibilityConstant(const CompatibilityConstant<Type>& cnst);
//- Construct and return a clone
virtual tmp<DataEntry<Type> > clone() const
{
return tmp<DataEntry<Type> >
(
new CompatibilityConstant<Type>(*this)
);
}
//- Destructor
virtual ~CompatibilityConstant();
// Member Functions
//- Return constant value
Type value(const scalar) const;
//- Integrate between two values
Type integrate(const scalar x1, const scalar x2) const;
// I/O
//- Ostream Operator
friend Ostream& operator<< <Type>
(
Ostream& os,
const CompatibilityConstant<Type>& cnst
);
//- Write in dictionary format
virtual void writeData(Ostream& os) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "CompatibilityConstant.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

69
src/OpenFOAM/primitives/functions/DataEntry/CompatibilityConstant/CompatibilityConstantIO.C Normal file
View File

@ -0,0 +1,69 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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.
OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "DataEntry.H"
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
template<class Type>
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const CompatibilityConstant<Type>& cnst
)
{
if (os.format() == IOstream::ASCII)
{
os << static_cast<const DataEntry<Type>& >(cnst)
<< token::SPACE << cnst.value_;
}
else
{
os << static_cast<const DataEntry<Type>& >(cnst);
os.write
(
reinterpret_cast<const char*>(&cnst.value_),
sizeof(cnst.value_)
);
}
// Check state of Ostream
os.check
(
"Ostream& operator<<(Ostream&, const CompatibilityConstant<Type>&)"
);
return os;
}
template<class Type>
void Foam::CompatibilityConstant<Type>::writeData(Ostream& os) const
{
os.writeKeyword(this->name_) << value_ << token::END_STATEMENT << nl;
}
// ************************************************************************* //

21
src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntryNew.C
View File

@ -35,17 +35,28 @@ Foam::autoPtr<Foam::DataEntry<Type> > Foam::DataEntry<Type>::New
)
{
Istream& is(dict.lookup(entryName));
word DataEntryType(is);
token firstToken(is);
word DataEntryType;
if (firstToken.isWord())
{
DataEntryType = firstToken.wordToken();
}
else
{
is.putBack(firstToken);
// DataEntryType = CompatibilityConstant<Type>::typeName;
DataEntryType = "CompatibilityConstant";
}
typename dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(DataEntryType);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorIn
(
"DataEntry<Type>::New(Istream&)"
) << "Unknown DataEntry type "
FatalErrorIn("DataEntry<Type>::New(const word&, const dictionary&)")
<< "Unknown DataEntry type "
<< DataEntryType << " for DataEntry "
<< entryName << nl << nl
<< "Valid DataEntry types are:" << nl

7
src/OpenFOAM/primitives/functions/DataEntry/makeDataEntries.C
View File

@ -23,6 +23,7 @@ License
\*---------------------------------------------------------------------------*/
#include "CompatibilityConstant.H"
#include "Constant.H"
#include "CSV.H"
#include "DataEntry.H"
@ -39,36 +40,42 @@ License
namespace Foam
{
makeDataEntry(label);
makeDataEntryType(CompatibilityConstant, label);
makeDataEntryType(Constant, label);
// makeDataEntryType(CSV, label);
makeDataEntryType(Table, label);
makeDataEntryType(TableFile, label);
makeDataEntry(scalar);
makeDataEntryType(CompatibilityConstant, scalar);
makeDataEntryType(Constant, scalar);
makeDataEntryType(CSV, scalar);
makeDataEntryType(Table, scalar);
makeDataEntryType(TableFile, scalar);
makeDataEntry(vector);
makeDataEntryType(CompatibilityConstant, vector);
makeDataEntryType(Constant, vector);
makeDataEntryType(CSV, vector);
makeDataEntryType(Table, vector);
makeDataEntryType(TableFile, vector);
makeDataEntry(sphericalTensor);
makeDataEntryType(CompatibilityConstant, sphericalTensor);
makeDataEntryType(Constant, sphericalTensor);
makeDataEntryType(CSV, sphericalTensor);
makeDataEntryType(Table, sphericalTensor);
makeDataEntryType(TableFile, sphericalTensor);
makeDataEntry(symmTensor);
makeDataEntryType(CompatibilityConstant, symmTensor);
makeDataEntryType(Constant, symmTensor);
makeDataEntryType(CSV, symmTensor);
makeDataEntryType(Table, symmTensor);
makeDataEntryType(TableFile, symmTensor);
makeDataEntry(tensor);
makeDataEntryType(CompatibilityConstant, tensor);
makeDataEntryType(Constant, tensor);
makeDataEntryType(CSV, tensor);
makeDataEntryType(Table, tensor);

12
src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
View File

@ -44,11 +44,15 @@ void Foam::MRFZone::relativeRhoFlux
const vector& origin = origin_.value();
const vector& Omega = Omega_.value();
const vectorField& Cfi = Cf.internalField();
const vectorField& Sfi = Sf.internalField();
scalarField& phii = phi.internalField();
// Internal faces
forAll(internalFaces_, i)
{
label facei = internalFaces_[i];
phi[facei] -= rho[facei]*(Omega ^ (Cf[facei] - origin)) & Sf[facei];
phii[facei] -= rho[facei]*(Omega ^ (Cfi[facei] - origin)) & Sfi[facei];
}
// Included patches
@ -91,11 +95,15 @@ void Foam::MRFZone::absoluteRhoFlux
const vector& origin = origin_.value();
const vector& Omega = Omega_.value();
const vectorField& Cfi = Cf.internalField();
const vectorField& Sfi = Sf.internalField();
scalarField& phii = phi.internalField();
// Internal faces
forAll(internalFaces_, i)
{
label facei = internalFaces_[i];
phi[facei] += rho[facei]*(Omega ^ (Cf[facei] - origin)) & Sf[facei];
phii[facei] += rho[facei]*(Omega ^ (Cfi[facei] - origin)) & Sfi[facei];
}
// Included patches

3
src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
View File

@ -196,9 +196,10 @@ bool Foam::pimpleControl::loop()
bool completed = false;
if (criteriaSatisfied())
{
Info<< algorithmName_ << ": converged in " << corr_ << " iterations"
Info<< algorithmName_ << ": converged in " << corr_ - 1 << " iterations"
<< endl;
completed = true;
corr_ = 0;
}
else
{

2
src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrad.H
View File

@ -65,7 +65,7 @@ class orthogonalSnGrad
public:
//- Runtime type information
TypeName("uncorrected");
TypeName("orthogonal");
// Constructors

12
src/finiteVolume/fvMesh/fvMesh.C
View File

@ -59,6 +59,8 @@ License
#include "upwindFECCellToFaceStencilObject.H"
#include "centredCFCFaceToCellStencilObject.H"
#include "meshSearchMeshObject.H"
#include "meshSearchFACECENTRETETSMeshObject.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -99,6 +101,11 @@ void Foam::fvMesh::clearGeom()
CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
skewCorrectionVectors::Delete(*this);
//quadraticFitSnGradData::Delete(*this);
// Note: should be in polyMesh::clearGeom but meshSearch not in OpenFOAM
// library
meshSearchMeshObject::Delete(*this);
meshSearchFACECENTRETETSMeshObject::Delete(*this);
}
@ -128,6 +135,11 @@ void Foam::fvMesh::clearAddressing()
upwindFECCellToFaceStencilObject::Delete(*this);
centredCFCFaceToCellStencilObject::Delete(*this);
// Note: should be in polyMesh::clearGeom but meshSearch not in OpenFOAM
// library
meshSearchMeshObject::Delete(*this);
meshSearchFACECENTRETETSMeshObject::Delete(*this);
}

8
src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C
View File

@ -83,8 +83,8 @@ void Foam::cellPointWeight::findTetrahedron
)
{
faceVertices_[0] = f[tetIs.faceBasePt()];
faceVertices_[1] = f[tetIs.facePtA()];;
faceVertices_[2] = f[tetIs.facePtB()];;
faceVertices_[1] = f[tetIs.facePtA()];
faceVertices_[2] = f[tetIs.facePtB()];
return;
}
@ -191,8 +191,8 @@ void Foam::cellPointWeight::findTriangle
weights_[3] = triWeights[2];
faceVertices_[0] = f[tetIs.faceBasePt()];
faceVertices_[1] = f[tetIs.facePtA()];;
faceVertices_[2] = f[tetIs.facePtB()];;
faceVertices_[1] = f[tetIs.facePtA()];
faceVertices_[2] = f[tetIs.facePtB()];
return;
}

3
src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C
View File

@ -90,7 +90,8 @@ Foam::COxidationDiffusionLimitedRate<CloudType>::COxidationDiffusionLimitedRate
O2GlobalId_(srm.O2GlobalId_),
CO2GlobalId_(srm.CO2GlobalId_),
WC_(srm.WC_),
WO2_(srm.WO2_)
WO2_(srm.WO2_),
HcCO2_(srm.HcCO2_)
{}

3
src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C
View File

@ -81,7 +81,8 @@ COxidationKineticDiffusionLimitedRate
O2GlobalId_(srm.O2GlobalId_),
CO2GlobalId_(srm.CO2GlobalId_),
WC_(srm.WC_),
WO2_(srm.WO2_)
WO2_(srm.WO2_),
HcCO2_(srm.HcCO2_)
{}

3
src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C
View File

@ -96,7 +96,8 @@ Foam::COxidationMurphyShaddix<CloudType>::COxidationMurphyShaddix
O2GlobalId_(srm.O2GlobalId_),
CO2GlobalId_(srm.CO2GlobalId_),
WC_(srm.WC_),
WO2_(srm.WO2_)
WO2_(srm.WO2_),
HcCO2_(srm.HcCO2_)
{}

38
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
View File

@ -230,6 +230,16 @@ void Foam::KinematicCloud<CloudType>::postEvolve()
functions_.postEvolve();
solution_.nextIter();
if (this->db().time().outputTime())
{
outputProperties_.writeObject
(
IOstream::ASCII,
IOstream::currentVersion,
this->db().time().writeCompression()
);
}
}
@ -281,6 +291,18 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
IOobject::NO_WRITE
)
),
outputProperties_
(
IOobject
(
cloudName + "OutputProperties",
mesh_.time().timeName(),
"uniform"/cloud::prefix/cloudName,
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
)
),
solution_(mesh_, particleProperties_.subDict("solution")),
constProps_(particleProperties_, solution_.active()),
subModelProperties_
@ -384,6 +406,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
cloudCopyPtr_(NULL),
mesh_(c.mesh_),
particleProperties_(c.particleProperties_),
outputProperties_(c.outputProperties_),
solution_(c.solution_),
constProps_(c.constProps_),
subModelProperties_(c.subModelProperties_),
@ -460,6 +483,19 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
false
)
),
outputProperties_
(
IOobject
(
name + "OutputProperties",
mesh_.time().timeName(),
"uniform"/cloud::prefix/name,
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
)
),
solution_(mesh),
constProps_(),
subModelProperties_(dictionary::null),
@ -639,7 +675,7 @@ void Foam::KinematicCloud<CloudType>::motion(TrackData& td)
template<class CloudType>
void Foam::KinematicCloud<CloudType>::info() const
void Foam::KinematicCloud<CloudType>::info()
{
vector linearMomentum = linearMomentumOfSystem();
reduce(linearMomentum, sumOp<vector>());

17
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
View File

@ -143,6 +143,9 @@ protected:
//- Dictionary of particle properties
IOdictionary particleProperties_;
//- Dictionary of output properties
IOdictionary outputProperties_;
//- Solution properties
cloudSolution solution_;
@ -324,6 +327,12 @@ public:
//- Return particle properties dictionary
inline const IOdictionary& particleProperties() const;
//- Return output properties dictionary
inline const IOdictionary& outputProperties() const;
//- Return non-const access to the output properties dictionary
inline IOdictionary& outputProperties();
//- Return const access to the solution properties
inline const cloudSolution& solution() const;
@ -468,6 +477,9 @@ public:
//- Mean diameter Dij
inline scalar Dij(const label i, const label j) const;
//- Max diameter
inline scalar Dmax() const;
// Fields
@ -543,8 +555,11 @@ public:
template<class TrackData>
void motion(TrackData& td);
// I-O
//- Print cloud information
void info() const;
void info();
};

31
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
View File

@ -50,6 +50,21 @@ Foam::KinematicCloud<CloudType>::particleProperties() const
}
template<class CloudType>
inline const Foam::IOdictionary&
Foam::KinematicCloud<CloudType>::outputProperties() const
{
return outputProperties_;
}
template<class CloudType>
inline Foam::IOdictionary& Foam::KinematicCloud<CloudType>::outputProperties()
{
return outputProperties_;
}
template<class CloudType>
inline const Foam::cloudSolution&
Foam::KinematicCloud<CloudType>::solution() const
@ -302,6 +317,22 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::Dij
}
template<class CloudType>
inline Foam::scalar Foam::KinematicCloud<CloudType>::Dmax() const
{
scalar d = -GREAT;
forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
{
const parcelType& p = iter();
d = max(d, p.d());
}
reduce(d, maxOp<scalar>());
return d;
}
template<class CloudType>
inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
(

11
src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
View File

@ -341,20 +341,13 @@ void Foam::ReactingCloud<CloudType>::evolve()
}
template<class CloudType>
void Foam::ReactingCloud<CloudType>::addToMassPhaseChange(const scalar dMass)
{
dMassPhaseChange_ += dMass;
}
template<class CloudType>
void Foam::ReactingCloud<CloudType>::info() const
void Foam::ReactingCloud<CloudType>::info()
{
CloudType::info();
Info<< " Mass transfer phase change = "
<< returnReduce(dMassPhaseChange_, sumOp<scalar>()) << nl;
this->phaseChange().info(Info);
}

20
src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
View File

@ -215,14 +215,18 @@ public:
// Sub-models
//- Return reference to reacting composition model
//- Return const access to reacting composition model
inline const CompositionModel<ReactingCloud<CloudType> >&
composition() const;
//- Return reference to reacting phase change model
//- Return const access to reacting phase change model
inline const PhaseChangeModel<ReactingCloud<CloudType> >&
phaseChange() const;
//- Return reference to reacting phase change model
inline PhaseChangeModel<ReactingCloud<CloudType> >&
phaseChange();
// Sources
@ -259,12 +263,6 @@ public:
inline tmp<fvScalarMatrix> Srho(volScalarField& rho) const;
// Check
//- Add to cumulative phase change mass transfer
void addToMassPhaseChange(const scalar dMass);
// Cloud evolution functions
//- Set parcel thermo properties
@ -300,12 +298,12 @@ public:
//- Evolve the cloud
void evolve();
//- Print cloud information
void info() const;
// I-O
//- Print cloud information
void info();
//- Write the field data for the cloud
virtual void writeFields() const;
};

8
src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
View File

@ -57,6 +57,14 @@ Foam::ReactingCloud<CloudType>::phaseChange() const
}
template<class CloudType>
inline Foam::PhaseChangeModel<Foam::ReactingCloud<CloudType> >&
Foam::ReactingCloud<CloudType>::phaseChange()
{
return phaseChangeModel_();
}
template<class CloudType>
inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ReactingCloud<CloudType>::rhoTrans(const label i)

29
src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
View File

@ -253,33 +253,12 @@ void Foam::ReactingMultiphaseCloud<CloudType>::evolve()
template<class CloudType>
void Foam::ReactingMultiphaseCloud<CloudType>::addToMassDevolatilisation
(
const scalar dMass
)
{
dMassDevolatilisation_ += dMass;
}
template<class CloudType>
void Foam::ReactingMultiphaseCloud<CloudType>::addToMassSurfaceReaction
(
const scalar dMass
)
{
dMassSurfaceReaction_ += dMass;
}
template<class CloudType>
void Foam::ReactingMultiphaseCloud<CloudType>::info() const
void Foam::ReactingMultiphaseCloud<CloudType>::info()
{
CloudType::info();
Info<< " Mass transfer devolatilisation = "
<< returnReduce(dMassDevolatilisation_, sumOp<scalar>()) << nl;
Info<< " Mass transfer surface reaction = "
<< returnReduce(dMassSurfaceReaction_, sumOp<scalar>()) << nl;
this->devolatilisation().info(Info);
this->surfaceReaction().info(Info);
}

31
src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.H
View File

@ -214,28 +214,33 @@ public:
// Sub-models
//- Return reference to devolatilisation model
//- Return const access to devolatilisation model
inline const DevolatilisationModel
<
ReactingMultiphaseCloud<CloudType>
>&
devolatilisation() const;
//- Return reference to reacting surface reaction model
//- Return reference to devolatilisation model
inline DevolatilisationModel
<
ReactingMultiphaseCloud<CloudType>
>&
devolatilisation();
//- Return const access to reacting surface reaction model
inline const SurfaceReactionModel
<
ReactingMultiphaseCloud<CloudType>
>&
surfaceReaction() const;
// Check
//- Add to cumulative volatilisation mass transfer
void addToMassDevolatilisation(const scalar dMass);
//- Add to cumulative surface reaction transfer
void addToMassSurfaceReaction(const scalar dMass);
//- Return reference to reacting surface reaction model
inline SurfaceReactionModel
<
ReactingMultiphaseCloud<CloudType>
>&
surfaceReaction();
// Cloud evolution functions
@ -267,12 +272,12 @@ public:
//- Evolve the cloud
void evolve();
//- Print cloud information
void info() const;
// I-O
//- Print cloud information
void info();
//- Write the field data for the cloud
virtual void writeFields() const;
};

22
src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloudI.H
View File

@ -52,6 +52,17 @@ Foam::ReactingMultiphaseCloud<CloudType>::devolatilisation() const
}
template<class CloudType>
inline Foam::DevolatilisationModel
<
Foam::ReactingMultiphaseCloud<CloudType>
>&
Foam::ReactingMultiphaseCloud<CloudType>::devolatilisation()
{
return devolatilisationModel_();
}
template<class CloudType>
inline const Foam::SurfaceReactionModel
<
@ -63,4 +74,15 @@ Foam::ReactingMultiphaseCloud<CloudType>::surfaceReaction() const
}
template<class CloudType>
inline Foam::SurfaceReactionModel
<
Foam::ReactingMultiphaseCloud<CloudType>
>&
Foam::ReactingMultiphaseCloud<CloudType>::surfaceReaction()
{
return surfaceReactionModel_();
}
// ************************************************************************* //

2
src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
View File

@ -334,7 +334,7 @@ void Foam::ThermoCloud<CloudType>::evolve()
template<class CloudType>
void Foam::ThermoCloud<CloudType>::info() const
void Foam::ThermoCloud<CloudType>::info()
{
CloudType::info();
}

2
src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
View File

@ -321,7 +321,7 @@ public:
// Check
//- Print cloud information
void info() const;
void info();
};

10
src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
View File

@ -240,6 +240,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
dt,
cellI,
Res,
Prs,
Ts,
mus/rhos,
d0,
@ -368,7 +369,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T0, idG, idL, idS);
td.cloud().addToMassPhaseChange(dm);
td.cloud().phaseChange().addToPhaseChangeMass(dm);
}
return;
@ -519,7 +520,10 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
scalar dMassTot = sum(dMassDV);
td.cloud().addToMassDevolatilisation(this->nParticle_*dMassTot);
td.cloud().devolatilisation().addToDevolatilisationMass
(
this->nParticle_*dMassTot
);
Sh -= dMassTot*td.cloud().constProps().LDevol()/dt;
@ -607,7 +611,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
dMassSRCarrier
);
td.cloud().addToMassSurfaceReaction
td.cloud().surfaceReaction().addToSurfaceReactionMass
(
this->nParticle_
*(sum(dMassSRGas) + sum(dMassSRLiquid) + sum(dMassSRSolid))

9
src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
View File

@ -329,6 +329,7 @@ void Foam::ReactingParcel<ParcelType>::calc
dt,
cellI,
Res,
Prs,
Ts,
mus/rhos,
d0,
@ -384,7 +385,7 @@ void Foam::ReactingParcel<ParcelType>::calc
}
td.cloud().UTrans()[cellI] += dm*U0;
td.cloud().addToMassPhaseChange(dm);
td.cloud().phaseChange().addToPhaseChangeMass(dm);
}
return;
@ -464,6 +465,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
const scalar dt,
const label cellI,
const scalar Re,
const scalar Pr,
const scalar Ts,
const scalar nus,
const scalar d,
@ -500,11 +502,14 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
dt,
cellI,
Re,
Pr,
d,
nus,
T,
Ts,
pc_,
this->Tc_,
YComponents,
dMassPC
);
@ -514,7 +519,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
const scalar dMassTot = sum(dMassPC);
// Add to cumulative phase change mass
td.cloud().addToMassPhaseChange(this->nParticle_*dMassTot);
td.cloud().phaseChange().addToPhaseChangeMass(this->nParticle_*dMassTot);
forAll(dMassPC, i)
{

1
src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
View File

@ -210,6 +210,7 @@ protected:
const scalar dt, // timestep
const label cellI, // owner cell
const scalar Re, // Reynolds number
const scalar Pr, // Prandtl number
const scalar Ts, // Surface temperature
const scalar nus, // Surface kinematic viscosity
const scalar d, // diameter

4
src/lagrangian/intermediate/parcels/include/makeReactingParcelPhaseChangeModels.H
View File

@ -30,6 +30,7 @@ License
#include "NoPhaseChange.H"
#include "LiquidEvaporation.H"
#include "LiquidEvaporationBoil.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -37,7 +38,8 @@ License
\
makePhaseChangeModel(CloudType); \
makePhaseChangeModelType(NoPhaseChange, CloudType); \
makePhaseChangeModelType(LiquidEvaporation, CloudType);
makePhaseChangeModelType(LiquidEvaporation, CloudType); \
makePhaseChangeModelType(LiquidEvaporationBoil, CloudType);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

2
src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.C
View File

@ -51,7 +51,7 @@ Foam::CloudFunctionObject<CloudType>::CloudFunctionObject
const word& type
)
:
SubModelBase<CloudType>(owner, dict, type, "")
SubModelBase<CloudType>(owner, dict, typeName, type, "")
{}

2
src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.H
View File

@ -64,7 +64,7 @@ class CloudFunctionObject
public:
//- Runtime type information
TypeName("postProcessingModel");
TypeName("cloudFunctionObject");
//- Declare runtime constructor selection table
declareRunTimeSelectionTable

2
src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/CollisionModel/CollisionModel.C
View File

@ -42,7 +42,7 @@ Foam::CollisionModel<CloudType>::CollisionModel
const word& type
)
:
SubModelBase<CloudType>(owner, dict, type)
SubModelBase<CloudType>(owner, dict, typeName, type)
{}

2
src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/DispersionModel/DispersionModel.C
View File

@ -42,7 +42,7 @@ Foam::DispersionModel<CloudType>::DispersionModel
const word& type
)
:
SubModelBase<CloudType>(owner, dict, type)
SubModelBase<CloudType>(owner, dict, typeName, type)
{}

117
src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
View File

@ -31,76 +31,6 @@ using namespace Foam::constant::mathematical;
// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
template<class CloudType>
void Foam::InjectionModel<CloudType>::readProps()
{
if (!this->owner().solution().transient())
{
return;
}
IOobject propsDictHeader
(
"injectionProperties",
this->owner().db().time().timeName(),
"uniform"/cloud::prefix/this->owner().name(),
this->owner().db(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
);
if (propsDictHeader.headerOk())
{
const IOdictionary propsDict(propsDictHeader);
propsDict.readIfPresent("massInjected", massInjected_);
propsDict.readIfPresent("nInjections", nInjections_);
propsDict.readIfPresent("parcelsAddedTotal", parcelsAddedTotal_);
propsDict.readIfPresent("timeStep0", timeStep0_);
}
}
template<class CloudType>
void Foam::InjectionModel<CloudType>::writeProps()
{
if (!this->owner().solution().transient())
{
return;
}
if (this->owner().db().time().outputTime())
{
IOdictionary propsDict
(
IOobject
(
"injectionProperties",
this->owner().db().time().timeName(),
"uniform"/cloud::prefix/this->owner().name(),
this->owner().db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
)
);
propsDict.add("massInjected", massInjected_);
propsDict.add("nInjections", nInjections_);
propsDict.add("parcelsAddedTotal", parcelsAddedTotal_);
propsDict.add("timeStep0", timeStep0_);
propsDict.writeObject
(
IOstream::ASCII,
IOstream::currentVersion,
this->owner().db().time().writeCompression()
);
}
}
template<class CloudType>
bool Foam::InjectionModel<CloudType>::validInjection(const label parcelI)
{
@ -328,9 +258,6 @@ void Foam::InjectionModel<CloudType>::postInjectCheck
// Increment number of injections
nInjections_++;
// Write current state to properties file
writeProps();
}
@ -343,16 +270,17 @@ Foam::InjectionModel<CloudType>::InjectionModel(CloudType& owner)
SOI_(0.0),
volumeTotal_(0.0),
massTotal_(0.0),
massInjected_(0.0),
nInjections_(0),
parcelsAddedTotal_(0),
massInjected_(this->template getBaseProperty<scalar>("massInjected")),
nInjections_(this->template getBaseProperty<scalar>("nInjections")),
parcelsAddedTotal_
(
this->template getBaseProperty<scalar>("parcelsAddedTotal")
),
parcelBasis_(pbNumber),
nParticleFixed_(0.0),
time0_(0.0),
timeStep0_(0.0)
{
readProps();
}
timeStep0_(this->template getBaseProperty<scalar>("timeStep0"))
{}
template<class CloudType>
@ -363,17 +291,20 @@ Foam::InjectionModel<CloudType>::InjectionModel
const word& type
)
:
SubModelBase<CloudType>(owner, dict, type),
SubModelBase<CloudType>(owner, dict, typeName, type),
SOI_(0.0),
volumeTotal_(0.0),
massTotal_(0.0),
massInjected_(0.0),
nInjections_(0),
parcelsAddedTotal_(0),
massInjected_(this->template getBaseProperty<scalar>("massInjected")),
nInjections_(this->template getBaseProperty<scalar>("nInjections")),
parcelsAddedTotal_
(
this->template getBaseProperty<scalar>("parcelsAddedTotal")
),
parcelBasis_(pbNumber),
nParticleFixed_(0.0),
time0_(owner.db().time().value()),
timeStep0_(0.0)
timeStep0_(this->template getBaseProperty<scalar>("timeStep0"))
{
// Provide some info
// - also serves to initialise mesh dimensions - needed for parallel runs
@ -424,8 +355,6 @@ Foam::InjectionModel<CloudType>::InjectionModel
)<< "parcelBasisType must be either 'number', 'mass' or 'fixed'" << nl
<< exit(FatalError);
}
readProps();
}
@ -803,10 +732,22 @@ bool Foam::InjectionModel<CloudType>::fullyDescribed() const
template<class CloudType>
void Foam::InjectionModel<CloudType>::info(Ostream& os) const
void Foam::InjectionModel<CloudType>::info(Ostream& os)
{
os << " Total number of parcels added = " << parcelsAddedTotal_ << nl
<< " Total mass introduced = " << massInjected_ << nl;
if
(
this->owner().solution().transient()
&& this->owner().db().time().outputTime()
)
{
this->setBaseProperty("massInjected", massInjected_);
this->setBaseProperty("nInjections", nInjections_);
this->setBaseProperty("parcelsAddedTotal", parcelsAddedTotal_);
this->setBaseProperty("timeStep0", timeStep0_);
}
}

13
src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
View File

@ -85,17 +85,6 @@ public:
};
private:
// Private Member Functions
//- Read injector properties from previous run (if applicable)
void readProps();
//- Write injector properties
void writeProps();
protected:
// Protected data
@ -328,7 +317,7 @@ public:
// I-O
//- Write injection info to stream
virtual void info(Ostream& os) const;
virtual void info(Ostream& os);
};

101
src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C
View File

@ -25,83 +25,6 @@ License
#include "LocalInteraction.H"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class CloudType>
void Foam::LocalInteraction<CloudType>::readProps()
{
if (!this->owner().solution().transient())
{
return;
}
IOobject propsDictHeader
(
"localInteractionProperties",
this->owner().db().time().timeName(),
"uniform"/cloud::prefix/this->owner().name(),
this->owner().db(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
);
if (propsDictHeader.headerOk())
{
const IOdictionary propsDict(propsDictHeader);
propsDict.readIfPresent("nEscape", nEscape0_);
propsDict.readIfPresent("massEscape", massEscape0_);
propsDict.readIfPresent("nStick", nStick0_);
propsDict.readIfPresent("massStick", massStick0_);
}
}
template<class CloudType>
void Foam::LocalInteraction<CloudType>::writeProps
(
const labelList& nEscape,
const scalarList& massEscape,
const labelList& nStick,
const scalarList& massStick
) const
{
if (!this->owner().solution().transient())
{
return;
}
if (this->owner().db().time().outputTime())
{
IOdictionary propsDict
(
IOobject
(
"localInteractionProperties",
this->owner().db().time().timeName(),
"uniform"/cloud::prefix/this->owner().name(),
this->owner().db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
)
);
propsDict.add("nEscape", nEscape);
propsDict.add("massEscape", massEscape);
propsDict.add("nStick", nStick);
propsDict.add("massStick", massStick);
propsDict.writeObject
(
IOstream::ASCII,
IOstream::currentVersion,
this->owner().db().time().writeCompression()
);
}
}
// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
template<class CloudType>
@ -122,6 +45,12 @@ Foam::LocalInteraction<CloudType>::LocalInteraction
nStick_(patchData_.size(), 0),
massStick_(patchData_.size(), 0.0)
{
// intialise starting counters
this->getModelProperty("nEscape", nEscape0_);
this->getModelProperty("massEscape", massEscape0_);
this->getModelProperty("nStick", nStick0_);
this->getModelProperty("massStick", massStick0_);
// check that interactions are valid/specified
forAll(patchData_, patchI)
{
@ -141,8 +70,6 @@ Foam::LocalInteraction<CloudType>::LocalInteraction
<< nl << exit(FatalError);
}
}
readProps();
}
@ -255,7 +182,7 @@ bool Foam::LocalInteraction<CloudType>::correct
"typename CloudType::parcelType&, "
"const polyPatch&, "
"bool&, "
"scalar&, "
"const scalar, "
"const tetIndices&"
") const"
) << "Unknown interaction type "
@ -275,7 +202,7 @@ bool Foam::LocalInteraction<CloudType>::correct
template<class CloudType>
void Foam::LocalInteraction<CloudType>::info(Ostream& os) const
void Foam::LocalInteraction<CloudType>::info(Ostream& os)
{
labelList npe(nEscape_);
Pstream::listCombineGather(npe, plusEqOp<label>());
@ -304,7 +231,17 @@ void Foam::LocalInteraction<CloudType>::info(Ostream& os) const
<< ", " << mps[i] << nl;
}
writeProps(npe, mpe, nps, mps);
if
(
this->owner().solution().transient()
&& this->owner().db().time().outputTime()
)
{
this->setModelProperty("nEscape", npe);
this->setModelProperty("massEscape", mpe);
this->setModelProperty("nStick", nps);
this->setModelProperty("massStick", mps);
}
}

17
src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.H
View File

@ -84,21 +84,6 @@ class LocalInteraction
List<scalar> massStick_;
// Private Member Functions
//- Read interaction properties from file
void readProps();
//- Write interaction properties to file
void writeProps
(
const labelList& nEscape,
const scalarList& massEscape,
const labelList& nStick,
const scalarList& massStick
) const;
public:
//- Runtime type information
@ -144,7 +129,7 @@ public:
// I-O
//- Write patch interaction info to stream
virtual void info(Ostream& os) const;
virtual void info(Ostream& os);
};

4
src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.C
View File

@ -122,7 +122,7 @@ Foam::PatchInteractionModel<CloudType>::PatchInteractionModel
const word& type
)
:
SubModelBase<CloudType>(owner, dict, type),
SubModelBase<CloudType>(owner, dict, typeName, type),
UName_(this->coeffDict().lookupOrDefault("UName", word("U")))
{}
@ -333,7 +333,7 @@ void Foam::PatchInteractionModel<CloudType>::patchData
template<class CloudType>
void Foam::PatchInteractionModel<CloudType>::info(Ostream& os) const
void Foam::PatchInteractionModel<CloudType>::info(Ostream& os)
{
// do nothing
}

2
src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.H
View File

@ -180,7 +180,7 @@ public:
// I-O
//- Write patch interaction info to stream
virtual void info(Ostream& os) const;
virtual void info(Ostream& os);
};

101
src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C
View File

@ -25,83 +25,6 @@ License
#include "StandardWallInteraction.H"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class CloudType>
void Foam::StandardWallInteraction<CloudType>::readProps()
{
if (!this->owner().solution().transient())
{
return;
}
IOobject propsDictHeader
(
"standardWallInteractionProperties",
this->owner().db().time().timeName(),
"uniform"/cloud::prefix/this->owner().name(),
this->owner().db(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
);
if (propsDictHeader.headerOk())
{
const IOdictionary propsDict(propsDictHeader);
propsDict.readIfPresent("nEscape", nEscape0_);
propsDict.readIfPresent("massEscape", massEscape0_);
propsDict.readIfPresent("nStick", nStick0_);
propsDict.readIfPresent("massStick", massStick0_);
}
}
template<class CloudType>
void Foam::StandardWallInteraction<CloudType>::writeProps
(
const label nEscape,
const scalar massEscape,
const label nStick,
const scalar massStick
) const
{
if (!this->owner().solution().transient())
{
return;
}
if (this->owner().db().time().outputTime())
{
IOdictionary propsDict
(
IOobject
(
"standardWallInteractionProperties",
this->owner().db().time().timeName(),
"uniform"/cloud::prefix/this->owner().name(),
this->owner().db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
)
);
propsDict.add("nEscape", nEscape);
propsDict.add("massEscape", massEscape);
propsDict.add("nStick", nStick);
propsDict.add("massStick", massStick);
propsDict.writeObject
(
IOstream::ASCII,
IOstream::currentVersion,
this->owner().db().time().writeCompression()
);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class CloudType>
@ -118,10 +41,10 @@ Foam::StandardWallInteraction<CloudType>::StandardWallInteraction
),
e_(0.0),
mu_(0.0),
nEscape0_(0),
massEscape0_(0.0),
nStick0_(0),
massStick0_(0.0),
nEscape0_(this->template getModelProperty<label>("nEscape")),
massEscape0_(this->template getModelProperty<scalar>("massEscape")),
nStick0_(this->template getModelProperty<label>("nStick")),
massStick0_(this->template getModelProperty<scalar>("massStick")),
nEscape_(0),
massEscape_(0.0),
nStick_(0),
@ -279,7 +202,7 @@ bool Foam::StandardWallInteraction<CloudType>::correct
template<class CloudType>
void Foam::StandardWallInteraction<CloudType>::info(Ostream& os) const
void Foam::StandardWallInteraction<CloudType>::info(Ostream& os)
{
label npe = returnReduce(nEscape_, sumOp<label>()) + nEscape0_;
scalar mpe = returnReduce(massEscape_, sumOp<scalar>()) + massEscape0_;
@ -287,11 +210,21 @@ void Foam::StandardWallInteraction<CloudType>::info(Ostream& os) const
label nps = returnReduce(nStick_, sumOp<label>()) + nStick0_;
scalar mps = returnReduce(massStick_, sumOp<scalar>()) + massStick0_;
os << " Parcel fates:" << nl
os << " Parcel fate (number, mass)" << nl
<< " - escape = " << npe << ", " << mpe << nl
<< " - stick = " << nps << ", " << mps << nl;
writeProps(npe, mpe, nps, mps);
if
(
this->owner().solution().transient()
&& this->owner().db().time().outputTime()
)
{
this->setModelProperty("nEscape", npe);
this->setModelProperty("massEscape", mpe);
this->setModelProperty("nStick", nps);
this->setModelProperty("massStick", mps);
}
}

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