zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: overset: various improvements in the framework

Browse Source 
The improvements include:

- Allowing overset patches to be displaced outside background domain.
  - The approach does not support overlapping of multiple inset meshes
    on top of background domain.
- Allowing fringe faces to walk away from hole cells in background domain.
  - The approach was not extensibly tested with overlapping patches.
- Improving mass conservation.
- Various experimental entries are removed: massFluxInterpolation, ddtCorr.
- New entries:
  - oversetAdjustPhi: adds a flux correction outside the pressure equation.
  - massCorrection: adds an implicit correction.
This commit is contained in:
sergio
2021-08-23 11:07:50 -07:00
committed by Kutalmis Bercin
parent 8c0679d25f
commit 2a406bbb25
80 changed files with 5118 additions and 2187 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/Make/options
View File

@ -3,6 +3,7 @@ EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/overset/lnInclude \
-I$(LIB_SRC)/overset/include/lnInclude
EXE_LIBS = \
-lfiniteVolume \

22
applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/createFields.H
View File

@ -31,3 +31,25 @@
Info<< "Reading diffusivity DT\n" << endl;
dimensionedScalar DT("DT", dimViscosity, transportProperties);
bool oversetPatchErrOutput =
simple.dict().getOrDefault("oversetPatchErrOutput", false);
// Dummy phi for oversetPatchErrOutput
tmp<surfaceScalarField> tdummyPhi;
if (oversetPatchErrOutput)
{
tdummyPhi = tmp<surfaceScalarField>::New
(
IOobject
(
"dummyPhi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar(dimless, Zero)
);
}

8
applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/overLaplacianDyMFoam.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2015 OpenFOAM Foundation
Copyright (C) 2016-2017 OpenCFD Ltd.
Copyright (C) 2016-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -58,6 +58,7 @@ Description
#include "fvOptions.H"
#include "simpleControl.H"
#include "dynamicFvMesh.H"
#include "oversetPatchPhiErr.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -99,6 +100,11 @@ int main(int argc, char *argv[])
fvOptions.constrain(TEqn);
TEqn.solve();
fvOptions.correct(T);
if (oversetPatchErrOutput)
{
oversetPatchPhiErr(TEqn, tdummyPhi.ref());
}
}
#include "write.H"

3
applications/solvers/basic/potentialFoam/overPotentialFoam/overPotentialFoam.C
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenCFD Ltd
Copyright (C) 2017-2022 OpenCFD Ltd
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -149,7 +149,6 @@ int main(int argc, char *argv[])
mesh.update();
surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
// Since solver contains no time loop it would never execute
// function objects so do it ourselves

4
applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createControls.H
View File

@ -1,4 +0,0 @@
bool ddtCorr
(
pimple.dict().getOrDefault("ddtCorr", true)
);

2
applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H
View File

@ -69,6 +69,8 @@ mesh.setFluxRequired(p.name());
// Mask field for zeroing out contributions on hole cells
#include "createCellMask.H"
// Create bool field with interpolated cells
#include "createInterpolatedCells.H"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence

50
applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/overRhoPimpleDyMFoam.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2015 OpenFOAM Foundation
Copyright (C) 2016-2017 OpenCFD Ltd.
Copyright (C) 2016-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -43,7 +43,6 @@ Description
#include "dynamicFvMesh.H"
#include "fluidThermo.H"
#include "turbulentFluidThermoModel.H"
#include "bound.H"
#include "pimpleControl.H"
#include "pressureControl.H"
#include "CorrectPhi.H"
@ -89,10 +88,8 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
// Store divrhoU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
@ -128,7 +125,6 @@ int main(int argc, char *argv[])
{
if (pimple.firstIter() || moveMeshOuterCorrectors)
{
// Do any mesh changes
mesh.update();
@ -137,52 +133,22 @@ int main(int argc, char *argv[])
MRF.update();
#include "setCellMask.H"
const surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
// Zero Uf on old faceMask (H-I)
rhoUf() *= faceMaskOld;
surfaceVectorField rhoUfint(fvc::interpolate(rho*U));
// Update Uf and phi on new C-I faces
rhoUf() += (1-faceMaskOld)*rhoUfint;
// Update Uf boundary
forAll(rhoUf().boundaryField(), patchI)
{
rhoUf().boundaryFieldRef()[patchI] =
rhoUfint.boundaryField()[patchI];
}
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & rhoUf();
#include "setInterpolatedCells.H"
#include "correctRhoPhiFaceMask.H"
if (correctPhi)
{
// Corrects flux on separated regions
#include "correctPhi.H"
}
// Zero phi on current H-I
const surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
phi *= faceMask;
U *= cellMask;
// Make the fluxes relative to the mesh-motion
fvc::makeRelative(phi, rho, U);
}
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}

15
applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H
View File

@ -25,17 +25,6 @@ surfaceScalarField phiHbyA
fvc::interpolate(rho)*fvc::flux(HbyA)
);
if (ddtCorr)
{
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
phiHbyA +=
faceMaskOld*MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf));
}
fvc::makeRelative(phiHbyA, rho, U);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
@ -134,8 +123,4 @@ if (thermo.dpdt())
}
}
surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
phi *= faceMask;

9
applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/readControls.H
View File

@ -1,9 +0,0 @@
#include "readTimeControls.H"
correctPhi = pimple.dict().getOrDefault("correctPhi", false);
checkMeshCourantNo =
pimple.dict().getOrDefault("checkMeshCourantNo", false);
ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);

3
applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/createFields.H
View File

@ -124,3 +124,6 @@ dimensionedScalar initialMass("initialMass", fvc::domainIntegrate(rho));
// Mask field for zeroing out contributions on hole cells
#include "createCellMask.H"
// Create bool field with interpolated cells
#include "createInterpolatedCells.H"

41
applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/overBuoyantPimpleDyMFoam.C
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -50,6 +50,7 @@ Description
#include "CorrectPhi.H"
#include "cellCellStencilObject.H"
#include "localMin.H"
#include "oversetAdjustPhi.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -86,9 +87,6 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readTimeControls.H"
#include "readControls.H"
#include "readDyMControls.H"
#include "compressibleCourantNo.H"
@ -128,45 +126,14 @@ int main(int argc, char *argv[])
MRF.update();
#include "setCellMask.H"
const surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
// Zero Uf on old faceMask (H-I)
rhoUf() *= faceMaskOld;
//fvc::correctRhoUf(rhoUfint, rho, U, phi);
surfaceVectorField rhoUfint(fvc::interpolate(rho*U));
// Update Uf and phi on new C-I faces
rhoUf() += (1-faceMaskOld)*rhoUfint;
// Update Uf boundary
forAll(rhoUf().boundaryField(), patchI)
{
rhoUf().boundaryFieldRef()[patchI] =
rhoUfint.boundaryField()[patchI];
}
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & rhoUf();
#include "setInterpolatedCells.H"
#include "correctRhoPhiFaceMask.H"
if (correctPhi)
{
#include "correctPhi.H"
}
// Zero phi on current H-I
const surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
phi *= faceMask;
U *= cellMask;
// Make the fluxes relative to the mesh-motion
fvc::makeRelative(phi, rho, U);
}

17
applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/pEqn.H
View File

@ -21,17 +21,14 @@ surfaceScalarField phiHbyA
fvc::flux(rho*HbyA) + phig
);
if (ddtCorr)
if (adjustFringe)
{
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
phiHbyA +=
faceMaskOld*MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi));
fvc::makeRelative(phiHbyA,rho, U);
oversetAdjustPhi(phiHbyA, U);
fvc::makeAbsolute(phiHbyA,rho, U);
}
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
@ -122,8 +119,4 @@ if (thermo.dpdt())
}
}
surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
phi *= faceMask;

9
applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/readControls.H
View File

@ -1,9 +0,0 @@
#include "readTimeControls.H"
correctPhi = pimple.dict().getOrDefault("correctPhi", false);
checkMeshCourantNo =
pimple.dict().getOrDefault("checkMeshCourantNo", false);
ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);

7
applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/CourantNo.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011 OpenFOAM Foundation
Copyright (C) 2016 OpenCFD Ltd.
Copyright (C) 2016-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -37,7 +37,10 @@ scalar meanCoNum = 0.0;
if (mesh.nInternalFaces())
{
surfaceScalarField phiMask(localMin<scalar>(mesh).interpolate(cellMask));
surfaceScalarField phiMask
(
localMin<scalar>(mesh).interpolate(cellMask + interpolatedCells)
);
scalarField sumPhi(fvc::surfaceSum(mag(phiMask*phi))().internalField());

1
applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/UEqn.H
View File

@ -1,5 +1,4 @@
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);
tmp<fvVectorMatrix> tUEqn

26
applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createControls.H
View File

@ -1,26 +0,0 @@
#include "createTimeControls.H"
bool correctPhi
(
pimple.dict().getOrDefault("correctPhi", false)
);
bool checkMeshCourantNo
(
pimple.dict().getOrDefault("checkMeshCourantNo", false)
);
bool massFluxInterpolation
(
pimple.dict().getOrDefault("massFluxInterpolation", false)
);
bool adjustFringe
(
pimple.dict().getOrDefault("oversetAdjustPhi", false)
);
bool ddtCorr
(
pimple.dict().getOrDefault("ddtCorr", true)
);

273
applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H
View File

@ -1,273 +0,0 @@
// Interpolation used
interpolationCellPoint<vector> UInterpolator(HbyA);
// Determine faces on outside of interpolated cells
bitSet isOwnerInterpolatedFace(mesh.nInternalFaces());
bitSet isNeiInterpolatedFace(mesh.nInternalFaces());
// Determine donor cells
labelListList donorCell(mesh.nInternalFaces());
scalarListList weightCellCells(mesh.nInternalFaces());
// Interpolated HbyA faces
vectorField UIntFaces(mesh.nInternalFaces(), Zero);
// Determine receptor neighbour cells
labelList receptorNeigCell(mesh.nInternalFaces(), -1);
{
const cellCellStencilObject& overlap = Stencil::New(mesh);
const labelList& cellTypes = overlap.cellTypes();
const labelIOList& zoneID = overlap.zoneID();
label nZones = gMax(zoneID)+1;
PtrList<fvMeshSubset> meshParts(nZones);
labelList nCellsPerZone(nZones, Zero);
// A mesh subset for each zone
forAll(meshParts, zonei)
{
meshParts.set
(
zonei,
// Select cells where the zoneID == zonei
new fvMeshSubset(mesh, zonei, zoneID)
);
}
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
label ownType = cellTypes[mesh.faceOwner()[faceI]];
label neiType = cellTypes[mesh.faceNeighbour()[faceI]];
if
(
ownType == cellCellStencil::INTERPOLATED
&& neiType == cellCellStencil::CALCULATED
)
{
isOwnerInterpolatedFace.set(faceI);
const vector& fc = mesh.faceCentres()[faceI];
for (label zoneI = 0; zoneI < nZones; zoneI++)
{
if (zoneI != zoneID[mesh.faceOwner()[faceI]])
{
const fvMesh& partMesh = meshParts[zoneI].subMesh();
const labelList& cellMap = meshParts[zoneI].cellMap();
label cellI = partMesh.findCell(fc);
if (cellI != -1)
{
// Determine weights
labelList stencil(partMesh.cellCells()[cellI]);
stencil.append(cellI);
label st = stencil.size();
donorCell[faceI].setSize(st);
weightCellCells[faceI].setSize(st);
scalarField weights(st);
forAll(stencil, i)
{
scalar d = mag
(
partMesh.cellCentres()[stencil[i]]
- fc
);
weights[i] = 1.0/d;
donorCell[faceI][i] = cellMap[stencil[i]];
}
weights /= sum(weights);
weightCellCells[faceI] = weights;
forAll(stencil, i)
{
UIntFaces[faceI] +=
weightCellCells[faceI][i]
*UInterpolator.interpolate
(
fc,
donorCell[faceI][i]
);
}
break;
}
}
}
receptorNeigCell[faceI] = mesh.faceNeighbour()[faceI];
}
else if
(
ownType == cellCellStencil::CALCULATED
&& neiType == cellCellStencil::INTERPOLATED
)
{
isNeiInterpolatedFace.set(faceI);
const vector& fc = mesh.faceCentres()[faceI];
for (label zoneI = 0; zoneI < nZones; zoneI++)
{
if (zoneI != zoneID[mesh.faceNeighbour()[faceI]])
{
const fvMesh& partMesh = meshParts[zoneI].subMesh();
const labelList& cellMap = meshParts[zoneI].cellMap();
label cellI = partMesh.findCell(fc);
if (cellI != -1)
{
// Determine weights
labelList stencil(partMesh.cellCells()[cellI]);
stencil.append(cellI);
label st = stencil.size();
donorCell[faceI].setSize(st);
weightCellCells[faceI].setSize(st);
scalarField weights(st);
forAll(stencil, i)
{
scalar d = mag
(
partMesh.cellCentres()[stencil[i]]
- fc
);
weights[i] = 1.0/d;
donorCell[faceI][i] = cellMap[stencil[i]];
}
weights /= sum(weights);
weightCellCells[faceI] = weights;
forAll(stencil, i)
{
UIntFaces[faceI] +=
weightCellCells[faceI][i]
*UInterpolator.interpolate
(
fc,
donorCell[faceI][i]
);
}
break;
}
}
}
receptorNeigCell[faceI] = mesh.faceOwner()[faceI];
}
}
}
// contravariant U
vectorField U1Contrav(mesh.nInternalFaces(), Zero);
surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf());
forAll(isNeiInterpolatedFace, faceI)
{
label cellId = -1;
if (isNeiInterpolatedFace.test(faceI))
{
cellId = mesh.faceNeighbour()[faceI];
}
else if (isOwnerInterpolatedFace.test(faceI))
{
cellId = mesh.faceOwner()[faceI];
}
if (cellId != -1)
{
const vector& n = faceNormals[faceI];
vector n1(Zero);
// 2-D cases
if (mesh.nSolutionD() == 2)
{
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
if (mesh.geometricD()[cmpt] == -1)
{
switch (cmpt)
{
case vector::X:
{
n1 = vector(0, n.z(), -n.y());
break;
}
case vector::Y:
{
n1 = vector(n.z(), 0, -n.x());
break;
}
case vector::Z:
{
n1 = vector(n.y(), -n.x(), 0);
break;
}
}
}
}
}
else if (mesh.nSolutionD() == 3)
{
//Determine which is the primary direction
if (mag(n.x()) > mag(n.y()) && mag(n.x()) > mag(n.z()))
{
n1 = vector(n.y(), -n.x(), 0);
}
else if (mag(n.y()) > mag(n.z()))
{
n1 = vector(0, n.z(), -n.y());
}
else
{
n1 = vector(-n.z(), 0, n.x());
}
}
n1.normalise();
const vector n2 = normalised(n ^ n1);
tensor rot =
tensor
(
n.x() ,n.y(), n.z(),
n1.x() ,n1.y(), n1.z(),
n2.x() ,n2.y(), n2.z()
);
// tensor rot =
// tensor
// (
// n & x ,n & y, n & z,
// n1 & x ,n1 & y, n1 & z,
// n2 & x ,n2 & y, n2 & z
// );
U1Contrav[faceI].x() =
2*transform(rot, UIntFaces[faceI]).x()
- transform(rot, HbyA[receptorNeigCell[faceI]]).x();
U1Contrav[faceI].y() = transform(rot, HbyA[cellId]).y();
U1Contrav[faceI].z() = transform(rot, HbyA[cellId]).z();
HbyA[cellId] = transform(inv(rot), U1Contrav[faceI]);
}
}

55
applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/overPimpleDyMFoam.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2016-2018 OpenCFD Ltd.
Copyright (C) 2016-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -46,12 +46,9 @@ Description
#include "fvOptions.H"
#include "cellCellStencilObject.H"
#include "zeroGradientFvPatchFields.H"
#include "localMin.H"
#include "interpolationCellPoint.H"
#include "transform.H"
#include "fvMeshSubset.H"
#include "oversetAdjustPhi.H"
#include "oversetPatchPhiErr.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -68,10 +65,9 @@ int main(int argc, char *argv[])
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createDyMControls.H"
#include "initContinuityErrs.H"
pimpleControl pimple(mesh);
#include "createFields.H"
#include "createUf.H"
#include "createMRF.H"
@ -88,7 +84,9 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
#include "readOversetDyMControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
@ -97,45 +95,20 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
bool changed = mesh.update();
mesh.update();
if (changed)
if (mesh.changing())
{
#include "setCellMask.H"
#include "setInterpolatedCells.H"
#include "correctPhiFaceMask.H"
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
fvc::makeRelative(phi, U);
// Zero Uf on old faceMask (H-I)
Uf *= faceMaskOld;
// Update Uf and phi on new C-I faces
Uf += (1-faceMaskOld)*fvc::interpolate(U);
phi = mesh.Sf() & Uf;
// Zero phi on current H-I
surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
phi *= faceMask;
}
if (mesh.changing() && correctPhi)
{
// Calculate absolute flux from the mapped surface velocity
#include "correctPhi.H"
}
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
if (mesh.changing() && checkMeshCourantNo)
{
#include "meshCourantNo.H"
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
// --- Pressure-velocity PIMPLE corrector loop

82
applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H
View File

@ -1,36 +1,11 @@
// Option 1: interpolate rAU, do not block out rAU on blocked cells
volScalarField rAU("rAU", 1.0/UEqn.A());
mesh.interpolate(rAU);
// Option 2: do not interpolate rAU but block out rAU
//surfaceScalarField rAUf("rAUf", fvc::interpolate(blockedCells*rAU));
// Option 3: do not interpolate rAU but zero out rAUf on faces on holes
// But what about:
//
// H
// H I C C C C
// H
//
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField H("H", UEqn.H());
volVectorField HbyA("HbyA", U);
HbyA = constrainHbyA(rAU*H, U, p);
if (massFluxInterpolation)
{
#include "interpolatedFaces.H"
}
if (runTime.outputTime())
{
H.write();
rAU.write();
HbyA.write();
}
if (pimple.nCorrPISO() <= 1)
{
tUEqn.clear();
@ -38,33 +13,16 @@ if (pimple.nCorrPISO() <= 1)
phiHbyA = fvc::flux(HbyA);
if (ddtCorr)
{
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
phiHbyA += rAUf*faceMaskOld*fvc::ddtCorr(U, Uf);
}
MRF.makeRelative(phiHbyA);
// WIP
if (p.needReference())
{
fvc::makeRelative(phiHbyA, U);
adjustPhi(phiHbyA, U, p);
fvc::makeAbsolute(phiHbyA, U);
}
if (adjustFringe)
{
fvc::makeRelative(phiHbyA, U);
oversetAdjustPhi(phiHbyA, U);
oversetAdjustPhi(phiHbyA, U, zoneIdMass);
fvc::makeAbsolute(phiHbyA, U);
}
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
@ -79,27 +37,26 @@ while (pimple.correctNonOrthogonal())
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
// option 2:
// rAUf*fvc::snGrad(p)*mesh.magSf();
pEqn.relax();
U =
cellMask*
(
HbyA
- rAU*fvc::reconstruct((pEqn.flux())/rAUf)
);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
if (oversetPatchErrOutput)
{
oversetPatchPhiErr(pEqn, phiHbyA);
}
}
// Excludes error in interpolated/hole cells
#include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
volVectorField gradP(fvc::grad(p));
// Option 2: zero out velocity on blocked out cells
//U = HbyA - rAU*cellMask*gradP;
// Option 3: zero out velocity on blocked out cells
// This is needed for the scalar Eq (k,epsilon, etc)
// which can use U as source term
U = cellMask*(HbyA - rAU*gradP);
U.correctBoundaryConditions();
fvOptions.correct(U);
{
Uf = fvc::interpolate(U);
@ -109,9 +66,4 @@ fvOptions.correct(U);
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
phi *= faceMask;

10
applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/readControls.H
View File

@ -1,10 +0,0 @@
#include "readTimeControls.H"
correctPhi = pimple.dict().getOrDefault("correctPhi", false);
checkMeshCourantNo = pimple.dict().getOrDefault("checkMeshCourantNo", false);
massFluxInterpolation =
pimple.dict().getOrDefault("massFluxInterpolation", false);
ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);

4
applications/solvers/incompressible/simpleFoam/overSimpleFoam/createOversetFields.H
View File

@ -24,7 +24,3 @@ bool adjustFringe
(
simple.dict().getOrDefault("oversetAdjustPhi", false)
);
bool massFluxInterpolation
(
simple.dict().getOrDefault("massFluxInterpolation", false)
);

8
applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H
View File

@ -1,18 +1,10 @@
{
surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
//mesh.interpolate(HbyA);
if (massFluxInterpolation)
{
#include "interpolatedFaces.H"
}
tUEqn.clear();
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));

5
applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/createFields.H
View File

@ -129,10 +129,5 @@ compressibleInterPhaseTransportModel turbulence
// Mask field for zeroing out contributions on hole cells
#include "createCellMask.H"
surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
// Create bool field with interpolated cells
#include "createInterpolatedCells.H"

41
applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/overCompressibleInterDyMFoam.C
View File

@ -80,9 +80,6 @@ int main(int argc, char *argv[])
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
pimpleControl pimple(mesh);
#include "createTimeControls.H"
#include "createDyMControls.H"
#include "createFields.H"
@ -109,7 +106,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
if (LTS)
{
@ -154,40 +151,10 @@ int main(int argc, char *argv[])
// Update cellMask field for blocking out hole cells
#include "setCellMask.H"
#include "setInterpolatedCells.H"
faceMask =
localMin<scalar>(mesh).interpolate(cellMask.oldTime());
// Zero Uf on old faceMask (H-I)
Uf *= faceMask;
const surfaceVectorField Uint(fvc::interpolate(U));
// Update Uf and phi on new C-I faces
Uf += (1-faceMask)*Uint;
// Update Uf boundary
forAll(Uf.boundaryField(), patchI)
{
Uf.boundaryFieldRef()[patchI] =
Uint.boundaryField()[patchI];
}
phi = mesh.Sf() & Uf;
// Correct phi on individual regions
if (correctPhi)
{
#include "correctPhi.H"
}
#include "correctPhiFaceMask.H"
mixture.correct();
// Zero phi on current H-I
faceMask = localMin<scalar>(mesh).interpolate(cellMask);
phi *= faceMask;
U *= cellMask;
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
@ -202,10 +169,6 @@ int main(int argc, char *argv[])
#include "alphaControls.H"
#include "compressibleAlphaEqnSubCycle.H"
const surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
rhoPhi *= faceMask;
turbulence.correctPhasePhi();

13
applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/pEqn.H
View File

@ -10,19 +10,6 @@
fvc::flux(HbyA)
);
if (ddtCorr)
{
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
phiHbyA +=
MRF.zeroFilter
(
fvc::interpolate(rho*rAU)*faceMaskOld*fvc::ddtCorr(U, Uf)
);
}
MRF.makeRelative(phiHbyA);
surfaceScalarField phig

7
applications/solvers/multiphase/interFoam/overInterDyMFoam/CourantNo.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011 OpenFOAM Foundation
Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -36,7 +36,10 @@ scalar meanCoNum = 0.0;
if (mesh.nInternalFaces())
{
surfaceScalarField phiMask(localMin<scalar>(mesh).interpolate(cellMask));
surfaceScalarField phiMask
(
localMin<scalar>(mesh).interpolate(cellMask + interpolatedCells)
);
scalarField sumPhi
(

2
applications/solvers/multiphase/interFoam/overInterDyMFoam/alphaSuSp.H
View File

@ -7,3 +7,5 @@ volScalarField::Internal divU
? fvc::div(phiCN() + mesh.phi())
: fvc::div(phiCN())
);
divU *= interpolatedCells*cellMask;

30
applications/solvers/multiphase/interFoam/overInterDyMFoam/createDyMControls.H
View File

@ -1,30 +0,0 @@
bool correctPhi
(
pimple.dict().getOrDefault("correctPhi", true)
);
bool checkMeshCourantNo
(
pimple.dict().getOrDefault("checkMeshCourantNo", false)
);
bool moveMeshOuterCorrectors
(
pimple.dict().getOrDefault("moveMeshOuterCorrectors", false)
);
bool massFluxInterpolation
(
pimple.dict().getOrDefault("massFluxInterpolation", false)
);
bool adjustFringe
(
pimple.dict().getOrDefault("oversetAdjustPhi", false)
);
bool ddtCorr
(
pimple.dict().getOrDefault("ddtCorr", true)
);

4
applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H
View File

@ -40,8 +40,12 @@ volVectorField U
nonInt.insert("HbyA");
nonInt.insert("grad(p_rgh)");
nonInt.insert("nHat");
nonInt.insert("surfaceIntegrate(nHatf)");
nonInt.insert("surfaceIntegrate(phi+meshPhi)");
nonInt.insert("surfaceIntegrate(phi)");
nonInt.insert("surfaceIntegrate(phiHbyA)");
nonInt.insert("surfaceSum(((S|magSf)*S)");
nonInt.insert("surfaceIntegrate(((rAUf*magSf)*snGradCorr(p_rgh)))");
nonInt.insert("cellMask");
nonInt.insert("cellDisplacement");
nonInt.insert("interpolatedCells");

65
applications/solvers/multiphase/interFoam/overInterDyMFoam/overInterDyMFoam.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2014 OpenFOAM Foundation
Copyright (C) 2016-2017 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -49,14 +49,11 @@ Description
#include "turbulentTransportModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
#include "fvcSmooth.H"
#include "cellCellStencilObject.H"
#include "localMin.H"
#include "interpolationCellPoint.H"
#include "transform.H"
#include "fvMeshSubset.H"
#include "oversetAdjustPhi.H"
#include "oversetPatchPhiErr.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -76,8 +73,7 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
pimpleControl pimple(mesh);
#include "createTimeControls.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
@ -97,11 +93,8 @@ int main(int argc, char *argv[])
dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0)
);
if (correctPhi)
{
#include "correctPhi.H"
}
#include "createUf.H"
#include "createControls.H"
#include "setCellMask.H"
#include "setInterpolatedCells.H"
@ -119,7 +112,8 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
#include "readOversetDyMControls.H"
if (LTS)
{
@ -158,50 +152,17 @@ int main(int argc, char *argv[])
talphaPhi1Corr0.clear();
}
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
// Update cellMask field for blocking out hole cells
#include "setCellMask.H"
#include "setInterpolatedCells.H"
#include "correctPhiFaceMask.H"
const surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
// Zero Uf on old faceMask (H-I)
Uf *= faceMaskOld;
const surfaceVectorField Uint(fvc::interpolate(U));
// Update Uf and phi on new C-I faces
Uf += (1-faceMaskOld)*Uint;
// Update Uf boundary
forAll(Uf.boundaryField(), patchI)
{
Uf.boundaryFieldRef()[patchI] =
Uint.boundaryField()[patchI];
}
phi = mesh.Sf() & Uf;
// Correct phi on individual regions
if (correctPhi)
{
#include "correctPhi.H"
}
mixture.correct();
// Zero phi on current H-I
const surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
phi *= faceMask;
U *= cellMask;
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
@ -213,14 +174,14 @@ int main(int argc, char *argv[])
}
}
if (adjustFringe)
{
oversetAdjustPhi(phi, U, zoneIdMass);
}
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
const surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
rhoPhi *= faceMask;
mixture.correct();

45
applications/solvers/multiphase/interFoam/overInterDyMFoam/pEqn.H
View File

@ -1,64 +1,31 @@
{
rAU = 1.0/UEqn.A();
//mesh.interpolate(rAU);
surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField H("H", UEqn.H());
volVectorField HbyA("HbyA", U);
//HbyA = rAU*UEqn.H();
HbyA = constrainHbyA(rAU*H, U, p_rgh);
if (massFluxInterpolation)
{
#include "interpolatedFaces.H"
}
HbyA = constrainHbyA(rAU*H, U, p_rgh);
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
if (ddtCorr)
{
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
phiHbyA +=
fvc::interpolate(rho*rAU)*faceMaskOld*fvc::ddtCorr(U, Uf);
}
MRF.makeRelative(phiHbyA);
if (p_rgh.needReference())
{
fvc::makeRelative(phiHbyA, U);
adjustPhi(phiHbyA, U, p_rgh);
fvc::makeAbsolute(phiHbyA, U);
}
if (adjustFringe)
{
fvc::makeRelative(phiHbyA, U);
oversetAdjustPhi(phiHbyA, U);
fvc::makeAbsolute(phiHbyA, U);
}
surfaceScalarField phig
(
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*faceMask*rAUf*mesh.magSf()
- ghf*fvc::snGrad(cellMask*rho)
)*rAUf*faceMask*mesh.magSf()
);
phiHbyA += phig;
if (adjustFringe)
{
fvc::makeRelative(phiHbyA, U);
oversetAdjustPhi(phiHbyA, U);
fvc::makeAbsolute(phiHbyA, U);
oversetAdjustPhi(phiHbyA, U, zoneIdMass);
}
// Update the pressure BCs to ensure flux consistency
@ -90,6 +57,10 @@
U.correctBoundaryConditions();
fvOptions.correct(U);
}
if (oversetPatchErrOutput)
{
oversetPatchPhiErr(p_rghEqn, phiHbyA);
}
}
#include "continuityErrs.H"

16
applications/solvers/multiphase/interFoam/overInterDyMFoam/readControls.H
View File

@ -1,16 +0,0 @@
#include "readTimeControls.H"
correctPhi = pimple.dict().getOrDefault("correctPhi", false);
checkMeshCourantNo =
pimple.dict().getOrDefault("checkMeshCourantNo", false);
moveMeshOuterCorrectors =
pimple.dict().getOrDefault("moveMeshOuterCorrectors", false);
massFluxInterpolation =
pimple.dict().getOrDefault("massFluxInterpolation", false);
ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);
adjustFringe = pimple.dict().getOrDefault("oversetAdjustPhi", false);

11
applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/correctPhi.H
View File

@ -1,11 +0,0 @@
CorrectPhi
(
U,
phi,
p_rgh,
surfaceScalarField("rAUf", fvc::interpolate(rAU)),
divU,
pimple
);
#include "continuityErrs.H"

5
applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/createFields.H
View File

@ -149,10 +149,5 @@ surfaceScalarField alphaPhi10
// Mask field for zeroing out contributions on hole cells
#include "createCellMask.H"
surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
// Create bool field with interpolated cells
#include "createInterpolatedCells.H"

40
applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/overInterPhaseChangeDyMFoam.C
View File

@ -82,9 +82,7 @@ int main(int argc, char *argv[])
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
pimpleControl pimple(mesh);
#include "createTimeControls.H"
#include "createDyMControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
@ -118,7 +116,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
// Store divU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
@ -153,40 +151,10 @@ int main(int argc, char *argv[])
// Update cellMask field for blocking out hole cells
#include "setCellMask.H"
#include "setInterpolatedCells.H"
faceMask =
localMin<scalar>(mesh).interpolate(cellMask.oldTime());
// Zero Uf on old faceMask (H-I)
Uf *= faceMask;
const surfaceVectorField Uint(fvc::interpolate(U));
// Update Uf and phi on new C-I faces
Uf += (1-faceMask)*Uint;
// Update Uf boundary
forAll(Uf.boundaryField(), patchI)
{
Uf.boundaryFieldRef()[patchI] =
Uint.boundaryField()[patchI];
}
phi = mesh.Sf() & Uf;
if (correctPhi)
{
#include "correctPhi.H"
}
#include "correctPhiFaceMask.H"
mixture->correct();
// Zero phi on current H-I
faceMask = localMin<scalar>(mesh).interpolate(cellMask);
phi *= faceMask;
U *= cellMask;
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
}
@ -214,10 +182,6 @@ int main(int argc, char *argv[])
mixture->correct();
#include "alphaEqnSubCycle.H"
const surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
rhoPhi *= faceMask;
interface.correct();

10
applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/pEqn.H
View File

@ -8,16 +8,6 @@
fvc::flux(HbyA)
);
if (ddtCorr)
{
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
phiHbyA += faceMaskOld*fvc::ddtCorr(U, Uf);
}
if (p_rgh.needReference())
{
fvc::makeRelative(phiHbyA, U);

6
src/fvOptions/Make/options
View File

@ -13,7 +13,8 @@ EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-I$(LIB_SRC)/overset/lnInclude
LIB_LIBS = \
-lfiniteVolume \
@ -26,4 +27,5 @@ LIB_LIBS = \
-lturbulenceModels \
-lincompressibleTurbulenceModels \
-lcompressibleTurbulenceModels \
-lreactionThermophysicalModels
-lreactionThermophysicalModels \
-loverset

22
src/fvOptions/cellSetOption/cellSetOption.C
View File

@ -30,6 +30,7 @@ License
#include "cellSet.H"
#include "cellBitSet.H"
#include "volFields.H"
#include "cellCellStencilObject.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -53,6 +54,7 @@ Foam::fv::cellSetOption::selectionModeTypeNames_
{ selectionModeType::smPoints, "points" },
{ selectionModeType::smCellSet, "cellSet" },
{ selectionModeType::smCellZone, "cellZone" },
{ selectionModeType::smCellType, "cellType" }
});
@ -97,6 +99,10 @@ void Foam::fv::cellSetOption::setSelection(const dictionary& dict)
}
break;
}
case smCellType:
{
break;
}
default:
{
FatalErrorInFunction
@ -234,6 +240,21 @@ void Foam::fv::cellSetOption::setCellSelection()
}
break;
}
case smCellType:
{
labelHashSet selectedCells;
const cellCellStencilObject& overlap = Stencil::New(mesh_);
const labelList& cellTypes = overlap.cellTypes();
forAll(cellTypes, celli)
{
if (cellTypes[celli] == cellCellStencil::POROUS)
{
selectedCells.insert(celli);
}
cells_ = selectedCells.sortedToc();
}
break;
}
default:
{
FatalErrorInFunction
@ -300,6 +321,7 @@ bool Foam::fv::cellSetOption::isActive()
(
selectionMode_ == smGeometric
|| selectionMode_ == smPoints
|| selectionMode_ == smCellType
)
{
// Geometric selection mode(s)

3
src/fvOptions/cellSetOption/cellSetOption.H
View File

@ -162,7 +162,8 @@ public:
smCellSet, //!< "cellSet"
smCellZone, //!< "cellZone"
smPoints, //!< "points"
smGeometric //!< "geometric"
smGeometric, //!< "geometric"
smCellType //!< "overset type cells"
};
//- List of selection mode type names

10
src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.C
View File

@ -30,6 +30,7 @@ License
#include "addToRunTimeSelectionTable.H"
#include "fvMatrix.H"
#include "volFields.H"
#include "cellCellStencilObject.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -67,14 +68,14 @@ void Foam::fv::velocityDampingConstraint::addDamping(fvMatrix<vector>& eqn)
label nDamped(0);
const label nTotCells(returnReduce(cells_.size(), sumOp<label>()));
for (label celli : cells_)
for (const label celli : cells_)
{
const scalar magU = mag(U[celli]);
if (magU > UMax_)
{
const scalar scale = sqr(Foam::cbrt(vol[celli]));
diag[celli] += scale*(magU-UMax_);
diag[celli] += C_*scale*(magU-UMax_);
++nDamped;
}
@ -105,7 +106,8 @@ Foam::fv::velocityDampingConstraint::velocityDampingConstraint
const fvMesh& mesh
)
:
fv::cellSetOption(name, modelType, dict, mesh)
fv::cellSetOption(name, modelType, dict, mesh),
C_(1)
{
read(dict);
}
@ -135,6 +137,8 @@ bool Foam::fv::velocityDampingConstraint::read(const dictionary& dict)
{
coeffs_.readEntry("UMax", UMax_);
coeffs_.readIfPresent("C", C_);
if (!coeffs_.readIfPresent("UNames", fieldNames_))
{
fieldNames_.resize(1);

7
src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2015-2020 OpenCFD Ltd.
Copyright (C) 2015-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -62,6 +62,7 @@ Usage
// Optional entries (runtime modifiable)
UNames (<Uname1> <Uname2> ... <UnameN>);
C <scalar>;
// Conditional optional entries (runtime modifiable)
@ -80,6 +81,7 @@ Usage
UMax | Maximum velocity magnitude | scalar | yes | -
UNames | Names of operand velocity fields | wordList | no | -
U | Name of operand velocity field | word | cndtnl | U
C | Model factor | scalar | no | 1
\endtable
The inherited entries are elaborated in:
@ -127,6 +129,9 @@ protected:
//- Maximum velocity magnitude
scalar UMax_;
//- Model factor
scalar C_;
// Protected Member Functions

14
src/overset/Make/files
View File

@ -8,7 +8,10 @@ cellCellStencil/trackingInverseDistance/voxelMeshSearch.C
cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
cellCellStencil/leastSquares/leastSquaresCellCellStencil.C
dynamicOversetFvMesh/dynamicOversetFvMesh.C
oversetFvMesh = oversetFvMesh
$(oversetFvMesh)/oversetFvMeshBase.C
$(oversetFvMesh)/dynamicOversetFvMesh/dynamicOversetFvMesh.C
$(oversetFvMesh)/staticOversetFvMesh/staticOversetFvMesh.C
fvMeshPrimitiveLduAddressing/fvMeshPrimitiveLduAddressing.C
@ -23,8 +26,15 @@ oversetAdjustPhi/oversetAdjustPhi.C
regionsToCell/regionsToCell.C
lduPrimitiveProcessorInterface/GAMG/calculatedProcessorGAMGInterface.C
lduPrimitiveProcessorInterface/GAMG/calculatedProcessorGAMGInterfaceField.C
oversetCoupledPolyPatch/oversetLduInterfaceField/oversetLduInterfaceField.C
oversetCoupledPolyPatch/oversetLduInterface/oversetLduInterface.C
oversetCoupledPolyPatch/oversetGAMGInterface/oversetGAMGInterface.C
oversetCoupledPolyPatch/oversetGAMGInterfaceField/oversetGAMGInterfaceField.C
oversetPatchPhiErr/oversetPatchPhiErr.C
LIB = $(FOAM_LIBBIN)/liboverset

4
src/overset/Make/options
View File

@ -1,7 +1,3 @@
/* Add to EXE_INC as required:
-DFULLDEBUG -O0 -g
*/
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/fileFormats/lnInclude \

527
src/overset/cellCellStencil/cellCellStencil/cellCellStencil.C
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017-2021 OpenCFD Ltd.
Copyright (C) 2017-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -30,6 +30,7 @@ License
#include "volFields.H"
#include "syncTools.H"
#include "globalIndex.H"
#include "oversetFvPatchFields.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -48,6 +49,8 @@ Foam::cellCellStencil::cellTypeNames_
{ cellType::CALCULATED, "calculated" },
{ cellType::INTERPOLATED, "interpolated" },
{ cellType::HOLE, "hole" },
{ cellType::SPECIAL, "special" },
{ cellType::POROUS, "porous" },
});
@ -96,6 +99,17 @@ Foam::cellCellStencil::~cellCellStencil()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::word Foam::cellCellStencil::baseName(const word& name)
{
if (name.ends_with("_0"))
{
return baseName(name.substr(0, name.size() - 2));
}
return name;
}
const Foam::labelIOList& Foam::cellCellStencil::zoneID(const fvMesh& mesh)
{
if (!mesh.foundObject<labelIOList>("zoneID"))
@ -128,9 +142,9 @@ const Foam::labelIOList& Foam::cellCellStencil::zoneID(const fvMesh& mesh)
),
mesh
);
forAll(volZoneID, cellI)
forAll(volZoneID, celli)
{
zoneID[cellI] = label(volZoneID[cellI]);
zoneID[celli] = label(volZoneID[celli]);
}
zoneIDPtr->store();
@ -296,5 +310,512 @@ void Foam::cellCellStencil::globalCellCells
}
}
void Foam::cellCellStencil::seedCell
(
const label celli,
const scalar wantedFraction,
bitSet& isFront,
scalarField& fraction
) const
{
const cell& cFaces = mesh_.cells()[celli];
forAll(cFaces, i)
{
const label nbrFacei = cFaces[i];
if (fraction[nbrFacei] < wantedFraction)
{
fraction[nbrFacei] = wantedFraction;
isFront.set(nbrFacei);
}
}
}
void Foam::cellCellStencil::setUpFront
(
const labelList& allCellTypes,
bitSet& isFront
) const
{
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
{
const label ownType = allCellTypes[own[facei]];
const label neiType = allCellTypes[nei[facei]];
if
(
(ownType == HOLE && neiType != HOLE)
|| (ownType != HOLE && neiType == HOLE)
)
{
isFront.set(facei);
}
}
labelList nbrCellTypes;
syncTools::swapBoundaryCellList(mesh_, allCellTypes, nbrCellTypes);
for
(
label facei = mesh_.nInternalFaces();
facei < mesh_.nFaces();
facei++
)
{
const label ownType = allCellTypes[own[facei]];
const label neiType = nbrCellTypes[facei-mesh_.nInternalFaces()];
if
(
(ownType == HOLE && neiType != HOLE)
|| (ownType != HOLE && neiType == HOLE)
)
{
isFront.set(facei);
}
}
}
void Foam::cellCellStencil::setUpFrontOnOversetPatch
(
const labelList& allCellTypes,
bitSet& isFront
) const
{
const fvBoundaryMesh& fvm = mesh_.boundary();
// 'overset' patches
forAll(fvm, patchi)
{
if (isA<oversetFvPatch>(fvm[patchi]))
{
const labelList& fc = fvm[patchi].faceCells();
forAll(fc, i)
{
const label celli = fc[i];
if (allCellTypes[celli] == INTERPOLATED)
{
// Note that acceptors might have been marked hole if
// there are no donors in which case we do not want to
// walk this out. This is an extreme situation.
isFront.set(fvm[patchi].start()+i);
}
}
}
}
}
void Foam::cellCellStencil::walkFront
(
const globalIndex& globalCells,
const label layerRelax,
const labelListList& allStencil,
labelList& allCellTypes,
scalarField& allWeight,
const scalarList& compactCellVol,
const labelListList& compactStencil,
const labelList& zoneID,
const label holeLayers,
const label useLayer
) const
{
if (useLayer > holeLayers)
{
FatalErrorInFunction<< "useLayer: " << useLayer
<< " is larger than : " << holeLayers
<< abort(FatalError);
}
if (useLayer == 0)
{
FatalErrorInFunction<< "useLayer: " << useLayer
<< " can not be zero."
<< abort(FatalError);
}
// Current front
bitSet isFront(mesh_.nFaces());
// List of cellTYpes
DynamicList<labelList> dataSet;
// List of cellTypes for increasing layers
DynamicList<scalarField> dAlllWeight;
// List of average volumen ration interpolated/donor
DynamicList<scalar> daverageVolRatio;
// Counting holes of different layers
DynamicList<label> dHoles;
const scalarField& V = mesh_.V();
// Current layer
labelField nLayer(allCellTypes.size(), 0);
for (label currLayer = 1; currLayer < holeLayers+1; ++currLayer)
{
// Set up original front
isFront = false;
setUpFront(allCellTypes, isFront);
labelList allCellTypesWork(allCellTypes);
bitSet isFrontWork(isFront);
label nCurrLayer = currLayer;
while (nCurrLayer > 1 && returnReduce(isFrontWork.any(), orOp<bool>()))
{
bitSet newIsFront(mesh_.nFaces());
forAll(isFrontWork, facei)
{
if (isFrontWork.test(facei))
{
const label own = mesh_.faceOwner()[facei];
if (allCellTypesWork[own] != HOLE)
{
allCellTypesWork[own] = HOLE;
newIsFront.set(mesh_.cells()[own]);
}
if (mesh_.isInternalFace(facei))
{
const label nei = mesh_.faceNeighbour()[facei];
if (allCellTypesWork[nei] != HOLE)
{
allCellTypesWork[nei] = HOLE;
newIsFront.set(mesh_.cells()[nei]);
}
}
}
}
syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
isFrontWork.transfer(newIsFront);
nCurrLayer--;
}
if ((debug&2) && (mesh_.time().outputTime()))
{
tmp<volScalarField> tfld
(
createField
(
mesh_,
"allCellTypesWork_Holes" + name(currLayer),
allCellTypesWork
)
);
tfld().write();
}
if (currLayer == 1)
{
setUpFrontOnOversetPatch(allCellTypes, isFront);
}
if (currLayer > 1)
{
isFront = false;
setUpFrontOnOversetPatch(allCellTypesWork, isFront);
setUpFront(allCellTypesWork, isFront);
}
// Current interpolation fraction
scalarField fraction(mesh_.nFaces(), Zero);
// Ratio between Inter/donor
scalarField volRatio(allCellTypes.size(), Zero);
forAll(isFront, facei)
{
if (isFront.test(facei))
{
fraction[facei] = 1.0;
}
}
scalarField allWeightWork(allCellTypes.size(), Zero);
bitSet nHoles(allCellTypes.size());
while (returnReduce(isFront.any(), orOp<bool>()))
{
// Interpolate cells on front
bitSet newIsFront(mesh_.nFaces());
scalarField newFraction(fraction);
forAll(isFront, facei)
{
if (isFront.test(facei))
{
const label own = mesh_.faceOwner()[facei];
if (allCellTypesWork[own] != HOLE)
{
if (allWeightWork[own] < fraction[facei])
{
// Cell wants to become interpolated (if sufficient
// stencil, otherwise becomes hole)
if (allStencil[own].size())
{
nLayer[own] = currLayer;
allWeightWork[own] = fraction[facei];
allCellTypesWork[own] = INTERPOLATED;
const label donorId = compactStencil[own][0];
volRatio[own] = V[own]/compactCellVol[donorId];
seedCell
(
own,
fraction[facei]-layerRelax,
newIsFront,
newFraction
);
nHoles[own] = true;
}
else
{
allWeightWork[own] = 0.0;
allCellTypesWork[own] = HOLE;
// Add faces of cell as new front
seedCell
(
own,
1.0,
newIsFront,
newFraction
);
}
}
}
if (mesh_.isInternalFace(facei))
{
label nei = mesh_.faceNeighbour()[facei];
if (allCellTypesWork[nei] != HOLE)
{
if (allWeightWork[nei] < fraction[facei])
{
if (allStencil[nei].size())
{
nLayer[nei] = currLayer;
allWeightWork[nei] = fraction[facei];
allCellTypesWork[nei] = INTERPOLATED;
const label donorId = compactStencil[nei][0];
volRatio[nei] = V[nei]/compactCellVol[donorId];
seedCell
(
nei,
fraction[facei]-layerRelax,
newIsFront,
newFraction
);
}
else
{
allWeightWork[nei] = 0.0;
allCellTypesWork[nei] = HOLE;
nHoles[nei] = true;
seedCell
(
nei,
1.0,
newIsFront,
newFraction
);
}
}
}
}
}
}
syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
syncTools::syncFaceList(mesh_, newFraction, maxEqOp<scalar>());
isFront.transfer(newIsFront);
fraction.transfer(newFraction);
}
if ((debug&2) && (mesh_.time().outputTime()))
{
tmp<volScalarField> tfld
(
createField
(
mesh_,
"allCellTypesWork_Layers" + name(currLayer),
allCellTypesWork
)
);
tfld().write();
}
dAlllWeight.append(allWeightWork);
dataSet.append(allCellTypesWork);
// Counting interpolated cells
label count = 1;
forAll (volRatio, i)
{
if (volRatio[i] > 0)
{
count++;
}
}
label nCount(count);
reduce(nCount, sumOp<label>());
const scalar aveVol = mag(gSum(volRatio));
daverageVolRatio.append(aveVol/nCount);
// Check holes number. A sudden increase occurs when the walk leaks
// out from the obstacle
label nTotalHoles(nHoles.count());
reduce(nTotalHoles, sumOp<label>());
dHoles.append(nTotalHoles);
// Check the increase between this layer and the last one
// if over 50% breaks the layer loop.
if
(
(currLayer > 1)
& (nTotalHoles > 2.0*dHoles[currLayer - 1])
)
{
break;
}
}
if ((debug&2) && (mesh_.time().outputTime()))
{
tmp<volScalarField> tfld
(
createField(mesh_, "walkFront_layers", nLayer)
);
tfld().write();
}
// Try to find the best averageVolRatio the further from the initial set
// As this one is next to HOLES
scalarList averageVolRatio;
averageVolRatio.transfer(daverageVolRatio);
label minVolId = findMin(averageVolRatio);
if (useLayer > -1)
{
Info<< nl << " Number of layers : " << averageVolRatio.size() << nl
<< " Average volumetric ratio : " << averageVolRatio << nl
<< " Number of holes cells : " << dHoles << nl
<< " Using layer : " << useLayer << nl << endl;
}
if (useLayer != -1)
{
minVolId = useLayer - 1;
}
allCellTypes.transfer(dataSet[minVolId]);
allWeight.transfer(dAlllWeight[minVolId]);
}
// * * * * * * * * * * * * * Ostream operator * * * * * * * * * * * * * * * //
template<>
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const InfoProxy<cellCellStencil>& ic
)
{
const cellCellStencil& e = ic.t_;
const labelUList& cellTypes = e.cellTypes();
const labelUList& interpolationCells = e.interpolationCells();
const labelListList& cellStencil = e.cellStencil();
labelList nCells(cellCellStencil::count(4, cellTypes));
label nInvalidInterpolated = 0;
label nLocal = 0;
label nMixed = 0;
label nRemote = 0;
forAll(interpolationCells, i)
{
label celli = interpolationCells[i];
const labelList& slots = cellStencil[celli];
if (slots.empty())
{
nInvalidInterpolated++;
}
bool hasLocal = false;
bool hasRemote = false;
forAll(slots, sloti)
{
if (slots[sloti] >= cellTypes.size())
{
hasRemote = true;
}
else
{
hasLocal = true;
}
}
if (hasRemote)
{
if (!hasLocal)
{
nRemote++;
}
else
{
nMixed++;
}
}
else if (hasLocal)
{
nLocal++;
}
}
reduce(nLocal, sumOp<label>());
reduce(nMixed, sumOp<label>());
reduce(nRemote, sumOp<label>());
reduce(nInvalidInterpolated, sumOp<label>());
Info<< "Overset analysis : nCells : "
<< returnReduce(cellTypes.size(), sumOp<label>()) << nl
<< incrIndent
<< indent << "calculated : " << nCells[cellCellStencil::CALCULATED]
<< nl
<< indent << "interpolated : " << nCells[cellCellStencil::INTERPOLATED]
<< " (interpolated from local:" << nLocal
<< " mixed local/remote:" << nMixed
<< " remote:" << nRemote
<< " special:" << nInvalidInterpolated << ")" << nl
<< indent << "hole : " << nCells[cellCellStencil::HOLE] << nl
<< indent << "special : " << nCells[cellCellStencil::SPECIAL] << nl
<< decrIndent << endl;
return os;
}
// ************************************************************************* //

