Skip to content
GitLab
Commit ba308083 authored by mattijs's avatar mattijs
Browse files

ENH: chtMultiRegionFoam: solve single h equation

parent 4b5049b0
Hide whitespace changes
Inline Side-by-side
applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
View file @ ba308083
......@@ -8,7 +8,9 @@ EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/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 = \
-lbasicThermophysicalModels \
......@@ -19,4 +21,6 @@ EXE_LIBS = \
-lcompressibleLESModels \
-lmeshTools \
-lfiniteVolume \
-lradiationModels
-lradiationModels \
-ldecompose \
-lreconstruct
applications/solvers/heatTransfer/chtMultiRegionFoam/allhEqn.H 0 → 100644
View file @ ba308083
// 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)));
}
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
View file @ ba308083
......@@ -39,6 +39,11 @@ Description
#include "solidRegionDiffNo.H"
#include "basicSolidThermo.H"
#include "radiationModel.H"
#include "fvFieldReconstructor.H"
#include "mixedFvPatchFields.H"
#include "fvFieldDecomposer.H"
#include "harmonic.H"
#include "rmap.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -49,12 +54,40 @@ int main(int argc, char *argv[])
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 "createSolidMeshes.H"
#include "createFluidFields.H"
#include "createSolidFields.H"
// Temperature solved on single mesh
#include "createAllMesh.H"
#include "createAllFields.H"
#include "initContinuityErrs.H"
#include "readTimeControls.H"
......@@ -81,9 +114,10 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
if (nOuterCorr != 1)
{
forAll(fluidRegions, i)
forAll(fluidToProc, i)
{
#include "setRegionFluidFields.H"
#include "storeOldFluidFields.H"
......@@ -96,22 +130,133 @@ int main(int argc, char *argv[])
{
bool finalIter = oCorr == nOuterCorr-1;
forAll(fluidRegions, i)
if (finalIter)
{
forAll(procMeshes, procI)
{
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 "readFluidMultiRegionPIMPLEControls.H"
if (oCorr == 0)
{
#include "rhoEqn.H"
}
if (temperatureCoupled)
{
// Map my properties to overall h equation
#include "rmapFluid.H"
}
else
{
#include "hEqn.H"
}
}
// Solve combined h equation
// ~~~~~~~~~~~~~~~~~~~~~~~~~
if (temperatureCoupled)
{
Info<< "\nSolving for fluid region "
<< fluidRegions[i].name() << endl;
Info<< "\nSolving single enthalpy for all regions"
<< 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"
#include "solveFluid.H"
thermo.correct();
rad.correct();
#include "solvePressureVelocityFluid.H"
}
forAll(solidRegions, i)
forAll(solidToProc, i)
{
Info<< "\nSolving for solid region "
<< solidRegions[i].name() << endl;
Info<< "\nUpdating thermo for solid region "
<< procMeshes[solidToProc[i]].name() << endl;
#include "setRegionSolidFields.H"
#include "readSolidMultiRegionPIMPLEControls.H"
#include "solveSolid.H"
thermo.correct();
}
if (finalIter)
{
forAll(procMeshes, procI)
{
procMeshes[procI].data::remove("finalIteration");
}
}
}
......
applications/solvers/heatTransfer/chtMultiRegionFoam/createAllFields.H 0 → 100644
View file @ ba308083
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)
);
}
applications/solvers/heatTransfer/chtMultiRegionFoam/createAllMesh.H 0 → 100644
View file @ ba308083
//
// 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]
)
);
}
}
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
View file @ ba308083
scalar CoNum = -GREAT;
forAll(fluidRegions, regionI)
forAll(fluidToProc, regionI)
{
CoNum = max
(
compressibleCourantNo
(
fluidRegions[regionI],
procMeshes[fluidToProc[regionI]],
runTime,
rhoFluid[regionI],
phiFluid[regionI]
......
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
View file @ ba308083
// Initialise fluid field pointer lists
PtrList<basicRhoThermo> thermoFluid(fluidRegions.size());
PtrList<volScalarField> rhoFluid(fluidRegions.size());
PtrList<volScalarField> KFluid(fluidRegions.size());
PtrList<volVectorField> UFluid(fluidRegions.size());
PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
PtrList<volScalarField> p_rghFluid(fluidRegions.size());
PtrList<volScalarField> ghFluid(fluidRegions.size());
PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
PtrList<radiation::radiationModel> radiation(fluidRegions.size());
PtrList<volScalarField> DpDtFluid(fluidRegions.size());
List<scalar> initialMassFluid(fluidRegions.size());
PtrList<basicRhoThermo> thermoFluid(rp.fluidRegionNames().size());
PtrList<volScalarField> rhoFluid(rp.fluidRegionNames().size());
PtrList<volScalarField> KFluid(rp.fluidRegionNames().size());
PtrList<volVectorField> UFluid(rp.fluidRegionNames().size());
PtrList<surfaceScalarField> phiFluid(rp.fluidRegionNames().size());
PtrList<uniformDimensionedVectorField> gFluid(rp.fluidRegionNames().size());
PtrList<compressible::turbulenceModel> turbulence
(
rp.fluidRegionNames().size()
);
PtrList<volScalarField> p_rghFluid(rp.fluidRegionNames().size());
PtrList<volScalarField> ghFluid(rp.fluidRegionNames().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());
// Populate fluid field pointer lists
forAll(fluidRegions, i)
forAll(rp.fluidRegionNames(), i)
{
Info<< "*** Reading fluid mesh thermophysical properties for region "
<< fluidRegions[i].name() << nl << endl;
<< rp.fluidRegionNames()[i] << nl << endl;
label procI = fluidToProc[i];
Info<< " Adding to thermoFluid\n" << endl;
thermoFluid.set
(
i,
basicRhoThermo::New(fluidRegions[i]).ptr()
basicRhoThermo::New(procMeshes[procI]).ptr()
);
Info<< " Adding to rhoFluid\n" << endl;
......@@ -37,7 +43,7 @@
(
"rho",
runTime.timeName(),
fluidRegions[i],
procMeshes[procI],
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
......@@ -55,7 +61,7 @@
(
"K",
runTime.timeName(),
fluidRegions[i],
procMeshes[procI],
IOobject::NO_READ,
IOobject::NO_WRITE
),
......@@ -73,11 +79,11 @@
(
"U",
runTime.timeName(),
fluidRegions[i],
procMeshes[procI],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
fluidRegions[i]
procMeshes[procI]
)
);
......@@ -91,12 +97,12 @@
(
"phi",
runTime.timeName(),
fluidRegions[i],
procMeshes[procI],
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(rhoFluid[i]*UFluid[i])
& fluidRegions[i].Sf()
& procMeshes[procI].Sf()
)
);
......@@ -110,7 +116,7 @@
(
"g",
runTime.constant(),
fluidRegions[i],
procMeshes[procI],
IOobject::MUST_READ,
IOobject::NO_WRITE
)
......@@ -137,14 +143,14 @@
ghFluid.set
(
i,
new volScalarField("gh", gFluid[i] & fluidRegions[i].C())
new volScalarField("gh", gFluid[i] & procMeshes[procI].C())
);
Info<< " Adding to ghfFluid\n" << endl;
ghfFluid.set
(
i,
new surfaceScalarField("ghf", gFluid[i] & fluidRegions[i].Cf())
new surfaceScalarField("ghf", gFluid[i] & procMeshes[procI].Cf())
);
p_rghFluid.set
......@@ -156,11 +162,11 @@
(
"p_rgh",
runTime.timeName(),
fluidRegions[i],
procMeshes[procI],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
fluidRegions[i]
procMeshes[procI]
)
);