From 142d9333c76334023383e871d433363f501f2b0a Mon Sep 17 00:00:00 2001 From: Henry <Henry> Date: Tue, 9 Apr 2013 14:28:54 +0100 Subject: [PATCH] CentredFitSnGrad: New fit-based snGrad --- src/finiteVolume/Make/files | 1 + .../CentredFitSnGrad/CentredFitSnGradData.C | 272 ++++++++++++++++++ .../CentredFitSnGrad/CentredFitSnGradData.H | 126 ++++++++ .../CentredFitSnGrad/CentredFitSnGradScheme.H | 182 ++++++++++++ .../quadraticFitSnGrad/quadraticFitSnGrads.C | 51 ++++ .../schemes/FitData/FitData.H | 32 ++- 6 files changed, 658 insertions(+), 6 deletions(-) create mode 100644 src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C create mode 100644 src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.H create mode 100644 src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradScheme.H create mode 100644 src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 1ce1f78949b..2a650f784f1 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -356,6 +356,7 @@ $(snGradSchemes)/faceCorrectedSnGrad/faceCorrectedSnGrads.C $(snGradSchemes)/limitedSnGrad/limitedSnGrads.C $(snGradSchemes)/uncorrectedSnGrad/uncorrectedSnGrads.C $(snGradSchemes)/orthogonalSnGrad/orthogonalSnGrads.C +$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGrads.C convectionSchemes = finiteVolume/convectionSchemes $(convectionSchemes)/convectionScheme/convectionSchemes.C diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C new file mode 100644 index 00000000000..135deac4bf9 --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C @@ -0,0 +1,272 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +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 "CentredFitSnGradData.H" +#include "surfaceFields.H" +#include "volFields.H" +#include "SVD.H" +#include "syncTools.H" +#include "extendedCentredCellToFaceStencil.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template<class Polynomial> +Foam::CentredFitSnGradData<Polynomial>::CentredFitSnGradData +( + const fvMesh& mesh, + const extendedCentredCellToFaceStencil& stencil, + const scalar linearLimitFactor, + const scalar centralWeight +) +: + FitData + < + CentredFitSnGradData<Polynomial>, + extendedCentredCellToFaceStencil, + Polynomial + > + ( + mesh, stencil, true, linearLimitFactor, centralWeight + ), + coeffs_(mesh.nFaces()) +{ + if (debug) + { + Info<< "Contructing CentredFitSnGradData<Polynomial>" << endl; + } + + calcFit(); + + if (debug) + { + Info<< "CentredFitSnGradData<Polynomial>::CentredFitSnGradData() :" + << "Finished constructing polynomialFit data" + << endl; + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class Polynomial> +void Foam::CentredFitSnGradData<Polynomial>::calcFit +( + scalarList& coeffsi, + const List<point>& C, + const scalar wLin, + const scalar deltaCoeff, + const label facei +) +{ + vector idir(1,0,0); + vector jdir(0,1,0); + vector kdir(0,0,1); + this->findFaceDirs(idir, jdir, kdir, facei); + + // Setup the point weights + scalarList wts(C.size(), scalar(1)); + wts[0] = this->centralWeight(); + wts[1] = this->centralWeight(); + + // Reference point + point p0 = this->mesh().faceCentres()[facei]; + + // p0 -> p vector in the face-local coordinate system + vector d; + + // Local coordinate scaling + scalar scale = 1; + + // Matrix of the polynomial components + scalarRectangularMatrix B(C.size(), this->minSize(), scalar(0)); + + forAll(C, ip) + { + const point& p = C[ip]; + const vector p0p = p - p0; + + d.x() = p0p & idir; + d.y() = p0p & jdir; + d.z() = p0p & kdir; + + if (ip == 0) + { + scale = cmptMax(cmptMag((d))); + } + + // Scale the radius vector + d /= scale; + + Polynomial::addCoeffs(B[ip], d, wts[ip], this->dim()); + } + + // Additional weighting for constant and linear terms + for (label i = 0; i < B.n(); i++) + { + B[i][0] *= wts[0]; + B[i][1] *= wts[0]; + } + + // Set the fit + label stencilSize = C.size(); + coeffsi.setSize(stencilSize); + + bool goodFit = false; + for (int iIt = 0; iIt < 8 && !goodFit; iIt++) + { + SVD svd(B, SMALL); + + for (label i=0; i<stencilSize; i++) + { + coeffsi[i] = wts[1]*wts[i]*svd.VSinvUt()[1][i]/scale; + } + + goodFit = + ( + mag(wts[0]*wts[0]*svd.VSinvUt()[0][0] - wLin) + < this->linearLimitFactor()*wLin) + && (mag(wts[0]*wts[1]*svd.VSinvUt()[0][1] - (1 - wLin) + ) < this->linearLimitFactor()*(1 - wLin)) + && coeffsi[0] < 0 && coeffsi[1] > 0 + && mag(coeffsi[0] + deltaCoeff) < 0.5*deltaCoeff + && mag(coeffsi[1] - deltaCoeff) < 0.5*deltaCoeff; + + if (!goodFit) + { + // (not good fit so increase weight in the centre and weight + // for constant and linear terms) + + WarningIn + ( + "CentredFitSnGradData<Polynomial>::calcFit" + "(const List<point>& C, const label facei" + ) << "Cannot fit face " << facei << " iteration " << iIt + << " with sum of weights " << sum(coeffsi) << nl + << " Weights " << coeffsi << nl + << " Linear weights " << wLin << " " << 1 - wLin << nl + << " deltaCoeff " << deltaCoeff << nl + << " sing vals " << svd.S() << nl + << "Components of goodFit:\n" + << " wts[0]*wts[0]*svd.VSinvUt()[0][0] = " + << wts[0]*wts[0]*svd.VSinvUt()[0][0] << nl + << " wts[0]*wts[1]*svd.VSinvUt()[0][1] = " + << wts[0]*wts[1]*svd.VSinvUt()[0][1] + << " dim = " << this->dim() << endl; + + wts[0] *= 10; + wts[1] *= 10; + + for (label j = 0; j < B.m(); j++) + { + B[0][j] *= 10; + B[1][j] *= 10; + } + + for (label i = 0; i < B.n(); i++) + { + B[i][0] *= 10; + B[i][1] *= 10; + } + } + } + + if (goodFit) + { + // Remove the uncorrected coefficients + coeffsi[0] += deltaCoeff; + coeffsi[1] -= deltaCoeff; + } + else + { + WarningIn + ( + "CentredFitSnGradData<Polynomial>::calcFit(..)" + ) << "Could not fit face " << facei + << " Coefficients = " << coeffsi + << ", reverting to uncorrected." << endl; + + coeffsi = 0; + } +} + + +template<class Polynomial> +void Foam::CentredFitSnGradData<Polynomial>::calcFit() +{ + const fvMesh& mesh = this->mesh(); + + // Get the cell/face centres in stencil order. + // Centred face stencils no good for triangles or tets. + // Need bigger stencils + List<List<point> > stencilPoints(mesh.nFaces()); + this->stencil().collectData(mesh.C(), stencilPoints); + + // find the fit coefficients for every face in the mesh + + const surfaceScalarField& w = mesh.surfaceInterpolation::weights(); + const surfaceScalarField& dC = mesh.deltaCoeffs(); + + for (label facei = 0; facei < mesh.nInternalFaces(); facei++) + { + calcFit + ( + coeffs_[facei], + stencilPoints[facei], + w[facei], + dC[facei], + facei + ); + } + + const surfaceScalarField::GeometricBoundaryField& bw = w.boundaryField(); + const surfaceScalarField::GeometricBoundaryField& bdC = dC.boundaryField(); + + forAll(bw, patchi) + { + const fvsPatchScalarField& pw = bw[patchi]; + const fvsPatchScalarField& pdC = bdC[patchi]; + + if (pw.coupled()) + { + label facei = pw.patch().start(); + + forAll(pw, i) + { + calcFit + ( + coeffs_[facei], + stencilPoints[facei], + pw[i], + pdC[i], + facei + ); + facei++; + } + } + } +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.H b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.H new file mode 100644 index 00000000000..c025a17dd5f --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.H @@ -0,0 +1,126 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +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::CentredFitSnGradData + +Description + Data for the quadratic fit correction interpolation scheme + +SourceFiles + CentredFitSnGradData.C + +\*---------------------------------------------------------------------------*/ + +#ifndef CentredFitSnGradData_H +#define CentredFitSnGradData_H + +#include "FitData.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class extendedCentredCellToFaceStencil; + +/*---------------------------------------------------------------------------*\ + Class CentredFitSnGradData Declaration +\*---------------------------------------------------------------------------*/ + +template<class Polynomial> +class CentredFitSnGradData +: + public FitData + < + CentredFitSnGradData<Polynomial>, + extendedCentredCellToFaceStencil, + Polynomial + > +{ + // Private data + + //- For each cell in the mesh store the values which multiply the + // values of the stencil to obtain the gradient for each direction + List<scalarList> coeffs_; + +public: + + TypeName("CentredFitSnGradData"); + + + // Constructors + + //- Construct from components + CentredFitSnGradData + ( + const fvMesh& mesh, + const extendedCentredCellToFaceStencil& stencil, + const scalar linearLimitFactor, + const scalar centralWeight + ); + + + //- Destructor + virtual ~CentredFitSnGradData() + {} + + + // Member functions + + //- Return reference to fit coefficients + const List<scalarList>& coeffs() const + { + return coeffs_; + } + + //- Calculate the fit for the specified face and set the coefficients + void calcFit + ( + scalarList& coeffsi, // coefficients to be set + const List<point>&, // Stencil points + const scalar wLin, // Weight for linear approximation (weights + // nearest neighbours) + const scalar deltaCoeff, // uncorrected delta coefficient + const label faci // Current face index + ); + + void calcFit(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "CentredFitSnGradData.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradScheme.H b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradScheme.H new file mode 100644 index 00000000000..ca307811e6f --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradScheme.H @@ -0,0 +1,182 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +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::CentredFitSnGradScheme + +Description + Centred fit snGrad scheme which applies an explicit correction to snGrad + +\*---------------------------------------------------------------------------*/ + +#ifndef CentredFitSnGradScheme_H +#define CentredFitSnGradScheme_H + +#include "CentredFitSnGradData.H" +#include "snGradScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class fvMesh; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ +/*---------------------------------------------------------------------------*\ + Class CentredFitSnGradSnGradScheme Declaration +\*---------------------------------------------------------------------------*/ + +template<class Type, class Polynomial, class Stencil> +class CentredFitSnGradScheme +: + public snGradScheme<Type> +{ + // Private Data + + //- Factor the fit is allowed to deviate from linear. + // This limits the amount of high-order correction and increases + // stability on bad meshes + const scalar linearLimitFactor_; + + //- Weights for central stencil + const scalar centralWeight_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + CentredFitSnGradScheme(const CentredFitSnGradScheme&); + + //- Disallow default bitwise assignment + void operator=(const CentredFitSnGradScheme&); + + +public: + + //- Runtime type information + TypeName("CentredFitSnGradScheme"); + + + // Constructors + + //- Construct from mesh and Istream + CentredFitSnGradScheme(const fvMesh& mesh, Istream& is) + : + snGradScheme<Type>(mesh), + linearLimitFactor_(readScalar(is)), + centralWeight_(1000) + {} + + + // Member Functions + + //- Return the interpolation weighting factors for the given field + virtual tmp<surfaceScalarField> deltaCoeffs + ( + const GeometricField<Type, fvPatchField, volMesh>& + ) const + { + return this->mesh().deltaCoeffs(); + } + + //- Return true if this scheme uses an explicit correction + virtual bool corrected() const + { + return true; + } + + //- Return the explicit correction to the face-interpolate + virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > + correction + ( + const GeometricField<Type, fvPatchField, volMesh>& vf + ) const + { + const fvMesh& mesh = this->mesh(); + + const extendedCentredCellToFaceStencil& stencil = Stencil::New + ( + mesh + ); + + const CentredFitSnGradData<Polynomial>& cfd = + CentredFitSnGradData<Polynomial>::New + ( + mesh, + stencil, + linearLimitFactor_, + centralWeight_ + ); + + const List<scalarList>& f = cfd.coeffs(); + + tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > sft + = stencil.weightedSum(vf, f); + + sft().dimensions() /= dimLength; + + return sft; + } +}; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// Add the patch constructor functions to the hash tables + +#define makeCentredFitSnGradTypeScheme(SS, POLYNOMIAL, STENCIL, TYPE) \ + \ +typedef CentredFitSnGradScheme<TYPE, POLYNOMIAL, STENCIL> \ + CentredFitSnGradScheme##TYPE##POLYNOMIAL##STENCIL##_; \ +defineTemplateTypeNameAndDebugWithName \ + (CentredFitSnGradScheme##TYPE##POLYNOMIAL##STENCIL##_, #SS, 0); \ + \ +snGradScheme<TYPE>::addMeshConstructorToTable \ +<CentredFitSnGradScheme<TYPE, POLYNOMIAL, STENCIL> > \ + add##SS##STENCIL##TYPE##MeshConstructorToTable_; + +#define makeCentredFitSnGradScheme(SS, POLYNOMIAL, STENCIL) \ + \ +makeCentredFitSnGradTypeScheme(SS,POLYNOMIAL,STENCIL,scalar) \ +makeCentredFitSnGradTypeScheme(SS,POLYNOMIAL,STENCIL,vector) \ +makeCentredFitSnGradTypeScheme(SS,POLYNOMIAL,STENCIL,sphericalTensor) \ +makeCentredFitSnGradTypeScheme(SS,POLYNOMIAL,STENCIL,symmTensor) \ +makeCentredFitSnGradTypeScheme(SS,POLYNOMIAL,STENCIL,tensor) + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C new file mode 100644 index 00000000000..63e4ea98fac --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C @@ -0,0 +1,51 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +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 "CentredFitSnGradScheme.H" +#include "quadraticFitPolynomial.H" +#include "centredCFCCellToFaceStencilObject.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + defineTemplateTypeNameAndDebug + ( + CentredFitSnGradData<quadraticFitPolynomial>, + 0 + ); + + namespace fv + { + makeCentredFitSnGradScheme + ( + quadraticFit, + quadraticFitPolynomial, + centredCFCCellToFaceStencilObject + ); + } +} + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H index cee373b6c4d..5faeed6eb19 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H @@ -29,7 +29,7 @@ Description The linearCorrection_ determines whether the fit is for a corrected linear scheme (first two coefficients are corrections for owner and neighbour) or a pure upwind scheme (first coefficient is correction for - owner ; weight on face taken as 1). + owner; weight on face taken as 1). SourceFiles FitData.C @@ -80,7 +80,7 @@ class FitData const label minSize_; - // Private Member Functions +protected: //- Find the normal direction (i) and j and k directions for face faci void findFaceDirs @@ -93,9 +93,6 @@ class FitData public: - //TypeName("FitData"); - - // Constructors //- Construct from components @@ -122,6 +119,30 @@ public: return stencil_; } + //- Factor the fit is allowed to deviate from the base scheme + scalar linearLimitFactor() const + { + return linearLimitFactor_; + } + + //- Return weight for central stencil + scalar centralWeight() const + { + return centralWeight_; + } + + //- Dimensionality of the geometry + label dim() const + { + return dim_; + } + + //- Minimum stencil size + label minSize() const + { + return minSize_; + } + bool linearCorrection() const { return linearCorrection_; @@ -140,7 +161,6 @@ public: //- Calculate the fit for all the faces virtual void calcFit() = 0; - //- Recalculate weights (but not stencil) when the mesh moves bool movePoints(); }; -- GitLab