99
src/overset/cellCellStencil/cellCellStencil/cellCellStencil.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017-2019 OpenCFD Ltd.
Copyright (C) 2017-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -74,7 +74,9 @@ public:
{
CALCULATED = 0, // normal operation
INTERPOLATED = 1, // interpolated
HOLE = 2 // hole
HOLE = 2, // hole
SPECIAL = 3, // hole that is donor
POROUS = 4 // tag for special treatment due to lack of donors
};
@ -91,12 +93,12 @@ protected:
//- Set of fields that should not be interpolated
wordHashSet nonInterpolatedFields_;
//- Dictionary of motion control parameters
const dictionary dict_;
// Protected Member Functions
//- Count occurrences (in parallel)
static labelList count(const label size, const labelUList& lst);
//- Helper: create volScalarField for postprocessing.
template<class Type>
static tmp<volScalarField> createField
@ -106,6 +108,9 @@ protected:
const UList<Type>&
);
//- Helper: strip off trailing _0
static word baseName(const word& name);
private:
@ -213,6 +218,9 @@ public:
return zoneID(mesh_);
}
//- Count occurrences (in parallel)
static labelList count(const label size, const labelUList& lst);
//- Helper: create cell-cell addressing in global numbering
static void globalCellCells
(
@ -223,8 +231,89 @@ public:
labelListList& cellCells,
pointListList& cellCellCentres
);
//- Interpolation of acceptor cells from donor cells
template<class T>
static void interpolate
(
Field<T>& psi,
const fvMesh& mesh,
const cellCellStencil& overlap,
const List<scalarList>& wghts
);
//- Explicit interpolation of acceptor cells from donor cells. Looks
//- up cellCellStencil.
template<class T>
void interpolate(const fvMesh& mesh, Field<T>& psi) const;
//- Explicit interpolation of acceptor cells from donor cells with
//- boundary condition handling
template<class GeoField>
void interpolate(GeoField& psi) const;
//- Explicit interpolation of all registered fields. No boundary
//- conditions. Excludes selected fields (and their old-time fields)
template<class GeoField>
void interpolate
(
const fvMesh& mesh,
const wordHashSet& suppressed
) const;
//- Version of correctBoundaryConditions that excludes 'overset' bcs
template<class GeoField, class SuppressBC>
static void correctBoundaryConditions(GeoField& psi);
//- Return info proxy.
//- Used to print stencil information to a stream
InfoProxy<cellCellStencil> info() const
{
return *this;
}
//- Surround holes with layer(s) of interpolated cells
void walkFront
(
const globalIndex& globalCells,
const label layerRelax,
const labelListList& allStencil,
labelList& allCellTypes,
scalarField& allWeight,
const scalarList& compactCellVol,
const labelListList& compactStencil,
const labelList& zoneID,
const label holeLayers,
const label useLayer
) const;
//- Set up front using allCellTypes
void setUpFront
(
const labelList& allCellTypes,
bitSet& isFront
) const;
//- Set up front on overset patches
void setUpFrontOnOversetPatch
(
const labelList& allCellTypes,
bitSet& isFront
) const;
//- Seed faces of cell with wantedFraction (if higher than current)
void seedCell
(
const label cellI,
const scalar wantedFraction,
bitSet& isFront,
scalarField& fraction
) const;
};
template<>
Ostream& operator<<(Ostream& os, const InfoProxy<cellCellStencil>& ic);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

184
src/overset/cellCellStencil/cellCellStencil/cellCellStencilTemplates.C
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2018-2019 OpenCFD Ltd.
Copyright (C) 2018-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -29,6 +29,126 @@ License
// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
template<class T>
void Foam::cellCellStencil::interpolate
(
Field<T>& psi,
const fvMesh& mesh,
const cellCellStencil& overlap,
const List<scalarList>& wghts
)
{
const labelListList& stencil = overlap.cellStencil();
if (stencil.size() != mesh.nCells())
{
return;
}
const mapDistribute& map = overlap.cellInterpolationMap();
const labelList& cellIDs = overlap.interpolationCells();
const scalarList& factor = overlap.cellInterpolationWeight();
Field<T> work(psi);
map.mapDistributeBase::distribute(work, UPstream::msgType()+1);
forAll(cellIDs, i)
{
const label celli = cellIDs[i];
const scalarList& w = wghts[celli];
const labelList& nbrs = stencil[celli];
const scalar f = factor[celli];
if (nbrs.size() == 0 && f != 0.0)
{
FatalErrorInFunction << "problem: cell:" << celli
<< " at:" << mesh.cellCentres()[celli]
<< " type:" << overlap.cellTypes()[celli]
<< " stencil:" << nbrs
<< " factor:" << f << exit(FatalError);
}
T s(pTraits<T>::zero);
forAll(nbrs, nbrI)
{
s += w[nbrI]*work[nbrs[nbrI]];
}
psi[celli] = (1.0-f)*psi[celli] + f*s;
}
}
template<class T>
void Foam::cellCellStencil::interpolate(const fvMesh& mesh, Field<T>& psi) const
{
const cellCellStencil& overlap = *this;
interpolate
(
psi,
mesh,
overlap,
overlap.cellInterpolationWeights()
);
}
template<class GeoField>
void Foam::cellCellStencil::interpolate(GeoField& psi) const
{
interpolate
(
psi.mesh(),
psi.primitiveFieldRef()
);
psi.correctBoundaryConditions();
}
template<class GeoField>
void Foam::cellCellStencil::interpolate
(
const fvMesh& mesh,
const wordHashSet& suppressed
) const
{
const cellCellStencil& overlap = *this;
auto flds(mesh.lookupClass<GeoField>());
for (auto fldPtr : flds)
{
const word& name = fldPtr->name();
if (!suppressed.found(baseName(name)))
{
if (debug)
{
Pout<< "cellCellStencil::interpolate: interpolating : "
<< name << endl;
}
auto& fld = const_cast<GeoField&>(*fldPtr);
interpolate
(
fld.primitiveFieldRef(),
mesh,
overlap,
overlap.cellInterpolationWeights()
);
}
else
{
if (debug)
{
Pout<< "cellCellStencil::interpolate: skipping : " << name
<< endl;
}
}
}
}
template<class Type>
Foam::tmp<Foam::volScalarField>
Foam::cellCellStencil::createField
@ -38,23 +158,20 @@ Foam::cellCellStencil::createField
const UList<Type>& psi
)
{
tmp<volScalarField> tfld
auto tfld = tmp<volScalarField>::New
(
new volScalarField
IOobject
(
IOobject
(
name,
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
name,
mesh.time().timeName(),
mesh,
dimensionedScalar(dimless, Zero),
zeroGradientFvPatchScalarField::typeName
)
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar(dimless, Zero),
zeroGradientFvPatchScalarField::typeName
);
volScalarField& fld = tfld.ref();
@ -66,4 +183,41 @@ Foam::cellCellStencil::createField
}
template<class GeoField, class SuppressBC>
void Foam::cellCellStencil::correctBoundaryConditions
(
GeoField& psi
)
{
// Version of GeoField::correctBoundaryConditions that exclude evaluation
// of oversetFvPatchFields
typename GeoField::Boundary& bfld = psi.boundaryFieldRef();
const label nReq = Pstream::nRequests();
forAll(bfld, patchi)
{
if (!isA<SuppressBC>(bfld[patchi]))
{
bfld[patchi].initEvaluate(Pstream::commsTypes::nonBlocking);
}
}
// Block for any outstanding requests
if (Pstream::parRun())
{
Pstream::waitRequests(nReq);
}
forAll(bfld, patchi)
{
if (!isA<SuppressBC>(bfld[patchi]))
{
bfld[patchi].evaluate(Pstream::commsTypes::nonBlocking);
}
}
}
// ************************************************************************* //

510
src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C
View File

@ -37,7 +37,8 @@ License
#include "oversetFvPatch.H"
#include "zeroGradientFvPatchFields.H"
#include "syncTools.H"
#include "dynamicOversetFvMesh.H"
#include "oversetFvMeshBase.H"
#include "oversetFvPatchField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -56,150 +57,6 @@ Foam::cellCellStencils::cellVolumeWeight::defaultOverlapTolerance_ = 1e-6;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::cellCellStencils::cellVolumeWeight::walkFront
(
const scalar layerRelax,
labelList& allCellTypes,
scalarField& allWeight
) const
{
// Current front
bitSet isFront(mesh_.nFaces());
// unused: bitSet doneCell(mesh_.nCells());
const fvBoundaryMesh& fvm = mesh_.boundary();
// 'overset' patches
forAll(fvm, patchI)
{
if (isA<oversetFvPatch>(fvm[patchI]))
{
//Pout<< "Storing faces on patch " << fvm[patchI].name() << endl;
forAll(fvm[patchI], i)
{
isFront.set(fvm[patchI].start()+i);
}
}
}
// Outside of 'hole' region
{
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
label ownType = allCellTypes[own[faceI]];
label neiType = allCellTypes[nei[faceI]];
if
(
(ownType == HOLE && neiType != HOLE)
|| (ownType != HOLE && neiType == HOLE)
)
{
//Pout<< "Front at face:" << faceI
// << " at:" << mesh_.faceCentres()[faceI] << endl;
isFront.set(faceI);
}
}
labelList nbrCellTypes;
syncTools::swapBoundaryCellList(mesh_, allCellTypes, nbrCellTypes);
for
(
label faceI = mesh_.nInternalFaces();
faceI < mesh_.nFaces();
faceI++
)
{
label ownType = allCellTypes[own[faceI]];
label neiType = nbrCellTypes[faceI-mesh_.nInternalFaces()];
if
(
(ownType == HOLE && neiType != HOLE)
|| (ownType != HOLE && neiType == HOLE)
)
{
//Pout<< "Front at coupled face:" << faceI
// << " at:" << mesh_.faceCentres()[faceI] << endl;
isFront.set(faceI);
}
}
}
// Current interpolation fraction
scalar fraction = 1.0;
while (fraction > SMALL && returnReduceOr(isFront.any()))
{
// Interpolate cells on front
Info<< "Front : fraction:" << fraction
<< " size:" << returnReduce(isFront.count(), sumOp<label>())
<< endl;
bitSet newIsFront(mesh_.nFaces());
forAll(isFront, faceI)
{
if (isFront.test(faceI))
{
label own = mesh_.faceOwner()[faceI];
if (allCellTypes[own] != HOLE)
{
if (allWeight[own] < fraction)
{
allWeight[own] = fraction;
//if (debug)
//{
// Pout<< " setting cell "
// << mesh_.cellCentres()[own]
// << " to " << fraction << endl;
//}
allCellTypes[own] = INTERPOLATED;
newIsFront.set(mesh_.cells()[own]);
}
}
if (mesh_.isInternalFace(faceI))
{
label nei = mesh_.faceNeighbour()[faceI];
if (allCellTypes[nei] != HOLE)
{
if (allWeight[nei] < fraction)
{
allWeight[nei] = fraction;
//if (debug)
//{
// Pout<< " setting cell "
// << mesh_.cellCentres()[nei]
// << " to " << fraction << endl;
//}
allCellTypes[nei] = INTERPOLATED;
newIsFront.set(mesh_.cells()[nei]);
}
}
}
}
}
syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
isFront.transfer(newIsFront);
fraction -= layerRelax;
}
}
void Foam::cellCellStencils::cellVolumeWeight::findHoles
(
const globalIndex& globalCells,
@ -213,13 +70,6 @@ void Foam::cellCellStencils::cellVolumeWeight::findHoles
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
// The input cellTypes will be
// - HOLE : cell part covered by other-mesh patch
// - INTERPOLATED : cell fully covered by other-mesh patch
// or next to 'overset' patch
// - CALCULATED : otherwise
//
// so we start a walk from our patches and any cell we cannot reach
// (because we walk is stopped by other-mesh patch) is a hole.
@ -503,6 +353,7 @@ void Foam::cellCellStencils::cellVolumeWeight::interpolatePatchTypes
{
// 'patch' overrides -1 and 'other'
result[cellI] = PATCH;
break;
}
else if (result[cellI] == -1)
{
@ -589,27 +440,26 @@ void Foam::cellCellStencils::cellVolumeWeight::combineCellTypes
// Patch-patch interaction... For now disable always
allCellTypes[cellI] = HOLE;
validDonors = false;
// Alternative is to look at the amount of overlap but this
// is not very robust
//if (allCellTypes[cellI] != HOLE)
//{
// scalar overlapVol = sum(weights[subCellI]);
// scalar v = mesh_.V()[cellI];
// if (overlapVol < (1.0-overlapTolerance_)*v)
// {
// //Pout<< "** Patch overlap:" << cellI
// // << " at:" << mesh_.cellCentres()[cellI] << endl;
// allCellTypes[cellI] = HOLE;
// validDonors = false;
// }
//}
// if (allCellTypes[cellI] != HOLE)
// {
// scalar overlapVol = sum(weights[subCellI]);
// scalar v = mesh_.V()[cellI];
// if (overlapVol < (1.0-overlapTolerance_)*v)
// {
// //Pout<< "** Patch overlap:" << cellI
// // << " at:" << mesh_.cellCentres()[cellI] << endl;
// allCellTypes[cellI] = HOLE;
// validDonors = false;
// }
// }
}
break;
case OVERSET:
{
validDonors = false;
validDonors = true;
}
break;
}
@ -685,6 +535,10 @@ Foam::cellCellStencils::cellVolumeWeight::cellVolumeWeight
mesh_,
dimensionedScalar(dimless, Zero),
zeroGradientFvPatchScalarField::typeName
),
allowInterpolatedDonors_
(
dict.getOrDefault("allowInterpolatedDonors", true)
)
{
// Protect local fields from interpolation
@ -701,6 +555,9 @@ Foam::cellCellStencils::cellVolumeWeight::cellVolumeWeight
// Read zoneID
this->zoneID();
overlapTolerance_ =
dict_.getOrDefault("overlapTolerance", defaultOverlapTolerance_);
// Read old-time cellTypes
IOobject io
(
@ -755,7 +612,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
}
Pstream::listCombineReduce(nCellsPerZone, plusEqOp<label>());
Info<< typeName << " : detected " << nZones
<< " mesh regions" << nl << endl;
@ -777,7 +633,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
Info<< decrIndent;
// Current best guess for cells. Includes best stencil. Weights should
// add up to volume.
labelList allCellTypes(mesh_.nCells(), CALCULATED);
@ -802,6 +657,15 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
markPatchCells(partMesh, partCellMap, allPatchTypes);
}
if ((debug&2) && (mesh_.time().outputTime()))
{
tmp<volScalarField> tfld
(
createField(mesh_, "allPatchTypes", allPatchTypes)
);
tfld().write();
}
labelList nCells(count(3, allPatchTypes));
Info<< nl
@ -815,7 +679,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
globalIndex globalCells(mesh_.nCells());
for (label srcI = 0; srcI < meshParts.size()-1; srcI++)
{
const fvMesh& srcMesh = meshParts[srcI].subMesh();
@ -837,7 +700,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
false // do not normalise
);
{
// Get tgt patch types on src mesh
labelList interpolatedTgtPatchTypes(srcMesh.nCells(), -1);
@ -863,6 +725,8 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
mapper.tgtMap()->distribute(tgtGlobalCells);
}
}
combineCellTypes
(
srcI,
@ -930,13 +794,20 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
}
if (debug)
if ((debug&2) && (mesh_.time().outputTime()))
{
tmp<volScalarField> tfld
(
createField(mesh_, "allCellTypes", allCellTypes)
);
tfld().write();
tmp<volScalarField> tdonors
(
createField(mesh_, "allDonorID", allDonorID)
);
tdonors().write();
}
@ -958,8 +829,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
}
else
{
//Pout<< "Holeing interpolated cell:" << cellI
// << " at:" << mesh_.cellCentres()[cellI] << endl;
allCellTypes[cellI] = HOLE;
allWeights[cellI].clear();
allStencil[cellI].clear();
@ -970,36 +839,97 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
}
}
/*
// Knock out cell with insufficient interpolation weights
forAll(allCellTypes, cellI)
{
if (allCellTypes[cellI] == INTERPOLATED)
{
scalar v = mesh_.V()[cellI];
scalar overlapVol = sum(allWeights[cellI]);
if (overlapVol < (1.0-overlapTolerance_)*v)
{
//Pout<< "Holeing cell:" << cellI
// << " at:" << mesh_.cellCentres()[cellI] << endl;
allCellTypes[cellI] = HOLE;
allWeights[cellI].clear();
allStencil[cellI].clear();
}
}
}
*/
if (debug)
if ((debug&2) && (mesh_.time().outputTime()))
{
tmp<volScalarField> tfld
(
createField(mesh_, "allCellTypes_patch", allCellTypes)
);
//tfld.ref().correctBoundaryConditions();
tfld().write();
labelList stencilSize(mesh_.nCells());
forAll(allStencil, celli)
{
stencilSize[celli] = allStencil[celli].size();
}
tmp<volScalarField> tfldStencil
(
createField(mesh_, "allStencil_patch", stencilSize)
);
tfldStencil().write();
}
// Mark unreachable bits
findHoles(globalCells, mesh_, zoneID, allStencil, allCellTypes);
if ((debug&2) && (mesh_.time().outputTime()))
{
tmp<volScalarField> tfld
(
createField(mesh_, "allCellTypes_hole", allCellTypes)
);
tfld().write();
labelList stencilSize(mesh_.nCells());
forAll(allStencil, celli)
{
stencilSize[celli] = allStencil[celli].size();
}
tmp<volScalarField> tfldStencil
(
createField(mesh_, "allStencil_hole", stencilSize)
);
tfldStencil().write();
}
// Add buffer interpolation layer around holes
scalarField allWeight(mesh_.nCells(), Zero);
labelListList compactStencil(allStencil);
List<Map<label>> compactStencilMap;
mapDistribute map(globalCells, compactStencil, compactStencilMap);
scalarList compactCellVol(mesh_.V());
map.distribute(compactCellVol);
walkFront
(
globalCells,
layerRelax,
allStencil,
allCellTypes,
allWeight,
compactCellVol,
compactStencil,
zoneID,
dict_.getOrDefault("holeLayers", 1),
dict_.getOrDefault("useLayer", -1)
);
if ((debug&2) && (mesh_.time().outputTime()))
{
tmp<volScalarField> tfld
(
createField(mesh_, "allCellTypes_front", allCellTypes)
);
tfld().write();
labelList stencilSize(mesh_.nCells());
forAll(allStencil, celli)
{
stencilSize[celli] = allStencil[celli].size();
}
tmp<volScalarField> tfldStencil
(
createField(mesh_, "allStencil_front", stencilSize)
);
tfldStencil().write();
}
// Check previous iteration cellTypes_ for any hole->calculated changes
// If so set the cell either to interpolated (if there are donors) or
// holes (if there are no donors). Note that any interpolated cell might
@ -1035,53 +965,29 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
}
// Mark unreachable bits
findHoles(globalCells, mesh_, zoneID, allStencil, allCellTypes);
labelList compactCellTypes(allCellTypes);
map.distribute(compactCellTypes);
if (debug)
{
tmp<volScalarField> tfld
(
createField(mesh_, "allCellTypes_hole", allCellTypes)
);
//tfld.ref().correctBoundaryConditions();
tfld().write();
}
// Add buffer interpolation layer around holes
scalarField allWeight(mesh_.nCells(), Zero);
walkFront(layerRelax, allCellTypes, allWeight);
if (debug)
{
tmp<volScalarField> tfld
(
createField(mesh_, "allCellTypes_front", allCellTypes)
);
//tfld.ref().correctBoundaryConditions();
tfld().write();
}
// Normalise weights, Clear storage
label nHoleDonors = 0;
forAll(allCellTypes, cellI)
{
if (allCellTypes[cellI] == INTERPOLATED)
{
if (allWeight[cellI] < SMALL || allStencil[cellI].size() == 0)
const labelList& slots = compactStencil[cellI];
forAll (slots, subCellI)
{
//Pout<< "Clearing cell:" << cellI
// << " at:" << mesh_.cellCentres()[cellI] << endl;
allWeights[cellI].clear();
allStencil[cellI].clear();
allWeight[cellI] = 0.0;
}
else
{
scalar s = sum(allWeights[cellI]);
forAll(allWeights[cellI], i)
if
(
compactCellTypes[slots[0]] == HOLE
||
(
!allowInterpolatedDonors_
&& compactCellTypes[slots[0]] == INTERPOLATED
)
)
{
allWeights[cellI][i] /= s;
allWeights[cellI][subCellI] = 0;
nHoleDonors++;
}
}
}
@ -1091,112 +997,52 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
allStencil[cellI].clear();
}
}
reduce(nHoleDonors, sumOp<label>());
// Normalize weights
forAll(allCellTypes, cellI)
{
if (allCellTypes[cellI] == INTERPOLATED)
{
const scalar s = sum(allWeights[cellI]);
if (s < SMALL)
{
allCellTypes[cellI] = POROUS;
allWeights[cellI].clear();
allStencil[cellI].clear();
}
else
{
forAll(allWeights[cellI], i)
{
allWeights[cellI][i] /= s;
}
}
}
}
// Write to volField for debugging
if (debug)
if ((debug&2) && (mesh_.time().outputTime()))
{
volScalarField patchTypes
(
IOobject
(
"patchTypes",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh_,
dimensionedScalar(dimless, Zero),
zeroGradientFvPatchScalarField::typeName
);
forAll(patchTypes.internalField(), cellI)
if ((debug&2) && (mesh_.time().outputTime()))
{
patchTypes[cellI] = allPatchTypes[cellI];
}
//patchTypes.correctBoundaryConditions();
dynamicOversetFvMesh::correctBoundaryConditions
<
volScalarField,
oversetFvPatchField<scalar>
>(patchTypes.boundaryFieldRef(), false);
patchTypes.write();
}
if (debug)
{
volScalarField volTypes
(
IOobject
tmp<volScalarField> tfld
(
"cellTypes",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh_,
dimensionedScalar(dimless, Zero),
zeroGradientFvPatchScalarField::typeName
);
forAll(volTypes.internalField(), cellI)
{
volTypes[cellI] = allCellTypes[cellI];
createField(mesh_, "allCellTypes_final", allCellTypes)
);
tfld().write();
}
//volTypes.correctBoundaryConditions();
dynamicOversetFvMesh::correctBoundaryConditions
<
volScalarField,
oversetFvPatchField<scalar>
>(volTypes.boundaryFieldRef(), false);
volTypes.write();
}
// // Check previous iteration cellTypes_ for any hole->calculated changes
// {
// label nCalculated = 0;
//
// forAll(cellTypes_, celli)
// {
// if (allCellTypes[celli] == CALCULATED && cellTypes_[celli] == HOLE)
// {
// if (allStencil[celli].size() == 0)
// {
// FatalErrorInFunction
// << "Cell:" << celli
// << " at:" << mesh_.cellCentres()[celli]
// << " zone:" << zoneID[celli]
// << " changed from hole to calculated"
// << " but there is no donor"
// << exit(FatalError);
// }
// else
// {
// allCellTypes[celli] = INTERPOLATED;
// nCalculated++;
// }
// }
// }
//
// if (debug)
// {
// Pout<< "Detected " << nCalculated << " cells changing from hole"
// << " to calculated. Changed these to interpolated"
// << endl;
// }
// }
cellTypes_.transfer(allCellTypes);
cellStencil_.transfer(allStencil);
cellInterpolationWeights_.transfer(allWeights);
cellInterpolationWeight_.transfer(allWeight);
//cellInterpolationWeight_.correctBoundaryConditions();
dynamicOversetFvMesh::correctBoundaryConditions
oversetFvMeshBase::correctBoundaryConditions
<
volScalarField,
oversetFvPatchField<scalar>
@ -1220,7 +1066,7 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
);
// Dump interpolation stencil
if (debug)
if ((debug&2) && (mesh_.time().outputTime()))
{
// Dump weight
cellInterpolationWeight_.instance() = mesh_.time().timeName();
@ -1247,15 +1093,16 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
forAll(slots, i)
{
const point& donorCc = cc[slots[i]];
const point& accCc = mesh_.cellCentres()[cellI];
str.writeLine(accCc, 0.1*accCc+0.9*donorCc);
if (cellInterpolationWeights_[cellI][slots[i]] > 0)
{
const point& donorCc = cc[slots[i]];
const point& accCc = mesh_.cellCentres()[cellI];
str.writeLine(accCc, 0.1*accCc+0.9*donorCc);
}
}
}
}
{
labelList nCells(count(3, cellTypes_));
Info<< "Overset analysis : nCells : "
@ -1267,7 +1114,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
<< decrIndent << endl;
}
// Tbd: detect if anything changed. Most likely it did!
return true;
}

10
src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2018 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -90,15 +90,11 @@ protected:
//- Amount of interpolation
volScalarField cellInterpolationWeight_;
//- Allow interpolared as donors
const bool allowInterpolatedDonors_;
// Protected Member Functions
void walkFront
(
const scalar layerRelax,
labelList& allCellTypes,
scalarField& allWeight
) const;
//- Find cells next to cells of type PATCH
void findHoles

911
src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
View File

File diff suppressed because it is too large Load Diff

55
src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017-2020 OpenCFD Ltd.
Copyright (C) 2017-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -69,6 +69,15 @@ class inverseDistance
:
public cellCellStencil
{
// Private Member Functions
//- No copy construct
inverseDistance(const inverseDistance&) = delete;
//- No copy assignment
void operator=(const inverseDistance&) = delete;
protected:
// Protected data
@ -76,6 +85,12 @@ protected:
//- Dictionary of motion control parameters
const dictionary dict_;
//- Allow holes as donors
const bool allowHoleDonors_;
//- Allow interpolared as donors
const bool allowInterpolatedDonors_;
//- Small perturbation vector for geometric tests
vector smallVec_;
@ -97,6 +112,13 @@ protected:
//- Amount of interpolation
volScalarField cellInterpolationWeight_;
//- Data for a set of interpolated/donor set
struct interpolatedDonorSet
{
scalar averVolRatio;
labelList cellTypes;
};
// Protected Member Functions
@ -217,6 +239,7 @@ protected:
const boolList& isBlockedFace,
labelList& cellRegion
) const;
autoPtr<globalIndex> compactedRegionSplit
(
const fvMesh& mesh,
@ -244,28 +267,18 @@ protected:
scalarField& fraction
) const;
//- Surround holes with layer(s) of interpolated cells
void walkFront
(
const scalar layerRelax,
const labelListList& allStencil,
labelList& allCellTypes,
scalarField& allWeight
) const;
//- Create stencil starting from the donor containing the acceptor
virtual void createStencil(const globalIndex&);
//- allowHoleDonors : allow donors of type HOLE (otherwise are filtered
//- out)
virtual void createStencil
(
const globalIndex&,
const bool allowHoleDonors
);
private:
// Private Member Functions
//- No copy construct
inverseDistance(const inverseDistance&) = delete;
//- No copy assignment
void operator=(const inverseDistance&) = delete;
//- Make holes next to live ones type SPECIAL and include in
//- interpolation stencil
void holeExtrapolationStencil(const globalIndex& globalCells);
public:

218
src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
View File

