Skip to content
GitLab
Commit 2a1f28a3 authored by henry's avatar henry
Browse files

Parallelized stencil-based fitting.

parent 809bbd5a
Hide whitespace changes
Inline Side-by-side
src/OpenFOAM/meshes/MeshObject/MeshObject.H
View file @ 2a1f28a3
......@@ -104,12 +104,12 @@ public:
);
// Destructor
static bool Delete(const Mesh& mesh);
// Destructors
virtual ~MeshObject();
static bool Delete(const Mesh& mesh);
// Member Functions
......
src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H
View file @ 2a1f28a3
......@@ -131,7 +131,6 @@ public:
const GeometricField<Type, fvPatchField, volMesh>& fld,
const List<List<scalar> >& stencilWeights
);
};
......
src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.H
View file @ 2a1f28a3
......@@ -70,10 +70,9 @@ public:
{}
// Destructor
virtual ~centredCFCStencilObject()
{}
//- Destructor
virtual ~centredCFCStencilObject()
{}
};
......
src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C
View file @ 2a1f28a3
......@@ -128,7 +128,6 @@ Foam::extendedStencil::weightedSum
// Boundaries. Either constrained or calculated so assign value
// directly (instead of nicely using operator==)
/*
typename GeometricField<Type, fvsPatchField, surfaceMesh>::
GeometricBoundaryField& bSfCorr = sf.boundaryField();
......@@ -136,22 +135,24 @@ Foam::extendedStencil::weightedSum
{
fvsPatchField<Type>& pSfCorr = bSfCorr[patchi];
label faceI = pSfCorr.patch().patch().start();
forAll(pSfCorr, i)
if (pSfCorr.coupled())
{
const List<Type>& stField = stencilFld[faceI];
const List<scalar>& stWeight = stencilWeights[faceI];
label faceI = pSfCorr.patch().patch().start();
forAll(stField, j)
forAll(pSfCorr, i)
{
pSfCorr[i] += stField[j]*stWeight[j];
}
const List<Type>& stField = stencilFld[faceI];
const List<scalar>& stWeight = stencilWeights[faceI];
faceI++;
forAll(stField, j)
{
pSfCorr[i] += stField[j]*stWeight[j];
}
faceI++;
}
}
}
*/
return tsfCorr;
}
......
src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C
View file @ 2a1f28a3
......@@ -102,32 +102,35 @@ Foam::extendedUpwindStencil::weightedSum
{
fvsPatchField<Type>& pSfCorr = bSfCorr[patchi];
label faceI = pSfCorr.patch().patch().start();
forAll(pSfCorr, i)
if (pSfCorr.coupled())
{
if (phi[faceI] > 0)
{
// Flux out of owner. Use upwind (= owner side) stencil.
const List<Type>& stField = ownFld[faceI];
const List<scalar>& stWeight = ownWeights[faceI];
label faceI = pSfCorr.patch().patch().start();
forAll(stField, j)
forAll(pSfCorr, i)
{
if (phi[faceI] > 0)
{
pSfCorr[i] += stField[j]*stWeight[j];
// Flux out of owner. Use upwind (= owner side) stencil.
const List<Type>& stField = ownFld[faceI];
const List<scalar>& stWeight = ownWeights[faceI];
forAll(stField, j)
{
pSfCorr[i] += stField[j]*stWeight[j];
}
}
}
else
{
const List<Type>& stField = neiFld[faceI];
const List<scalar>& stWeight = neiWeights[faceI];
forAll(stField, j)
else
{
pSfCorr[i] += stField[j]*stWeight[j];
const List<Type>& stField = neiFld[faceI];
const List<scalar>& stWeight = neiWeights[faceI];
forAll(stField, j)
{
pSfCorr[i] += stField[j]*stWeight[j];
}
}
faceI++;
}
faceI++;
}
}
......
src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.H
View file @ 2a1f28a3
......@@ -68,7 +68,6 @@ public:
//- Construct from mesh
explicit cellFaceCellStencil(const polyMesh& mesh);
};
......
src/finiteVolume/fvMesh/fvMesh.C
View file @ 2a1f28a3
......@@ -42,9 +42,9 @@ License
#include "extendedLeastSquaresVectors.H"
#include "extendedLeastSquaresVectors.H"
#include "leastSquaresVectors.H"
//#include "linearFitData.H"
#include "quadraticFitData.H"
//#include "quadraticFitSnGradData.H"
#include "CentredFitData.H"
#include "linearFitPolynomial.H"
#include "quadraticLinearFitPolynomial.H"
#include "skewCorrectionVectors.H"
#include "centredCECStencilObject.H"
......@@ -89,15 +89,13 @@ void Foam::fvMesh::clearGeom()
// Mesh motion flux cannot be deleted here because the old-time flux
// needs to be saved.
// Things geometry dependent that are not updated.
volPointInterpolation::Delete(*this);
extendedLeastSquaresVectors::Delete(*this);
extendedLeastSquaresVectors::Delete(*this);
leastSquaresVectors::Delete(*this);
//linearFitData::Delete(*this);
quadraticFitData::Delete(*this);
//quadraticFitSnGradData::Delete(*this);
CentredFitData<linearFitPolynomial>::Delete(*this);
CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
skewCorrectionVectors::Delete(*this);
}
......@@ -112,9 +110,8 @@ void Foam::fvMesh::clearAddressing()
extendedLeastSquaresVectors::Delete(*this);
extendedLeastSquaresVectors::Delete(*this);
leastSquaresVectors::Delete(*this);
//linearFitData::Delete(*this);
quadraticFitData::Delete(*this);
//quadraticFitSnGradData::Delete(*this);
CentredFitData<linearFitPolynomial>::Delete(*this);
CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
skewCorrectionVectors::Delete(*this);
centredCECStencilObject::Delete(*this);
......@@ -364,7 +361,6 @@ Foam::fvMesh::~fvMesh()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Helper function for construction from pieces
void Foam::fvMesh::addFvPatches
(
const List<polyPatch*> & p,
......@@ -552,6 +548,30 @@ void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
}
// Temporary helper function to call move points on
// MeshObjects
template<class Type>
void MeshObjectMovePoints(const Foam::fvMesh& mesh)
{
if
(
mesh.db().objectRegistry::foundObject<Type>
(
Type::typeName
)
)
{
const_cast<Type&>
(
mesh.db().objectRegistry::lookupObject<Type>
(
Type::typeName
)
).movePoints();
}
}
Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
{
// Grab old time volumes if the time has been incremented
......@@ -642,133 +662,12 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
// Hack until proper callbacks. Below are all the fvMesh MeshObjects with a
// movePoints function.
// volPointInterpolation
if
(
db().objectRegistry::foundObject<volPointInterpolation>
(
volPointInterpolation::typeName
)
)
{
const_cast<volPointInterpolation&>
(
db().objectRegistry::lookupObject<volPointInterpolation>
(
volPointInterpolation::typeName
)
).movePoints();
}
// extendedLeastSquaresVectors
if
(
db().objectRegistry::foundObject<extendedLeastSquaresVectors>
(
extendedLeastSquaresVectors::typeName
)
)
{
const_cast<extendedLeastSquaresVectors&>
(
db().objectRegistry::lookupObject<extendedLeastSquaresVectors>
(
extendedLeastSquaresVectors::typeName
)
).movePoints();
}
// leastSquaresVectors
if
(
db().objectRegistry::foundObject<leastSquaresVectors>
(
leastSquaresVectors::typeName
)
)
{
const_cast<leastSquaresVectors&>
(
db().objectRegistry::lookupObject<leastSquaresVectors>
(
leastSquaresVectors::typeName
)
).movePoints();
}
//// linearFitData
//if
//(
// db().objectRegistry::foundObject<linearFitData>
// (
// linearFitData::typeName
// )
//)
//{
// const_cast<linearFitData&>
// (
// db().objectRegistry::lookupObject<linearFitData>
// (
// linearFitData::typeName
// )
// ).movePoints();
//}
// quadraticFitData
if
(
db().objectRegistry::foundObject<quadraticFitData>
(
quadraticFitData::typeName
)
)
{
const_cast<quadraticFitData&>
(
db().objectRegistry::lookupObject<quadraticFitData>
(
quadraticFitData::typeName
)
).movePoints();
}
//// quadraticFitSnGradData
//if
//(
// db().objectRegistry::foundObject<quadraticFitSnGradData>
// (
// quadraticFitSnGradData::typeName
// )
//)
//{
// const_cast<quadraticFitSnGradData&>
// (
// db().objectRegistry::lookupObject<quadraticFitSnGradData>
// (
// quadraticFitSnGradData::typeName
// )
// ).movePoints();
//}
// skewCorrectionVectors
if
(
db().objectRegistry::foundObject<skewCorrectionVectors>
(
skewCorrectionVectors::typeName
)
)
{
const_cast<skewCorrectionVectors&>
(
db().objectRegistry::lookupObject<skewCorrectionVectors>
(
skewCorrectionVectors::typeName
)
).movePoints();
}
MeshObjectMovePoints<volPointInterpolation>(*this);
MeshObjectMovePoints<extendedLeastSquaresVectors>(*this);
MeshObjectMovePoints<leastSquaresVectors>(*this);
MeshObjectMovePoints<CentredFitData<linearFitPolynomial> >(*this);
MeshObjectMovePoints<CentredFitData<quadraticLinearFitPolynomial> >(*this);
MeshObjectMovePoints<skewCorrectionVectors>(*this);
return tsweptVols;
}
......
src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C
View file @ 2a1f28a3
......@@ -43,6 +43,7 @@ Foam::CentredFitData<Polynomial>::CentredFitData
)
:
MeshObject<fvMesh, CentredFitData<Polynomial> >(mesh),
stencil_(stencil),
linearLimitFactor_(linearLimitFactor),
centralWeight_(centralWeight),
# ifdef SPHERICAL_GEOMETRY
......@@ -51,62 +52,29 @@ Foam::CentredFitData<Polynomial>::CentredFitData
dim_(mesh.nGeometricD()),
# endif
minSize_(Polynomial::nTerms(dim_)),
coeffs_(mesh.nInternalFaces())
coeffs_(mesh.nFaces())
{
if (debug)
{
Info<< "Contructing CentredFitData<Polynomial>" << endl;
}
// check input
if (centralWeight_ < 1 - SMALL)
// Check input
if (linearLimitFactor > 1)
{
FatalErrorIn("CentredFitData<Polynomial>::CentredFitData")
<< "centralWeight requested = " << centralWeight_
<< "linearLimitFactor requested = " << linearLimitFactor
<< " should not be less than one"
<< exit(FatalError);
}
if (minSize_ == 0)
{
FatalErrorIn("CentredFitSnGradData")
<< " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError);
}
// store the polynomial size for each cell to write out
surfaceScalarField interpPolySize
(
IOobject
(
"CentredFitInterpPolySize",
"constant",
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("CentredFitInterpPolySize", dimless, scalar(0))
);
// Get the cell/face centres in stencil order.
// Centred face stencils no good for triangles of tets. Need bigger stencils
List<List<point> > stencilPoints(mesh.nFaces());
stencil.collectData(mesh.C(), stencilPoints);
// find the fit coefficients for every face in the mesh
for(label faci = 0; faci < mesh.nInternalFaces(); faci++)
{
interpPolySize[faci] = calcFit(stencilPoints[faci], faci);
}
calcFit();
if (debug)
{
Info<< "CentredFitData<Polynomial>::CentredFitData() :"
<< "Finished constructing polynomialFit data"
<< endl;
interpPolySize.write();
}
}
......@@ -120,15 +88,14 @@ void Foam::CentredFitData<Polynomial>::findFaceDirs
vector& jdir, // value changed in return
vector& kdir, // value changed in return
const fvMesh& mesh,
const label faci
const label facei
)
{
idir = mesh.Sf()[faci];
//idir = mesh.C()[mesh.neighbour()[faci]] - mesh.C()[mesh.owner()[faci]];
idir = mesh.faceAreas()[facei];
idir /= mag(idir);
# ifndef SPHERICAL_GEOMETRY
if (mesh.nGeometricD() <= 2) // find the normal direcion
if (mesh.nGeometricD() <= 2) // find the normal direction
{
if (mesh.directions()[0] == -1)
{
......@@ -143,14 +110,14 @@ void Foam::CentredFitData<Polynomial>::findFaceDirs
kdir = vector(0, 0, 1);
}
}
else // 3D so find a direction in the place of the face
else // 3D so find a direction in the plane of the face
{
const face& f = mesh.faces()[faci];
kdir = mesh.points()[f[0]] - mesh.points()[f[1]];
const face& f = mesh.faces()[facei];
kdir = mesh.points()[f[0]] - mesh.faceCentres()[facei];
}
# else
// Spherical geometry so kdir is the radial direction
kdir = mesh.Cf()[faci];
kdir = mesh.faceCentres()[facei];
# endif
if (mesh.nGeometricD() == 3)
......@@ -175,17 +142,58 @@ void Foam::CentredFitData<Polynomial>::findFaceDirs
}
template<class Polynomial>
void Foam::CentredFitData<Polynomial>::calcFit()
{
const fvMesh& mesh = this->mesh();
// Get the cell/face centres in stencil order.
// Centred face stencils no good for triangles of tets.
// Need bigger stencils
List<List<point> > stencilPoints(mesh.nFaces());
stencil_.collectData(mesh.C(), stencilPoints);
// find the fit coefficients for every face in the mesh
const surfaceScalarField& w = this->mesh().surfaceInterpolation::weights();
for(label facei = 0; facei < mesh.nInternalFaces(); facei++)
{
calcFit(stencilPoints[facei], w[facei], facei);
}
const surfaceScalarField::GeometricBoundaryField& bw = w.boundaryField();
forAll(bw, patchi)
{
const fvsPatchScalarField& pw = bw[patchi];
if (pw.coupled())
{
label facei = pw.patch().patch().start();
forAll(pw, i)
{
calcFit(stencilPoints[facei], pw[i], facei);
facei++;
}
}
}
}
template<class Polynomial>
Foam::label Foam::CentredFitData<Polynomial>::calcFit
(
const List<point>& C,
const label faci
const scalar wLin,
const label facei
)
{
vector idir(1,0,0);
vector jdir(0,1,0);
vector kdir(0,0,1);
findFaceDirs(idir, jdir, kdir, this->mesh(), faci);
findFaceDirs(idir, jdir, kdir, this->mesh(), facei);
// Setup the point weights
scalarList wts(C.size(), scalar(1));
......@@ -193,7 +201,7 @@ Foam::label Foam::CentredFitData<Polynomial>::calcFit
wts[1] = centralWeight_;
// Reference point
point p0 = this->mesh().faceCentres()[faci];
point p0 = this->mesh().faceCentres()[facei];
// p0 -> p vector in the face-local coordinate system
vector d;
......@@ -235,13 +243,10 @@ Foam::label Foam::CentredFitData<Polynomial>::calcFit
// Set the fit
label stencilSize = C.size();
coeffs_[faci].setSize(stencilSize);
coeffs_[facei].setSize(stencilSize);
scalarList singVals(minSize_);
label nSVDzeros = 0;
const GeometricField<scalar, fvsPatchField, surfaceMesh>& w =
this->mesh().surfaceInterpolation::weights();
bool goodFit = false;
for(int iIt = 0; iIt < 10 && !goodFit; iIt++)
{
......@@ -251,11 +256,11 @@ Foam::label Foam::CentredFitData<Polynomial>::calcFit
scalar fit1 = wts[1]*svd.VSinvUt()[0][1];
goodFit =
(mag(fit0 - w[faci]) < linearLimitFactor_*w[faci])
&& (mag(fit1 - (1 - w[faci])) < linearLimitFactor_*(1 - w[faci]));
(mag(fit0 - wLin) < linearLimitFactor_*wLin)
&& (mag(fit1 - (1 - wLin)) < linearLimitFactor_*(1 - wLin));
//scalar w0Err = fit0/w[faci];
//scalar w1Err = fit1/(1 - w[faci]);
//scalar w0Err = fit0/wLin;
//scalar w1Err = fit1/(1 - wLin);
//goodFit =
// (w0Err > linearLimitFactor_ && w0Err < (1 + linearLimitFactor_))
......@@ -263,12 +268,12 @@ Foam::label Foam::CentredFitData<Polynomial>::calcFit
if (goodFit)
{
coeffs_[faci][0] = fit0;
coeffs_[faci][1] = fit1;
coeffs_[facei][0] = fit0;
coeffs_[facei][1] = fit1;
for(label i=2; i<stencilSize; i++)
{
coeffs_[faci][i] = wts[i]*svd.VSinvUt()[0][i];
coeffs_[facei][i] = wts[i]*svd.VSinvUt()[0][i];
}
singVals = svd.S();
......@@ -300,22 +305,24 @@ Foam::label Foam::CentredFitData<Polynomial>::calcFit
// (
// min