Commit ddde3308 authored by mattijs's avatar mattijs
Browse files

ENH: overset: new solvers, new stencil

parent 74b557d5
......@@ -7,8 +7,8 @@ surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", constrainHbyA(rAU*UEqn.H(), U, p));
//mesh.interpolate(HbyA);
volVectorField HbyA("HbyA", U);
HbyA = constrainHbyA(rAU*UEqn.H(), U, p);
if (pimple.nCorrPISO() <= 1)
{
......
......@@ -16,6 +16,9 @@
fvOptions.constrain(UEqn);
solve(UEqn == -cellMask*fvc::grad(p));
if (simple.momentumPredictor())
{
solve(UEqn == -cellMask*fvc::grad(p));
}
fvOptions.correct(U);
......@@ -5,7 +5,7 @@
surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = constrainHbyA(rAU*UEqn.H(), U, p);
HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
//mesh.interpolate(HbyA);
if (massFluxInterpolation)
......
......@@ -138,3 +138,4 @@ volScalarField pDivU
mesh,
dimensionedScalar("pDivU", p.dimensions()/dimTime, 0)
);
#include "readTimeControls.H"
correctPhi = pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
correctPhi = pimple.dict().lookupOrDefault<Switch>("correctPhi", false);
checkMeshCourantNo =
pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);
......
......@@ -32,6 +32,7 @@ bodies
type sphere;
mass 1;
radius 0.05;
centreOfMass (0 10 0);
transform (1 0 0 0 1 0 0 0 1) (0 -1 0);
mergeWith hinge;
}
......
......@@ -14,6 +14,7 @@ bodies
parent root;
mass 6e4;
radius 0.01;
centreOfMass (0 0 0);
transform (1 0 0 0 1 0 0 0 1) (0 0 0);
joint
{
......@@ -37,6 +38,7 @@ bodies
parent M;
mass 6e3;
radius 0.01;
centreOfMass (0 0 0);
transform (1 0 0 0 1 0 0 0 1) (0 -5 0);
mergeWith M;
}
......
......@@ -233,10 +233,10 @@ public:
void reorder(const labelUList&, const bool validBoundary);
//- writeData member function required by regIOobject
bool writeData(Ostream&) const;
virtual bool writeData(Ostream&) const;
//- Write using given format, version and form uncompression
bool writeObject
virtual bool writeObject
(
IOstream::streamFormat fmt,
IOstream::versionNumber ver,
......
......@@ -101,7 +101,18 @@ bool Foam::patchDistMethods::Poisson::correct
volVectorField gradyPsi(fvc::grad(yPsi));
volScalarField magGradyPsi(mag(gradyPsi));
y = sqrt(magSqr(gradyPsi) + 2*yPsi) - magGradyPsi;
// Need to stabilise the y for overset meshes since the holed cells
// keep the initial value (0.0) so the gradient of that will be
// zero as well. Turbulence models do not like zero wall distance.
y = max
(
sqrt(magSqr(gradyPsi) + 2*yPsi) - magGradyPsi,
dimensionedScalar("smallY", dimLength, SMALL)
);
// For overset: enforce smooth y field (yPsi is smooth, magGradyPsi is
// not)
mesh_.interpolate(y);
// Need to stabilise the y for overset meshes since the holed cells
// keep the initial value (0.0) so the gradient of that will be
......
......@@ -6,6 +6,7 @@ cellCellStencil/inverseDistance/waveMethod.C
cellCellStencil/inverseDistance/meshToMeshData.C
cellCellStencil/trackingInverseDistance/voxelMeshSearch.C
cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
cellCellStencil/leastSquares/leastSquaresCellCellStencil.C
dynamicOversetFvMesh/dynamicOversetFvMesh.C
......@@ -21,6 +22,9 @@ oversetPolyPatch/oversetGAMGInterfaceField.C
oversetPolyPatch/oversetPointPatch.C
oversetPolyPatch/oversetPointPatchFields.C
oversetPolyPatch/semiImplicitOversetFvPatchFields.C
oversetPolyPatch/semiImplicitOversetGAMGInterfaceField.C
oversetAdjustPhi/oversetAdjustPhi.C
regionsToCell/regionsToCell.C
......
......@@ -94,22 +94,22 @@ Foam::cellCellStencil::~cellCellStencil()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::labelIOList& Foam::cellCellStencil::zoneID() const
const Foam::labelIOList& Foam::cellCellStencil::zoneID(const fvMesh& mesh)
{
if (!mesh_.foundObject<labelIOList>("zoneID"))
if (!mesh.foundObject<labelIOList>("zoneID"))
{
labelIOList* zoneIDPtr = new labelIOList
(
IOobject
(
"zoneID",
mesh_.facesInstance(),
mesh.facesInstance(),
polyMesh::meshSubDir,
mesh_,
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_.nCells()
mesh.nCells()
);
labelIOList& zoneID = *zoneIDPtr;
......@@ -118,13 +118,13 @@ const Foam::labelIOList& Foam::cellCellStencil::zoneID() const
IOobject
(
"zoneID",
mesh_.time().findInstance(mesh_.dbDir(), "zoneID"),
mesh_,
mesh.time().findInstance(mesh.dbDir(), "zoneID"),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
),
mesh_
mesh
);
forAll(volZoneID, cellI)
{
......@@ -133,7 +133,7 @@ const Foam::labelIOList& Foam::cellCellStencil::zoneID() const
zoneIDPtr->store();
}
return mesh_.lookupObject<labelIOList>("zoneID");
return mesh.lookupObject<labelIOList>("zoneID");
}
......
......@@ -168,12 +168,27 @@ public:
// calculated, 1 = use interpolated)
virtual const scalarList& cellInterpolationWeight() const = 0;
//- Calculate weights for a single acceptor
virtual void stencilWeights
(
const point& sample,
const pointList& donorCcs,
scalarList& weights
) const = 0;
//- Helper: is stencil fully local
bool localStencil(const labelUList&) const;
//- Helper: get reference to registered zoneID. Loads volScalarField
// if not registered.
const labelIOList& zoneID() const;
static const labelIOList& zoneID(const fvMesh&);
//- Helper: get reference to registered zoneID. Loads volScalarField
// if not registered.
const labelIOList& zoneID() const
{
return zoneID(mesh_);
}
//- Helper: create cell-cell addressing in global numbering
static void globalCellCells
......
......@@ -154,6 +154,17 @@ public:
{
return stencilPtr_().cellInterpolationWeight();
}
//- Calculate weights for a single acceptor
virtual void stencilWeights
(
const point& sample,
const pointList& donorCcs,
scalarList& weights
) const
{
stencilPtr_().stencilWeights(sample, donorCcs, weights);
}
};
......
......@@ -1100,4 +1100,39 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
}
void Foam::cellCellStencils::cellVolumeWeight::stencilWeights
(
const point& sample,
const pointList& donorCcs,
scalarList& weights
) const
{
// Inverse-distance weighting
weights.setSize(donorCcs.size());
scalar sum = 0.0;
forAll(donorCcs, i)
{
scalar d = mag(sample-donorCcs[i]);
if (d > ROOTVSMALL)
{
weights[i] = 1.0/d;
sum += weights[i];
}
else
{
// Short circuit
weights = 0.0;
weights[i] = 1.0;
return;
}
}
forAll(weights, i)
{
weights[i] /= sum;
}
}
// ************************************************************************* //
......@@ -227,6 +227,15 @@ public:
{
return cellInterpolationWeight_;
}
//- Calculate inverse distance weights for a single acceptor. Revert
// to inverse distance (so not consistent with volume overlap!)
virtual void stencilWeights
(
const point& sample,
const pointList& donorCcs,
scalarList& weights
) const;
};
......
......@@ -155,6 +155,7 @@ void Foam::cellCellStencils::inverseDistance::fill
void Foam::cellCellStencils::inverseDistance::markBoundaries
(
const fvMesh& mesh,
const vector& smallVec,
const boundBox& bb,
const labelVector& nDivs,
......@@ -188,6 +189,9 @@ void Foam::cellCellStencils::inverseDistance::markBoundaries
// Mark in voxel mesh
boundBox faceBb(pp.points(), pp[i], false);
faceBb.min() -= smallVec;
faceBb.max() += smallVec;
if (bb.overlaps(faceBb))
{
fill(patchTypes, bb, nDivs, faceBb, patchCellType::PATCH);
......@@ -213,6 +217,9 @@ void Foam::cellCellStencils::inverseDistance::markBoundaries
// Mark in voxel mesh
boundBox faceBb(pp.points(), pp[i], false);
faceBb.min() -= smallVec;
faceBb.max() += smallVec;
if (bb.overlaps(faceBb))
{
fill(patchTypes, bb, nDivs, faceBb, patchCellType::OVERSET);
......@@ -333,6 +340,10 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
forAll(tgtCellMap, tgtCelli)
{
label celli = tgtCellMap[tgtCelli];
treeBoundBox cBb(cellBb(mesh_, celli));
cBb.min() -= smallVec_;
cBb.max() += smallVec_;
if
(
overlaps
......@@ -340,7 +351,7 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
srcPatchBb,
zoneDivs,
srcPatchTypes,
cellBb(mesh_, celli),
cBb,
patchCellType::PATCH
)
)
......@@ -399,6 +410,9 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
forAll(tgtCellMap, tgtCelli)
{
label celli = tgtCellMap[tgtCelli];
treeBoundBox cBb(cellBb(mesh_, celli));
cBb.min() -= smallVec_;
cBb.max() += smallVec_;
if
(
overlaps
......@@ -406,7 +420,7 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
srcPatchBb,
zoneDivs,
srcPatchTypes,
cellBb(mesh_, celli),
cBb,
patchCellType::PATCH
)
)
......@@ -501,7 +515,9 @@ void Foam::cellCellStencils::inverseDistance::markDonors
label celli = tgtCellMap[tgtCelli];
if (allStencil[celli].empty())
{
const treeBoundBox subBb(cellBb(mesh_, celli));
treeBoundBox subBb(cellBb(mesh_, celli));
subBb.min() -= smallVec_;
subBb.max() += smallVec_;
forAll(srcOverlapProcs, i)
{
......@@ -1194,9 +1210,17 @@ void Foam::cellCellStencils::inverseDistance::walkFront
{
if (isA<oversetFvPatch>(fvm[patchI]))
{
forAll(fvm[patchI], i)
const labelList& fc = fvm[patchI].faceCells();
forAll(fc, i)
{
isFront[fvm[patchI].start()+i] = true;
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[fvm[patchI].start()+i] = true;
}
}
}
}
......@@ -1353,12 +1377,12 @@ void Foam::cellCellStencils::inverseDistance::walkFront
}
void Foam::cellCellStencils::inverseDistance::calcStencilWeights
void Foam::cellCellStencils::inverseDistance::stencilWeights
(
const point& sample,
const pointList& donorCcs,
scalarList& weights
)
) const
{
// Inverse-distance weighting
......@@ -1500,7 +1524,7 @@ void Foam::cellCellStencils::inverseDistance::createStencil
{
label cellI = donorCells[i];
const pointList& donorCentres = donorCellCentres[cellI];
calcStencilWeights
stencilWeights
(
samples[cellI],
donorCentres,
......@@ -1565,6 +1589,7 @@ Foam::cellCellStencils::inverseDistance::inverseDistance
:
cellCellStencil(mesh),
dict_(dict),
smallVec_(vector::zero),
cellTypes_(labelList(mesh.nCells(), CALCULATED)),
interpolationCells_(0),
cellInterpolationMap_(),
......@@ -1605,6 +1630,9 @@ bool Foam::cellCellStencils::inverseDistance::update()
{
scalar layerRelax(dict_.lookupOrDefault("layerRelax", 1.0));
scalar tol = dict_.lookupOrDefault("tolerance", 1e-10);
smallVec_ = mesh_.bounds().span()*tol;
const labelIOList& zoneID = this->zoneID();
label nZones = gMax(zoneID)+1;
......@@ -1745,6 +1773,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
markBoundaries
(
meshParts[zoneI].subMesh(),
smallVec_,
patchBb[zoneI][Pstream::myProcNo()],
patchDivisions[zoneI],
......@@ -1758,7 +1787,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
// Print a bit
{
Info<< typeName << " : detected " << nZones
Info<< type() << " : detected " << nZones
<< " mesh regions" << endl;
Info<< incrIndent;
forAll(nCellsPerZone, zoneI)
......@@ -1862,6 +1891,16 @@ bool Foam::cellCellStencils::inverseDistance::update()
}
else
{
// If there are no donors we can either set the
// acceptor cell to 'hole' i.e. nothing gets calculated
// if the acceptor cells go outside the donor meshes or
// we can set it to 'calculated' to have something
// like an acmi type behaviour where only covered
// acceptor cell use the interpolation and non-covered
// become calculated. Uncomment below line. Note: this
// should be switchable?
//allCellTypes[cellI] = CALCULATED;
allCellTypes[cellI] = HOLE;
}
}
......@@ -1915,7 +1954,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
// Dump stencil
mkDir(mesh_.time().timePath());
OBJstream str(mesh_.time().timePath()/"injectionStencil.obj");
Pout<< typeName << " : dumping injectionStencil to "
Pout<< type() << " : dumping injectionStencil to "
<< str.name() << endl;
pointField cc(mesh_.cellCentres());
cellInterpolationMap_.distribute(cc);
......@@ -1973,7 +2012,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
// Dump stencil
mkDir(mesh_.time().timePath());
OBJstream str(mesh_.time().timePath()/"stencil.obj");
Pout<< typeName << " : dumping to " << str.name() << endl;
Pout<< type() << " : dumping to " << str.name() << endl;
pointField cc(mesh_.cellCentres());
cellInterpolationMap_.distribute(cc);
......
......@@ -74,6 +74,9 @@ protected:
//- Dictionary of motion control parameters
const dictionary dict_;
//- Small perturbation vector for geometric tests
vector smallVec_;
//- Per cell the cell type
labelList cellTypes_;
......@@ -143,6 +146,8 @@ protected:
static void markBoundaries
(
const fvMesh& mesh,
const vector& smallVec,
const boundBox& bb,
const labelVector& nDivs,
PackedList<2>& patchTypes,
......@@ -237,15 +242,6 @@ protected:
scalarField& allWeight
) const;
//- Calculate inverse distance weights
static void calcStencilWeights
(
const point& sample,
const pointList& donorCcs,
scalarList& weights
);
//- Create stencil starting from the donor containing the acceptor
virtual void createStencil(const globalIndex&);
......@@ -319,6 +315,14 @@ public:
{
return cellInterpolationWeight_;
}
//- Calculate inverse distance weights for a single acceptor
virtual void stencilWeights
(
const point& sample,
const pointList& donorCcs,
scalarList& weights
) const;
};
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "leastSquaresCellCellStencil.H"
#include "addToRunTimeSelectionTable.H"
#include "SVD.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace cellCellStencils
{
defineTypeNameAndDebug(leastSquares, 0);
addToRunTimeSelectionTable(cellCellStencil, leastSquares, mesh);
}
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::cellCellStencils::leastSquares::stencilWeights
(
const point& sample,
const pointList& donorCcs,
scalarList& weights
) const
{
// Implicit least squares weighting
// Number of donors
label nD = donorCcs.size();
weights.setSize(nD);
// List for distance vectors and LSQ weights
List<vector> d(nD);
scalarList LSQw(nD);
// Sum of weights
scalar W = 0;
// Sum of weighted distance vectors
vector dw = vector::zero;
RectangularMatrix<scalar> A(nD, 3);
bool shortC = false;
// Compute distance vectors and fill rectangular matrix
forAll(donorCcs, j)
{
// Neighbour weights
d[j] = donorCcs[j] - sample;