@ -37,6 +37,7 @@ License
#include "treeBoundBoxList.H"
#include "voxelMeshSearch.H"
#include "dynamicOversetFvMesh.H"
#include "waveMethod.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -337,17 +338,24 @@ void Foam::cellCellStencils::trackingInverseDistance::markDonors
const fvMesh& srcMesh = meshParts_[srcI].subMesh();
const labelList& srcCellMap = meshParts_[srcI].cellMap();
const voxelMeshSearch& meshSearch = meshSearches[srcI];
const fvMesh& tgtMesh = meshParts_[tgtI].subMesh();
const pointField& tgtCc = tgtMesh.cellCentres();
const labelList& tgtCellMap = meshParts_[tgtI].cellMap();
// 1. do processor-local src/tgt overlap
{
labelList tgtToSrcAddr;
waveMethod::calculate(tgtMesh, srcMesh, tgtToSrcAddr);
forAll(tgtCellMap, tgtCelli)
{
label srcCelli = meshSearch.findCell(tgtCc[tgtCelli]);
if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
const label srcCelli = tgtToSrcAddr[tgtCelli];
if
(
srcCelli != -1
&& allCellTypes[tgtCellMap[tgtCelli]] != HOLE
)
{
label celli = tgtCellMap[tgtCelli];
@ -444,8 +452,9 @@ void Foam::cellCellStencils::trackingInverseDistance::markDonors
labelList donors(samples.size(), -1);
forAll(samples, sampleI)
{
label srcCelli = meshSearch.findCell(samples[sampleI]);
if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
const point& sample = samples[sampleI];
label srcCelli = srcMesh.findCell(sample, polyMesh::CELL_TETS);
if (srcCelli != -1)
{
donors[sampleI] = globalCells_.toGlobal(srcCellMap[srcCelli]);
}
@ -736,7 +745,7 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
}
if (debug)
if (debug & 2)
{
forAll(patchParts, zonei)
{
@ -832,13 +841,19 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
DebugInfo<< FUNCTION_NAME << " : Determined holes and donor-acceptors"
<< endl;
if (debug)
if ((debug & 2) && (mesh_.time().outputTime()))
{
tmp<volScalarField> tfld
(
createField(mesh_, "allCellTypes", allCellTypes)
);
tfld().write();
tmp<volScalarField> tallDonorID
(
createField(mesh_, "allDonorID", allDonorID)
);
tallDonorID().write();
}
@ -867,15 +882,86 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
}
DebugInfo<< FUNCTION_NAME << " : Removed bad donors" << endl;
if (debug)
if ((debug & 2) && (mesh_.time().outputTime()))
{
tmp<volScalarField> tfld
(
createField(mesh_, "allCellTypes_patch", allCellTypes)
);
tfld().write();
tmp<volScalarField> tfldOld
(
createField(mesh_, "allCellTypes_old", cellTypes_)
);
tfldOld().write();
}
// Mark unreachable bits
findHoles(globalCells_, mesh_, zoneID, allStencil, allCellTypes);
DebugInfo<< FUNCTION_NAME << " : Flood-filled holes" << endl;
if ((debug & 2) && (mesh_.time().outputTime()))
{
tmp<volScalarField> tfld
(
createField(mesh_, "allCellTypes_hole", allCellTypes)
);
tfld().write();
labelList stencilSize(mesh_.nCells());
forAll(allStencil, celli)
{
stencilSize[celli] = allStencil[celli].size();
}
tmp<volScalarField> tfldsten
(
createField(mesh_, "allStencil_hole", stencilSize)
);
tfldsten().write();
}
// Update allStencil with new fill HOLES
forAll(allCellTypes, celli)
{
if (allCellTypes[celli] == HOLE && allStencil[celli].size())
{
allStencil[celli].clear();
}
}
// Add buffer interpolation layer(s) around holes
scalarField allWeight(mesh_.nCells(), Zero);
labelListList compactStencil(allStencil);
List<Map<label>> compactStencilMap;
mapDistribute map(globalCells, compactStencil, compactStencilMap);
scalarList compactCellVol(mesh_.V());
map.distribute(compactCellVol);
walkFront
(
globalCells,
layerRelax,
allStencil,
allCellTypes,
allWeight,
compactCellVol,
compactStencil,
zoneID,
dict_.getOrDefault("holeLayers", 1),
dict_.getOrDefault("useLayer", -1)
);
if ((debug & 2) && (mesh_.time().outputTime()))
{
tmp<volScalarField> tfld
(
createField(mesh_, "allCellTypes_front", allCellTypes)
);
tfld().write();
}
// Check previous iteration cellTypes_ for any hole->calculated changes
// If so set the cell either to interpolated (if there are donors) or
@ -911,48 +997,47 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
}
// Mark unreachable bits
findHoles(globalCells_, mesh_, zoneID, allStencil, allCellTypes);
DebugInfo<< FUNCTION_NAME << " : Flood-filled holes" << endl;
if (debug)
{
tmp<volScalarField> tfld
(
createField(mesh_, "allCellTypes_hole", allCellTypes)
);
tfld().write();
}
// Add buffer interpolation layer(s) around holes
scalarField allWeight(mesh_.nCells(), Zero);
walkFront(layerRelax, allStencil, allCellTypes, allWeight);
DebugInfo<< FUNCTION_NAME << " : Implemented layer relaxation" << endl;
if (debug)
{
tmp<volScalarField> tfld
(
createField(mesh_, "allCellTypes_front", allCellTypes)
);
tfld().write();
}
// Convert cell-cell addressing to stencil in compact notation
cellTypes_.transfer(allCellTypes);
cellStencil_.setSize(mesh_.nCells());
cellInterpolationWeights_.setSize(mesh_.nCells());
DynamicList<label> interpolationCells;
labelList compactCellTypes(cellTypes_);
map.distribute(compactCellTypes);
label nHoleDonors = 0;
forAll(cellTypes_, celli)
{
if (cellTypes_[celli] == INTERPOLATED)
{
cellStencil_[celli].transfer(allStencil[celli]);
cellInterpolationWeights_[celli].setSize(1);
cellInterpolationWeights_[celli][0] = 1.0;
interpolationCells.append(celli);
const labelList& slots = compactStencil[celli];
if (slots.size())
{
if
(
compactCellTypes[slots[0]] == HOLE
||
(
!allowInterpolatedDonors_
&& compactCellTypes[slots[0]] == INTERPOLATED
)
)
{
cellTypes_[celli] = POROUS;
cellStencil_[celli].clear();
cellInterpolationWeights_[celli].clear();
nHoleDonors++;
}
else
{
cellStencil_[celli].transfer(allStencil[celli]);
cellInterpolationWeights_[celli].setSize(1);
cellInterpolationWeights_[celli][0] = 1.0;
interpolationCells.append(celli);
}
}
}
else
{
@ -960,6 +1045,11 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
cellInterpolationWeights_[celli].clear();
}
}
reduce(nHoleDonors, sumOp<label>());
DebugInfo<< FUNCTION_NAME << "nHole Donors : " << nHoleDonors << endl;
interpolationCells_.transfer(interpolationCells);
List<Map<label>> compactMap;
@ -973,16 +1063,18 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
)
);
cellInterpolationWeight_.transfer(allWeight);
//cellInterpolationWeight_.correctBoundaryConditions();
dynamicOversetFvMesh::correctBoundaryConditions
oversetFvMeshBase::correctBoundaryConditions
<
volScalarField,
oversetFvPatchField<scalar>
>(cellInterpolationWeight_.boundaryFieldRef(), false);
if (debug & 2)
if ((debug & 2) && (mesh_.time().outputTime()))
{
// Dump mesh
mesh_.time().write();
// Dump stencil
mkDir(mesh_.time().timePath());
OBJstream str(mesh_.time().timePath()/"injectionStencil.obj");
@ -1010,45 +1102,23 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
// Extend stencil to get inverse distance weighted neighbours
createStencil(globalCells);
createStencil(globalCells, allowHoleDonors_);
DebugInfo<< FUNCTION_NAME << " : Extended stencil" << endl;
// Optional: convert hole cells next to non-hole cells into
// interpolate-from-neighbours (of cell type SPECIAL)
if (allowHoleDonors_)
{
holeExtrapolationStencil(globalCells_);
}
if (debug & 2)
if ((debug & 2) && (mesh_.time().outputTime()))
{
// Dump weight
cellInterpolationWeight_.instance() = mesh_.time().timeName();
cellInterpolationWeight_.write();
// Dump cell types
volScalarField volTypes
(
IOobject
(
"cellTypes",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh_,
dimensionedScalar(dimless, Zero),
zeroGradientFvPatchScalarField::typeName
);
forAll(volTypes.internalField(), celli)
{
volTypes[celli] = cellTypes_[celli];
}
//volTypes.correctBoundaryConditions();
dynamicOversetFvMesh::correctBoundaryConditions
<
volScalarField,
oversetFvPatchField<scalar>
>(volTypes.boundaryFieldRef(), false);
volTypes.write();
// Dump stencil
mkDir(mesh_.time().timePath());
OBJstream str(mesh_.time().timePath()/"stencil.obj");

44
src/overset/include/correctPhiFaceMask.H Normal file
View File

@ -0,0 +1,44 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
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/>.
Description
Correct phi on new faces C-I faces
\*---------------------------------------------------------------------------*/
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
// Zero Uf on old faceMask (H-I)
Uf *= faceMaskOld;
// Update Uf and phi on new C-I faces
Uf += (1-faceMaskOld)*fvc::interpolate(U);
phi = mesh.Sf() & Uf;
phi *= faceMask;
// ************************************************************************* //

55
src/overset/include/correctRhoPhiFaceMask.H Normal file
View File

@ -0,0 +1,55 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
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/>.
Description
Correct phi on new faces C-I faces
\*---------------------------------------------------------------------------*/
const surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
// Zero Uf on old faceMask (H-I)
rhoUf() *= faceMaskOld;
surfaceVectorField rhoUfint(fvc::interpolate(rho*U));
// Update Uf and phi on new C-I faces
rhoUf() += (1-faceMaskOld)*rhoUfint;
// Update Uf boundary
forAll(rhoUf().boundaryField(), patchI)
{
rhoUf().boundaryFieldRef()[patchI] = rhoUfint.boundaryField()[patchI];
}
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & rhoUf();
phi *= faceMask;
// ************************************************************************* //

7
src/overset/include/createCellMask.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 OpenCFD Ltd.
Copyright (C) 2016-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -50,6 +50,11 @@ volScalarField cellMask
zeroGradientFvPatchScalarField::typeName
);
surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
#include "setCellMask.H"
// ************************************************************************* //

9
src/overset/include/createControls.H Normal file
View File

@ -0,0 +1,9 @@
bool adjustFringe
(
pimple.dict().getOrDefault("oversetAdjustPhi", false)
);
bool oversetPatchErrOutput
(
pimple.dict().getOrDefault("oversetPatchErrOutput", false)
);

4
src/overset/include/readOversetDyMControls.H Normal file
View File

@ -0,0 +1,4 @@
adjustFringe = pimple.dict().getOrDefault("oversetAdjustPhi", false);
label zoneIdMass = pimple.dict().getOrDefault("adjustZone", -1);
oversetPatchErrOutput =
pimple.dict().getOrDefault("oversetPatchErrOutput", false);

11
src/overset/include/setCellMask.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 OpenCFD Ltd.
Copyright (C) 2016-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -35,12 +35,19 @@ Description
cellMask.primitiveFieldRef() = 1.0;
forAll(cellMask, cellI)
{
if (cellTypes[cellI] == cellCellStencil::HOLE)
if
(
cellTypes[cellI] == cellCellStencil::HOLE
|| cellTypes[cellI] == cellCellStencil::SPECIAL
)
{
cellMask[cellI] = 0.0;
}
}
cellMask.correctBoundaryConditions();
faceMask = localMin<scalar>(mesh).interpolate(cellMask);
}
// ************************************************************************* //

289
src/overset/oversetAdjustPhi/oversetAdjustPhi.C
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017-2018 OpenCFD Ltd.
Copyright (C) 2017-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -37,64 +37,62 @@ License
bool Foam::oversetAdjustPhi
(
surfaceScalarField& phi,
const volVectorField& U
const volVectorField& U,
const label zoneId
)
{
const fvMesh& mesh = U.mesh();
const cellCellStencilObject& overlap = Stencil::New(mesh);
const labelList& cellTypes = overlap.cellTypes();
const labelList& zoneID = overlap.zoneID();
label nZones = gMax(zoneID)+1;
// Pass1: accumulate all fluxes, calculate correction factor
scalarField massIn(nZones, Zero);
scalarField adjustableMassOut(nZones, Zero);
surfaceScalarField::Boundary& bphi =
phi.boundaryFieldRef();
scalar massIn = 0;
scalar massOut = 0;
surfaceScalarField::Boundary& bphi = phi.boundaryFieldRef();
// Check all faces on the outside of interpolated cells
const labelUList& own = mesh.owner();
const labelUList& nei = mesh.neighbour();
forAll(own, facei)
{
forAll(own, facei)
const label zonei = zoneID[own[facei]];
const label ownType = cellTypes[own[facei]];
const label neiType = cellTypes[nei[facei]];
const bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
const bool neiCalc =
(ownType == cellCellStencil::INTERPOLATED)
&& (neiType == cellCellStencil::CALCULATED);
const bool ownOrCalc = (ownCalc || neiCalc);
if
(
(ownOrCalc && (zonei == zoneId))
|| (ownOrCalc && (zoneId == -1))
)
{
label zonei = zoneID[own[facei]]; // note:own and nei in same zone
// Calculate flux w.r.t. calculated cell
scalar flux = phi[facei];
label ownType = cellTypes[own[facei]];
label neiType = cellTypes[nei[facei]];
bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
bool neiCalc =
(ownType == cellCellStencil::INTERPOLATED)
&& (neiType == cellCellStencil::CALCULATED);
if (ownCalc || neiCalc)
if (ownCalc)
{
// Calculate flux w.r.t. calculated cell
scalar flux = phi[facei];
if (neiCalc)
{
flux = -flux;
}
flux = -flux;
}
if (flux < 0.0)
{
massIn[zonei] -= flux;
}
else
{
adjustableMassOut[zonei] += flux;
}
if (flux < 0.0)
{
massIn -= flux;
}
else
{
massOut += flux;
}
}
}
@ -103,131 +101,95 @@ bool Foam::oversetAdjustPhi
// Check all coupled faces on the outside of interpolated cells
labelList neiCellTypes;
syncTools::swapBoundaryCellList(mesh, cellTypes, neiCellTypes);
forAll(bphi, patchi)
{
forAll(bphi, patchi)
const fvPatchVectorField& Up = U.boundaryField()[patchi];
const fvsPatchScalarField& phip = bphi[patchi];
const labelUList& fc = Up.patch().faceCells();
const label start = Up.patch().start();
forAll(fc, i)
{
const fvPatchVectorField& Up = U.boundaryField()[patchi];
const fvsPatchScalarField& phip = bphi[patchi];
const labelUList& fc = Up.patch().faceCells();
const label facei = start + i;
const label celli = fc[i];
const label ownType = cellTypes[celli];
const label neiType = neiCellTypes[facei - mesh.nInternalFaces()];
label start = Up.patch().start();
const label zonei = zoneID[celli];
forAll(fc, i)
const bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
if (ownCalc && (zonei == zoneId))
{
label facei = start+i;
label celli = fc[i];
label ownType = cellTypes[celli];
label neiType = neiCellTypes[facei-mesh.nInternalFaces()];
// Calculate flux w.r.t. calculated cell
scalar flux = phip[i];
bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
if (ownCalc)
if (flux < 0.0)
{
// Calculate flux w.r.t. calculated cell
scalar flux = phip[i];
if (flux < 0.0)
{
massIn[zoneID[celli]] -= flux;
}
else
{
adjustableMassOut[zoneID[celli]] += flux;
}
massIn -= flux;
}
else
{
massOut += flux;
}
}
}
}
reduce(massIn, sumOp<scalar>());
reduce(massOut, sumOp<scalar>());
// Calculate the total flux in the domain, used for normalisation
scalar totalFlux = VSMALL + sum(mag(phi)).value();
const scalar massCorr = massIn/(massOut + SMALL);
forAll(massIn, zonei)
if (fv::debug)
{
reduce(massIn[zonei], sumOp<scalar>());
reduce(adjustableMassOut[zonei], sumOp<scalar>());
Info<< "Zone : " << zoneId << nl
<< "mass outflow : " << massOut << nl
<< "mass inflow : " << massIn << nl
<< "correction factor : " << massCorr << endl;
}
scalarField massCorr(nZones, 1.0);
forAll(massIn, zonei)
{
scalar magAdjustableMassOut = mag(adjustableMassOut[zonei]);
if
(
magAdjustableMassOut > VSMALL
&& magAdjustableMassOut/totalFlux > SMALL
)
{
massCorr[zonei] = massIn[zonei]/adjustableMassOut[zonei];
}
else if (mag(massIn[zonei])/totalFlux > 1e-8)
{
WarningInFunction
<< "Continuity error cannot be removed by adjusting the"
" flow at fringe faces.\n Please check the cell types"
<< " from the overset analysis."
<< nl
<< "Zone : " << zonei << nl
<< "Total flux : " << totalFlux << nl
<< "Specified mass inflow : " << massIn[zonei] << nl
<< "Adjustable mass outflow : " << adjustableMassOut[zonei]
<< nl << endl;
}
if (fv::debug)
{
Info<< "Zone : " << zonei << nl
<< "Total flux : " << totalFlux << nl
<< "Specified mass inflow : " << massIn[zonei] << nl
<< "Adjustable mass outflow : " << adjustableMassOut[zonei]
<< nl
<< "Correction factor : " << massCorr[zonei] << nl
<< endl;
}
}
// Pass2: adjust fluxes
forAll(own, facei)
{
label zonei = zoneID[own[facei]]; // note:own and nei in same zone
const label zonei = zoneID[own[facei]]; // own and nei in same zone
label ownType = cellTypes[own[facei]];
label neiType = cellTypes[nei[facei]];
const label ownType = cellTypes[own[facei]];
const label neiType = cellTypes[nei[facei]];
bool ownCalc =
const bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
bool neiCalc =
const bool neiCalc =
(ownType == cellCellStencil::INTERPOLATED)
&& (neiType == cellCellStencil::CALCULATED);
if (ownCalc || neiCalc)
const bool ownOrCalc = (ownCalc || neiCalc);
if
(
(ownOrCalc && (zonei == zoneId)) || (ownOrCalc && (zoneId == -1))
)
{
// Calculate flux w.r.t. calculated cell
scalar flux = phi[facei];
if (neiCalc)
if (ownCalc)
{
flux = -flux;
}
if (flux < 0.0)
{
phi[facei] /= Foam::sqrt(massCorr[zonei]);
phi[facei] /= Foam::sqrt(massCorr);
}
else
{
phi[facei] *= Foam::sqrt(massCorr[zonei]);
phi[facei] *= Foam::sqrt(massCorr);
}
}
}
@ -238,28 +200,31 @@ bool Foam::oversetAdjustPhi
fvsPatchScalarField& phip = bphi[patchi];
const labelUList& fc = Up.patch().faceCells();
label start = Up.patch().start();
const label start = Up.patch().start();
forAll(fc, i)
{
label facei = start+i;
label celli = fc[i];
label zonei = zoneID[celli]; // note:own and nei in same zone
label ownType = cellTypes[celli];
label neiType = neiCellTypes[facei-mesh.nInternalFaces()];
const label facei = start + i;
const label celli = fc[i];
const label zonei = zoneID[celli]; // note:own and nei in same zone
const label ownType = cellTypes[celli];
const label neiType = neiCellTypes[facei - mesh.nInternalFaces()];
bool ownCalc =
const bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
bool neiCalc =
const bool neiCalc =
(ownType == cellCellStencil::INTERPOLATED)
&& (neiType == cellCellStencil::CALCULATED);
if (ownCalc || neiCalc)
const bool ownOrCalc = (ownCalc || neiCalc);
if (ownOrCalc && (zonei == zoneId))
{
// Calculate flux w.r.t. calculated cell
scalar flux = phip[i];
if (neiCalc)
{
flux = -flux;
@ -267,18 +232,64 @@ bool Foam::oversetAdjustPhi
if (flux < 0.0)
{
phip[i] /= Foam::sqrt(massCorr[zonei]);
// phip[i] /= Foam::sqrt(massCorr[zonei]);
}
else
{
phip[i] *= Foam::sqrt(massCorr[zonei]);
phip[i] *= massCorr;
}
}
}
}
return sum(mag(massIn))/totalFlux < SMALL
&& sum(mag(adjustableMassOut))/totalFlux < SMALL;
// Check correction
if (fv::debug)
{
scalar massOutCheck = 0;
scalar massInCheck = 0;
forAll(own, facei)
{
const label zonei = zoneID[own[facei]];
const label ownType = cellTypes[own[facei]];
const label neiType = cellTypes[nei[facei]];
const bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
const bool neiCalc =
(ownType == cellCellStencil::INTERPOLATED)
&& (neiType == cellCellStencil::CALCULATED);
const bool ownOrCalc = (ownCalc || neiCalc);
if
(
(ownOrCalc && (zonei == zoneId))||(ownOrCalc && (zoneId == -1))
)
{
scalar flux = phi[facei];
if (ownCalc)
{
flux = -flux;
}
if (flux < 0.0)
{
massInCheck -= flux;
}
else
{
massOutCheck += flux;
}
}
}
}
return true;
}

5
src/overset/oversetAdjustPhi/oversetAdjustPhi.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenCFD Ltd.
Copyright (C) 2017-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -51,7 +51,8 @@ namespace Foam
bool oversetAdjustPhi
(
surfaceScalarField& phi,
const volVectorField& U
const volVectorField& U,
const label zoneId = -1
);
} // End namespace Foam

85
src/overset/oversetCoupledPolyPatch/oversetGAMGInterface/oversetGAMGInterface.C Normal file
View File

@ -0,0 +1,85 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "oversetGAMGInterface.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(oversetGAMGInterface, 0);
addToRunTimeSelectionTable
(
GAMGInterface,
oversetGAMGInterface,
lduInterface
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::oversetGAMGInterface::oversetGAMGInterface
(
const label index,
const lduInterfacePtrsList& coarseInterfaces,
const lduInterface& fineInterface,
const labelField& localRestrictAddressing,
const labelField& neighbourRestrictAddressing,
const label fineLevelIndex,
const label coarseComm
)
:
GAMGInterface
(
index,
coarseInterfaces
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::oversetGAMGInterface::~oversetGAMGInterface()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::labelField> Foam::oversetGAMGInterface::internalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& iF
) const
{
return tmp<labelField>::New(iF);
}
// ************************************************************************* //

122
src/overset/oversetCoupledPolyPatch/oversetGAMGInterface/oversetGAMGInterface.H Normal file
View File

@ -0,0 +1,122 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenFOAM Foundation
-------------------------------------------------------------------------------
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::oversetGAMGInterface
Description
GAMG agglomerated cyclic AMI interface.
SourceFiles
oversetGAMGInterface.C
\*---------------------------------------------------------------------------*/
#ifndef oversetGAMGInterface_H
#define oversetGAMGInterface_H
#include "GAMGInterface.H"
#include "oversetLduInterface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class oversetGAMGInterface Declaration
\*---------------------------------------------------------------------------*/
class oversetGAMGInterface
:
public GAMGInterface,
virtual public oversetLduInterface
{
// Private Member Functions
//- No copy construct
oversetGAMGInterface(const oversetGAMGInterface&) = delete;
//- No copy assignment
void operator=(const oversetGAMGInterface&) = delete;
public:
//- Runtime type information
TypeName("overset");
// Constructors
//- Construct from fine level interface,
//- local and neighbour restrict addressing
oversetGAMGInterface
(
const label index,
const lduInterfacePtrsList& coarseInterfaces,
const lduInterface& fineInterface,
const labelField& restrictAddressing,
const labelField& neighbourRestrictAddressing,
const label fineLevelIndex,
const label coarseComm
);
//- Destructor
virtual ~oversetGAMGInterface();
// Member Functions
// Interface transfer functions
//- Transfer and return internal field adjacent to the interface
virtual tmp<labelField> internalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& iF
) const;
// I/O
//- Write to stream
virtual void write(Ostream&) const
{
NotImplemented;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

97
src/overset/oversetCoupledPolyPatch/oversetGAMGInterfaceField/oversetGAMGInterfaceField.C Normal file
View File

@ -0,0 +1,97 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "oversetGAMGInterfaceField.H"
#include "addToRunTimeSelectionTable.H"
#include "lduMatrix.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(oversetGAMGInterfaceField, 0);
addToRunTimeSelectionTable
(
GAMGInterfaceField,
oversetGAMGInterfaceField,
lduInterface
);
addToRunTimeSelectionTable
(
GAMGInterfaceField,
oversetGAMGInterfaceField,
lduInterfaceField
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::oversetGAMGInterfaceField::oversetGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const lduInterfaceField& fineInterface
)
:
GAMGInterfaceField(GAMGCp, fineInterface)
{}
Foam::oversetGAMGInterfaceField::oversetGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
)
:
GAMGInterfaceField(GAMGCp, doTransform, rank)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::oversetGAMGInterfaceField::~oversetGAMGInterfaceField()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::oversetGAMGInterfaceField::updateInterfaceMatrix
(
solveScalarField& result,
const bool add,
const lduAddressing& lduAddr,
const label patchId,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
) const
{}
// ************************************************************************* //

125
src/overset/oversetCoupledPolyPatch/oversetGAMGInterfaceField/oversetGAMGInterfaceField.H Normal file
View File

@ -0,0 +1,125 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
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::oversetGAMGInterfaceField
Description
GAMG agglomerated
SourceFiles
oversetGAMGInterfaceField.C
\*---------------------------------------------------------------------------*/
#ifndef oversetGAMGInterfaceField_H
#define oversetGAMGInterfaceField_H
#include "GAMGInterfaceField.H"
#include "oversetGAMGInterface.H"
#include "oversetLduInterfaceField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class oversetGAMGInterfaceField Declaration
\*---------------------------------------------------------------------------*/
class oversetGAMGInterfaceField
:
public GAMGInterfaceField,
virtual public oversetLduInterfaceField
{
// Private Member Functions
//- No copy construct
oversetGAMGInterfaceField
(
const oversetGAMGInterfaceField&
) = delete;
//- No copy assignment
void operator=(const oversetGAMGInterfaceField&) = delete;
public:
//- Runtime type information
TypeName("overset");
// Constructors
//- Construct from GAMG interface and fine level interface field
oversetGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const lduInterfaceField& fineInterfaceField
);
//- Construct from GAMG interface and fine level interface field
oversetGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
);
//- Destructor
virtual ~oversetGAMGInterfaceField();
// Member Functions
// Interface matrix update
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
solveScalarField& result,
const bool add,
const lduAddressing& lduAddr,
const label patchId,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

44
src/overset/oversetCoupledPolyPatch/oversetLduInterface/oversetLduInterface.C Normal file
View File

@ -0,0 +1,44 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "oversetLduInterface.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(oversetLduInterface, 0);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::oversetLduInterface::~oversetLduInterface()
{}
// ************************************************************************* //

79
src/overset/oversetCoupledPolyPatch/oversetLduInterface/oversetLduInterface.H Normal file
View File

@ -0,0 +1,79 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
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::oversetLduInterface
Description
An abstract base class for overset coupled interfaces.
SourceFiles
oversetLduInterface.C
\*---------------------------------------------------------------------------*/
#ifndef oversetLduInterface_H
#define oversetLduInterface_H
#include "typeInfo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class oversetLduInterface Declaration
\*---------------------------------------------------------------------------*/
class oversetLduInterface
{
public:
//- Runtime type information
TypeName("oversetLduInterface");
// Constructors
//- Construct null
oversetLduInterface()
{}
//- Destructor
virtual ~oversetLduInterface();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

45
src/overset/oversetCoupledPolyPatch/oversetLduInterfaceField/oversetLduInterfaceField.C Normal file
View File

@ -0,0 +1,45 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "oversetLduInterfaceField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(oversetLduInterfaceField, 0);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::oversetLduInterfaceField::~oversetLduInterfaceField()
{}
// ************************************************************************* //

80
src/overset/oversetCoupledPolyPatch/oversetLduInterfaceField/oversetLduInterfaceField.H Normal file
View File

@ -0,0 +1,80 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
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::oversetLduInterfaceField
Description
Abstract base class for overset coupled interface fields.
SourceFiles
oversetLduInterfaceField.C
\*---------------------------------------------------------------------------*/
#ifndef oversetLduInterfaceField_H
#define oversetLduInterfaceField_H
#include "primitiveFieldsFwd.H"
#include "typeInfo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class oversetLduInterfaceField Declaration
\*---------------------------------------------------------------------------*/
class oversetLduInterfaceField
{
public:
//- Runtime type information
TypeName("oversetLduInterfaceField");
// Constructors
//- Construct given coupled patch
oversetLduInterfaceField()
{}
//- Destructor
virtual ~oversetLduInterfaceField();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

86
src/overset/oversetFvMesh/dynamicOversetFvMesh/dynamicOversetFvMesh.C Normal file
View File

@ -0,0 +1,86 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify i
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 "dynamicOversetFvMesh.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(dynamicOversetFvMesh, 0);
addToRunTimeSelectionTable(dynamicFvMesh, dynamicOversetFvMesh, IOobject);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::dynamicOversetFvMesh::dynamicOversetFvMesh
(
const IOobject& io,
const bool doInit
)
:
dynamicMotionSolverListFvMesh(io, doInit),
oversetFvMeshBase(static_cast<const fvMesh&>(*this))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::dynamicOversetFvMesh::~dynamicOversetFvMesh()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::dynamicOversetFvMesh::update()
{
if (!dynamicMotionSolverListFvMesh::update())
{
return false;
}
oversetFvMeshBase::update();
return true;
}
bool Foam::dynamicOversetFvMesh::writeObject
(
IOstreamOption streamOpt,
const bool valid
) const
{
bool ok = dynamicMotionSolverListFvMesh::writeObject(streamOpt, valid);
ok = oversetFvMeshBase::writeObject(streamOpt, valid) && ok;
return ok;
}
// ************************************************************************* //

226
src/overset/oversetFvMesh/dynamicOversetFvMesh/dynamicOversetFvMesh.H Normal file
View File

@ -0,0 +1,226 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
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::dynamicOversetFvMesh
Description
dynamicFvMesh with support for overset meshes.
SourceFiles
dynamicOversetFvMesh.C
\*---------------------------------------------------------------------------*/
#ifndef dynamicOversetFvMesh_H
#define dynamicOversetFvMesh_H
#include "dynamicMotionSolverListFvMesh.H"
#include "oversetFvMeshBase.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class dynamicOversetFvMesh Declaration
\*---------------------------------------------------------------------------*/
class dynamicOversetFvMesh
:
public dynamicMotionSolverListFvMesh,
public oversetFvMeshBase
{
// Private Member Functions
//- No copy construct
dynamicOversetFvMesh(const dynamicOversetFvMesh&) = delete;
//- No copy assignment
void operator=(const dynamicOversetFvMesh&) = delete;
public:
//- Runtime type information
TypeName("dynamicOversetFvMesh");
// Constructors
//- Construct from IOobject
dynamicOversetFvMesh(const IOobject& io, const bool doInit=true);
//- Destructor
virtual ~dynamicOversetFvMesh();
// Member Functions
//- Override ldu addressing
virtual const lduAddressing& lduAddr() const
{
return oversetFvMeshBase::lduAddr();
}
//- Override ldu addressing
virtual lduInterfacePtrsList interfaces() const
{
return oversetFvMeshBase::interfaces();
}
// Overset
// Explicit interpolation
virtual void interpolate(scalarField& psi) const
{
Stencil::New(*this).interpolate<scalar>(*this, psi);
}
virtual void interpolate(vectorField& psi) const
{
Stencil::New(*this).interpolate<vector>(*this, psi);
}
virtual void interpolate(sphericalTensorField& psi) const
{
Stencil::New
(
*this
).interpolate<sphericalTensor>(*this, psi);
}
virtual void interpolate(symmTensorField& psi) const
{
Stencil::New(*this).interpolate<symmTensor>(*this, psi);
}
virtual void interpolate(tensorField& psi) const
{
Stencil::New(*this).interpolate<tensor>(*this, psi);
}
virtual void interpolate(volScalarField& psi) const
{
Stencil::New(*this).interpolate<volScalarField>(psi);
}
virtual void interpolate(volVectorField& psi) const
{
Stencil::New(*this).interpolate<volVectorField>(psi);
}
virtual void interpolate(volSphericalTensorField& psi) const
{
Stencil::New(*this).interpolate
<
volSphericalTensorField
>(psi);
}
virtual void interpolate(volSymmTensorField& psi) const
{
Stencil::New(*this).interpolate<volSymmTensorField>(psi);
}
virtual void interpolate(volTensorField& psi) const
{
Stencil::New(*this).interpolate<volTensorField>(psi);
}
// Implicit interpolation (matrix manipulation)
//- Solve returning the solution statistics given convergence
//- tolerance. Use the given solver controls
virtual SolverPerformance<scalar> solve
(
fvMatrix<scalar>& m,
const dictionary& dict
) const
{
return oversetFvMeshBase::solveOverset<scalar>(m, dict);
}
//- Solve returning the solution statistics given convergence
//- tolerance. Use the given solver controls
virtual SolverPerformance<vector> solve
(
fvMatrix<vector>& m,
const dictionary& dict
) const
{
return oversetFvMeshBase::solveOverset<vector>(m, dict);
}
//- Solve returning the solution statistics given convergence
//- tolerance. Use the given solver controls
virtual SolverPerformance<symmTensor> solve
(
fvMatrix<symmTensor>& m,
const dictionary& dict
) const
{
return oversetFvMeshBase::solveOverset<symmTensor>
(m, dict);
}
//- Solve returning the solution statistics given convergence
//- tolerance. Use the given solver controls
virtual SolverPerformance<tensor> solve
(
fvMatrix<tensor>& m,
const dictionary& dict
) const
{
return oversetFvMeshBase::solveOverset<tensor>(m, dict);
}
//- Update the mesh for both mesh motion and topology change
virtual bool update();
//- Write using given format, version and compression
virtual bool writeObject
(
IOstreamOption streamOpt,
const bool valid
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

200
src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.C → src/overset/oversetFvMesh/oversetFvMeshBase.C
View File

@ -25,28 +25,27 @@ License
\*---------------------------------------------------------------------------*/
#include "dynamicOversetFvMesh.H"
#include "oversetFvMeshBase.H"
#include "addToRunTimeSelectionTable.H"
#include "cellCellStencilObject.H"
#include "zeroGradientFvPatchFields.H"
#include "lduPrimitiveProcessorInterface.H"
#include "globalIndex.H"
#include "GAMGAgglomeration.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(dynamicOversetFvMesh, 0);
addToRunTimeSelectionTable(dynamicFvMesh, dynamicOversetFvMesh, IOobject);
addToRunTimeSelectionTable(dynamicFvMesh, dynamicOversetFvMesh, doInit);
defineTypeNameAndDebug(oversetFvMeshBase, 0);
}
// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
bool Foam::dynamicOversetFvMesh::updateAddressing() const
bool Foam::oversetFvMeshBase::updateAddressing() const
{
const cellCellStencilObject& overlap = Stencil::New(*this);
const cellCellStencilObject& overlap = Stencil::New(mesh_);
// The (processor-local part of the) stencil determines the local
// faces to add to the matrix. tbd: parallel
@ -54,7 +53,7 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
// Get the base addressing
//const lduAddressing& baseAddr = dynamicMotionSolverFvMesh::lduAddr();
const lduAddressing& baseAddr = dynamicFvMesh::lduAddr();
const lduAddressing& baseAddr = mesh_.fvMesh::lduAddr();
// Add to the base addressing
labelList lowerAddr;
@ -89,13 +88,12 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
if (debug)
{
Pout<< "dynamicOversetFvMesh::update() : extended addressing from"
Pout<< "oversetFvMeshBase::update() : extended addressing from"
<< " nFaces:" << baseAddr.lowerAddr().size()
<< " to nFaces:" << lowerAddr.size()
<< " nExtraFaces:" << nExtraFaces << endl;
}
// Now for the tricky bits. We want to hand out processor faces according
// to the localFaceCells/remoteFaceCells. Ultimately we need
// per entry in stencil:
@ -162,7 +160,6 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
// proc1. We first want the pair : proc1 receive and proc2 send
// and then the pair : proc1 send, proc2 receive
labelList procToInterface(Pstream::nProcs(), -1);
forAll(procOwner, proci)
@ -266,32 +263,31 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
{
if (interface != -1)
{
interface = procToInterface[interface]+boundary().size();
interface = procToInterface[interface]+mesh_.boundary().size();
}
}
}
// Get addressing and interfaces of all interfaces
UPtrList<const labelUList> patchAddr;
{
const fvBoundaryMesh& fvp = boundary();
const fvBoundaryMesh& fvp = mesh_.boundary();
patchAddr.setSize(fvp.size() + remoteStencilInterfaces_.size());
//allInterfaces_ = dynamicMotionSolverFvMesh::interfaces();
allInterfaces_ = dynamicFvMesh::interfaces();
allInterfaces_ = mesh_.fvMesh::interfaces();
allInterfaces_.setSize(patchAddr.size());
forAll(fvp, patchi)
{
patchAddr.set(patchi, &fvp[patchi].faceCells());
}
forAll(remoteStencilInterfaces_, i)
{
const label patchi = fvp.size()+i;
label patchi = fvp.size()+i;
const lduPrimitiveProcessorInterface& pp =
remoteStencilInterfaces_[i];
@ -305,6 +301,7 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
allInterfaces_.set(patchi, &pp);
}
}
const lduSchedule ps
(
lduPrimitiveMesh::nonBlockingSchedule<processorLduInterface>
@ -317,7 +314,7 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
(
new fvMeshPrimitiveLduAddressing
(
nCells(),
mesh_.nCells(),
std::move(lowerAddr),
std::move(upperAddr),
patchAddr,
@ -336,22 +333,22 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
<< " upper:" << addr.upperAddr().size() << endl;
// Using lduAddressing::patch
forAll(patchAddr, patchI)
forAll(patchAddr, patchi)
{
Pout<< " " << patchI << "\tpatchAddr:"
<< addr.patchAddr(patchI).size()
Pout<< " " << patchi << "\tpatchAddr:"
<< addr.patchAddr(patchi).size()
<< endl;
}
// Using interfaces
const lduInterfacePtrsList& iFaces = allInterfaces_;
const lduInterfacePtrsList& iFaces = mesh_.interfaces();
Pout<< "Adapted interFaces:" << iFaces.size() << endl;
forAll(iFaces, patchI)
forAll(iFaces, patchi)
{
if (iFaces.set(patchI))
if (iFaces.set(patchi))
{
Pout<< " " << patchI << "\tinterface:"
<< iFaces[patchI].type() << endl;
Pout<< " " << patchi << "\tinterface:"
<< iFaces[patchi].type() << endl;
}
}
}
@ -360,7 +357,7 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
}
Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
Foam::scalar Foam::oversetFvMeshBase::cellAverage
(
const labelList& types,
const labelList& nbrTypes,
@ -370,16 +367,16 @@ Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
bitSet& isFront
) const
{
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const cell& cFaces = cells()[celli];
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
const cell& cFaces = mesh_.cells()[celli];
scalar avg = 0.0;
label n = 0;
label nFront = 0;
for (const label facei : cFaces)
{
if (isInternalFace(facei))
if (mesh_.isInternalFace(facei))
{
label nbrCelli = (own[facei] == celli ? nei[facei] : own[facei]);
if (norm[nbrCelli] == -GREAT)
@ -399,7 +396,7 @@ Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
}
else
{
if (nbrNorm[facei-nInternalFaces()] == -GREAT)
if (nbrNorm[facei-mesh_.nInternalFaces()] == -GREAT)
{
if (isFront.set(facei))
{
@ -408,7 +405,7 @@ Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
}
else
{
avg += nbrNorm[facei-nInternalFaces()];
avg += nbrNorm[facei-mesh_.nInternalFaces()];
n++;
}
}
@ -425,13 +422,13 @@ Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
}
void Foam::dynamicOversetFvMesh::writeAgglomeration
void Foam::oversetFvMeshBase::writeAgglomeration
(
const GAMGAgglomeration& agglom
) const
{
labelList cellToCoarse(identity(nCells()));
labelListList coarseToCell(invertOneToMany(nCells(), cellToCoarse));
labelList cellToCoarse(identity(mesh_.nCells()));
labelListList coarseToCell(invertOneToMany(mesh_.nCells(), cellToCoarse));
// Write initial agglomeration
{
@ -440,13 +437,13 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
IOobject
(
"agglomeration",
this->time().timeName(),
*this,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
*this,
mesh_,
dimensionedScalar(dimless, Zero)
);
scalarField& fld = scalarAgglomeration.primitiveFieldRef();
@ -463,7 +460,7 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
scalarAgglomeration.write();
Info<< "Writing initial cell distribution to "
<< this->time().timeName() << endl;
<< mesh_.time().timeName() << endl;
}
@ -495,13 +492,13 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
IOobject
(
"agglomeration_" + Foam::name(level),
this->time().timeName(),
*this,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
*this,
mesh_,
dimensionedScalar(dimless, Zero)
);
scalarField& fld = scalarAgglomeration.primitiveFieldRef();
@ -526,52 +523,30 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::dynamicOversetFvMesh::dynamicOversetFvMesh
(
const IOobject& io,
const bool doInit
)
Foam::oversetFvMeshBase::oversetFvMeshBase(const fvMesh& mesh, bool doInit)
:
dynamicMotionSolverListFvMesh(io, doInit)
mesh_(mesh),
active_(false)
{
if (doInit)
{
init(false); // do not initialise lower levels
}
}
bool Foam::dynamicOversetFvMesh::init(const bool doInit)
{
if (doInit)
{
dynamicMotionSolverListFvMesh::init(doInit);
}
active_ = false;
// Load stencil (but do not update)
(void)Stencil::New(*this, false);
// Assume something changed
return true;
(void)Stencil::New(mesh_, false);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::dynamicOversetFvMesh::~dynamicOversetFvMesh()
Foam::oversetFvMeshBase::~oversetFvMeshBase()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::lduAddressing& Foam::dynamicOversetFvMesh::lduAddr() const
const Foam::lduAddressing& Foam::oversetFvMeshBase::lduAddr() const
{
if (!active_)
{
//return dynamicMotionSolverFvMesh::lduAddr();
return dynamicFvMesh::lduAddr();
return mesh_.fvMesh::lduAddr();
}
if (!lduPtr_)
{
@ -582,12 +557,11 @@ const Foam::lduAddressing& Foam::dynamicOversetFvMesh::lduAddr() const
}
Foam::lduInterfacePtrsList Foam::dynamicOversetFvMesh::interfaces() const
Foam::lduInterfacePtrsList Foam::oversetFvMeshBase::interfaces() const
{
if (!active_)
{
//return dynamicMotionSolverFvMesh::interfaces();
return dynamicFvMesh::interfaces();
return mesh_.fvMesh::interfaces();
}
if (!lduPtr_)
{
@ -599,7 +573,7 @@ Foam::lduInterfacePtrsList Foam::dynamicOversetFvMesh::interfaces() const
const Foam::fvMeshPrimitiveLduAddressing&
Foam::dynamicOversetFvMesh::primitiveLduAddr() const
Foam::oversetFvMeshBase::primitiveLduAddr() const
{
if (!lduPtr_)
{
@ -611,10 +585,8 @@ Foam::dynamicOversetFvMesh::primitiveLduAddr() const
}
bool Foam::dynamicOversetFvMesh::update()
bool Foam::oversetFvMeshBase::update()
{
//if (dynamicMotionSolverFvMesh::update())
if (dynamicMotionSolverListFvMesh::update())
{
// Calculate the local extra faces for the interpolation. Note: could
// let demand-driven lduAddr() trigger it but just to make sure.
@ -631,55 +603,43 @@ bool Foam::dynamicOversetFvMesh::update()
}
Foam::word Foam::dynamicOversetFvMesh::baseName(const word& name)
bool Foam::oversetFvMeshBase::interpolateFields()
{
if (name.ends_with("_0"))
{
return baseName(name.substr(0, name.size()-2));
}
const cellCellStencilObject& overlap = Stencil::New(mesh_);
return name;
}
bool Foam::dynamicOversetFvMesh::interpolateFields()
{
// Add the stencil suppression list
wordHashSet suppressed(Stencil::New(*this).nonInterpolatedFields());
wordHashSet suppressed(overlap.nonInterpolatedFields());
// Use whatever the solver has set up as suppression list
const dictionary* dictPtr
(
this->schemesDict().findDict("oversetInterpolationSuppressed")
mesh_.schemesDict().findDict("oversetInterpolationSuppressed")
);
if (dictPtr)
{
suppressed.insert(dictPtr->toc());
}
interpolate<volScalarField>(suppressed);
interpolate<volVectorField>(suppressed);
interpolate<volSphericalTensorField>(suppressed);
interpolate<volSymmTensorField>(suppressed);
interpolate<volTensorField>(suppressed);
overlap.interpolate<volScalarField>(mesh_, suppressed);
overlap.interpolate<volVectorField>(mesh_, suppressed);
overlap.interpolate<volSphericalTensorField>(mesh_, suppressed);
overlap.interpolate<volSymmTensorField>(mesh_, suppressed);
overlap.interpolate<volTensorField>(mesh_, suppressed);
return true;
}
bool Foam::dynamicOversetFvMesh::writeObject
bool Foam::oversetFvMeshBase::writeObject
(
IOstreamOption streamOpt,
const bool valid
) const
{
//bool ok = dynamicMotionSolverFvMesh::writeObject(streamOpt, valid);
bool ok = dynamicMotionSolverListFvMesh::writeObject(streamOpt, valid);
// For postprocessing : write cellTypes and zoneID
bool ok = false;
{
const cellCellStencilObject& overlap = Stencil::New(*this);
const cellCellStencilObject& overlap = Stencil::New(mesh_);
const labelUList& cellTypes = overlap.cellTypes();
@ -688,13 +648,13 @@ bool Foam::dynamicOversetFvMesh::writeObject
IOobject
(
"cellTypes",
this->time().timeName(),
*this,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
*this,
mesh_,
dimensionedScalar(dimless, Zero),
zeroGradientFvPatchScalarField::typeName
);
@ -712,18 +672,18 @@ bool Foam::dynamicOversetFvMesh::writeObject
IOobject
(
"zoneID",
this->time().timeName(),
*this,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
*this,
mesh_,
dimensionedScalar(dimless, Zero),
zeroGradientFvPatchScalarField::typeName
);
const cellCellStencilObject& overlap = Stencil::New(*this);
const cellCellStencilObject& overlap = Stencil::New(mesh_);
const labelIOList& zoneID = overlap.zoneID();
forAll(zoneID, cellI)
@ -735,7 +695,7 @@ bool Foam::dynamicOversetFvMesh::writeObject
}
if (debug)
{
const cellCellStencilObject& overlap = Stencil::New(*this);
const cellCellStencilObject& overlap = Stencil::New(mesh_);
const labelIOList& zoneID = overlap.zoneID();
const labelListList& cellStencil = overlap.cellStencil();
@ -744,7 +704,7 @@ bool Foam::dynamicOversetFvMesh::writeObject
overlap.cellInterpolationMap().distribute(donorZoneID);
// Get remote cellCentres
pointField cc(C());
pointField cc(mesh_.C());
overlap.cellInterpolationMap().distribute(cc);
volScalarField volDonorZoneID
@ -752,13 +712,13 @@ bool Foam::dynamicOversetFvMesh::writeObject
IOobject
(
"donorZoneID",
this->time().timeName(),
*this,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
*this,
mesh_,
dimensionedScalar("minOne", dimless, scalar(-1)),
zeroGradientFvPatchScalarField::typeName
);
@ -774,7 +734,7 @@ bool Foam::dynamicOversetFvMesh::writeObject
if (donorZoneID[stencil[i]] != volDonorZoneID[cellI])
{
WarningInFunction << "Mixed donor meshes for cell "
<< cellI << " at " << C()[cellI]
<< cellI << " at " << mesh_.C()[cellI]
<< " donors:" << UIndirectList<point>(cc, stencil)
<< endl;
volDonorZoneID[cellI] = -2;
@ -784,7 +744,15 @@ bool Foam::dynamicOversetFvMesh::writeObject
}
//- Do not correctBoundaryConditions since re-interpolates!
//volDonorZoneID.correctBoundaryConditions();
volDonorZoneID.writeObject(streamOpt, valid);
cellCellStencil::correctBoundaryConditions
<
volScalarField,
oversetFvPatchField<scalar>
>
(
volDonorZoneID
);
ok = volDonorZoneID.writeObject(streamOpt, valid);
}
return ok;

227
src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.H → src/overset/oversetFvMesh/oversetFvMeshBase.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015-2020 OpenCFD Ltd.
Copyright (C) 2015-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -24,45 +24,53 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::dynamicOversetFvMesh
Foam::oversetFvMeshBase
Description
dynamicFvMesh with support for overset meshes.
Support for overset functionality.
SourceFiles
dynamicOversetFvMesh.C
oversetFvMeshBase.C
\*---------------------------------------------------------------------------*/
#ifndef dynamicOversetFvMesh_H
#define dynamicOversetFvMesh_H
#ifndef oversetFvMeshBase_H
#define oversetFvMeshBase_H
#include "dynamicMotionSolverListFvMesh.H"
#include "labelIOList.H"
#include "fvMeshPrimitiveLduAddressing.H"
#include "fvMeshPrimitiveLduAddressing.H"
#include "lduInterfaceFieldPtrsList.H"
#include "volFields.H"
#include "volFieldsFwd.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class mapDistribute;
class lduPrimitiveProcessorInterface;
class GAMGAgglomeration;
/*---------------------------------------------------------------------------*\
Class dynamicOversetFvMesh Declaration
Class oversetFvMeshBase Declaration
\*---------------------------------------------------------------------------*/
class dynamicOversetFvMesh
:
public dynamicMotionSolverListFvMesh
class oversetFvMeshBase
{
// Private Member Functions
//- No copy construct
oversetFvMeshBase(const oversetFvMeshBase&) = delete;
//- No copy assignment
void operator=(const oversetFvMeshBase&) = delete;
protected:
// Protected data
// Protected Data
//- Reference to mesh
const fvMesh& mesh_;
//- Select base addressing (false) or locally stored extended
// lduAddressing (true)
@ -72,12 +80,11 @@ protected:
mutable autoPtr<fvMeshPrimitiveLduAddressing> lduPtr_;
//- Added (processor)lduInterfaces for remote bits of stencil.
//PtrList<const lduInterface> remoteStencilInterfaces_;
mutable PtrList<const lduPrimitiveProcessorInterface>
remoteStencilInterfaces_;
//- Interfaces for above mesh. Contains both original and
// above added processorLduInterfaces
//- above added processorLduInterfaces
mutable lduInterfacePtrsList allInterfaces_;
//- Corresponding faces (in above lduPtr) to the stencil
@ -105,44 +112,29 @@ protected:
const lduInterfacePtrsList&
) const;
//- Explicit interpolation of acceptor cells from donor cells
template<class T>
void interpolate(Field<T>& psi) const;
//- Explicit interpolation of acceptor cells from donor cells with
// boundary condition handling
template<class GeoField>
void interpolate(GeoField& psi) const;
//- Helper: strip off trailing _0
static word baseName(const word& name);
//- Explicit interpolation of all registered fields. Excludes
// selected fields (and their old-time fields)
template<class GeoField>
void interpolate(const wordHashSet& suppressed);
//- Freeze values at holes
//template<class Type>
//void freezeHoles(fvMatrix<Type>&) const;
//- Get scalar interfaces of certain type
//template<class GeoField, class PatchType>
//lduInterfaceFieldPtrsList scalarInterfaces(const GeoField& psi) const;
//- Determine normalisation for interpolation. This equals the
// original diagonal except stabilised for zero diagonals (possible
// in hole cells)
//- Scale coefficient depending on cell type
template<class Type>
tmp<scalarField> normalisation(const fvMatrix<Type>& m) const;
//- Add interpolation to matrix (coefficients)
template<class Type>
void addInterpolation(fvMatrix<Type>&, const scalarField& norm) const;
void scaleConnection
(
Field<Type>& coeffs,
const labelUList& types,
const scalarList& factor,
const bool setHoleCellValue,
const label celli,
const label facei
) const;
//- Solve given dictionary with settings
template<class Type>
SolverPerformance<Type> solve(fvMatrix<Type>&, const dictionary&) const;
SolverPerformance<Type> solveOverset
(
fvMatrix<Type>& m,
const dictionary&
) const;
//- Debug: correct coupled bc
template<class GeoField>
@ -166,36 +158,24 @@ protected:
) const;
private:
// Private Member Functions
//- No copy construct
dynamicOversetFvMesh(const dynamicOversetFvMesh&) = delete;
//- No copy assignment
void operator=(const dynamicOversetFvMesh&) = delete;
public:
//- Runtime type information
TypeName("dynamicOversetFvMesh");
TypeName("oversetFvMeshBase");
// Constructors
//- Construct from IOobject
dynamicOversetFvMesh(const IOobject& io, const bool doInit=true);
oversetFvMeshBase(const fvMesh& mesh, const bool doInit=true);
//- Destructor
virtual ~dynamicOversetFvMesh();
virtual ~oversetFvMeshBase();
// Member Functions
// Extended addressing
//- Return extended ldu addressing
@ -231,12 +211,14 @@ public:
{
DebugInfo<< "Switching to extended addressing with nFaces:"
<< primitiveLduAddr().lowerAddr().size()
<< " nInterfaces:" << allInterfaces_.size()
<< endl;
}
else
{
DebugInfo<< "Switching to base addressing with nFaces:"
<< fvMesh::lduAddr().lowerAddr().size()
<< mesh_.fvMesh::lduAddr().lowerAddr().size()
<< " nInterfaces:" << mesh_.fvMesh::interfaces().size()
<< endl;
}
}
@ -244,108 +226,17 @@ public:
// Overset
// Explicit interpolation
//- Manipulate the matrix to add the interpolation/set hole
// values
template<class Type>
void addInterpolation
(
fvMatrix<Type>& m,
const scalarField& normalisation,
const bool setHoleCellValue,
const Type& holeCellValue
) const;
virtual void interpolate(scalarField& psi) const
{
interpolate<scalar>(psi);
}
virtual void interpolate(vectorField& psi) const
{
interpolate<vector>(psi);
}
virtual void interpolate(sphericalTensorField& psi) const
{
interpolate<sphericalTensor>(psi);
}
virtual void interpolate(symmTensorField& psi) const
{
interpolate<symmTensor>(psi);
}
virtual void interpolate(tensorField& psi) const
{
interpolate<tensor>(psi);
}
virtual void interpolate(volScalarField& psi) const
{
interpolate<volScalarField>(psi);
}
virtual void interpolate(volVectorField& psi) const
{
interpolate<volVectorField>(psi);
}
virtual void interpolate(volSphericalTensorField& psi) const
{
interpolate<volSphericalTensorField>(psi);
}
virtual void interpolate(volSymmTensorField& psi) const
{
interpolate<volSymmTensorField>(psi);
}
virtual void interpolate(volTensorField& psi) const
{
interpolate<volTensorField>(psi);
}
// Implicit interpolation (matrix manipulation)
//- Solve returning the solution statistics given convergence
// tolerance. Use the given solver controls
virtual SolverPerformance<scalar> solve
(
fvMatrix<scalar>& m,
const dictionary& dict
) const
{
return solve<scalar>(m, dict);
}
//- Solve returning the solution statistics given convergence
// tolerance. Use the given solver controls
virtual SolverPerformance<vector> solve
(
fvMatrix<vector>& m,
const dictionary& dict
) const
{
return solve<vector>(m, dict);
}
//- Solve returning the solution statistics given convergence
// tolerance. Use the given solver controls
virtual SolverPerformance<symmTensor> solve
(
fvMatrix<symmTensor>& m,
const dictionary& dict
) const
{
return solve<symmTensor>(m, dict);
}
//- Solve returning the solution statistics given convergence
// tolerance. Use the given solver controls
virtual SolverPerformance<tensor> solve
(
fvMatrix<tensor>& m,
const dictionary& dict
) const
{
return solve<tensor>(m, dict);
}
//- Initialise all non-demand-driven data
virtual bool init(const bool doInit);
//- Update the mesh for both mesh motion and topology change
virtual bool update();
@ -372,6 +263,12 @@ public:
typename GeoField::Boundary& bfld,
const bool typeOnly
);
//- Determine normalisation for interpolation. This equals the
//- original diagonal except stabilised for zero diagonals (possible
//- in hole cells)
template<class Type>
tmp<scalarField> normalisation(const fvMatrix<Type>& m) const;
};
@ -382,7 +279,7 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "dynamicOversetFvMeshTemplates.C"
# include "oversetFvMeshBaseTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

324
src/overset/dynamicOversetFvMesh/dynamicOversetFvMeshTemplates.C → src/overset/oversetFvMesh/oversetFvMeshBaseTemplates.C
View File

@ -27,7 +27,6 @@ License
#include "volFields.H"
#include "fvMatrix.H"
#include "cellCellStencilObject.H"
#include "oversetFvPatchField.H"
#include "calculatedProcessorFvPatchField.H"
#include "lduInterfaceFieldPtrsList.H"
@ -36,84 +35,46 @@ License
// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
template<class T>
void Foam::dynamicOversetFvMesh::interpolate(Field<T>& psi) const
template<class Type>
void Foam::oversetFvMeshBase::scaleConnection
(
Field<Type>& coeffs,
const labelUList& types,
const scalarList& factor,
const bool setHoleCellValue,
const label celli,
const label facei
) const
{
const cellCellStencil& overlap = Stencil::New(*this);
const labelListList& stencil = overlap.cellStencil();
const label cType = types[celli];
const scalar f = factor[celli];
if (stencil.size() != nCells())
if (cType == cellCellStencil::INTERPOLATED)
{
return;
coeffs[facei] *= 1.0-f;
}
const mapDistribute& map = overlap.cellInterpolationMap();
const List<scalarList>& wghts = overlap.cellInterpolationWeights();
const labelList& cellIDs = overlap.interpolationCells();
const scalarList& factor = overlap.cellInterpolationWeight();
Field<T> work(psi);
map.mapDistributeBase::distribute(work, UPstream::msgType()+1);
forAll(cellIDs, i)
else if (cType == cellCellStencil::HOLE)
{
label celli = cellIDs[i];
const scalarList& w = wghts[celli];
const labelList& nbrs = stencil[celli];
const scalar f = factor[celli];
T s(pTraits<T>::zero);
forAll(nbrs, nbrI)
{
s += w[nbrI]*work[nbrs[nbrI]];
}
//Pout<< "Interpolated value:" << s << endl;
//T oldPsi = psi[celli];
psi[celli] = (1.0-f)*psi[celli] + f*s;
//Pout<< "psi was:" << oldPsi << " now:" << psi[celli] << endl;
// Disconnect hole cell from influence of neighbour
coeffs[facei] = pTraits<Type>::zero;
}
}
template<class GeoField>
void Foam::dynamicOversetFvMesh::interpolate(GeoField& psi) const
{
interpolate(psi.primitiveFieldRef());
psi.correctBoundaryConditions();
}
template<class GeoField>
void Foam::dynamicOversetFvMesh::interpolate(const wordHashSet& suppressed)
{
auto flds(this->lookupClass<GeoField>());
for (auto fldPtr : flds)
else if (cType == cellCellStencil::SPECIAL)
{
const word& name = fldPtr->name();
if (!suppressed.found(baseName(name)))
if (setHoleCellValue)
{
if (debug)
{
Pout<< "dynamicOversetFvMesh::interpolate: interpolating : "
<< name << endl;
}
interpolate(fldPtr->primitiveFieldRef());
// Behave like hole
coeffs[facei] = pTraits<Type>::zero;
}
else
{
if (debug)
{
Pout<< "dynamicOversetFvMesh::interpolate: skipping : " << name
<< endl;
}
// Behave like interpolated
coeffs[facei] *= 1.0-f;
}
}
}
template<class GeoField, class PatchType>
void Foam::dynamicOversetFvMesh::correctBoundaryConditions
void Foam::oversetFvMeshBase::correctBoundaryConditions
(
typename GeoField::Boundary& bfld,
const bool typeOnly
@ -133,7 +94,7 @@ void Foam::dynamicOversetFvMesh::correctBoundaryConditions
// Wait for outstanding requests
if (commsType == UPstream::commsTypes::nonBlocking)
{
UPstream::waitRequests(startOfRequests);
Pstream::waitRequests(startOfRequests);
}
forAll(bfld, patchi)
@ -147,7 +108,7 @@ void Foam::dynamicOversetFvMesh::correctBoundaryConditions
template<class Type>
Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
Foam::tmp<Foam::scalarField> Foam::oversetFvMeshBase::normalisation
(
const fvMatrix<Type>& m
) const
@ -165,7 +126,7 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
{
forAll(internalCoeffs, patchi)
{
const labelUList& fc = lduAddr().patchAddr(patchi);
const labelUList& fc = mesh_.lduAddr().patchAddr(patchi);
const Field<Type>& intCoeffs = internalCoeffs[patchi];
const scalarField cmptCoeffs(intCoeffs.component(cmpt));
forAll(fc, i)
@ -199,9 +160,9 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
{
// Walk out the norm across hole cells
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const cellCellStencilObject& overlap = Stencil::New(*this);
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
const cellCellStencilObject& overlap = Stencil::New(mesh_);
const labelUList& types = overlap.cellTypes();
label nHoles = 0;
@ -215,8 +176,8 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
}
}
bitSet isFront(nFaces());
for (label facei = 0; facei < nInternalFaces(); facei++)
bitSet isFront(mesh_.nFaces());
for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
{
label ownType = types[own[facei]];
label neiType = types[nei[facei]];
@ -230,11 +191,16 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
}
}
labelList nbrTypes;
syncTools::swapBoundaryCellList(*this, types, nbrTypes);
for (label facei = nInternalFaces(); facei < nFaces(); facei++)
syncTools::swapBoundaryCellList(mesh_, types, nbrTypes);
for
(
label facei = mesh_.nInternalFaces();
facei < mesh_.nFaces();
++facei
)
{
label ownType = types[own[facei]];
label neiType = nbrTypes[facei-nInternalFaces()];
const label ownType = types[own[facei]];
const label neiType = nbrTypes[facei-mesh_.nInternalFaces()];
if
(
(ownType == cellCellStencil::HOLE)
@ -249,9 +215,9 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
while (true)
{
scalarField nbrNorm;
syncTools::swapBoundaryCellList(*this, extrapolatedNorm, nbrNorm);
syncTools::swapBoundaryCellList(mesh_, extrapolatedNorm, nbrNorm);
bitSet newIsFront(nFaces());
bitSet newIsFront(mesh_.nFaces());
scalarField newNorm(extrapolatedNorm);
label nChanged = 0;
@ -273,7 +239,7 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
}
if
(
isInternalFace(facei)
mesh_.isInternalFace(facei)
&& extrapolatedNorm[nei[facei]] == -GREAT
)
{
@ -299,7 +265,7 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
// Transfer new front
extrapolatedNorm.transfer(newNorm);
isFront.transfer(newIsFront);
syncTools::syncFaceList(*this, isFront, maxEqOp<unsigned int>());
syncTools::syncFaceList(mesh_, isFront, maxEqOp<unsigned int>());
}
@ -324,16 +290,17 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
template<class Type>
void Foam::dynamicOversetFvMesh::addInterpolation
void Foam::oversetFvMeshBase::addInterpolation
(
fvMatrix<Type>& m,
const scalarField& normalisation
const scalarField& normalisation,
const bool setHoleCellValue,
const Type& holeCellValue
) const
{
const cellCellStencilObject& overlap = Stencil::New(*this);
const cellCellStencilObject& overlap = Stencil::New(mesh_);
const List<scalarList>& wghts = overlap.cellInterpolationWeights();
const labelListList& stencil = overlap.cellStencil();
const labelList& cellIDs = overlap.interpolationCells();
const scalarList& factor = overlap.cellInterpolationWeight();
const labelUList& types = overlap.cellTypes();
@ -350,8 +317,6 @@ void Foam::dynamicOversetFvMesh::addInterpolation
const lduAddressing& addr = lduAddr();
const labelUList& upperAddr = addr.upperAddr();
const labelUList& lowerAddr = addr.lowerAddr();
const labelUList& ownerStartAddr = addr.ownerStartAddr();
const labelUList& losortAddr = addr.losortAddr();
const lduInterfacePtrsList& interfaces = allInterfaces_;
if (!isA<fvMeshPrimitiveLduAddressing>(addr))
@ -362,7 +327,6 @@ void Foam::dynamicOversetFvMesh::addInterpolation
}
// 1. Adapt lduMatrix for additional faces and new ordering
upper.setSize(upperAddr.size(), 0.0);
inplaceReorder(reverseFaceMap_, upper);
@ -371,7 +335,7 @@ void Foam::dynamicOversetFvMesh::addInterpolation
//const label nOldInterfaces = dynamicMotionSolverFvMesh::interfaces().size();
const label nOldInterfaces = dynamicFvMesh::interfaces().size();
const label nOldInterfaces = mesh_.fvMesh::interfaces().size();
if (interfaces.size() > nOldInterfaces)
@ -443,7 +407,6 @@ void Foam::dynamicOversetFvMesh::addInterpolation
// fvMatrix:correction, both of which are outside the linear solver.
// Clear out existing connections on cells to be interpolated
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Note: could avoid doing the zeroing of the new faces since these
@ -451,6 +414,11 @@ void Foam::dynamicOversetFvMesh::addInterpolation
forAll(upperAddr, facei)
{
const label l = lowerAddr[facei];
scaleConnection(upper, types, factor, setHoleCellValue, l, facei);
const label u = upperAddr[facei];
scaleConnection(lower, types, factor, setHoleCellValue, u, facei);
/*
if (types[upperAddr[facei]] == cellCellStencil::INTERPOLATED)
{
// Disconnect upper from lower
@ -463,10 +431,21 @@ void Foam::dynamicOversetFvMesh::addInterpolation
label celli = lowerAddr[facei];
upper[facei] *= 1.0-factor[celli];
}
*/
}
for (label patchi = 0; patchi < nOldInterfaces; ++patchi)
{
const labelUList& fc = addr.patchAddr(patchi);
Field<Type>& bCoeffs = m.boundaryCoeffs()[patchi];
Field<Type>& iCoeffs = m.internalCoeffs()[patchi];
forAll(fc, i)
{
scaleConnection(bCoeffs, types, factor, setHoleCellValue, fc[i], i);
scaleConnection(iCoeffs, types, factor, setHoleCellValue, fc[i], i);
}
/*
const labelUList& fc = addr.patchAddr(patchi);
Field<Type>& intCoeffs = m.internalCoeffs()[patchi];
Field<Type>& bouCoeffs = m.boundaryCoeffs()[patchi];
@ -487,6 +466,7 @@ void Foam::dynamicOversetFvMesh::addInterpolation
}
}
}
*/
}
@ -497,8 +477,88 @@ void Foam::dynamicOversetFvMesh::addInterpolation
// Do hole cells. Note: maybe put into interpolationCells() loop above?
forAll(types, celli)
{
if (types[celli] == cellCellStencil::HOLE)
if
(
types[celli] == cellCellStencil::HOLE
|| (setHoleCellValue && types[celli] == cellCellStencil::SPECIAL)
)
{
const Type wantedValue
(
setHoleCellValue
? holeCellValue
: m.psi()[celli]
);
diag[celli] = normalisation[celli];
source[celli] = normalisation[celli]*wantedValue;
}
else if
(
types[celli] == cellCellStencil::INTERPOLATED
|| (!setHoleCellValue && types[celli] == cellCellStencil::SPECIAL)
)
{
const scalar f = factor[celli];
const scalarList& w = wghts[celli];
const labelList& nbrs = stencil[celli];
const labelList& nbrFaces = stencilFaces_[celli];
const labelList& nbrPatches = stencilPatches_[celli];
diag[celli] *= (1.0-f);
source[celli] *= (1.0-f);
forAll(nbrs, nbri)
{
const label patchi = nbrPatches[nbri];
const label facei = nbrFaces[nbri];
if (patchi == -1)
{
const label nbrCelli = nbrs[nbri];
// Add the coefficients
const scalar s = normalisation[celli]*f*w[nbri];
scalar& u = upper[facei];
scalar& l = lower[facei];
if (celli < nbrCelli)
{
diag[celli] += s;
u += -s;
}
else
{
diag[celli] += s;
l += -s;
}
}
else
{
// Patch face. Store in boundaryCoeffs. Note sign change.
//const label globalCelli = globalCellIDs[nbrs[nbri]];
//const label proci =
// globalNumbering.whichProcID(globalCelli);
//const label remoteCelli =
// globalNumbering.toLocal(proci, globalCelli);
//
//Pout<< "for cell:" << celli
// << " need weight from remote slot:" << nbrs[nbri]
// << " proc:" << proci << " remote cell:" << remoteCelli
// << " patch:" << patchi
// << " patchFace:" << facei
// << " weight:" << w[nbri]
// << endl;
const scalar s = normalisation[celli]*f*w[nbri];
m.boundaryCoeffs()[patchi][facei] += pTraits<Type>::one*s;
m.internalCoeffs()[patchi][facei] += pTraits<Type>::one*s;
// Note: do NOT add to diagonal - this is in the
// internalCoeffs and gets added to the diagonal
// inside fvMatrix::solve
}
}
}
/*
label startLabel = ownerStartAddr[celli];
label endLabel = ownerStartAddr[celli + 1];
@ -518,7 +578,7 @@ void Foam::dynamicOversetFvMesh::addInterpolation
diag[celli] = normalisation[celli];
source[celli] = normalisation[celli]*m.psi()[celli];
}
*/
}
@ -530,7 +590,7 @@ void Foam::dynamicOversetFvMesh::addInterpolation
//}
//overlap.cellInterpolationMap().distribute(globalCellIDs);
/*
forAll(cellIDs, i)
{
label celli = cellIDs[i];
@ -636,11 +696,12 @@ void Foam::dynamicOversetFvMesh::addInterpolation
//}
}
}
*/
}
template<class Type>
Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
Foam::SolverPerformance<Type> Foam::oversetFvMeshBase::solveOverset
(
fvMatrix<Type>& m,
const dictionary& dict
@ -665,17 +726,17 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
{
if (debug)
{
Pout<< "dynamicOversetFvMesh::solve() :"
Pout<< "oversetFvMeshBase::solveOverset() :"
<< " bypassing matrix adjustment for field " << m.psi().name()
<< endl;
}
//return dynamicMotionSolverFvMesh::solve(m, dict);
return dynamicFvMesh::solve(m, dict);
return mesh_.fvMesh::solve(m, dict);
}
if (debug)
{
Pout<< "dynamicOversetFvMesh::solve() :"
Pout<< "oversetFvMeshBase::solveOverset() :"
<< " adjusting matrix for interpolation for field "
<< m.psi().name() << endl;
}
@ -683,20 +744,20 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
// Calculate stabilised diagonal as normalisation for interpolation
const scalarField norm(normalisation(m));
if (debug)
if (debug && mesh_.time().outputTime())
{
volScalarField scale
(
IOobject
(
m.psi().name() + "_normalisation",
this->time().timeName(),
*this,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
*this,
mesh_,
dimensionedScalar(dimless, Zero)
);
scale.ref().field() = norm;
@ -709,7 +770,7 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
if (debug)
{
Pout<< "dynamicOversetFvMesh::solve() :"
Pout<< "oversetFvMeshBase::solveOverset() :"
<< " writing matrix normalisation for field " << m.psi().name()
<< " to " << scale.name() << endl;
}
@ -725,9 +786,15 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
scalarField oldLower(m.lower());
FieldField<Field, Type> oldInt(m.internalCoeffs());
FieldField<Field, Type> oldBou(m.boundaryCoeffs());
Field<Type> oldSource(m.source());
scalarField oldDiag(m.diag());
const label oldSize = bpsi.size();
addInterpolation(m, norm);
// Insert the interpolation into the matrix (done inside
// oversetFvPatchField<Type>::manipulateMatrix)
m.boundaryManipulate(bpsi);
// Swap psi values so added patches have patchNeighbourField
correctBoundaryConditions<GeoField, calculatedProcessorFvPatchField<Type>>
@ -736,18 +803,9 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
true
);
// Print a bit
//write(Pout, m, lduAddr(), interfaces());
//{
// const fvSolution& sol = static_cast<const fvSolution&>(*this);
// const dictionary& pDict = sol.subDict("solvers").subDict("p");
// writeAgglomeration(GAMGAgglomeration::New(m, pDict));
//}
// Use lower level solver
//SolverPerformance<Type> s(dynamicMotionSolverFvMesh::solve(m, dict));
SolverPerformance<Type> s(dynamicFvMesh::solve(m, dict));
SolverPerformance<Type> s(mesh_.fvMesh::solve(m, dict));
// Restore boundary field
bpsi.setSize(oldSize);
@ -755,9 +813,45 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
// Restore matrix
m.upper().transfer(oldUpper);
m.lower().transfer(oldLower);
m.source().transfer(oldSource);
m.diag().transfer(oldDiag);
m.internalCoeffs().transfer(oldInt);
m.boundaryCoeffs().transfer(oldBou);
const cellCellStencilObject& overlap = Stencil::New(mesh_);
const labelUList& types = overlap.cellTypes();
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
auto& psi =
const_cast<GeometricField<Type, fvPatchField, volMesh>&>(m.psi());
// Mirror hole cell values next to calculated cells
for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
{
const label ownType = types[own[facei]];
const label neiType = types[nei[facei]];
if
(
ownType == cellCellStencil::HOLE
&& neiType == cellCellStencil::CALCULATED)
{
psi[own[facei]] = psi[nei[facei]];
}
else if
(
ownType == cellCellStencil::CALCULATED
&& neiType == cellCellStencil::HOLE
)
{
psi[nei[facei]] = psi[own[facei]];
}
}
// Switch to original addressing
active(false);
@ -766,7 +860,7 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
template<class Type>
void Foam::dynamicOversetFvMesh::write
void Foam::oversetFvMeshBase::write
(
Ostream& os,
const fvMatrix<Type>& m,
@ -893,7 +987,7 @@ void Foam::dynamicOversetFvMesh::write
template<class GeoField>
void Foam::dynamicOversetFvMesh::correctCoupledBoundaryConditions(GeoField& fld)
void Foam::oversetFvMeshBase::correctCoupledBoundaryConditions(GeoField& fld)
{
typename GeoField::Boundary& bfld = fld.boundaryFieldRef();
@ -912,7 +1006,7 @@ void Foam::dynamicOversetFvMesh::correctCoupledBoundaryConditions(GeoField& fld)
// Wait for outstanding requests
if (commsType == UPstream::commsTypes::nonBlocking)
{
UPstream::waitRequests(startOfRequests);
Pstream::waitRequests(startOfRequests);
}
forAll(bfld, patchi)
@ -927,7 +1021,7 @@ void Foam::dynamicOversetFvMesh::correctCoupledBoundaryConditions(GeoField& fld)
template<class GeoField>
void Foam::dynamicOversetFvMesh::checkCoupledBC(const GeoField& fld)
void Foam::oversetFvMeshBase::checkCoupledBC(const GeoField& fld)
{
Pout<< "** starting checking of " << fld.name() << endl;

84
src/overset/oversetFvMesh/staticOversetFvMesh/staticOversetFvMesh.C Normal file
View File

@ -0,0 +1,84 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify i
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 "staticOversetFvMesh.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(staticOversetFvMesh, 0);
addToRunTimeSelectionTable(dynamicFvMesh, staticOversetFvMesh, IOobject);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::staticOversetFvMesh::staticOversetFvMesh(const IOobject& io)
:
staticFvMesh(io),
oversetFvMeshBase(static_cast<const fvMesh&>(*this))
{
oversetFvMeshBase::update();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::staticOversetFvMesh::~staticOversetFvMesh()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::staticOversetFvMesh::update()
{
if (!staticFvMesh::update())
{
return false;
}
oversetFvMeshBase::update();
return true;
}
bool Foam::staticOversetFvMesh::writeObject
(
IOstreamOption streamOpt,
const bool valid
) const
{
bool ok = staticFvMesh::writeObject(streamOpt, valid);
ok = oversetFvMeshBase::writeObject(streamOpt, valid) && ok;
return ok;
}
// ************************************************************************* //

225
src/overset/oversetFvMesh/staticOversetFvMesh/staticOversetFvMesh.H Normal file
View File

@ -0,0 +1,225 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
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::staticOversetFvMesh
Description
fvMesh with support for overset meshes.
SourceFiles
staticOversetFvMesh.C
\*---------------------------------------------------------------------------*/
#ifndef staticOversetFvMesh_H
#define staticOversetFvMesh_H
#include "staticFvMesh.H"
#include "oversetFvMeshBase.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class staticOversetFvMesh Declaration
\*---------------------------------------------------------------------------*/
class staticOversetFvMesh
:
public staticFvMesh,
public oversetFvMeshBase
{
// Private Member Functions
//- No copy construct
staticOversetFvMesh(const staticOversetFvMesh&) = delete;
//- No copy assignment
void operator=(const staticOversetFvMesh&) = delete;
public:
//- Runtime type information
TypeName("staticOversetFvMesh");
// Constructors
//- Construct from IOobject
staticOversetFvMesh(const IOobject& io);
//- Destructor
virtual ~staticOversetFvMesh();
// Member Functions
//- Override ldu addressing
virtual const lduAddressing& lduAddr() const
{
return oversetFvMeshBase::lduAddr();
}
//- Override ldu addressing
virtual lduInterfacePtrsList interfaces() const
{
return oversetFvMeshBase::interfaces();
}
// Overset
// Explicit interpolation
virtual void interpolate(scalarField& psi) const
{
Stencil::New(*this).interpolate<scalar>(*this, psi);
}
virtual void interpolate(vectorField& psi) const
{
Stencil::New(*this).interpolate<vector>(*this, psi);
}
virtual void interpolate(sphericalTensorField& psi) const
{
Stencil::New
(
*this
).interpolate<sphericalTensor>(*this, psi);
}
virtual void interpolate(symmTensorField& psi) const
{
Stencil::New(*this).interpolate<symmTensor>(*this, psi);
}
virtual void interpolate(tensorField& psi) const
{
Stencil::New(*this).interpolate<tensor>(*this, psi);
}
virtual void interpolate(volScalarField& psi) const
{
Stencil::New(*this).interpolate<volScalarField>(psi);
}
virtual void interpolate(volVectorField& psi) const
{
Stencil::New(*this).interpolate<volVectorField>(psi);
}
virtual void interpolate(volSphericalTensorField& psi) const
{
Stencil::New(*this).interpolate
<
volSphericalTensorField
>(psi);
}
virtual void interpolate(volSymmTensorField& psi) const
{
Stencil::New(*this).interpolate<volSymmTensorField>(psi);
}
virtual void interpolate(volTensorField& psi) const
{
Stencil::New(*this).interpolate<volTensorField>(psi);
}
// Implicit interpolation (matrix manipulation)
//- Solve returning the solution statistics given convergence
//- tolerance. Use the given solver controls
virtual SolverPerformance<scalar> solve
(
fvMatrix<scalar>& m,
const dictionary& dict
) const
{
return oversetFvMeshBase::solveOverset<scalar>(m, dict);
}
//- Solve returning the solution statistics given convergence
//- tolerance. Use the given solver controls
virtual SolverPerformance<vector> solve
(
fvMatrix<vector>& m,
const dictionary& dict
) const
{
return oversetFvMeshBase::solveOverset<vector>(m, dict);
}
//- Solve returning the solution statistics given convergence
//- tolerance. Use the given solver controls
virtual SolverPerformance<symmTensor> solve
(
fvMatrix<symmTensor>& m,
const dictionary& dict
) const
{
return oversetFvMeshBase::solveOverset<symmTensor>(m, dict);
}
//- Solve returning the solution statistics given convergence
//- tolerance. Use the given solver controls
virtual SolverPerformance<tensor> solve
(
fvMatrix<tensor>& m,
const dictionary& dict
) const
{
return oversetFvMeshBase::solveOverset<tensor>(m, dict);
}
//- Update the mesh for both mesh motion and topology change
virtual bool update();
//- Write using given format, version and compression
virtual bool writeObject
(
IOstreamOption streamOpt,
const bool valid
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

58
src/overset/oversetPatchPhiErr/oversetPatchPhiErr.C Normal file
View File

@ -0,0 +1,58 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "oversetPatchPhiErr.H"
#include "oversetFvPatchField.H"
#include "volFields.H"
#include "volFieldsFwd.H"
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
void Foam::oversetPatchPhiErr
(
const fvScalarMatrix& m,
const surfaceScalarField& phi
)
{
const volScalarField::Boundary& bm = m.psi().boundaryField();
forAll(bm, patchi)
{
const auto& fvp = bm[patchi];
if (isA<oversetFvPatchField<scalar>>(fvp))
{
const oversetFvPatchField<scalar>& op =
refCast<const oversetFvPatchField<scalar>>(fvp);
op.fringeFlux(m, phi);
}
}
}
// ************************************************************************* //

62
src/overset/oversetPatchPhiErr/oversetPatchPhiErr.H Normal file
View File

@ -0,0 +1,62 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
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/>.
InNamespace
Foam
Description
oversetPatchPhiErr
SourceFiles
oversetPatchPhiErr.C
\*---------------------------------------------------------------------------*/
#ifndef oversetPatchPhiErr_H
#define oversetPatchPhiErr_H
#include "volFieldsFwd.H"
#include "surfaceFieldsFwd.H"
#include "fvScalarMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
//- Prints out sum of m.flux and phi over the fringe faces
void oversetPatchPhiErr
(
const fvScalarMatrix& m,
const surfaceScalarField& phi
);
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

63
src/overset/oversetPolyPatch/oversetFvPatch.C
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2019 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -41,4 +41,65 @@ namespace Foam
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::tmp<Foam::labelField> Foam::oversetFvPatch::interfaceInternalField
(
const labelUList& internalData
) const
{
// Return all internal values
return patchInternalField(internalData);
}
Foam::tmp<Foam::labelField> Foam::oversetFvPatch::interfaceInternalField
(
const labelUList& internalData,
const labelUList& faceCells
) const
{
return patchInternalField(internalData, faceCells);
}
void Foam::oversetFvPatch::initInternalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& iF
) const
{}
void Foam::oversetFvPatch::initInternalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& iF,
const labelUList& faceCells
) const
{}
Foam::tmp<Foam::labelField> Foam::oversetFvPatch::internalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& iF
) const
{
return tmp<labelField>::New(iF);
}
Foam::tmp<Foam::labelField> Foam::oversetFvPatch::internalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& iF,
const labelUList& nbrFaceCells
) const
{
NotImplemented
return tmp<labelField>::New(iF);
}
// ************************************************************************* //

66
src/overset/oversetPolyPatch/oversetFvPatch.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2019 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -37,6 +37,9 @@ SourceFiles
#ifndef oversetFvPatch_H
#define oversetFvPatch_H
#include "coupledFvPatch.H"
#include "oversetLduInterface.H"
#include "lduInterface.H"
#include "fvPatch.H"
#include "oversetPolyPatch.H"
#include "fvBoundaryMesh.H"
@ -52,9 +55,11 @@ namespace Foam
class oversetFvPatch
:
public fvPatch
public fvPatch,
public lduInterface,
public oversetLduInterface
{
// Private data
// Private Data
const oversetPolyPatch& oversetPolyPatch_;
@ -75,7 +80,7 @@ public:
{}
// Member functions
// Member Functions
// Access
@ -85,9 +90,6 @@ public:
return oversetPolyPatch_;
}
// Access
//- Return faceCell addressing
virtual const labelUList& faceCells() const
{
@ -99,6 +101,56 @@ public:
{
return oversetPolyPatch_.master();
}
// Interface Transfer Functions
//- Return the values of the given internal data adjacent to
//- the interface as a field
virtual tmp<labelField> interfaceInternalField
(
const labelUList& internalData
) const;
//- Return the values of the given internal data adjacent to
//- the interface as a field using faceCell mapping
virtual tmp<labelField> interfaceInternalField
(
const labelUList& internalData,
const labelUList& faceCells
) const;
//- Initialise transfer of internal field adjacent to the interface
virtual void initInternalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& iF
) const;
//- Initialise transfer of internal field adjacent to the interface
//- using faceCells mapping
virtual void initInternalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& iF,
const labelUList& faceCells
) const;
//- Transfer and return internal field adjacent to the interface
virtual tmp<labelField> internalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& iF
) const;
//- Transfer and return internal field adjacent to the interface
//- using nbrFaceCells mapping
virtual tmp<labelField> internalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& iF,
const labelUList& nbrFaceCells
) const;
};

885
src/overset/oversetPolyPatch/oversetFvPatchField.C
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2018 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -28,7 +28,8 @@ License
#include "volFields.H"
#include "cellCellStencil.H"
#include "cellCellStencilObject.H"
#include "dynamicOversetFvMesh.H"
#include "oversetFvMeshBase.H"
#include "syncTools.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -39,8 +40,16 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
const DimensionedField<Type, volMesh>& iF
)
:
zeroGradientFvPatchField<Type>(p, iF),
oversetPatch_(refCast<const oversetFvPatch>(p))
coupledFvPatchField<Type>(p, iF),
oversetPatch_(refCast<const oversetFvPatch>(p)),
setHoleCellValue_(false),
massCorrection_(false),
interpolateHoleCellValue_(false),
holeCellValue_(pTraits<Type>::min),
fringeUpperCoeffs_(),
fringeLowerCoeffs_(),
fringeFaces_(),
zoneId_(-1)
{}
@ -53,8 +62,16 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
const fvPatchFieldMapper& mapper
)
:
zeroGradientFvPatchField<Type>(ptf, p, iF, mapper),
oversetPatch_(refCast<const oversetFvPatch>(p))
coupledFvPatchField<Type>(ptf, p, iF, mapper),
oversetPatch_(refCast<const oversetFvPatch>(p)),
setHoleCellValue_(ptf.setHoleCellValue_),
massCorrection_(ptf.massCorrection_),
interpolateHoleCellValue_(ptf.interpolateHoleCellValue_),
holeCellValue_(ptf.holeCellValue_),
fringeUpperCoeffs_(ptf.fringeUpperCoeffs_),
fringeLowerCoeffs_(ptf.fringeLowerCoeffs_),
fringeFaces_(ptf.fringeFaces_),
zoneId_(ptf.zoneId_)
{}
@ -66,9 +83,37 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
const dictionary& dict
)
:
zeroGradientFvPatchField<Type>(p, iF, dict),
oversetPatch_(refCast<const oversetFvPatch>(p, dict))
{}
coupledFvPatchField<Type>(p, iF, dict, false),
oversetPatch_(refCast<const oversetFvPatch>(p, dict)),
setHoleCellValue_(dict.getOrDefault("setHoleCellValue", false)),
massCorrection_(dict.getOrDefault("massCorrection", false)),
interpolateHoleCellValue_
(
dict.getOrDefault("interpolateHoleCellValue", false)
),
holeCellValue_
(
setHoleCellValue_
? dict.get<Type>("holeCellValue")
: pTraits<Type>::min
),
fringeUpperCoeffs_(),
fringeLowerCoeffs_(),
fringeFaces_(),
zoneId_(dict.getOrDefault<label>("zone", -1))
{
if (dict.found("value"))
{
Field<Type>::operator=
(
Field<Type>("value", dict, p.size())
);
}
else
{
Field<Type>::operator=(this->patchInternalField());
}
}
template<class Type>
@ -77,8 +122,16 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
const oversetFvPatchField<Type>& ptf
)
:
zeroGradientFvPatchField<Type>(ptf),
oversetPatch_(ptf.oversetPatch_)
coupledFvPatchField<Type>(ptf),
oversetPatch_(ptf.oversetPatch_),
setHoleCellValue_(ptf.setHoleCellValue_),
massCorrection_(ptf.massCorrection_),
interpolateHoleCellValue_(ptf.interpolateHoleCellValue_),
holeCellValue_(ptf.holeCellValue_),
fringeUpperCoeffs_(ptf.fringeUpperCoeffs_),
fringeLowerCoeffs_(ptf.fringeLowerCoeffs_),
fringeFaces_(ptf.fringeFaces_),
zoneId_(ptf.zoneId_)
{}
@ -89,13 +142,454 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
const DimensionedField<Type, volMesh>& iF
)
:
zeroGradientFvPatchField<Type>(ptf, iF),
oversetPatch_(ptf.oversetPatch_)
coupledFvPatchField<Type>(ptf, iF),
oversetPatch_(ptf.oversetPatch_),
setHoleCellValue_(ptf.setHoleCellValue_),
massCorrection_(ptf.massCorrection_),
interpolateHoleCellValue_(ptf.interpolateHoleCellValue_),
holeCellValue_(ptf.holeCellValue_),
fringeUpperCoeffs_(ptf.fringeUpperCoeffs_),
fringeLowerCoeffs_(ptf.fringeLowerCoeffs_),
fringeFaces_(ptf.fringeFaces_),
zoneId_(ptf.zoneId_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::oversetFvPatchField<Type>::storeFringeCoefficients
(
const fvMatrix<Type>& matrix
)
{
const fvMesh& mesh = this->internalField().mesh();
const cellCellStencilObject& overlap = Stencil::New(mesh);
const labelList& cellTypes = overlap.cellTypes();
const labelList& zoneID = overlap.zoneID();
const labelUList& own = mesh.owner();
const labelUList& nei = mesh.neighbour();
const polyBoundaryMesh& boundaryMesh = mesh.boundaryMesh();
label fringesFaces = 0;
forAll(own, facei)
{
const label zonei = zoneID[own[facei]];
const label ownType = cellTypes[own[facei]];
const label neiType = cellTypes[nei[facei]];
const bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
const bool neiCalc =
(ownType == cellCellStencil::INTERPOLATED)
&& (neiType == cellCellStencil::CALCULATED);
const bool ownNei = (ownCalc || neiCalc);
if
(
(ownNei && (zonei == zoneId_)) || (ownNei && (zoneId_ == -1))
)
{
fringesFaces++;
}
}
const fvPatchList& patches = mesh.boundary();
labelList neiCellTypes;
syncTools::swapBoundaryCellList(mesh, cellTypes, neiCellTypes);
{
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
const labelUList& fc = curPatch.faceCells();
const label start = curPatch.start();
forAll(fc, i)
{
const label facei = start + i;
const label celli = fc[i];
const label ownType = cellTypes[celli];
const label neiType = neiCellTypes[facei-mesh.nInternalFaces()];
const label zonei = zoneID[celli];
const bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
if (ownCalc && (zonei == zoneId_))
{
fringesFaces++;
}
}
}
}
fringeUpperCoeffs_.setSize(fringesFaces, Zero);
fringeLowerCoeffs_.setSize(fringesFaces, Zero);
fringeFaces_.setSize(fringesFaces, -1);
const scalarField& upper = matrix.upper();
const scalarField& lower = matrix.lower();
fringesFaces = 0;
forAll(own, facei)
{
const label zonei = zoneID[own[facei]];
const label ownType = cellTypes[own[facei]];
const label neiType = cellTypes[nei[facei]];
const bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
const bool neiCalc =
(ownType == cellCellStencil::INTERPOLATED)
&& (neiType == cellCellStencil::CALCULATED);
const bool ownNei = (ownCalc || neiCalc);
if
(
(ownNei && (zonei == zoneId_)) || (ownNei && (zoneId_ == -1))
)
{
fringeUpperCoeffs_[fringesFaces] = upper[facei];
fringeLowerCoeffs_[fringesFaces] = lower[facei];
fringeFaces_[fringesFaces] = facei;
fringesFaces++;
}
}
forAll(boundaryMesh, patchi)
{
const polyPatch& p = boundaryMesh[patchi];
if (isA<coupledPolyPatch>(p))
{
const labelUList& fc = p.faceCells();
const label start = p.start();
forAll(fc, i)
{
const label facei = start + i;
const label celli = fc[i];
const label ownType = cellTypes[celli];
const label neiType = neiCellTypes[facei-mesh.nInternalFaces()];
const label zonei = zoneID[celli];
const bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
const bool neiCalc =
(neiType == cellCellStencil::CALCULATED)
&& (ownType == cellCellStencil::INTERPOLATED);
if ((ownCalc||neiCalc) && (zonei == zoneId_))
{
fringeLowerCoeffs_[fringesFaces] =
component
(
matrix.internalCoeffs()[patchi][facei],
0
);
fringeUpperCoeffs_[fringesFaces] =
component
(
matrix.boundaryCoeffs()[patchi][facei],
0
);
fringeFaces_[fringesFaces] = facei;
fringesFaces++;
}
}
}
}
}
template<class Type>
void Foam::oversetFvPatchField<Type>::fringeFlux
(
const fvMatrix<Type>& matrix,
const surfaceScalarField& phi
) const
{
scalar massIn = 0;
scalar phiIn = 0;
const Field<Type>& psi = matrix.psi();
const scalarField& upper = matrix.upper();
const scalarField& lower = matrix.lower();
if (this->oversetPatch_.master())
{
const fvMesh& mesh = this->internalField().mesh();
const cellCellStencilObject& overlap = Stencil::New(mesh);
const labelList& cellTypes = overlap.cellTypes();
const labelList& zoneID = overlap.zoneID();
// Check all faces on the outside of interpolated cells
const labelUList& own = mesh.owner();
const labelUList& nei = mesh.neighbour();
label fringesFaces = 0;
forAll(own, facei)
{
const label zonei = zoneID[own[facei]];
const label ownType = cellTypes[own[facei]];
const label neiType = cellTypes[nei[facei]];
const bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
const bool neiCalc =
(ownType == cellCellStencil::INTERPOLATED)
&& (neiType == cellCellStencil::CALCULATED);
const bool ownNei = (ownCalc || neiCalc);
if
(
(ownNei && (zonei == zoneId_)) || (ownNei && (zoneId_ == -1))
)
{
const label fringei = fringeFaces_[fringesFaces];
// Get fringe upper/lower coeffs
//const scalar& ufc = fringeUpperCoeffs_[fringesFaces];
//const scalar& lfc = fringeLowerCoeffs_[fringesFaces];
const scalar& ufc = upper[fringei];
const scalar& lfc = lower[fringei];
const Type curFlux =
ufc*psi[nei[fringei]] - lfc*psi[own[fringei]];
if (neiCalc)
{
phiIn -= phi[fringei];
massIn -= curFlux;
}
else
{
phiIn += phi[fringei];
massIn += curFlux;
}
fringesFaces++;
}
}
}
reduce(massIn, sumOp<scalar>());
reduce(phiIn, sumOp<scalar>());
Info << " gSum(phi) on fringes " << phiIn << endl;
Info << " gSum(p.flux) on fringes " << massIn << endl;
}
template<class Type>
void Foam::oversetFvPatchField<Type>::adjustPsi
(
Field<scalar>& psi,
const lduAddressing& lduAddr,
solveScalarField& result
) const
{
const fvMesh& mesh = this->internalField().mesh();
const cellCellStencilObject& overlap = Stencil::New(mesh);
const labelList& cellTypes = overlap.cellTypes();
const labelList& zoneID = overlap.zoneID();
// Pass-1: accumulate all fluxes, calculate correction factor
scalarField interpolatedCoeffs(fringeUpperCoeffs_.size(), Zero);
// Options for scaling corrections
scalar massIn = 0;
scalar offDiagCoeffs = 0;
labelField facePerCell(cellTypes.size(), 0);
// Check all faces on the outside of interpolated cells
const labelUList& own = mesh.owner();
const labelUList& nei = mesh.neighbour();
label fringesFaces = 0;
{
forAll(own, facei)
{
const label zonei = zoneID[own[facei]];
const label ownType = cellTypes[own[facei]];
const label neiType = cellTypes[nei[facei]];
const bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
const bool neiCalc =
(ownType == cellCellStencil::INTERPOLATED)
&& (neiType == cellCellStencil::CALCULATED);
const bool ownNei = (ownCalc || neiCalc);
if
(
(ownNei && (zonei == zoneId_)) || (ownNei && (zoneId_ == -1))
)
{
// Get fringe upper/lower coeffs
const scalar& ufc = fringeUpperCoeffs_[fringesFaces];
const scalar& lfc = fringeLowerCoeffs_[fringesFaces];
const scalar curFlux =
ufc*psi[nei[facei]] - lfc*psi[own[facei]];
if (neiCalc) // interpolated is owner
{
massIn -= curFlux;
offDiagCoeffs += lfc;
facePerCell[own[facei]]++;
}
else
{
massIn += curFlux;
offDiagCoeffs += ufc;
facePerCell[nei[facei]]++;
}
fringesFaces++;
}
}
}
scalarField weights(facePerCell.size(), scalar(1));
forAll (weights, celli)
{
if (facePerCell[celli] > 1)
{
weights[celli] = scalar(1)/facePerCell[celli];
}
}
// Check all coupled faces on the outside of interpolated cells
labelList neiCellTypes;
syncTools::swapBoundaryCellList(mesh, cellTypes, neiCellTypes);
const fvPatchList& boundaryMesh = mesh.boundary();
forAll(boundaryMesh, patchi)
{
const polyPatch& p = mesh.boundaryMesh()[patchi];
if (isA<coupledPolyPatch>(p))
{
const auto& coupledPatch = refCast<const coupledPolyPatch>(p);
const labelUList& fc = p.faceCells();
const label start = p.start();
forAll(fc, i)
{
const label facei = start + i;
const label celli = fc[i];
const label ownType = cellTypes[celli];
const label neiType = neiCellTypes[facei-mesh.nInternalFaces()];
const label zonei = zoneID[celli];
const bool ownCalc =
(ownType == cellCellStencil::CALCULATED)
&& (neiType == cellCellStencil::INTERPOLATED);
const bool neiCalc =
(neiType == cellCellStencil::CALCULATED)
&& (ownType == cellCellStencil::INTERPOLATED);
if ((ownCalc||neiCalc) && (zonei == zoneId_))
{
const scalar psiOwn = psi[celli];
const scalar& lfc = fringeLowerCoeffs_[fringesFaces];
const scalar curFlux = lfc*psiOwn;
if (ownCalc)
{
massIn -= curFlux;
if (coupledPatch.owner())
{
offDiagCoeffs -= lfc;
}
fringesFaces++;
}
else
{
massIn += curFlux;
if (coupledPatch.owner())
{
offDiagCoeffs -= lfc;
}
fringesFaces++;
}
}
}
}
}
reduce(massIn, sumOp<scalar>());
reduce(offDiagCoeffs, sumOp<scalar>());
scalar psiCorr = -massIn/offDiagCoeffs;
forAll (cellTypes, celli)
{
const bool bInter = (cellTypes[celli] == cellCellStencil::INTERPOLATED);
const label zonei = zoneID[celli];
if
(
(bInter && (zonei == zoneId_)) ||(bInter && (zoneId_ == -1))
)
{
psi[celli] += psiCorr;
}
}
}
template<class Type>
Foam::tmp<Foam::Field<Type>> Foam::oversetFvPatchField<Type>::
patchNeighbourField() const
{
return tmp<Field<Type>>
(
new Field<Type>(this->size(), Zero)
);
}
template<class Type>
void Foam::oversetFvPatchField<Type>::initEvaluate
(
@ -205,25 +699,374 @@ void Foam::oversetFvPatchField<Type>::initEvaluate
Info<< "Interpolating non-suppressed field " << fldName
<< endl;
}
mesh.interpolate
// Interpolate without bc update (since would trigger infinite
// recursion back to oversetFvPatchField<Type>::evaluate)
// The only problem is bcs that use the new cell values
// (e.g. zeroGradient, processor). These need to appear -after-
// the 'overset' bc.
const cellCellStencil& overlap = Stencil::New(mesh);
Field<Type>& fld =
const_cast<Field<Type>&>(this->primitiveField());
// tbd: different weights for different variables ...
cellCellStencil::interpolate
(
const_cast<Field<Type>&>
(
this->primitiveField()
)
fld,
mesh,
overlap,
overlap.cellInterpolationWeights()
);
if (this->setHoleCellValue_)
{
const labelUList& types = overlap.cellTypes();
label nConstrained = 0;
forAll(types, celli)
{
const label cType = types[celli];
if
(
cType == cellCellStencil::HOLE
|| cType == cellCellStencil::SPECIAL
)
{
fld[celli] = this->holeCellValue_;
nConstrained++;
}
}
if (debug)
{
Pout<< FUNCTION_NAME << " field:" << fldName
<< " patch:" << this->oversetPatch_.name()
<< " set:" << nConstrained << " cells to:"
<< this->holeCellValue_ << endl;
}
}
}
}
}
coupledFvPatchField<Type>::initEvaluate(commsType);
}
template<class Type>
void Foam::oversetFvPatchField<Type>::manipulateMatrix
(
fvMatrix<Type>& matrix
)
{
if (this->manipulatedMatrix())
{
return;
}
const oversetFvPatch& ovp = this->oversetPatch_;
if (ovp.master())
{
if (massCorrection_ || (debug & 2))
{
storeFringeCoefficients(matrix);
}
// Trigger interpolation
const fvMesh& mesh = this->internalField().mesh();
const word& fldName = this->internalField().name();
// Try to find out if the solve routine comes from the unadapted mesh
// TBD. This should be cleaner.
if (&mesh.lduAddr() == &mesh.fvMesh::lduAddr())
{
return;
}
if (debug)
{
Pout<< FUNCTION_NAME << " field:" << fldName
<< " patch:" << ovp.name() << endl;
}
// Calculate stabilised diagonal as normalisation for interpolation
const scalarField norm
(
dynamic_cast<const oversetFvMeshBase&>(mesh).normalisation(matrix)
);
const cellCellStencil& overlap = Stencil::New(mesh);
const labelUList& types = overlap.cellTypes();
const labelListList& stencil = overlap.cellStencil();
dynamic_cast<const oversetFvMeshBase&>(mesh).addInterpolation
(
matrix,
norm,
this->setHoleCellValue_,
this->holeCellValue_
);
if (debug)
{
pointField allCs(mesh.cellCentres());
const mapDistribute& map = overlap.cellInterpolationMap();
map.distribute(allCs, false, UPstream::msgType()+1);
// Make sure we don't interpolate from a hole
scalarField marker(this->primitiveField().size(), 0);
forAll(types, celli)
{
if (types[celli] == cellCellStencil::HOLE)
{
marker[celli] = 1.0;
}
}
cellCellStencil::interpolate
(
marker,
mesh,
overlap,
overlap.cellInterpolationWeights()
);
forAll(marker, celli)
{
if
(
types[celli] == cellCellStencil::INTERPOLATED
&& marker[celli] > SMALL
)
{
WarningInFunction
<< " field:" << fldName
<< " patch:" << ovp.name()
<< " found:" << celli
<< " at:" << mesh.cellCentres()[celli]
<< " donorSlots:" << stencil[celli]
<< " at:"
<< UIndirectList<point>(allCs, stencil[celli])
<< " amount-of-hole:" << marker[celli]
<< endl;
}
}
// Make sure we don't have matrix coefficients for interpolated
// or hole cells
const lduAddressing& addr = mesh.lduAddr();
const labelUList& upperAddr = addr.upperAddr();
const labelUList& lowerAddr = addr.lowerAddr();
const scalarField& lower = matrix.lower();
const scalarField& upper = matrix.upper();
forAll(lowerAddr, facei)
{
const label l = lowerAddr[facei];
const bool lHole = (types[l] == cellCellStencil::HOLE);
const label u = upperAddr[facei];
const bool uHole = (types[u] == cellCellStencil::HOLE);
if
(
(lHole && upper[facei] != 0.0)
|| (uHole && lower[facei] != 0.0)
)
{
FatalErrorInFunction
<< "Hole-neighbouring face:" << facei
<< " lower:" << l
<< " type:" << types[l]
<< " coeff:" << lower[facei]
<< " upper:" << upperAddr[facei]
<< " type:" << types[u]
<< " coeff:" << upper[facei]
<< exit(FatalError);
}
// Empty donor list: treat like hole but still allow to
// influence neighbouring domains
const bool lEmpty =
(
types[l] == cellCellStencil::INTERPOLATED
&& stencil[l].empty()
);
const bool uEmpty =
(
types[u] == cellCellStencil::INTERPOLATED
&& stencil[u].empty()
);
if
(
(lEmpty && upper[facei] != 0.0)
|| (uEmpty && lower[facei] != 0.0)
)
{
FatalErrorInFunction
<< "Still connected face:" << facei << " lower:" << l
<< " type:" << types[l]
<< " coeff:" << lower[facei]
<< " upper:" << u
<< " type:" << types[u]
<< " coeff:" << upper[facei]
<< exit(FatalError);
}
}
forAll(matrix.internalCoeffs(), patchi)
{
const labelUList& fc = addr.patchAddr(patchi);
const Field<Type>& bouCoeffs = matrix.boundaryCoeffs()[patchi];
forAll(fc, i)
{
const label celli = fc[i];
const bool lHole = (types[celli] == cellCellStencil::HOLE);
if (lHole && bouCoeffs[i] != pTraits<Type>::zero)
{
FatalErrorInFunction
<< "Patch:" << patchi
<< " patchFace:" << i
<< " lower:" << celli
<< " type:" << types[celli]
<< " bouCoeff:" << bouCoeffs[i]
<< exit(FatalError);
}
// Check whether this is influenced by neighbouring domains
const bool lEmpty =
(
types[celli] == cellCellStencil::INTERPOLATED
&& stencil[celli].empty()
);
if (lEmpty && bouCoeffs[i] != pTraits<Type>::zero)
{
FatalErrorInFunction
<< "Patch:" << patchi
<< " patchFace:" << i
<< " lower:" << celli
<< " type:" << types[celli]
<< " bouCoeff:" << bouCoeffs[i]
<< exit(FatalError);
}
}
}
// Make sure that diagonal is non-zero. Note: should add
// boundaryCoeff ...
const FieldField<Field, Type>& internalCoeffs =
matrix.internalCoeffs();
for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; ++cmpt)
{
// Replacement for m.addBoundaryDiag(norm, cmpt);
scalarField diag(matrix.diag());
forAll(internalCoeffs, patchi)
{
const labelUList& fc = addr.patchAddr(patchi);
const Field<Type>& intCoeffs = internalCoeffs[patchi];
const scalarField cmptCoeffs(intCoeffs.component(cmpt));
forAll(fc, i)
{
diag[fc[i]] += cmptCoeffs[i];
}
}
forAll(diag, celli)
{
if (mag(diag[celli]) < SMALL)
{
FatalErrorInFunction
<< "Patch:" << ovp.name()
<< " cell:" << celli
<< " at:" << mesh.cellCentres()[celli]
<< " diag:" << diag[celli]
<< exit(FatalError);
}
}
}
}
}
fvPatchField<Type>::manipulateMatrix(matrix);
}
template<class Type>
void Foam::oversetFvPatchField<Type>::initInterfaceMatrixUpdate
(
solveScalarField& result,
const bool add,
const lduAddressing& lduAddr,
const label interfacei,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const
{
}
template<class Type>
void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
(
solveScalarField& result,
const bool add,
const lduAddressing& lduAddr,
const label patchId,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
) const
{
scalarField& psi = const_cast<scalarField&>(psiInternal);
if (massCorrection_ && this->oversetPatch_.master())
{
adjustPsi(psi, lduAddr, result);
}
}
template<class Type>
void Foam::oversetFvPatchField<Type>::write(Ostream& os) const
{
zeroGradientFvPatchField<Type>::write(os);
// Make sure to write the value for ease of postprocessing.
this->writeEntry("value", os);
coupledFvPatchField<Type>::write(os);
if (this->setHoleCellValue_)
{
os.writeEntry("setHoleCellValue", setHoleCellValue_);
os.writeEntry("holeCellValue", holeCellValue_);
os.writeEntryIfDifferent
(
"interpolateHoleCellValue",
false,
interpolateHoleCellValue_
);
}
os.writeEntryIfDifferent
(
"massCorrection",
false,
massCorrection_
);
os.writeEntryIfDifferent
(
"zone",
-1,
zoneId_
);
}

174
src/overset/oversetPolyPatch/oversetFvPatchField.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -42,7 +42,8 @@ SourceFiles
#define oversetFvPatchField_H
#include "oversetFvPatch.H"
#include "zeroGradientFvPatchField.H"
#include "coupledFvPatchFields.H"
#include "oversetLduInterfaceField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -56,11 +57,12 @@ namespace Foam
template<class Type>
class oversetFvPatchField
:
public zeroGradientFvPatchField<Type>
public oversetLduInterfaceField,
public coupledFvPatchField<Type>
{
protected:
// Protected data
// Protected Data
//- Local reference cast into the overset patch
const oversetFvPatch& oversetPatch_;
@ -68,6 +70,32 @@ protected:
//- Master patch ID
mutable label masterPatchID_;
// Hole cell controls
//- Set hole cell value
const bool setHoleCellValue_;
//- Correct mass
const bool massCorrection_;
//- Interpolate hole cell value (from nearby non-hole cell)
const bool interpolateHoleCellValue_;
//- Hole cell value
const Type holeCellValue_;
//- Fringe upper coefficients
mutable scalarField fringeUpperCoeffs_;
//- Fringe lower coefficients
mutable scalarField fringeLowerCoeffs_;
//- Fringe faces
mutable labelField fringeFaces_;
//- Zone to sum flux for mass conservation
label zoneId_;
public:
@ -133,13 +161,147 @@ public:
}
// Member functions
// Member Functions
// Evaluation functions
// Coupled and adjust flux
//- Return neighbour field. Dummy
virtual tmp<Field<Type> > patchNeighbourField() const;
//- Adjust psi for mass correction. Requires storeFringeCoefficients
// to have been called before
void adjustPsi
(
Field<scalar>& psi,
const lduAddressing& lduAddr,
solveScalarField& result
) const;
//- Store fringe coefficients and faces
void storeFringeCoefficients(const fvMatrix<Type>& matrix);
//- Calculate patch flux (helper function). Requires
// storeFringeCoefficients to have been called before
void fringeFlux
(
const fvMatrix<Type>& m,
const surfaceScalarField& phi
) const;
// Evaluation
//- Initialise the evaluation of the patch field
virtual void initEvaluate(const Pstream::commsTypes commsType);
//- Return the matrix diagonal coefficients corresponding to the
//- evaluation of the value of this patchField with given weights
virtual tmp<Field<Type> > valueInternalCoeffs
(
const tmp<scalarField>&
) const
{
return tmp<Field<Type>>
(
new Field<Type>(this->size(), Zero)
);
}
//- Return the matrix source coefficients corresponding to the
//- evaluation of the value of this patchField with given weights
virtual tmp<Field<Type> > valueBoundaryCoeffs
(
const tmp<scalarField>&
) const
{
return tmp<Field<Type>>
(
new Field<Type>(this->size(), Zero)
);
}
//- Return the matrix diagonal coefficients corresponding to the
//- evaluation of the gradient of this patchField
tmp<Field<Type> > gradientInternalCoeffs() const
{
return tmp<Field<Type>>
(
new Field<Type>(this->size(), Zero)
);
}
//- Return the matrix source coefficients corresponding to the
//- evaluation of the gradient of this patchField
tmp<Field<Type> > gradientBoundaryCoeffs() const
{
return tmp<Field<Type>>
(
new Field<Type>(this->size(), Zero)
);
}
//- Manipulate matrix
virtual void manipulateMatrix(fvMatrix<Type>& matrix);
// Coupled interface functionality
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
solveScalarField& result,
const bool add,
const lduAddressing& lduAddr,
const label patchId,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
) const;
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
Field<Type>&,
const bool add,
const lduAddressing& lduAddr,
const label patchId,
const Field<Type>&,
const scalarField&,
const Pstream::commsTypes commsType
) const
{
NotImplemented;
}
//- Initialise neighbour matrix update
virtual void initInterfaceMatrixUpdate
(
Field<Type>&,
const bool add,
const lduAddressing&,
const label interfacei,
const Field<Type>&,
const scalarField&,
const Pstream::commsTypes commsType
) const
{
NotImplemented;
}
virtual void initInterfaceMatrixUpdate
(
solveScalarField& result,
const bool add,
const lduAddressing&,
const label interfacei,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const;
// I-O
//- Write