exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Feature: block coupled updates

Browse source
This commit is contained in:
Hrvoje Jasak 2015-04-27 18:09:57 +01:00 committed by Dominik Christ
parent 76b941e7ff
commit 95fb99d2a6
2 changed files with 44 additions and 265 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

273
src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C
View file

@ -30,7 +30,7 @@ License
template<class Type> template<class Type>
template<class fieldType> template<class fieldType>
void fvBlockMatrix<Type>::insertSolutionVector void Foam::fvBlockMatrix<Type>::insertSolutionVector
( (
const direction dir, const direction dir,
const Field<fieldType>& xSingle const Field<fieldType>& xSingle
@ -59,7 +59,7 @@ void fvBlockMatrix<Type>::insertSolutionVector
template<class Type> template<class Type>
template<class matrixType> template<class matrixType>
void fvBlockMatrix<Type>::insertDiagSource void Foam::fvBlockMatrix<Type>::insertDiagSource
( (
const direction dir, const direction dir,
fvMatrix<matrixType>& matrix fvMatrix<matrixType>& matrix
@ -162,7 +162,7 @@ void fvBlockMatrix<Type>::insertDiagSource
template<class Type> template<class Type>
template<class matrixType> template<class matrixType>
void fvBlockMatrix<Type>::insertUpperLower void Foam::fvBlockMatrix<Type>::insertUpperLower
( (
const direction dir, const direction dir,
const fvMatrix<matrixType>& matrix const fvMatrix<matrixType>& matrix
@ -290,7 +290,7 @@ void fvBlockMatrix<Type>::insertUpperLower
template<class Type> template<class Type>
template<class matrixType> template<class matrixType>
void fvBlockMatrix<Type>::updateCouplingCoeffs void Foam::fvBlockMatrix<Type>::updateCouplingCoeffs
( (
const direction dir, const direction dir,
const fvMatrix<matrixType>& matrix const fvMatrix<matrixType>& matrix
@ -299,7 +299,9 @@ void fvBlockMatrix<Type>::updateCouplingCoeffs
const direction nCmpts = pTraits<matrixType>::nComponents; const direction nCmpts = pTraits<matrixType>::nComponents;
direction localDir = dir; direction localDir = dir;
const GeometricField<matrixType, fvPatchField, volMesh>& psi = matrix.psi(); const GeometricField<matrixType, fvPatchField, volMesh>& psi =
matrix.psi();
forAll(psi.boundaryField(), patchI) forAll(psi.boundaryField(), patchI)
{ {
const fvPatchField<matrixType>& pf = psi.boundaryField()[patchI]; const fvPatchField<matrixType>& pf = psi.boundaryField()[patchI];
@ -373,7 +375,7 @@ void fvBlockMatrix<Type>::updateCouplingCoeffs
template<class Type> template<class Type>
template<class blockType, class fieldType> template<class blockType, class fieldType>
void fvBlockMatrix<Type>::insertBlock void Foam::fvBlockMatrix<Type>::insertBlock
( (
const direction dirI, const direction dirI,
const direction dirJ, const direction dirJ,
@ -471,7 +473,7 @@ void fvBlockMatrix<Type>::insertBlock
template<class Type> template<class Type>
template<class blockType, class fieldType> template<class blockType, class fieldType>
void fvBlockMatrix<Type>::insertBoundaryContributions void Foam::fvBlockMatrix<Type>::insertBoundaryContributions
( (
const direction dirI, const direction dirI,
const direction dirJ, const direction dirJ,
@ -561,39 +563,30 @@ void fvBlockMatrix<Type>::insertBoundaryContributions
template<class Type> template<class Type>
void fvBlockMatrix<Type>::insertCouplingDiagSource void Foam::fvBlockMatrix<Type>::insertCouplingDiag
( (
const direction dirI, const direction dirI,
const direction dirJ, 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 // Get reference to block diagonal of the block system
typename CoeffField<Type>::squareTypeField& blockDiag = typename CoeffField<Type>::squareTypeField& blockDiag =
this->diag().asSquare(); this->diag().asSquare();
// Get reference to this source field of the block system
Field<Type>& b = this->source();
// Set off-diagonal coefficient // Set off-diagonal coefficient
forAll(diag, cellI) forAll(coeffIJ, cellI)
{ {
blockDiag[cellI](dirI, dirJ) += diag[cellI]; blockDiag[cellI](dirI, dirJ) += coeffIJ[cellI];
b[cellI](dirI) += source[cellI];
} }
// Source compensation is done in function updateSourceCoupling()
// after all coupling terms are added. HJ, 27/Apr/2015
} }
template<class Type> template<class Type>
void fvBlockMatrix<Type>::insertCouplingUpperLower void Foam::fvBlockMatrix<Type>::insertCouplingUpperLower
( (
const direction dirI, const direction dirI,
const direction dirJ, const direction dirJ,
@ -707,7 +700,7 @@ Foam::fvBlockMatrix<Type>::fvBlockMatrix
template<class Type> template<class Type>
template<class fieldType> template<class fieldType>
void fvBlockMatrix<Type>::retrieveSolution void Foam::fvBlockMatrix<Type>::retrieveSolution
( (
const direction dir, const direction dir,
Field<fieldType>& xSingle Field<fieldType>& xSingle
@ -736,7 +729,7 @@ void fvBlockMatrix<Type>::retrieveSolution
template<class Type> template<class Type>
template<class matrixType> template<class matrixType>
void fvBlockMatrix<Type>::insertEquation void Foam::fvBlockMatrix<Type>::insertEquation
( (
const direction dir, const direction dir,
fvMatrix<matrixType>& matrix fvMatrix<matrixType>& matrix
@ -751,7 +744,7 @@ void fvBlockMatrix<Type>::insertEquation
template<class Type> template<class Type>
template<class blockType, class fieldType> template<class blockType, class fieldType>
void fvBlockMatrix<Type>::insertBlockCoupling void Foam::fvBlockMatrix<Type>::insertBlockCoupling
( (
const direction dirI, const direction dirI,
const direction dirJ, const direction dirJ,
@ -765,20 +758,22 @@ void fvBlockMatrix<Type>::insertBlockCoupling
template<class Type> template<class Type>
void fvBlockMatrix<Type>::insertEquationCoupling void Foam::fvBlockMatrix<Type>::insertEquationCoupling
( (
const direction dirI, const direction dirI,
const direction dirJ, const direction dirJ,
const fvScalarMatrix& matrix const scalarField& coeffIJ
) )
{ {
insertCouplingDiagSource(dirI, dirJ, matrix); // Multiply coefficients by volume
insertCouplingUpperLower(dirI, dirJ, matrix); scalarField coeffIJVol = coeffIJ*psi_.mesh().V();
insertCouplingDiag(dirI, dirJ, coeffIJVol);
} }
template<class Type> template<class Type>
void fvBlockMatrix<Type>::blockAdd void Foam::fvBlockMatrix<Type>::blockAdd
( (
const direction dir, const direction dir,
const scalarField& xSingle, const scalarField& xSingle,
@ -793,7 +788,7 @@ void fvBlockMatrix<Type>::blockAdd
template<class Type> template<class Type>
void fvBlockMatrix<Type>::updateSourceCoupling() void Foam::fvBlockMatrix<Type>::updateSourceCoupling()
{ {
// Eliminate off-diagonal block coefficients from the square diagonal // Eliminate off-diagonal block coefficients from the square diagonal
// With this change, the segregated matrix can be assembled with complete // With this change, the segregated matrix can be assembled with complete
@ -827,215 +822,7 @@ void fvBlockMatrix<Type>::updateSourceCoupling()
template<class Type> template<class Type>
void fvBlockMatrix<Type>::insertPicardTensor Foam::BlockSolverPerformance<Type> Foam::fvBlockMatrix<Type>::solve
(
const direction UEqnDir,
const volVectorField& U,
const surfaceScalarField& phi
)
{
const direction nCmpts = pTraits<vector>::nComponents;
direction localDirI = UEqnDir;
direction localDirJ = UEqnDir;
// Sanity check
{
const direction blockMatrixSize = pTraits<Type>::nComponents;
if (nCmpts > blockMatrixSize)
{
FatalErrorIn
(
"void fvBlockMatrix<Type>::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<surfaceInterpolationScheme<vector> > tinterpScheme =
fvc::scheme<vector>(mesh, U.name());
tmp<surfaceScalarField> 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<vector, vector> piSystem(mesh);
// Get references to ldu fields of pi block system always as square
typename CoeffField<vector>::squareTypeField& piDiag =
piSystem.diag().asSquare();
typename CoeffField<vector>::squareTypeField& piUpper =
piSystem.upper().asSquare();
typename CoeffField<vector>::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<vector>::squareTypeField& pipCoupleUpper =
piSystem.coupleUpper()[patchI].asSquare();
typename CoeffField<vector>::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<Type>::squareTypeField& blockDiag =
this->diag().asSquare();
typename CoeffField<Type>::squareTypeField& blockUpper =
this->upper().asSquare();
typename CoeffField<Type>::squareTypeField& blockLower =
this->lower().asSquare();
Field<Type>& 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<Type>::squareTypeField& pcoupleUpper =
this->coupleUpper()[patchI].asSquare();
typename CoeffField<Type>::squareTypeField& pcoupleLower =
this->coupleLower()[patchI].asSquare();
const typename CoeffField<vector>::squareTypeField& pipcu =
piSystem.coupleUpper()[patchI].asSquare();
const typename CoeffField<vector>::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<class Type>
Foam::BlockSolverPerformance<Type> fvBlockMatrix<Type>::solve
( (
const dictionary& solverControls const dictionary& solverControls
) )
@ -1057,7 +844,7 @@ Foam::BlockSolverPerformance<Type> fvBlockMatrix<Type>::solve
template<class Type> template<class Type>
Foam::BlockSolverPerformance<Type> fvBlockMatrix<Type>::solve() Foam::BlockSolverPerformance<Type> Foam::fvBlockMatrix<Type>::solve()
{ {
return solve(psi_.mesh().solutionDict().solverDict(psi_.name())); return solve(psi_.mesh().solutionDict().solverDict(psi_.name()));
} }

36
src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H
View file

@ -26,11 +26,12 @@ Class
Description Description
fvBlockMatrix is an extension of fvMatrix for block coupled types. It holds 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. functions for p-U coupling.
Author Author
Vuko Vukcevic, FMENA Zagreb. Vuko Vukcevic, FMENA Zagreb.
Update by Hrvoje Jasak
SourceFiles SourceFiles
fvBlockMatrix.C fvBlockMatrix.C
@ -112,14 +113,13 @@ class fvBlockMatrix
// Insertion functions for fvScalarMatrix into off-diagonal positions // Insertion functions for fvScalarMatrix into off-diagonal positions
// (coupling matrices) // (coupling matrices)
// These two functions are dubious. Reconsider. HJ, 17/July/2014.
//- Insert coupling matrix diag and source into this fvBlockMatrix //- Insert coupling matrix diag element into this fvBlockMatrix
void insertCouplingDiagSource void insertCouplingDiag
( (
const direction dirI, const direction dirI,
const direction dirJ, const direction dirJ,
const fvScalarMatrix& matrix const scalarField& coeffIJ
); );
//- Insert coupling matrix lower and upper into this fvBlockMatrix //- Insert coupling matrix lower and upper into this fvBlockMatrix
@ -144,8 +144,8 @@ class fvBlockMatrix
const bool incrementColumnDir const bool incrementColumnDir
); );
//- Insert source and coupling coeffs of BlockLduSystem (obtained by //- Insert source and coupling coeffs of BlockLduSystem
// implicit grad/div operator) into this fvBlockMatrix // (eg. obtained by implicit grad/div operator)
template<class blockType, class fieldType> template<class blockType, class fieldType>
void insertBoundaryContributions void insertBoundaryContributions
( (
@ -192,7 +192,7 @@ public:
// Insertion and retrieval public tools // Insertion and retrieval public tools
//- Retrieve part of internal (solved) field from this fvBlockMatrix //- Retrieve part of internal field from this fvBlockMatrix
template<class fieldType> template<class fieldType>
void retrieveSolution void retrieveSolution
( (
@ -220,12 +220,13 @@ public:
); );
//- Insert scalar equation coupling into this fvBlockMatrix //- 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 void insertEquationCoupling
( (
const direction dirI, const direction dirI,
const direction dirJ, const direction dirJ,
const fvScalarMatrix& matrix const scalarField& coeffIJ
); );
//- Add field into block field //- Add field into block field
@ -237,20 +238,11 @@ public:
); );
//- Update coupling of block system. //- Update coupling of block system.
// Subtracts the block-coefficient coupling as specified by the user // Subtracts the block-coefficient coupling as specified by the
// from the source, leaving the implicit update given by // user from the source, leaving the implicit update given by
// linearisation // linearisation
void updateSourceCoupling(); 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 // Solver calls for fvBlockMatrix