From 2a1f28a3a8dfeffaf2971bff7533ea26a988fa06 Mon Sep 17 00:00:00 2001 From: henry Date: Thu, 23 Oct 2008 22:44:27 +0100 Subject: [PATCH] Parallelized stencil-based fitting. --- src/OpenFOAM/meshes/MeshObject/MeshObject.H | 6 +- .../fvMesh/extendedStencil/extendedStencil.H | 1 - .../centredCFCStencilObject.H | 7 +- .../extendedStencilTemplates.C | 23 +-- .../extendedUpwindStencilTemplates.C | 43 +++-- .../faceStencil/cellFaceCellStencil.H | 1 - src/finiteVolume/fvMesh/fvMesh.C | 175 ++++-------------- .../schemes/CentredFitScheme/CentredFitData.C | 154 +++++++-------- .../schemes/CentredFitScheme/CentredFitData.H | 15 +- 9 files changed, 172 insertions(+), 253 deletions(-) diff --git a/src/OpenFOAM/meshes/MeshObject/MeshObject.H b/src/OpenFOAM/meshes/MeshObject/MeshObject.H index 3bf11b5f6a..eac6096c79 100644 --- a/src/OpenFOAM/meshes/MeshObject/MeshObject.H +++ b/src/OpenFOAM/meshes/MeshObject/MeshObject.H @@ -104,12 +104,12 @@ public: ); - // Destructor - - static bool Delete(const Mesh& mesh); + // Destructors virtual ~MeshObject(); + static bool Delete(const Mesh& mesh); + // Member Functions diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H index 12d3179e62..25ba241f1e 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H @@ -131,7 +131,6 @@ public: const GeometricField& fld, const List >& stencilWeights ); - }; diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.H b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.H index 704f9a92d4..e36be07993 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.H +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.H @@ -70,10 +70,9 @@ public: {} - // Destructor - - virtual ~centredCFCStencilObject() - {} + //- Destructor + virtual ~centredCFCStencilObject() + {} }; diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C index 4aac55df75..ec024775d3 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C @@ -128,7 +128,6 @@ Foam::extendedStencil::weightedSum // Boundaries. Either constrained or calculated so assign value // directly (instead of nicely using operator==) - /* typename GeometricField:: GeometricBoundaryField& bSfCorr = sf.boundaryField(); @@ -136,22 +135,24 @@ Foam::extendedStencil::weightedSum { fvsPatchField& pSfCorr = bSfCorr[patchi]; - label faceI = pSfCorr.patch().patch().start(); - - forAll(pSfCorr, i) + if (pSfCorr.coupled()) { - const List& stField = stencilFld[faceI]; - const List& stWeight = stencilWeights[faceI]; + label faceI = pSfCorr.patch().patch().start(); - forAll(stField, j) + forAll(pSfCorr, i) { - pSfCorr[i] += stField[j]*stWeight[j]; - } + const List& stField = stencilFld[faceI]; + const List& stWeight = stencilWeights[faceI]; - faceI++; + forAll(stField, j) + { + pSfCorr[i] += stField[j]*stWeight[j]; + } + + faceI++; + } } } - */ return tsfCorr; } diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C b/src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C index 1cfa79eaa5..00cd58e17d 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C @@ -102,32 +102,35 @@ Foam::extendedUpwindStencil::weightedSum { fvsPatchField& 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& stField = ownFld[faceI]; - const List& stWeight = ownWeights[faceI]; + label faceI = pSfCorr.patch().patch().start(); - forAll(stField, j) - { - pSfCorr[i] += stField[j]*stWeight[j]; - } - } - else + forAll(pSfCorr, i) { - const List& stField = neiFld[faceI]; - const List& stWeight = neiWeights[faceI]; - - forAll(stField, j) + if (phi[faceI] > 0) { - pSfCorr[i] += stField[j]*stWeight[j]; + // Flux out of owner. Use upwind (= owner side) stencil. + const List& stField = ownFld[faceI]; + const List& stWeight = ownWeights[faceI]; + + forAll(stField, j) + { + pSfCorr[i] += stField[j]*stWeight[j]; + } } + else + { + const List& stField = neiFld[faceI]; + const List& stWeight = neiWeights[faceI]; + + forAll(stField, j) + { + pSfCorr[i] += stField[j]*stWeight[j]; + } + } + faceI++; } - faceI++; } } diff --git a/src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.H b/src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.H index 93e9d9117a..67179e4607 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.H +++ b/src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.H @@ -68,7 +68,6 @@ public: //- Construct from mesh explicit cellFaceCellStencil(const polyMesh& mesh); - }; diff --git a/src/finiteVolume/fvMesh/fvMesh.C b/src/finiteVolume/fvMesh/fvMesh.C index 8021b09f55..3438e54bfe 100644 --- a/src/finiteVolume/fvMesh/fvMesh.C +++ b/src/finiteVolume/fvMesh/fvMesh.C @@ -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::Delete(*this); + CentredFitData::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::Delete(*this); + CentredFitData::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 & p, @@ -552,6 +548,30 @@ void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap) } +// Temporary helper function to call move points on +// MeshObjects +template +void MeshObjectMovePoints(const Foam::fvMesh& mesh) +{ + if + ( + mesh.db().objectRegistry::foundObject + ( + Type::typeName + ) + ) + { + const_cast + ( + mesh.db().objectRegistry::lookupObject + ( + Type::typeName + ) + ).movePoints(); + } +} + + Foam::tmp Foam::fvMesh::movePoints(const pointField& p) { // Grab old time volumes if the time has been incremented @@ -642,133 +662,12 @@ Foam::tmp 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::typeName - ) - ) - { - const_cast - ( - db().objectRegistry::lookupObject - ( - volPointInterpolation::typeName - ) - ).movePoints(); - } - - // extendedLeastSquaresVectors - if - ( - db().objectRegistry::foundObject - ( - extendedLeastSquaresVectors::typeName - ) - ) - { - const_cast - ( - db().objectRegistry::lookupObject - ( - extendedLeastSquaresVectors::typeName - ) - ).movePoints(); - } - - // leastSquaresVectors - if - ( - db().objectRegistry::foundObject - ( - leastSquaresVectors::typeName - ) - ) - { - const_cast - ( - db().objectRegistry::lookupObject - ( - leastSquaresVectors::typeName - ) - ).movePoints(); - } - - //// linearFitData - //if - //( - // db().objectRegistry::foundObject - // ( - // linearFitData::typeName - // ) - //) - //{ - // const_cast - // ( - // db().objectRegistry::lookupObject - // ( - // linearFitData::typeName - // ) - // ).movePoints(); - //} - - // quadraticFitData - if - ( - db().objectRegistry::foundObject - ( - quadraticFitData::typeName - ) - ) - { - const_cast - ( - db().objectRegistry::lookupObject - ( - quadraticFitData::typeName - ) - ).movePoints(); - } - - //// quadraticFitSnGradData - //if - //( - // db().objectRegistry::foundObject - // ( - // quadraticFitSnGradData::typeName - // ) - //) - //{ - // const_cast - // ( - // db().objectRegistry::lookupObject - // ( - // quadraticFitSnGradData::typeName - // ) - // ).movePoints(); - //} - - // skewCorrectionVectors - if - ( - db().objectRegistry::foundObject - ( - skewCorrectionVectors::typeName - ) - ) - { - const_cast - ( - db().objectRegistry::lookupObject - ( - skewCorrectionVectors::typeName - ) - ).movePoints(); - } - + MeshObjectMovePoints(*this); + MeshObjectMovePoints(*this); + MeshObjectMovePoints(*this); + MeshObjectMovePoints >(*this); + MeshObjectMovePoints >(*this); + MeshObjectMovePoints(*this); return tsweptVols; } diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C index dc8720ac89..dea1b41143 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C @@ -43,6 +43,7 @@ Foam::CentredFitData::CentredFitData ) : MeshObject >(mesh), + stencil_(stencil), linearLimitFactor_(linearLimitFactor), centralWeight_(centralWeight), # ifdef SPHERICAL_GEOMETRY @@ -51,62 +52,29 @@ Foam::CentredFitData::CentredFitData dim_(mesh.nGeometricD()), # endif minSize_(Polynomial::nTerms(dim_)), - coeffs_(mesh.nInternalFaces()) + coeffs_(mesh.nFaces()) { if (debug) { Info<< "Contructing CentredFitData" << endl; } - // check input - if (centralWeight_ < 1 - SMALL) + // Check input + if (linearLimitFactor > 1) { FatalErrorIn("CentredFitData::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 > 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::CentredFitData() :" << "Finished constructing polynomialFit data" << endl; - - interpPolySize.write(); } } @@ -120,15 +88,14 @@ void Foam::CentredFitData::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::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::findFaceDirs } +template +void Foam::CentredFitData::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 > 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 Foam::label Foam::CentredFitData::calcFit ( const List& 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::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::calcFit // Set the fit label stencilSize = C.size(); - coeffs_[faci].setSize(stencilSize); + coeffs_[facei].setSize(stencilSize); scalarList singVals(minSize_); label nSVDzeros = 0; - const GeometricField& w = - this->mesh().surfaceInterpolation::weights(); - bool goodFit = false; for(int iIt = 0; iIt < 10 && !goodFit; iIt++) { @@ -251,11 +256,11 @@ Foam::label Foam::CentredFitData::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::calcFit if (goodFit) { - coeffs_[faci][0] = fit0; - coeffs_[faci][1] = fit1; + coeffs_[facei][0] = fit0; + coeffs_[facei][1] = fit1; for(label i=2; i::calcFit // ( // min // ( - // min(alpha - beta*mag(coeffs_[faci][0] - w[faci])/w[faci], 1), - // min(alpha - beta*mag(coeffs_[faci][1] - (1 - w[faci]))/(1 - w[faci]), 1) + // min(alpha - beta*mag(coeffs_[facei][0] - wLin)/wLin, 1), + // min(alpha - beta*mag(coeffs_[facei][1] - (1 - wLin)) + // /(1 - wLin), 1) // ), 0 // ); - //Info<< w[faci] << " " << coeffs_[faci][0] << " " << (1 - w[faci]) << " " << coeffs_[faci][1] << endl; + //Info<< wLin << " " << coeffs_[facei][0] + // << " " << (1 - wLin) << " " << coeffs_[facei][1] << endl; - // Remove the uncorrected linear coefficients - coeffs_[faci][0] -= w[faci]; - coeffs_[faci][1] -= 1 - w[faci]; + // Remove the uncorrected linear ocefficients + coeffs_[facei][0] -= wLin; + coeffs_[facei][1] -= 1 - wLin; // if (limiter < 0.99) // { // for(label i = 0; i < stencilSize; i++) // { - // coeffs_[faci][i] *= limiter; + // coeffs_[facei][i] *= limiter; // } // } } @@ -326,16 +333,16 @@ Foam::label Foam::CentredFitData::calcFit WarningIn ( "CentredFitData::calcFit" - "(const List& C, const label faci" - ) << "Could not fit face " << faci + "(const List& C, const label facei" + ) << "Could not fit face " << facei << ", reverting to linear." << nl << " Weights " - << coeffs_[faci][0] << " " << w[faci] << nl + << coeffs_[facei][0] << " " << wLin << nl << " Linear weights " - << coeffs_[faci][1] << " " << 1 - w[faci] << endl; + << coeffs_[facei][1] << " " << 1 - wLin << endl; } - coeffs_[faci] = 0; + coeffs_[facei] = 0; } return minSize_ - nSVDzeros; @@ -345,8 +352,7 @@ Foam::label Foam::CentredFitData::calcFit template bool Foam::CentredFitData::movePoints() { - notImplemented("CentredFitData::movePoints()"); - + calcFit(); return true; } diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H index c225642964..1630e06a92 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H @@ -57,6 +57,9 @@ class CentredFitData { // Private data + //- The stencil the fit is based on + const extendedCentredStencil& stencil_; + //- Factor the fit is allowed to deviate from linear. // This limits the amount of high-order correction and increases // stability on bad meshes @@ -88,7 +91,17 @@ class CentredFitData const label faci ); - label calcFit(const List&, const label faci); + //- Calculate the fit for the all the mesh faces + // and set the coefficients + void calcFit(); + + //- Calculate the fit for the specified face and set the coefficients + label calcFit + ( + const List&, // Stencil points + const scalar wLin, // Linear weight + const label faci // Current face index + ); public: