diff --git a/src/fvMeshStitchers/moving/fvMeshStitchersMoving.C b/src/fvMeshStitchers/moving/fvMeshStitchersMoving.C index 8b67c4e14a..309f888e72 100644 --- a/src/fvMeshStitchers/moving/fvMeshStitchersMoving.C +++ b/src/fvMeshStitchers/moving/fvMeshStitchersMoving.C @@ -24,8 +24,11 @@ License \*---------------------------------------------------------------------------*/ #include "fvMeshStitchersMoving.H" +#include "FvFaceCellWave.H" #include "fvMeshSubset.H" #include "fvmLaplacian.H" +#include "meshPhiCorrectInfo.H" +#include "meshPhiPreCorrectInfo.H" #include "movingWallVelocityFvPatchVectorField.H" #include "movingWallSlipVelocityFvPatchVectorField.H" #include "regionSplit.H" @@ -537,18 +540,34 @@ void Foam::fvMeshStitchers::moving::internalFaceCorrectMeshPhi // across the couplings, as the mismatch will be fixed later when area and // mesh flux is added to the error faces. So, here, we correct the // correction (!) by propagating the error from the cells furthest from the - // interface back to the interface. This process is wave-like, but the - // sub-mesh is thin, so it should only require 1 or 2 iterations. + // interface back to the interface. This is done using waves. - // Initialise layers and weights - List subFaceLayerAndWeight(subMesh.nFaces(), {-1, NaN}); - List subCellLayerAndWeight(subMesh.nCells(), {-1, Zero}); + // Dynamic memory for wave initialisation + DynamicList subChangedPatchAndFaces; + DynamicList subChangedFacePci; + DynamicList subChangedFaceCi; + + // Allocate data for the pre-correction wave + List subInternalFacePci(subMesh.nInternalFaces()); + List> subPatchFacePci + ( + FvFaceCellWave::template + sizesListList>> + ( + FvFaceCellWave::template + listListSizes(subMesh.boundary()), + meshPhiPreCorrectInfo() + ) + ); + List subCellPci(subMesh.nCells()); + + // Initialisation for the pre-correction wave + subChangedPatchAndFaces.clear(); + subChangedFacePci.clear(); forAll(subMesh.boundary(), subPatchi) { const fvPatch& subFvp = subMesh.boundary()[subPatchi]; - if (subFvp.coupled()) continue; - forAll(subFvp, subPatchFacei) { const label subFacei = subFvp.start() + subPatchFacei; @@ -557,171 +576,127 @@ void Foam::fvMeshStitchers::moving::internalFaceCorrectMeshPhi if (bFacei >= 0 && bFaceIsOwnerOrig[bFacei]) { - subFaceLayerAndWeight[subFacei] = {0, bFaceNccMagSf[bFacei]}; + subChangedPatchAndFaces.append({subPatchi, subPatchFacei}); + subChangedFacePci.append + ( + meshPhiPreCorrectInfo(0, bFaceNccMagSf[bFacei]) + ); } } } - // Create layers and weights by waving from owner-orig faces - while (gMin(subFaceLayerAndWeight).layer < 0) + // Pre-correction wave + FvFaceCellWave preWave + ( + subMesh, + subInternalFacePci, + subPatchFacePci, + subCellPci + ); + preWave.setFaceInfo(subChangedPatchAndFaces, subChangedFacePci); + const label nWaveLayers = + preWave.iterate(subMesh.globalData().nTotalCells() + 1); + + // Allocate data for the correction wave + List subInternalFaceCi(subMesh.nInternalFaces()); + List> subPatchFaceCi + ( + FvFaceCellWave::template + sizesListList>> + ( + FvFaceCellWave::template + listListSizes(subMesh.boundary()), + meshPhiCorrectInfo() + ) + ); + List subCellCi(subMesh.nCells()); + + // Calculate the current error in the rate of change of volume + const volScalarField::Internal dVdtError + ( + (subV - subV0)/subMesh.time().deltaT() + - fvc::surfaceIntegrate(subPhi + subDeltaPhi)()*subV + ); + + // Construct track data for the correction wave + meshPhiCorrectInfo::trackData td + ( + subInternalFacePci, + subPatchFacePci, + subCellPci, + dVdtError + ); + + // Wave backwards through the layers to generate the corrections. Note that + // this has to be done in stages, so that later layers complete in their + // entirety before the earlier layers begin. Otherwise they interfere. + for (label waveLayeri = nWaveLayers - 1; waveLayeri >= 0; waveLayeri --) { - // Face to cell - List subCellLayerAndWeight0(subCellLayerAndWeight); - forAll(subMesh.faces(), subFacei) + // The layer indices on the faces that we want to wave from + const label faceLayeri = (waveLayeri + 1)*2; + + // Initialisation for the correction wave + subChangedPatchAndFaces.clear(); + subChangedFaceCi.clear(); + forAll(subInternalFacePci, subFacei) { - if (subFaceLayerAndWeight[subFacei].layer != -1) + if (subInternalFacePci[subFacei].layer() == faceLayeri) { - const label subOwnCelli = subMesh.faceOwner()[subFacei]; - if (subCellLayerAndWeight0[subOwnCelli].layer == -1) - { - subCellLayerAndWeight[subOwnCelli].layer = - subFaceLayerAndWeight[subFacei].layer + 1; - subCellLayerAndWeight[subOwnCelli].weight += - subFaceLayerAndWeight[subFacei].weight; - } - - if (!subMesh.isInternalFace(subFacei)) continue; - - const label subNbrCelli = subMesh.faceNeighbour()[subFacei]; - if (subCellLayerAndWeight0[subNbrCelli].layer == -1) - { - subCellLayerAndWeight[subNbrCelli].layer = - subFaceLayerAndWeight[subFacei].layer + 1; - subCellLayerAndWeight[subNbrCelli].weight += - subFaceLayerAndWeight[subFacei].weight; - } + subChangedPatchAndFaces.append({-1, subFacei}); + subChangedFaceCi.append + ( + subInternalFaceCi[subFacei].valid(td) + ? subInternalFaceCi[subFacei] + : meshPhiCorrectInfo(meshPhiCorrectInfo::shape::face) + ); } } - - // Cell to face - List subFaceLayerAndWeight0(subFaceLayerAndWeight); - forAll(subMesh.cells(), subCelli) + forAll(subPatchFacePci, subPatchi) { - if (subCellLayerAndWeight[subCelli].layer != -1) + forAll(subPatchFacePci[subPatchi], subPatchFacei) { - forAll(subMesh.cells()[subCelli], subCellFacei) - { - const label subFacei = - subMesh.cells()[subCelli][subCellFacei]; - if (subFaceLayerAndWeight0[subFacei].layer == -1) - { - subFaceLayerAndWeight[subFacei].layer = - subCellLayerAndWeight[subCelli].layer + 1; - subFaceLayerAndWeight[subFacei].weight = - subCellLayerAndWeight[subCelli].weight; - } - } - } - } - - // Synchronise - SubList subBFaceLayerAndWeight - ( - subFaceLayerAndWeight, - subMesh.nFaces() - subMesh.nInternalFaces(), - subMesh.nInternalFaces() - ); - syncTools::syncBoundaryFaceList - ( - subMesh, - subBFaceLayerAndWeight, - maxEqOp(), - [](const coupledPolyPatch&, List&){} - ); - } - - // Function to extract a value from a surface field given a face index - auto subFaceValue = [&subMesh] - ( - surfaceScalarField& sf, - const label subFacei - ) -> scalar& - { - const label subBFacei = - subMesh.isInternalFace(subFacei) - ? -1 : subFacei - subMesh.nInternalFaces(); - const label subPatchi = - subMesh.isInternalFace(subFacei) - ? -1 : subMesh.boundaryMesh().patchID()[subBFacei]; - const label subPatchFacei = - subMesh.isInternalFace(subFacei) - ? -1 : subFacei - subMesh.boundaryMesh()[subPatchi].start(); - - return - subMesh.isInternalFace(subFacei) - ? sf[subFacei] - : sf.boundaryFieldRef()[subPatchi][subPatchFacei]; - }; - - // Correct volume conservation error by waving back to the owner-orig faces - for - ( - label cellLayeri = gMax(subCellLayerAndWeight).layer; - cellLayeri > 0; - cellLayeri -= 2 - ) - { - const surfaceScalarField subDeltaPhi0(subDeltaPhi); - - forAll(subMesh.cells(), subCelli) - { - if (subCellLayerAndWeight[subCelli].layer != cellLayeri) continue; - - // Compute the error for this cell - scalar deltaV = subV[subCelli] - subV0[subCelli]; - forAll(subV.mesh().cells()[subCelli], subCellFacei) - { - const label subFacei = - subV.mesh().cells()[subCelli][subCellFacei]; - const bool o = subV.mesh().faceOwner()[subFacei] == subCelli; - - deltaV -= - ( - subFaceValue(subPhi, subFacei) - + subFaceValue(subDeltaPhi, subFacei) - ) - *(o ? +1 : -1); - } - - // Push the error onto the lower-layered faces - forAll(subMesh.cells()[subCelli], subCellFacei) - { - const label subFacei = - subMesh.cells()[subCelli][subCellFacei]; - const bool o = subMesh.faceOwner()[subFacei] == subCelli; - const scalar w = - subFaceLayerAndWeight[subFacei].weight - /subCellLayerAndWeight[subCelli].weight; - if ( - subCellLayerAndWeight[subCelli].layer - > subFaceLayerAndWeight[subFacei].layer + subPatchFacePci[subPatchi][subPatchFacei].layer() + == faceLayeri ) { - subFaceValue(subDeltaPhi, subFacei) += - deltaV*w*(o ? +1 : -1); + subChangedPatchAndFaces.append({subPatchi, subPatchFacei}); + subChangedFaceCi.append + ( + subPatchFaceCi[subPatchi][subPatchFacei].valid(td) + ? subPatchFaceCi[subPatchi][subPatchFacei] + : meshPhiCorrectInfo(meshPhiCorrectInfo::shape::face) + ); } } } - // Synchronise - const surfaceScalarField::Boundary dSubDeltaPhiNbrBf + // Correction wave + FvFaceCellWave wave ( - surfaceScalarField::Internal::null(), - (subDeltaPhi - subDeltaPhi0)() - .boundaryField() - .boundaryNeighbourField() + subMesh, + subInternalFaceCi, + subPatchFaceCi, + subCellCi, + td ); - forAll(subMesh.boundary(), subPatchi) - { - const fvPatch& subFvp = subMesh.boundary()[subPatchi]; + wave.setFaceInfo(subChangedPatchAndFaces, subChangedFaceCi); + wave.iterate(1); + } - if (subFvp.coupled()) - { - subDeltaPhi.boundaryFieldRef()[subPatchi] -= - dSubDeltaPhiNbrBf[subPatchi]; - } + // Apply corrections + forAll(subInternalFaceCi, subFacei) + { + subDeltaPhi.primitiveFieldRef()[subFacei] += + subInternalFaceCi[subFacei].deltaPhi(); + } + forAll(subPatchFacePci, subPatchi) + { + forAll(subPatchFacePci[subPatchi], subPatchFacei) + { + subDeltaPhi.boundaryFieldRef()[subPatchi][subPatchFacei] += + subPatchFaceCi[subPatchi][subPatchFacei].deltaPhi(); } } diff --git a/src/fvMeshStitchers/moving/meshPhiCorrectInfo.H b/src/fvMeshStitchers/moving/meshPhiCorrectInfo.H new file mode 100644 index 0000000000..61b6bfddca --- /dev/null +++ b/src/fvMeshStitchers/moving/meshPhiCorrectInfo.H @@ -0,0 +1,243 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2022 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 . + +Class + Foam::meshPhiCorrectInfo + +SourceFiles + meshPhiCorrectInfoI.H + +\*---------------------------------------------------------------------------*/ + +#ifndef meshPhiCorrectInfo_H +#define meshPhiCorrectInfo_H + +#include "meshPhiPreCorrectInfo.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class fvPatch; +class fvMesh; +class transformer; + +// Forward declaration of friend functions and operators +class meshPhiCorrectInfo; +Ostream& operator<<(Ostream&, const meshPhiCorrectInfo&); +Istream& operator>>(Istream&, meshPhiCorrectInfo&); + +/*---------------------------------------------------------------------------*\ + Class meshPhiCorrectInfo Declaration +\*---------------------------------------------------------------------------*/ + +class meshPhiCorrectInfo +{ +public: + + // Public Classes + + //- Enumeration to define the mesh shape the info applies to + enum class shape + { + invalid, + face, + cell + }; + + + //- Tracking data. Mostly just a collection of references to the + // pre-correct info. + class trackData + { + public: + + const List& internalFacePci_; + const List>& patchFacePci_; + const List& cellPci_; + const DimensionedField& dVdtError_; + + trackData + ( + const List& internalFacePci, + const List>& patchFacePci, + const List& cellPci, + const DimensionedField& dVdtError + ) + : + internalFacePci_(internalFacePci), + patchFacePci_(patchFacePci), + cellPci_(cellPci), + dVdtError_(dVdtError) + {} + }; + + +private: + + // Private Data + + //- The primitive shape that this info applies to + shape shape_; + + //- Flux correction (if on a face) or volume conservation correction + // (if in a cell) + scalar delta_; + + + // Private Member Functions + + //- Throw an error if this is not a face + void validateFace() const; + + //- Throw an error if this is not a cell + void validateCell() const; + + +public: + + // Constructors + + //- Construct null + inline meshPhiCorrectInfo(); + + //- Construct for given shape + inline meshPhiCorrectInfo(const shape& s); + + + // Member Functions + + // Access + + //- Return the flux correction + inline scalar deltaPhi() const; + + //- Access the flux correction + inline scalar& deltaPhi(); + + //- Return the volume change rate error correction + inline scalar deltaDVdtError() const; + + //- Access the volume change rate error correction + inline scalar& deltaDVdtError(); + + + // Needed by FaceCellWave + + //- Check whether the meshPhiCorrectInfo has been changed at all + // or still contains original (invalid) value. + template + inline bool valid(TrackingData& td) const; + + //- Check for identical geometrical data. Used for checking + // consistency across cyclics. + template + inline bool sameGeometry + ( + const fvMesh& mesh, + const meshPhiCorrectInfo&, + const scalar, + TrackingData& td + ) const; + + //- Transform across an interface + template + inline void transform + ( + const fvPatch& patch, + const label patchFacei, + const transformer& transform, + TrackingData& td + ); + + //- Influence of neighbouring face + template + inline bool updateCell + ( + const fvMesh& mesh, + const label thisCelli, + const labelPair& neighbourPatchAndFacei, + const meshPhiCorrectInfo& neighbourInfo, + const scalar tol, + TrackingData& td + ); + + //- Influence of neighbouring cell + template + inline bool updateFace + ( + const fvMesh& mesh, + const labelPair& thisPatchAndFacei, + const label neighbourCelli, + const meshPhiCorrectInfo& neighbourInfo, + const scalar tol, + TrackingData& td + ); + + //- Influence of different value on same face + template + inline bool updateFace + ( + const fvMesh& mesh, + const labelPair& thisPatchAndFacei, + const meshPhiCorrectInfo& neighbourInfo, + const scalar tol, + TrackingData& td + ); + + //- Test equality + template + inline bool equal + ( + const meshPhiCorrectInfo&, + TrackingData& td + ) const; + + + // Member Operators + + inline bool operator==(const meshPhiCorrectInfo&) const; + inline bool operator!=(const meshPhiCorrectInfo&) const; + + + // IOstream Operators + + friend Ostream& operator<<(Ostream&, const meshPhiCorrectInfo&); + friend Istream& operator>>(Istream&, meshPhiCorrectInfo&); +}; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "meshPhiCorrectInfoI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/fvMeshStitchers/moving/meshPhiCorrectInfoI.H b/src/fvMeshStitchers/moving/meshPhiCorrectInfoI.H new file mode 100644 index 0000000000..9d31bcc0ee --- /dev/null +++ b/src/fvMeshStitchers/moving/meshPhiCorrectInfoI.H @@ -0,0 +1,292 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2022 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 . + +\*---------------------------------------------------------------------------*/ + +#include "meshPhiCorrectInfo.H" +#include "fvMesh.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +inline void Foam::meshPhiCorrectInfo::validateFace() const +{ + if (shape_ != shape::face) + { + FatalError + << "Face data requested from a non-face correction info" + << exit(FatalError); + } +} + + +inline void Foam::meshPhiCorrectInfo::validateCell() const +{ + if (shape_ != shape::cell) + { + FatalError + << "Cell data requested from a non-cell correction info" + << exit(FatalError); + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +inline Foam::meshPhiCorrectInfo::meshPhiCorrectInfo() +: + shape_(shape::invalid), + delta_(NaN) +{} + + +inline Foam::meshPhiCorrectInfo::meshPhiCorrectInfo(const shape& s) +: + shape_(s), + delta_(s == shape::invalid ? NaN : 0) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline Foam::scalar Foam::meshPhiCorrectInfo::deltaPhi() const +{ + validateFace(); + return delta_; +} + + +inline Foam::scalar& Foam::meshPhiCorrectInfo::deltaPhi() +{ + validateFace(); + return delta_; +} + + +inline Foam::scalar Foam::meshPhiCorrectInfo::deltaDVdtError() const +{ + validateCell(); + return delta_; +} + + +inline Foam::scalar& Foam::meshPhiCorrectInfo::deltaDVdtError() +{ + validateCell(); + return delta_; +} + + +template +inline bool Foam::meshPhiCorrectInfo::valid(TrackingData& td) const +{ + return shape_ != shape::invalid; +} + + +template +inline bool Foam::meshPhiCorrectInfo::sameGeometry +( + const fvMesh&, + const meshPhiCorrectInfo& l, + const scalar tol, + TrackingData& td +) const +{ + return true; +} + + +template +inline void Foam::meshPhiCorrectInfo::transform +( + const fvPatch& patch, + const label patchFacei, + const transformer& transform, + TrackingData& td +) +{} + + +template +inline bool Foam::meshPhiCorrectInfo::updateCell +( + const fvMesh& mesh, + const label thisCelli, + const labelPair& neighbourPatchAndFacei, + const meshPhiCorrectInfo& neighbourInfo, + const scalar tol, + TrackingData& td +) +{ + const label neighbourPatchi = neighbourPatchAndFacei.first(); + const label neighbourFacei = neighbourPatchAndFacei.second(); + + const meshPhiPreCorrectInfo& thisPci = + td.cellPci_[thisCelli]; + const meshPhiPreCorrectInfo& neighbourPci = + neighbourPatchi == -1 + ? td.internalFacePci_[neighbourFacei] + : td.patchFacePci_[neighbourPatchi][neighbourFacei]; + + if (!valid(td)) + { + shape_ = shape::cell; + deltaDVdtError() = 0; + } + + if (thisPci.layer() < neighbourPci.layer()) + { + const label s = + neighbourPatchi == -1 + ? (mesh.owner()[neighbourFacei] == thisCelli ? +1 : -1) + : +1; + + deltaDVdtError() -= s*neighbourInfo.deltaPhi(); + + return true; + } + else + { + return false; + } +} + + +template +inline bool Foam::meshPhiCorrectInfo::updateFace +( + const fvMesh& mesh, + const labelPair& thisPatchAndFacei, + const label neighbourCelli, + const meshPhiCorrectInfo& neighbourInfo, + const scalar tol, + TrackingData& td +) +{ + const label thisPatchi = thisPatchAndFacei.first(); + const label thisFacei = thisPatchAndFacei.second(); + + const meshPhiPreCorrectInfo& thisPci = + thisPatchi == -1 + ? td.internalFacePci_[thisFacei] + : td.patchFacePci_[thisPatchi][thisFacei]; + const meshPhiPreCorrectInfo& neighbourPci = + td.cellPci_[neighbourCelli]; + + if (!valid(td)) + { + shape_ = shape::face; + deltaPhi() = 0; + } + + if (thisPci.layer() < neighbourPci.layer()) + { + const label s = + thisPatchi == -1 + ? (mesh.owner()[thisFacei] == neighbourCelli ? +1 : -1) + : +1; + + deltaPhi() = + s*thisPci.weight()/neighbourPci.weight() + *(td.dVdtError_[neighbourCelli] + neighbourInfo.deltaDVdtError()); + + return true; + } + else + { + return false; + } +} + + +template +inline bool Foam::meshPhiCorrectInfo::updateFace +( + const fvMesh& mesh, + const labelPair& thisPatchAndFacei, + const meshPhiCorrectInfo& neighbourInfo, + const scalar tol, + TrackingData& td +) +{ + if (!valid(td)) + { + shape_ = shape::face; + deltaPhi() = 0; + } + + deltaPhi() = - neighbourInfo.deltaPhi(); + + return true; +} + + +template +inline bool Foam::meshPhiCorrectInfo::equal +( + const meshPhiCorrectInfo& rhs, + TrackingData& td +) const +{ + return operator==(rhs); +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +inline bool Foam::meshPhiCorrectInfo::operator== +( + const Foam::meshPhiCorrectInfo& rhs +) const +{ + return shape_ == rhs.shape_ && delta_ == rhs.delta_; +} + + +inline bool Foam::meshPhiCorrectInfo::operator!= +( + const Foam::meshPhiCorrectInfo& rhs +) const +{ + return !(*this == rhs); +} + + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + +Foam::Ostream& Foam::operator<<(Ostream& os, const meshPhiCorrectInfo& l) +{ + return os << static_cast