From 132cbd97de069126a008efbd22c18c081f3479f8 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Thu, 7 Apr 2016 09:03:05 +0200 Subject: [PATCH] Updates in fvBlockMatrix --- .../fvMatrices/fvBlockMatrix/fvBlockMatrix.C | 303 ++++++++++++++++-- .../fvMatrices/fvBlockMatrix/fvBlockMatrix.H | 27 +- 2 files changed, 305 insertions(+), 25 deletions(-) diff --git a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C index 294935d87..9051b50b7 100644 --- a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C +++ b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C @@ -26,16 +26,11 @@ License #include "fvBlockMatrix.H" #include "IOstreams.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template template -void fvBlockMatrix::insertSolutionVector +void Foam::fvBlockMatrix::insertSolutionVector ( const direction dir, const Field& xSingle @@ -64,7 +59,7 @@ void fvBlockMatrix::insertSolutionVector template template -void fvBlockMatrix::insertDiagSource +void Foam::fvBlockMatrix::insertDiagSource ( const direction dir, fvMatrix& matrix @@ -167,7 +162,7 @@ void fvBlockMatrix::insertDiagSource template template -void fvBlockMatrix::insertUpperLower +void Foam::fvBlockMatrix::insertUpperLower ( const direction dir, const fvMatrix& matrix @@ -295,7 +290,7 @@ void fvBlockMatrix::insertUpperLower template template -void fvBlockMatrix::updateCouplingCoeffs +void Foam::fvBlockMatrix::updateCouplingCoeffs ( const direction dir, const fvMatrix& matrix @@ -568,7 +563,7 @@ void Foam::fvBlockMatrix::insertBoundaryContributions template -void fvBlockMatrix::insertCouplingDiag +void Foam::fvBlockMatrix::insertCouplingDiag ( const direction dirI, const direction dirJ, @@ -591,7 +586,39 @@ void fvBlockMatrix::insertCouplingDiag template -void fvBlockMatrix::insertCouplingUpperLower +void Foam::fvBlockMatrix::insertCouplingDiagSource +( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& matrix +) +{ + // Prepare the diagonal and source + scalarField diag = matrix.diag(); + scalarField source = matrix.source(); + + // Add boundary source contribution + matrix.addBoundaryDiag(diag, 0); + matrix.addBoundarySource(source, false); + + // Get reference to block diagonal of the block system + typename CoeffField::squareTypeField& blockDiag = + this->diag().asSquare(); + + // Get reference to this source field of the block system + Field& b = this->source(); + + // Set off-diagonal coefficient + forAll(diag, cellI) + { + blockDiag[cellI](dirI, dirJ) += diag[cellI]; + b[cellI](dirI) += source[cellI]; + } +} + + +template +void Foam::fvBlockMatrix::insertCouplingUpperLower ( const direction dirI, const direction dirJ, @@ -678,7 +705,7 @@ void fvBlockMatrix::insertCouplingUpperLower // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -fvBlockMatrix::fvBlockMatrix +Foam::fvBlockMatrix::fvBlockMatrix ( GeometricField& psi ) @@ -691,7 +718,7 @@ fvBlockMatrix::fvBlockMatrix template -fvBlockMatrix::fvBlockMatrix +Foam::fvBlockMatrix::fvBlockMatrix ( const fvBlockMatrix& bxs ) @@ -705,7 +732,7 @@ fvBlockMatrix::fvBlockMatrix template template -void fvBlockMatrix::retrieveSolution +void Foam::fvBlockMatrix::retrieveSolution ( const direction dir, Field& xSingle @@ -734,7 +761,7 @@ void fvBlockMatrix::retrieveSolution template template -void fvBlockMatrix::insertEquation +void Foam::fvBlockMatrix::insertEquation ( const direction dir, fvMatrix& matrix @@ -749,7 +776,7 @@ void fvBlockMatrix::insertEquation template template -void fvBlockMatrix::insertBlockCoupling +void Foam::fvBlockMatrix::insertBlockCoupling ( const direction dirI, const direction dirJ, @@ -778,7 +805,20 @@ void Foam::fvBlockMatrix::insertEquationCoupling template -void fvBlockMatrix::blockAdd +void Foam::fvBlockMatrix::insertEquationCoupling +( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& matrix +) +{ + insertCouplingDiagSource(dirI, dirJ, matrix); + insertCouplingUpperLower(dirI, dirJ, matrix); +} + + +template +void Foam::fvBlockMatrix::blockAdd ( const direction dir, const scalarField& xSingle, @@ -827,7 +867,226 @@ void Foam::fvBlockMatrix::updateSourceCoupling() template -BlockSolverPerformance fvBlockMatrix::solve +void Foam::fvBlockMatrix::insertPicardTensor +( + const direction UEqnDir, + const volVectorField& U, + const surfaceScalarField& phi +) +{ + const direction nCmpts = pTraits::nComponents; + direction localDirI = UEqnDir; + direction localDirJ = UEqnDir; + + // Sanity check + { + const direction blockMatrixSize = pTraits::nComponents; + + if (nCmpts > blockMatrixSize) + { + FatalErrorIn + ( + "void fvBlockMatrix::insertPicardTensor\n" + "(\n" + " const direction UEqnDir,\n" + " const volVectorField& U,\n" + " const surfaceScalarField& phi\n" + ")" + ) << "Trying to insert Picard tensor term into smaller " + << "fvBlockMatrix. Do you have momentum equation?" + << abort(FatalError); + } + } + + // Get weights for U which makes the implicit flux part + const fvMesh& mesh = U.mesh(); + + tmp > tinterpScheme = + fvc::scheme(mesh, U.name()); + + tmp tweights = tinterpScheme().weights(U); + const scalarField& wIn = tweights().internalField(); + + // Calculate the Pi tensor. Consider hard coding the interpolation scheme to + // correspond to the div(phi, U) interpolation scheme for consistency. + // VV, 21/July/2014. + const surfaceTensorField pi + ( + fvc::interpolate(U, -phi, "Uconvection")*mesh.Sf() + ); + const tensorField& piIn = pi.internalField(); + + BlockLduSystem piSystem(mesh); + + // Get references to ldu fields of pi block system always as square + typename CoeffField::squareTypeField& piDiag = + piSystem.diag().asSquare(); + typename CoeffField::squareTypeField& piUpper = + piSystem.upper().asSquare(); + typename CoeffField::squareTypeField& piLower = + piSystem.lower().asSquare(); + + vectorField& piSource = piSystem.source(); + + piLower = -wIn*piIn; + piUpper = piLower + piIn; + piSystem.negSumDiag(); + + // Boundary contributions - hard coded or explicit because of the problem + // with inconsistent return type of internalCoeffs. VV, 21/July/2014. + forAll(U.boundaryField(), patchI) + { + const fvPatchVectorField& Ub = U.boundaryField()[patchI]; + const fvPatch& patch = Ub.patch(); + const fvsPatchTensorField& pib = pi.boundaryField()[patchI]; + const fvsPatchScalarField& wb = tweights().boundaryField()[patchI]; + const unallocLabelList& fc = patch.faceCells(); + + // Check for empty patches. Needed since the boundary conditions are + // hard coded. VV. 18/Sep/2014. + if (patch.type() == "empty") + { + continue; + } + + // Hard coded zeroGradient if patch does not fix value + if (!Ub.fixesValue()) + { + forAll(Ub, faceI) + { + piDiag[fc[faceI]] += pib[faceI]; + } + } + // Coupled patches + else if (patch.coupled()) + { + typename CoeffField::squareTypeField& pipCoupleUpper = + piSystem.coupleUpper()[patchI].asSquare(); + typename CoeffField::squareTypeField& pipCoupleLower = + piSystem.coupleLower()[patchI].asSquare(); + + const tensorField pcl = -wb*pib; + const tensorField pcu = pcl + pib; + + // Coupling contributions + pipCoupleLower -= pcl; + pipCoupleUpper -= pcu; + } + else if (Ub.fixesValue()) + { + const vectorField boundaryCoeffs(Ub.valueBoundaryCoeffs(wb)); + + // Boundary contribution + forAll(Ub, faceI) + { + piSource[fc[faceI]] -= pib[faceI] & boundaryCoeffs[faceI]; + } + } + else + { + FatalErrorIn + ( + "void fvBlockMatrix::insertPicardTensor\n" + "(\n" + " const direction UEqnDir,\n" + " const volVectorField& U,\n" + " const surfaceScalarField& phi\n" + ")" + ) << "Patch does not fix value, nor doesn't fix value nor is" + << " coupled." + << abort(FatalError); + } + } + + // Consider chucking the above code into fvm::picardTensor operator. + // VV, 21/July/2014. + + // Finally add Picard piSystem into this fvBlockMatrix + typename CoeffField::squareTypeField& blockDiag = + this->diag().asSquare(); + typename CoeffField::squareTypeField& blockUpper = + this->upper().asSquare(); + typename CoeffField::squareTypeField& blockLower = + this->lower().asSquare(); + + Field& blockSource = this->source(); + + // Add diagonal, source, lower and upper + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + for (direction cmptJ = 0; cmptJ < nCmpts; cmptJ++) + { + forAll(piDiag, cellI) + { + blockDiag[cellI](localDirI, localDirJ) += + piDiag[cellI](cmptI, cmptJ); + } + + forAll(piUpper, faceI) + { + blockUpper[faceI](localDirI, localDirJ) += + piUpper[faceI](cmptI, cmptJ); + blockLower[faceI](localDirI, localDirJ) += + piLower[faceI](cmptI, cmptJ); + } + + localDirJ++; + } + + forAll(piSource, cellI) + { + blockSource[cellI](localDirI) += piSource[cellI](cmptI); + } + + localDirI++; + + // Reset localDirJ + localDirJ = UEqnDir; + } + + // Reset local direction for coupling contributions + localDirI = UEqnDir; + localDirJ = UEqnDir; + + // Add coupling contributions + forAll(U.boundaryField(), patchI) + { + if (U.boundaryField()[patchI].patch().coupled()) + { + typename CoeffField::squareTypeField& pcoupleUpper = + this->coupleUpper()[patchI].asSquare(); + typename CoeffField::squareTypeField& pcoupleLower = + this->coupleLower()[patchI].asSquare(); + + const typename CoeffField::squareTypeField& pipcu = + piSystem.coupleUpper()[patchI].asSquare(); + const typename CoeffField::squareTypeField& pipcl = + piSystem.coupleLower()[patchI].asSquare(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + for (direction cmptJ = 0; cmptJ < nCmpts; cmptJ++) + { + forAll(pipcu, faceI) + { + pcoupleUpper[faceI](localDirI, localDirJ) += + pipcu[faceI](cmptI, cmptJ); + pcoupleLower[faceI](localDirI, localDirJ) += + pipcl[faceI](cmptI, cmptJ); + } + } + } + + // Reset local directions for other patches + localDirI = UEqnDir; + localDirJ = UEqnDir; + } + } +} + + +template +Foam::BlockSolverPerformance Foam::fvBlockMatrix::solve ( const dictionary& solverControls ) @@ -849,7 +1108,7 @@ BlockSolverPerformance fvBlockMatrix::solve template -BlockSolverPerformance fvBlockMatrix::solve() +Foam::BlockSolverPerformance Foam::fvBlockMatrix::solve() { return solve(psi_.mesh().solutionDict().solverDict(psi_.name())); } @@ -858,7 +1117,7 @@ BlockSolverPerformance fvBlockMatrix::solve() // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // template -Ostream& operator<< +Foam::Ostream& Foam::operator<< ( Ostream& os, const fvBlockMatrix& bxs @@ -871,8 +1130,4 @@ Ostream& operator<< } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - // ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H index 4b56c8d33..c567e7eae 100644 --- a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H +++ b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H @@ -122,6 +122,14 @@ class fvBlockMatrix const scalarField& coeffIJ ); + //- Insert coupling matrix diag and source into this fvBlockMatrix + void insertCouplingDiagSource + ( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& matrix + ); + //- Insert coupling matrix lower and upper into this fvBlockMatrix void insertCouplingUpperLower ( @@ -219,7 +227,7 @@ public: const bool incrementColumnDir ); - //- Insert scalar equation coupling into this fvBlockMatrix + //- Insert diagonal only equation coupling into this fvBlockMatrix // Source compensation is done in function updateSourceCoupling() // after all coupling terms are added. HJ, 27/Apr/2015 void insertEquationCoupling @@ -229,6 +237,14 @@ public: const scalarField& coeffIJ ); + //- Insert scalar equation coupling into this fvBlockMatrix + void insertEquationCoupling + ( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& matrix + ); + //- Add field into block field void blockAdd ( @@ -243,6 +259,15 @@ public: // linearisation void updateSourceCoupling(); + //- Insert Picard tensor term that comes from Picard linearisation + // of convection term in momentum equation. VV, 21/July/2014. + void insertPicardTensor + ( + const direction UEqnDir, + const volVectorField& U, + const surfaceScalarField& phi + ); + // Solver calls for fvBlockMatrix