ENH: chtMultiRegionFoam - reverted calculation of global hEqn

This commit is contained in:
andy
2011-11-23 10:31:25 +00:00
parent 0a69856c5d
commit 971d1f1449
23 changed files with 75 additions and 1139 deletions

View File

@ -8,9 +8,7 @@ EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude
-I$(LIB_SRC)/parallel/decompose/decompose/lnInclude \
-I$(LIB_SRC)/parallel/reconstruct/reconstruct/lnInclude
EXE_LIBS = \ EXE_LIBS = \
-lbasicThermophysicalModels \ -lbasicThermophysicalModels \
@ -21,6 +19,4 @@ EXE_LIBS = \
-lcompressibleLESModels \ -lcompressibleLESModels \
-lmeshTools \ -lmeshTools \
-lfiniteVolume \ -lfiniteVolume \
-lradiationModels \ -lradiationModels
-ldecompose \
-lreconstruct

View File

@ -1,84 +0,0 @@
// Get mapped alpha (surfaceScalarField)
tmp<surfaceScalarField> tallAlpha = procToAllMapper().reconstructFvSurfaceField
(
IOobject
(
"alpha",
allMesh().time().timeName(),
allMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
procAlpha
);
// Get alpha from harmonic interpolation of vol quantities
// (Note: really only needed at internal faces originating patches
// inbetween regions)
tmp<surfaceScalarField> allHarmonicAlpha
(
harmonic(allMesh()).interpolate(allVolAlpha())
);
// Loop over all fluid and solid regions to transfer
// allHarmonicAlpha to allAlpha
surfaceScalarField& allAlpha = tallAlpha();
forAll(boundaryProcAddressing, procI)
{
forAll(boundaryProcAddressing[procI], patchI)
{
if (boundaryProcAddressing[procI][patchI] == -1)
{
// Interface patch
const labelList::subList cp =
procMeshes[procI].boundary()[patchI].patchSlice
(
faceProcAddressing[procI]
);
forAll(cp, faceI)
{
label curF = mag(cp[faceI])-1;
if (curF < allMesh().nInternalFaces())
{
allAlpha[curF] = allHarmonicAlpha()[curF];
}
}
}
}
}
tmp<surfaceScalarField> allPhi
(
procToAllMapper().reconstructFvSurfaceField
(
IOobject
(
"phi",
allMesh().time().timeName(),
allMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
procPhi
)
);
// So we have nNonOrthCorr
//#include "readSolidMultiRegionPIMPLEControls.H"
//for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix hEqn
(
fvm::ddt(allRho(), allh())
+ fvm::div(allPhi(), allh())
- fvm::laplacian(allAlpha, allh())
==
allSource()
);
hEqn.relax();
hEqn.solve(allMesh().solver(allh().select(finalIter)));
}

View File

@ -39,11 +39,6 @@ Description
#include "solidRegionDiffNo.H" #include "solidRegionDiffNo.H"
#include "basicSolidThermo.H" #include "basicSolidThermo.H"
#include "radiationModel.H" #include "radiationModel.H"
#include "fvFieldReconstructor.H"
#include "mixedFvPatchFields.H"
#include "fvFieldDecomposer.H"
#include "harmonic.H"
#include "rmap.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -54,40 +49,12 @@ int main(int argc, char *argv[])
regionProperties rp(runTime); regionProperties rp(runTime);
const label nAllRegions =
rp.fluidRegionNames().size()
+ rp.solidRegionNames().size();
PtrList<labelIOList> cellProcAddressing(nAllRegions);
PtrList<labelIOList> faceProcAddressing(nAllRegions);
PtrList<labelIOList> boundaryProcAddressing(nAllRegions);
PtrList<fvMesh> procMeshes(nAllRegions);
// Load meshes, fluid first
labelList fluidToProc(identity(rp.fluidRegionNames().size()));
labelList solidToProc(rp.solidRegionNames().size());
forAll(solidToProc, i)
{
solidToProc[i] = fluidToProc.size()+i;
}
// Get the coupled solution flag
#include "readPIMPLEControls.H"
if (temperatureCoupled)
{
Info<< "Solving single enthalpy for all equations" << nl << endl;
}
#include "createFluidMeshes.H" #include "createFluidMeshes.H"
#include "createSolidMeshes.H" #include "createSolidMeshes.H"
#include "createFluidFields.H" #include "createFluidFields.H"
#include "createSolidFields.H" #include "createSolidFields.H"
// Temperature solved on single mesh
#include "createAllMesh.H"
#include "createAllFields.H"
#include "initContinuityErrs.H" #include "initContinuityErrs.H"
#include "readTimeControls.H" #include "readTimeControls.H"
@ -114,10 +81,9 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl; Info<< "Time = " << runTime.timeName() << nl << endl;
if (nOuterCorr != 1) if (nOuterCorr != 1)
{ {
forAll(fluidToProc, i) forAll(fluidRegions, i)
{ {
#include "setRegionFluidFields.H" #include "setRegionFluidFields.H"
#include "storeOldFluidFields.H" #include "storeOldFluidFields.H"
@ -130,133 +96,22 @@ int main(int argc, char *argv[])
{ {
bool finalIter = oCorr == nOuterCorr-1; bool finalIter = oCorr == nOuterCorr-1;
if (finalIter) forAll(fluidRegions, i)
{ {
forAll(procMeshes, procI) Info<< "\nSolving for fluid region "
{ << fluidRegions[i].name() << endl;
procMeshes[procI].data::add("finalIteration", true);
}
}
PtrList<surfaceScalarField> procPhi(nAllRegions);
PtrList<surfaceScalarField> procAlpha(nAllRegions);
// Solve (uncoupled) or set up (coupled) the temperature equation
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
forAll(solidToProc, i)
{
label procI = solidToProc[i];
Info<< "\nSolving temperature for solid region "
<< procMeshes[procI].name() << endl;
#include "setRegionSolidFields.H"
#include "readSolidMultiRegionPIMPLEControls.H"
if (temperatureCoupled)
{
// Map my properties to overall h equation
#include "rmapSolid.H"
}
else
{
#include "solveSolid.H"
}
}
forAll(fluidToProc, i)
{
label procI = fluidToProc[i];
Info<< "\nSolving temperature for fluid region "
<< procMeshes[procI].name() << endl;
#include "setRegionFluidFields.H" #include "setRegionFluidFields.H"
#include "readFluidMultiRegionPIMPLEControls.H" #include "readFluidMultiRegionPIMPLEControls.H"
#include "solveFluid.H"
if (oCorr == 0)
{
#include "rhoEqn.H"
}
if (temperatureCoupled)
{
// Map my properties to overall h equation
#include "rmapFluid.H"
}
else
{
#include "hEqn.H"
}
} }
forAll(solidRegions, i)
// Solve combined h equation
// ~~~~~~~~~~~~~~~~~~~~~~~~~
if (temperatureCoupled)
{ {
Info<< "\nSolving single enthalpy for all regions" Info<< "\nSolving for solid region "
<< endl; << solidRegions[i].name() << endl;
// Solve combined h
#include "allhEqn.H"
forAll(solidToProc, i)
{
label procI = solidToProc[i];
#include "setRegionSolidFields.H"
#include "readSolidMultiRegionPIMPLEControls.H"
#include "mapSolid.H"
}
forAll(fluidToProc, i)
{
label procI = fluidToProc[i];
#include "setRegionFluidFields.H"
#include "readFluidMultiRegionPIMPLEControls.H"
#include "mapFluid.H"
}
}
// Update thermos
// ~~~~~~~~~~~~~~
forAll(fluidToProc, i)
{
Info<< "\nUpdating thermo for fluid region "
<< procMeshes[fluidToProc[i]].name() << endl;
#include "setRegionFluidFields.H"
#include "readFluidMultiRegionPIMPLEControls.H"
thermo.correct();
rad.correct();
#include "solvePressureVelocityFluid.H"
}
forAll(solidToProc, i)
{
Info<< "\nUpdating thermo for solid region "
<< procMeshes[solidToProc[i]].name() << endl;
#include "setRegionSolidFields.H" #include "setRegionSolidFields.H"
#include "readSolidMultiRegionPIMPLEControls.H" #include "readSolidMultiRegionPIMPLEControls.H"
#include "solveSolid.H"
thermo.correct();
}
if (finalIter)
{
forAll(procMeshes, procI)
{
procMeshes[procI].data::remove("finalIteration");
}
} }
} }

View File

@ -1,64 +0,0 @@
autoPtr<volScalarField> allRho;
autoPtr<volScalarField> allh;
autoPtr<volScalarField> allVolAlpha;
autoPtr<fvScalarMatrix> allSource;
if (temperatureCoupled)
{
allRho.reset
(
new volScalarField
(
IOobject
(
"rho",
runTime.timeName(),
allMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
allMesh(),
dimensionedScalar("rho", dimDensity, 0.0)
)
);
allh.reset
(
new volScalarField
(
IOobject
(
"h",
runTime.timeName(),
allMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
allMesh(),
dimensionedScalar("h", dimEnergy/dimMass, 0.0),
mixedFvPatchScalarField::typeName
)
);
allVolAlpha.reset
(
new volScalarField
(
IOobject
(
"volAlpha",
runTime.timeName(),
allMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
allMesh(),
dimensionedScalar("volAlpha", dimMass/dimLength/dimTime, 0.0)
)
);
allSource.reset
(
new fvMatrix<scalar>(allh(), allh().dimensions()*dimMass/dimTime)
);
}

View File

@ -1,63 +0,0 @@
//
// createAllMesh.H
// ~~~~~~~~~~~~~~~
autoPtr<fvMesh> allMesh;
autoPtr<fvFieldReconstructor> procToAllMapper;
PtrList<fvFieldDecomposer> allToProcMappers;
if (temperatureCoupled)
{
Foam::Info
<< "Create mesh for time = "
<< runTime.timeName() << Foam::nl << Foam::endl;
allMesh.reset
(
new Foam::fvMesh
(
Foam::IOobject
(
Foam::fvMesh::defaultRegion,
runTime.timeName(),
runTime,
Foam::IOobject::MUST_READ
)
)
);
procToAllMapper.reset
(
new fvFieldReconstructor
(
allMesh(),
procMeshes,
faceProcAddressing,
cellProcAddressing,
boundaryProcAddressing
)
);
allToProcMappers.setSize
(
rp.fluidRegionNames().size()
+ rp.solidRegionNames().size()
);
forAll(allToProcMappers, i)
{
allToProcMappers.set
(
i,
new fvFieldDecomposer
(
allMesh(),
procMeshes[i],
faceProcAddressing[i],
cellProcAddressing[i],
boundaryProcAddressing[i]
)
);
}
}

View File

@ -1,12 +1,12 @@
scalar CoNum = -GREAT; scalar CoNum = -GREAT;
forAll(fluidToProc, regionI) forAll(fluidRegions, regionI)
{ {
CoNum = max CoNum = max
( (
compressibleCourantNo compressibleCourantNo
( (
procMeshes[fluidToProc[regionI]], fluidRegions[regionI],
runTime, runTime,
rhoFluid[regionI], rhoFluid[regionI],
phiFluid[regionI] phiFluid[regionI]

View File

@ -1,36 +1,30 @@
// Initialise fluid field pointer lists // Initialise fluid field pointer lists
PtrList<basicRhoThermo> thermoFluid(rp.fluidRegionNames().size()); PtrList<basicRhoThermo> thermoFluid(fluidRegions.size());
PtrList<volScalarField> rhoFluid(rp.fluidRegionNames().size()); PtrList<volScalarField> rhoFluid(fluidRegions.size());
PtrList<volScalarField> KFluid(rp.fluidRegionNames().size()); PtrList<volScalarField> KFluid(fluidRegions.size());
PtrList<volVectorField> UFluid(rp.fluidRegionNames().size()); PtrList<volVectorField> UFluid(fluidRegions.size());
PtrList<surfaceScalarField> phiFluid(rp.fluidRegionNames().size()); PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
PtrList<uniformDimensionedVectorField> gFluid(rp.fluidRegionNames().size()); PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
PtrList<compressible::turbulenceModel> turbulence PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
( PtrList<volScalarField> p_rghFluid(fluidRegions.size());
rp.fluidRegionNames().size() PtrList<volScalarField> ghFluid(fluidRegions.size());
); PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
PtrList<volScalarField> p_rghFluid(rp.fluidRegionNames().size()); PtrList<radiation::radiationModel> radiation(fluidRegions.size());
PtrList<volScalarField> ghFluid(rp.fluidRegionNames().size()); PtrList<volScalarField> DpDtFluid(fluidRegions.size());
PtrList<surfaceScalarField> ghfFluid(rp.fluidRegionNames().size());
PtrList<radiation::radiationModel> radiation(rp.fluidRegionNames().size());
PtrList<volScalarField> DpDtFluid(rp.fluidRegionNames().size());
List<scalar> initialMassFluid(rp.fluidRegionNames().size()); List<scalar> initialMassFluid(fluidRegions.size());
// Populate fluid field pointer lists // Populate fluid field pointer lists
forAll(rp.fluidRegionNames(), i) forAll(fluidRegions, i)
{ {
Info<< "*** Reading fluid mesh thermophysical properties for region " Info<< "*** Reading fluid mesh thermophysical properties for region "
<< rp.fluidRegionNames()[i] << nl << endl; << fluidRegions[i].name() << nl << endl;
label procI = fluidToProc[i];
Info<< " Adding to thermoFluid\n" << endl; Info<< " Adding to thermoFluid\n" << endl;
thermoFluid.set thermoFluid.set
( (
i, i,
basicRhoThermo::New(procMeshes[procI]).ptr() basicRhoThermo::New(fluidRegions[i]).ptr()
); );
Info<< " Adding to rhoFluid\n" << endl; Info<< " Adding to rhoFluid\n" << endl;
@ -43,7 +37,7 @@
( (
"rho", "rho",
runTime.timeName(), runTime.timeName(),
procMeshes[procI], fluidRegions[i],
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
@ -61,7 +55,7 @@
( (
"K", "K",
runTime.timeName(), runTime.timeName(),
procMeshes[procI], fluidRegions[i],
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
@ -79,11 +73,11 @@
( (
"U", "U",
runTime.timeName(), runTime.timeName(),
procMeshes[procI], fluidRegions[i],
IOobject::MUST_READ, IOobject::MUST_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
procMeshes[procI] fluidRegions[i]
) )
); );
@ -97,12 +91,12 @@
( (
"phi", "phi",
runTime.timeName(), runTime.timeName(),
procMeshes[procI], fluidRegions[i],
IOobject::READ_IF_PRESENT, IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
linearInterpolate(rhoFluid[i]*UFluid[i]) linearInterpolate(rhoFluid[i]*UFluid[i])
& procMeshes[procI].Sf() & fluidRegions[i].Sf()
) )
); );
@ -116,7 +110,7 @@
( (
"g", "g",
runTime.constant(), runTime.constant(),
procMeshes[procI], fluidRegions[i],
IOobject::MUST_READ, IOobject::MUST_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
) )
@ -143,14 +137,14 @@
ghFluid.set ghFluid.set
( (
i, i,
new volScalarField("gh", gFluid[i] & procMeshes[procI].C()) new volScalarField("gh", gFluid[i] & fluidRegions[i].C())
); );
Info<< " Adding to ghfFluid\n" << endl; Info<< " Adding to ghfFluid\n" << endl;
ghfFluid.set ghfFluid.set
( (
i, i,
new surfaceScalarField("ghf", gFluid[i] & procMeshes[procI].Cf()) new surfaceScalarField("ghf", gFluid[i] & fluidRegions[i].Cf())
); );
p_rghFluid.set p_rghFluid.set
@ -162,11 +156,11 @@
( (
"p_rgh", "p_rgh",
runTime.timeName(), runTime.timeName(),
procMeshes[procI], fluidRegions[i],
IOobject::MUST_READ, IOobject::MUST_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
procMeshes[procI] fluidRegions[i]
) )
); );

View File

@ -1,13 +1,13 @@
PtrList<fvMesh> fluidRegions(rp.fluidRegionNames().size());
forAll(rp.fluidRegionNames(), i) forAll(rp.fluidRegionNames(), i)
{ {
Info<< "Create fluid mesh for region " << rp.fluidRegionNames()[i] Info<< "Create fluid mesh for region " << rp.fluidRegionNames()[i]
<< " for time = " << runTime.timeName() << nl << endl; << " for time = " << runTime.timeName() << nl << endl;
label procI = fluidToProc[i]; fluidRegions.set
procMeshes.set
( (
procI, i,
new fvMesh new fvMesh
( (
IOobject IOobject
@ -19,58 +19,4 @@
) )
) )
); );
if (temperatureCoupled)
{
cellProcAddressing.set
(
procI,
new labelIOList
(
IOobject
(
"cellRegionAddressing",
procMeshes[procI].facesInstance(),
procMeshes[procI].meshSubDir,
procMeshes[procI],
IOobject::MUST_READ,
IOobject::NO_WRITE
)
)
);
faceProcAddressing.set
(
procI,
new labelIOList
(
IOobject
(
"faceRegionAddressing",
procMeshes[procI].facesInstance(),
procMeshes[procI].meshSubDir,
procMeshes[procI],
IOobject::MUST_READ,
IOobject::NO_WRITE
)
)
);
boundaryProcAddressing.set
(
procI,
new labelIOList
(
IOobject
(
"boundaryRegionAddressing",
procMeshes[procI].facesInstance(),
procMeshes[procI].meshSubDir,
procMeshes[procI],
IOobject::MUST_READ,
IOobject::NO_WRITE
)
)
);
}
} }

View File

@ -1 +1 @@
List<scalar> cumulativeContErr(fluidToProc.size(), 0.0); List<scalar> cumulativeContErr(fluidRegions.size(), 0.0);

View File

@ -1,3 +0,0 @@
h = allToProcMappers[procI].decomposeField(allh(), true);
h.oldTime().timeIndex() = allh().oldTime().timeIndex();
h.correctBoundaryConditions();

View File

@ -1,56 +0,0 @@
// Note:Map rho and rho.oldTime() since fluid rho assigned to at
// end of iteration.
rmap
(
allRho(),
rho,
faceProcAddressing[procI],
cellProcAddressing[procI],
boundaryProcAddressing[procI]
);
rmap
(
allRho().oldTime(),
rho.oldTime(),
faceProcAddressing[procI],
cellProcAddressing[procI],
boundaryProcAddressing[procI]
);
// Necessary? Probably only for boundary values since bcs on
// h are not the same as those on allh
rmap
(
allh(),
h,
faceProcAddressing[procI],
cellProcAddressing[procI],
boundaryProcAddressing[procI]
);
procAlpha.set(procI, fvc::interpolate(turb.alphaEff()));
rmap
(
allVolAlpha(),
turb.alphaEff()(),
faceProcAddressing[procI],
cellProcAddressing[procI],
boundaryProcAddressing[procI]
);
rmap
(
allSource(),
(DpDt + rad.Sh(thermo))(),
faceProcAddressing[procI],
cellProcAddressing[procI],
boundaryProcAddressing[procI]
);
procPhi.set
(
procI,
new surfaceScalarField(phi)
);

View File

@ -1,4 +1,4 @@
fvMesh& mesh = procMeshes[fluidToProc[i]]; fvMesh& mesh = fluidRegions[i];
basicRhoThermo& thermo = thermoFluid[i]; basicRhoThermo& thermo = thermoFluid[i];
volScalarField& rho = rhoFluid[i]; volScalarField& rho = rhoFluid[i];

View File

@ -1,11 +0,0 @@
#include "UEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}
turb.correct();
rho = thermo.rho();

View File

@ -6,6 +6,3 @@
const int nOuterCorr = const int nOuterCorr =
pimple.lookupOrDefault<int>("nOuterCorrectors", 1); pimple.lookupOrDefault<int>("nOuterCorrectors", 1);
const Switch temperatureCoupled =
pimple.lookupOrDefault<Switch>("temperatureCoupled", false);

View File

@ -1,96 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
Global
rmap
Description
map field on subset of mesh onto overall field. Rewrites boundary
conditions to be mixed.
The source fields can have different patch types for the same destination
patch. To work around this it attempts to convert all patch fields into
mixed type since this can accomodate anything from fixedValue to
fixedGradient.
SourceFiles
rmap.C
\*---------------------------------------------------------------------------*/
#ifndef rmap_H
#define rmap_H
#include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Map patchField
template<class Type>
static void rmap
(
fvPatchField<Type>& destBC,
const labelList& reverseAddressing,
const fvPatchField<Type>& sourceBC
);
//- Map volField
template<class Type>
static void rmap
(
GeometricField<Type, fvPatchField, volMesh>& dest,
const GeometricField<Type, fvPatchField, volMesh>& source,
const labelList& faceProcAddressing,
const labelList& cellProcAddressing,
const labelList& boundaryProcAddressing
);
//- Map fvMatrix
template<class Type>
static void rmap
(
fvMatrix<Type>& dest,
const fvMatrix<Type>& source,
const labelList& faceProcAddressing,
const labelList& cellProcAddressing,
const labelList& boundaryProcAddressing
);
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "rmapTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,309 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "rmap.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
void Foam::rmap
(
fvPatchField<Type>& destBC,
const labelList& reverseAddressing,
const fvPatchField<Type>& sourceBC
)
{
// Assign value
destBC.Field<Type>::rmap(sourceBC, reverseAddressing);
// Assign other properties
if (isA<mixedFvPatchField<Type> >(destBC))
{
mixedFvPatchField<Type>& mp =
refCast<mixedFvPatchField<Type> >(destBC);
if (isA<mixedFvPatchField<Type> >(sourceBC))
{
const mixedFvPatchField<Type>& Tp =
refCast<const mixedFvPatchField<Type> >(sourceBC);
mp.refValue().rmap(Tp.refValue(), reverseAddressing);
mp.refGrad().rmap(Tp.refGrad(), reverseAddressing);
mp.valueFraction().rmap(Tp.valueFraction(), reverseAddressing);
}
else if (isA<fixedGradientFvPatchField<Type> >(sourceBC))
{
const fixedGradientFvPatchField<Type>& Tp =
refCast<const fixedGradientFvPatchField<Type> >
(
sourceBC
);
// Make pure fixedGradient
mp.refValue().rmap(Tp, reverseAddressing); // unused
mp.refGrad().rmap(Tp.gradient(), reverseAddressing);
mp.valueFraction().rmap
(
Field<Type>(reverseAddressing.size(), 0.0),
reverseAddressing
);
}
else if (isA<zeroGradientFvPatchField<Type> >(sourceBC))
{
// Make pure fixedGradient with gradient = 0
mp.refValue().rmap(sourceBC, reverseAddressing); // unused
mp.refGrad().rmap
(
Field<Type>(reverseAddressing.size(), 0.0),
reverseAddressing
);
mp.valueFraction().rmap
(
Field<Type>(reverseAddressing.size(), 0.0),
reverseAddressing
);
}
else if (isA<fixedValueFvPatchField<Type> >(sourceBC))
{
// Make pure fixedValue
mp.refValue().rmap(sourceBC, reverseAddressing);
mp.refGrad().rmap
(
Field<Type>(reverseAddressing.size(), 0.0),
reverseAddressing
); // unused
mp.valueFraction().rmap
(
Field<Type>(reverseAddressing.size(), 1.0),
reverseAddressing
);
}
else if (isA<calculatedFvPatchField<Type> >(sourceBC))
{
// Make pure fixedValue
mp.refValue().rmap(sourceBC, reverseAddressing);
mp.refGrad().rmap
(
Field<Type>(reverseAddressing.size(), 0.0),
reverseAddressing
); // unused
mp.valueFraction().rmap
(
Field<Type>(reverseAddressing.size(), 1.0),
reverseAddressing
);
}
else
{
FatalErrorIn("rmap(..)")
<< "Don't know how to map source bc "
<< sourceBC.type()
<< " into a mixed boundary condition at "
<< destBC.patch().name()
<< exit(FatalError);
}
}
else if (isA<fixedGradientFvPatchField<Type> >(destBC))
{
fixedGradientFvPatchField<Type>& mp =
refCast<fixedGradientFvPatchField<Type> >(destBC);
if (isA<fixedGradientFvPatchField<Type> >(sourceBC))
{
const fixedGradientFvPatchField<Type>& Tp =
refCast<const fixedGradientFvPatchField<Type> >
(
sourceBC
);
mp.gradient().rmap(Tp.gradient(), reverseAddressing);
}
else if (isA<mixedFvPatchField<Type> >(sourceBC))
{
const mixedFvPatchField<Type>& Tp =
refCast<const mixedFvPatchField<Type> >(sourceBC);
mp.gradient().rmap(Tp.snGrad(), reverseAddressing);
}
else if (isA<zeroGradientFvPatchField<Type> >(sourceBC))
{
mp.gradient().rmap
(
Field<Type>(reverseAddressing.size(), 0.0),
reverseAddressing
);
}
else
{
FatalErrorIn("rmap(..)")
<< "Don't know how to map source bc "
<< sourceBC.type()
<< " into a fixedGradient boundary condition at "
<< destBC.patch().name()
<< exit(FatalError);
}
}
}
template<class Type>
void Foam::rmap
(
GeometricField<Type, fvPatchField, volMesh>& dest,
const GeometricField<Type, fvPatchField, volMesh>& source,
const labelList& faceProcAddressing,
const labelList& cellProcAddressing,
const labelList& boundaryProcAddressing
)
{
if (dest.dimensions() != source.dimensions())
{
FatalErrorIn("rmap(..)")
<< "Different dimensions for = for fields " << dest.name()
<< " and " << source.name() << endl
<< " dimensions : " << dest.dimensions()
<< " = " << source.dimensions() << endl
<< exit(FatalError);
}
// Copy internal field
dest.internalField().rmap(source.internalField(), cellProcAddressing);
// Copy boundary properties as mixed
forAll(source.boundaryField(), patchI)
{
label curBPatch = boundaryProcAddressing[patchI];
if (curBPatch == -1)
{
// Unknown patch. Do not change any values.
}
else
{
// Get addressing slice for this patch
const labelList::subList cp =
source.mesh().boundary()[patchI].patchSlice
(
faceProcAddressing
);
const label curPatchStart =
dest.mesh().boundaryMesh()[curBPatch].start();
labelList reverseAddressing(cp.size());
forAll(cp, faceI)
{
// Subtract one to take into account offsets for
// face direction.
if (cp[faceI] <= 0)
{
FatalErrorIn("rmap(..)")
<< "Problem:"
<< " patch:" << source.mesh().boundary()[patchI].name()
<< " field:" << source.name()
<< " local face:" << faceI
<< " mapped to:" << cp[faceI] << exit(FatalError);
}
reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
}
// Map curBPatch from source patch. Is like rmap but also
// copies non-value properties from alike patchFields.
rmap
(
dest.boundaryField()[curBPatch],
reverseAddressing,
source.boundaryField()[patchI]
);
}
}
// Copy timeIndex
dest.timeIndex() = source.timeIndex();
}
template<class Type>
void Foam::rmap
(
fvMatrix<Type>& dest,
const fvMatrix<Type>& source,
const labelList& faceProcAddressing,
const labelList& cellProcAddressing,
const labelList& boundaryProcAddressing
)
{
dest.source().rmap(source.source(), cellProcAddressing);
FieldField<Field, Type>& sourceInternal =
const_cast<fvMatrix<Type>&>(source).internalCoeffs();
FieldField<Field, Type>& sourceBoundary =
const_cast<fvMatrix<Type>&>(source).boundaryCoeffs();
forAll(sourceInternal, patchI)
{
label curBPatch = boundaryProcAddressing[patchI];
if (curBPatch == -1)
{
// Unknown patch. Do not change any values.
}
else
{
// Get addressing slice for this patch
const fvMesh& sourceMesh = source.psi().mesh();
const labelList::subList cp =
sourceMesh.boundary()[patchI].patchSlice
(
faceProcAddressing
);
const label curPatchStart =
dest.psi().mesh().boundaryMesh()[curBPatch].start();
labelList reverseAddressing(cp.size());
forAll(cp, faceI)
{
// Subtract one to take into account offsets for
// face direction.
reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
}
dest.internalCoeffs()[curBPatch].rmap
(
sourceInternal[patchI],
reverseAddressing
);
dest.boundaryCoeffs()[curBPatch].rmap
(
sourceBoundary[patchI],
reverseAddressing
);
}
}
}
// ************************************************************************* //

View File

@ -1,36 +1,12 @@
// Initialise solid field pointer lists // Initialise solid field pointer lists
PtrList<basicSolidThermo> thermoSolid(rp.solidRegionNames().size()); PtrList<basicSolidThermo> thermos(solidRegions.size());
PtrList<volScalarField> hSolid(rp.solidRegionNames().size());
// Populate solid field pointer lists // Populate solid field pointer lists
forAll(rp.solidRegionNames(), i) forAll(solidRegions, i)
{ {
Info<< "*** Reading solid mesh thermophysical properties for region " Info<< "*** Reading solid mesh thermophysical properties for region "
<< rp.solidRegionNames()[i] << nl << endl; << solidRegions[i].name() << nl << endl;
label procI = solidToProc[i]; Info<< " Adding to thermos\n" << endl;
thermos.set(i, basicSolidThermo::New(solidRegions[i]));
Info<< " Adding to thermoSolid\n" << endl;
thermoSolid.set(i, basicSolidThermo::New(procMeshes[procI]));
if (temperatureCoupled)
{
Info<< " Adding to hSolid\n" << endl;
hSolid.set
(
i,
new volScalarField
(
IOobject
(
"h",
runTime.timeName(),
procMeshes[procI],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
procMeshes[procI]
)
);
}
} }

View File

@ -1,13 +1,13 @@
PtrList<fvMesh> solidRegions(rp.solidRegionNames().size());
forAll(rp.solidRegionNames(), i) forAll(rp.solidRegionNames(), i)
{ {
Info<< "Create solid mesh for region " << rp.solidRegionNames()[i] Info<< "Create solid mesh for region " << rp.solidRegionNames()[i]
<< " for time = " << runTime.timeName() << nl << endl; << " for time = " << runTime.timeName() << nl << endl;
label procI = solidToProc[i]; solidRegions.set
procMeshes.set
( (
procI, i,
new fvMesh new fvMesh
( (
IOobject IOobject
@ -20,57 +20,8 @@
) )
); );
if (temperatureCoupled) // Force calculation of geometric properties to prevent it being done
{ // later in e.g. some boundary evaluation
cellProcAddressing.set //(void)solidRegions[i].weights();
( //(void)solidRegions[i].deltaCoeffs();
procI,
new labelIOList
(
IOobject
(
"cellRegionAddressing",
procMeshes[procI].facesInstance(),
procMeshes[procI].meshSubDir,
procMeshes[procI],
IOobject::MUST_READ,
IOobject::NO_WRITE
)
)
);
faceProcAddressing.set
(
procI,
new labelIOList
(
IOobject
(
"faceRegionAddressing",
procMeshes[procI].facesInstance(),
procMeshes[procI].meshSubDir,
procMeshes[procI],
IOobject::MUST_READ,
IOobject::NO_WRITE
)
)
);
boundaryProcAddressing.set
(
procI,
new labelIOList
(
IOobject
(
"boundaryRegionAddressing",
procMeshes[procI].facesInstance(),
procMeshes[procI].meshSubDir,
procMeshes[procI],
IOobject::MUST_READ,
IOobject::NO_WRITE
)
)
);
}
} }

View File

@ -1,15 +0,0 @@
{
volScalarField& h = hSolid[i];
h = allToProcMappers[procI].decomposeField(allh(), true);
h.oldTime().timeIndex() = allh().oldTime().timeIndex();
T += (h-h.oldTime())/cp;
// Correct T boundary conditions and update h boundary
// conditions accordingly.
volScalarField::GeometricBoundaryField Told = T.boundaryField();
T.correctBoundaryConditions();
h.boundaryField() +=
cp.boundaryField()
* (T.boundaryField()-Told);
}

View File

@ -1,90 +0,0 @@
{
volScalarField& h = hSolid[i];
procPhi.setSize(nAllRegions);
procAlpha.setSize(nAllRegions);
rmap
(
allRho(),
rho,
faceProcAddressing[procI],
cellProcAddressing[procI],
boundaryProcAddressing[procI]
);
// Necessary? Probably only for boundary values since bcs on
// h are not the same as those on allh
rmap
(
allh(),
h,
faceProcAddressing[procI],
cellProcAddressing[procI],
boundaryProcAddressing[procI]
);
tmp<volScalarField> Kcp(K/cp);
rmap
(
allVolAlpha(),
Kcp(),
faceProcAddressing[procI],
cellProcAddressing[procI],
boundaryProcAddressing[procI]
);
procAlpha.set(procI, fvc::interpolate(Kcp));
// allSource is initialised to zero already
//rmap
//(
// allSource(),
// volScalarField
// (
// IOobject
// (
// "procSource",
// runTime.timeName(),
// mesh,
// IOobject::NO_READ,
// IOobject::AUTO_WRITE
// ),
// mesh,
// dimensionedScalar
// (
// "procSource",
// allh().dimensions()*dimDensity/dimTime,
// 0.0
// )
// ),
// faceProcAddressing[procI],
// cellProcAddressing[procI],
// boundaryProcAddressing[procI]
//);
procPhi.set
(
procI,
new surfaceScalarField
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar
(
"phi",
dimDensity*dimVelocity*dimArea,
0.0
)
)
);
}

View File

@ -1,5 +1,5 @@
fvMesh& mesh = procMeshes[solidToProc[i]]; fvMesh& mesh = solidRegions[i];
basicSolidThermo& thermo = thermoSolid[i]; basicSolidThermo& thermo = thermos[i];
tmp<volScalarField> trho = thermo.rho(); tmp<volScalarField> trho = thermo.rho();
const volScalarField& rho = trho(); const volScalarField& rho = trho();

View File

@ -1,6 +1,6 @@
scalar DiNum = -GREAT; scalar DiNum = -GREAT;
forAll(solidToProc, i) forAll(solidRegions, i)
{ {
# include "setRegionSolidFields.H" # include "setRegionSolidFields.H"
@ -8,7 +8,7 @@
( (
solidRegionDiffNo solidRegionDiffNo
( (
procMeshes[solidToProc[i]], solidRegions[i],
runTime, runTime,
rho*cp, rho*cp,
K K

View File

@ -1,3 +1,8 @@
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
{ {
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{ {
@ -12,3 +17,10 @@
Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl; Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl;
} }
thermo.correct();
if (finalIter)
{
mesh.data::remove("finalIteration");
}