diff --git a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C index f2697f212..18dbee6ff 100644 --- a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C +++ b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C @@ -30,7 +30,7 @@ License template template -void fvBlockMatrix::insertSolutionVector +void Foam::fvBlockMatrix::insertSolutionVector ( const direction dir, const Field& xSingle @@ -59,7 +59,7 @@ void fvBlockMatrix::insertSolutionVector template template -void fvBlockMatrix::insertDiagSource +void Foam::fvBlockMatrix::insertDiagSource ( const direction dir, fvMatrix& matrix @@ -162,7 +162,7 @@ void fvBlockMatrix::insertDiagSource template template -void fvBlockMatrix::insertUpperLower +void Foam::fvBlockMatrix::insertUpperLower ( const direction dir, const fvMatrix& matrix @@ -290,7 +290,7 @@ void fvBlockMatrix::insertUpperLower template template -void fvBlockMatrix::updateCouplingCoeffs +void Foam::fvBlockMatrix::updateCouplingCoeffs ( const direction dir, const fvMatrix& matrix @@ -299,7 +299,9 @@ void fvBlockMatrix::updateCouplingCoeffs const direction nCmpts = pTraits::nComponents; direction localDir = dir; - const GeometricField& psi = matrix.psi(); + const GeometricField& psi = + matrix.psi(); + forAll(psi.boundaryField(), patchI) { const fvPatchField& pf = psi.boundaryField()[patchI]; @@ -373,7 +375,7 @@ void fvBlockMatrix::updateCouplingCoeffs template template -void fvBlockMatrix::insertBlock +void Foam::fvBlockMatrix::insertBlock ( const direction dirI, const direction dirJ, @@ -471,7 +473,7 @@ void fvBlockMatrix::insertBlock template template -void fvBlockMatrix::insertBoundaryContributions +void Foam::fvBlockMatrix::insertBoundaryContributions ( const direction dirI, const direction dirJ, @@ -561,39 +563,30 @@ void fvBlockMatrix::insertBoundaryContributions template -void fvBlockMatrix::insertCouplingDiagSource +void Foam::fvBlockMatrix::insertCouplingDiag ( const direction dirI, const direction dirJ, - const fvScalarMatrix& matrix + const scalarField& coeffIJ ) { - // 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) + forAll(coeffIJ, cellI) { - blockDiag[cellI](dirI, dirJ) += diag[cellI]; - b[cellI](dirI) += source[cellI]; + blockDiag[cellI](dirI, dirJ) += coeffIJ[cellI]; } + + // Source compensation is done in function updateSourceCoupling() + // after all coupling terms are added. HJ, 27/Apr/2015 } template -void fvBlockMatrix::insertCouplingUpperLower +void Foam::fvBlockMatrix::insertCouplingUpperLower ( const direction dirI, const direction dirJ, @@ -707,7 +700,7 @@ Foam::fvBlockMatrix::fvBlockMatrix template template -void fvBlockMatrix::retrieveSolution +void Foam::fvBlockMatrix::retrieveSolution ( const direction dir, Field& xSingle @@ -736,7 +729,7 @@ void fvBlockMatrix::retrieveSolution template template -void fvBlockMatrix::insertEquation +void Foam::fvBlockMatrix::insertEquation ( const direction dir, fvMatrix& matrix @@ -751,7 +744,7 @@ void fvBlockMatrix::insertEquation template template -void fvBlockMatrix::insertBlockCoupling +void Foam::fvBlockMatrix::insertBlockCoupling ( const direction dirI, const direction dirJ, @@ -765,20 +758,22 @@ void fvBlockMatrix::insertBlockCoupling template -void fvBlockMatrix::insertEquationCoupling +void Foam::fvBlockMatrix::insertEquationCoupling ( const direction dirI, const direction dirJ, - const fvScalarMatrix& matrix + const scalarField& coeffIJ ) { - insertCouplingDiagSource(dirI, dirJ, matrix); - insertCouplingUpperLower(dirI, dirJ, matrix); + // Multiply coefficients by volume + scalarField coeffIJVol = coeffIJ*psi_.mesh().V(); + + insertCouplingDiag(dirI, dirJ, coeffIJVol); } template -void fvBlockMatrix::blockAdd +void Foam::fvBlockMatrix::blockAdd ( const direction dir, const scalarField& xSingle, @@ -793,7 +788,7 @@ void fvBlockMatrix::blockAdd template -void fvBlockMatrix::updateSourceCoupling() +void Foam::fvBlockMatrix::updateSourceCoupling() { // Eliminate off-diagonal block coefficients from the square diagonal // With this change, the segregated matrix can be assembled with complete @@ -827,215 +822,7 @@ void fvBlockMatrix::updateSourceCoupling() template -void 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(U, phi) interpolation scheme for consistency. - // VV, 21/July/2014. - const surfaceTensorField pi(mesh.Sf()*fvc::interpolate(U, phi, "(phi,U)")); -// const surfaceTensorField pi(fvc::interpolate(U, phi, "(phi,U)")*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(); - - Info << "Diag coeff: " << piDiag[125] << nl - << "Lower coeff: " << piLower[125] << nl - << "Upper coeff: " << piUpper[125] << endl; - - // Boundary contributions - treated explicitly 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 tensorField& pib = pi.boundaryField()[patchI]; - const fvsPatchScalarField& wb = tweights().boundaryField()[patchI]; - const unallocLabelList& fc = patch.faceCells(); - -// const vectorField internalCoeffs(Ub.valueInternalCoeffs(wb)); - - // Explicit diag boundary contribution -// forAll(Ub, faceI) -// { -// piSource[fc[faceI]] -= pib[faceI] & Ub[faceI]; -// piSource[fc[faceI]] -= Ub[faceI] & pib[faceI]; -// } - - // Hard coded zeroGradient if patch does not fix value - if (!U.boundaryField()[patchI].fixesValue()) - { - forAll(Ub, faceI) - { - piDiag[fc[faceI]] += pib[faceI]; - } - } - // Coupled patches treated implicitly - 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 - { - const vectorField boundaryCoeffs(Ub.valueBoundaryCoeffs(wb)); - - // Boundary contribution - forAll(Ub, faceI) - { - piSource[fc[faceI]] -= pib[faceI] & boundaryCoeffs[faceI]; - } - } - } - - // 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); -// blockSource[cellI](localDirI, localDirJ) += -// piSource[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 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 fvBlockMatrix::solve +Foam::BlockSolverPerformance Foam::fvBlockMatrix::solve ( const dictionary& solverControls ) @@ -1057,7 +844,7 @@ Foam::BlockSolverPerformance fvBlockMatrix::solve template -Foam::BlockSolverPerformance fvBlockMatrix::solve() +Foam::BlockSolverPerformance Foam::fvBlockMatrix::solve() { return solve(psi_.mesh().solutionDict().solverDict(psi_.name())); } diff --git a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H index d212a7253..972db44d8 100644 --- a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H +++ b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H @@ -26,11 +26,12 @@ Class Description fvBlockMatrix is an extension of fvMatrix for block coupled types. It holds - general insertion and retrieval tools for block coupling along with specific + general insertion and retrieval tools for block coupling and specific functions for p-U coupling. Author Vuko Vukcevic, FMENA Zagreb. + Update by Hrvoje Jasak SourceFiles fvBlockMatrix.C @@ -112,14 +113,13 @@ class fvBlockMatrix // Insertion functions for fvScalarMatrix into off-diagonal positions // (coupling matrices) - // These two functions are dubious. Reconsider. HJ, 17/July/2014. - //- Insert coupling matrix diag and source into this fvBlockMatrix - void insertCouplingDiagSource + //- Insert coupling matrix diag element into this fvBlockMatrix + void insertCouplingDiag ( const direction dirI, const direction dirJ, - const fvScalarMatrix& matrix + const scalarField& coeffIJ ); //- Insert coupling matrix lower and upper into this fvBlockMatrix @@ -144,8 +144,8 @@ class fvBlockMatrix const bool incrementColumnDir ); - //- Insert source and coupling coeffs of BlockLduSystem (obtained by - // implicit grad/div operator) into this fvBlockMatrix + //- Insert source and coupling coeffs of BlockLduSystem + // (eg. obtained by implicit grad/div operator) template void insertBoundaryContributions ( @@ -192,7 +192,7 @@ public: // Insertion and retrieval public tools - //- Retrieve part of internal (solved) field from this fvBlockMatrix + //- Retrieve part of internal field from this fvBlockMatrix template void retrieveSolution ( @@ -220,12 +220,13 @@ public: ); //- Insert scalar equation coupling into this fvBlockMatrix - // Not tested: VV, 10/July/2014 + // Source compensation is done in function updateSourceCoupling() + // after all coupling terms are added. HJ, 27/Apr/2015 void insertEquationCoupling ( const direction dirI, const direction dirJ, - const fvScalarMatrix& matrix + const scalarField& coeffIJ ); //- Add field into block field @@ -237,20 +238,11 @@ public: ); //- Update coupling of block system. - // Subtracts the block-coefficient coupling as specified by the user - // from the source, leaving the implicit update given by - // linearisation + // Subtracts the block-coefficient coupling as specified by the + // user from the source, leaving the implicit update given by + // 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