From 0c31cec90dd9f085087395715681fb4008a25f84 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 18 Aug 2011 17:27:12 +0100 Subject: [PATCH 1/6] Version update --- .../conjugateHeatFoam/conjugateCavity/system/fvSolution | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tutorials/coupled/conjugateHeatFoam/conjugateCavity/system/fvSolution b/tutorials/coupled/conjugateHeatFoam/conjugateCavity/system/fvSolution index f4149aab5..379da69cf 100644 --- a/tutorials/coupled/conjugateHeatFoam/conjugateCavity/system/fvSolution +++ b/tutorials/coupled/conjugateHeatFoam/conjugateCavity/system/fvSolution @@ -22,7 +22,7 @@ solvers preconditioner DIC; tolerance 1e-06; relTol 0; - }; + } U { @@ -30,7 +30,7 @@ solvers preconditioner DILU; tolerance 1e-05; relTol 0; - }; + } T+T { @@ -41,7 +41,7 @@ solvers maxIter 1000; tolerance 1e-6; relTol 0.0; - }; + } } PISO From b1a9acbf012c22bd96f7fd5fdef9112887ed9bc5 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 18 Aug 2011 17:28:49 +0100 Subject: [PATCH 2/6] Version update --- .../coupled/conjugateHeatFoam/heatedBlock/system/fvSolution | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tutorials/coupled/conjugateHeatFoam/heatedBlock/system/fvSolution b/tutorials/coupled/conjugateHeatFoam/heatedBlock/system/fvSolution index 3f610bfb4..204128d05 100644 --- a/tutorials/coupled/conjugateHeatFoam/heatedBlock/system/fvSolution +++ b/tutorials/coupled/conjugateHeatFoam/heatedBlock/system/fvSolution @@ -15,9 +15,7 @@ FoamFile // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers -{ - T ICCG 1e-06 0; -} +{} SIMPLE { From a63606d91eba481748ee1f34745f01362ee8a782 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 19 Aug 2011 10:12:03 +0100 Subject: [PATCH 3/6] Added debug option for uncovered GGI faces --- .../polyPatches/constraint/ggi/ggiPolyPatch.C | 51 ++++++++++++++++++- 1 file changed, 50 insertions(+), 1 deletion(-) diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C index 3e2424f87..102d462af 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C @@ -37,6 +37,8 @@ Contributor #include "polyPatchID.H" #include "ZoneIDs.H" #include "SubField.H" +#include "Time.H" +#include "indirectPrimitivePatch.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -854,6 +856,52 @@ void Foam::ggiPolyPatch::calcTransforms() << "Master: " << name() << " Slave: " << shadowName() << endl; + if (patchToPatch().uncoveredMasterFaces().size() > 0) + { + // Write uncovered master faces + Info<< "Writing uncovered master faces for patch " + << name() << " as VTK." << endl; + + const polyMesh& mesh = boundaryMesh().mesh(); + + fileName fvPath(mesh.time().path()/"VTK"); + mkDir(fvPath); + + indirectPrimitivePatch::writeVTK + ( + fvPath/fileName("uncoveredGgiFaces" + name()), + IndirectList + ( + localFaces(), + patchToPatch().uncoveredMasterFaces() + ), + localPoints() + ); + } + + if (patchToPatch().uncoveredSlaveFaces().size() > 0) + { + // Write uncovered master faces + Info<< "Writing uncovered shadow faces for patch " + << shadowName() << " as VTK." << endl; + + const polyMesh& mesh = boundaryMesh().mesh(); + + fileName fvPath(mesh.time().path()/"VTK"); + mkDir(fvPath); + + indirectPrimitivePatch::writeVTK + ( + fvPath/fileName("uncoveredGgiFaces" + shadowName()), + IndirectList + ( + shadow().localFaces(), + patchToPatch().uncoveredSlaveFaces() + ), + shadow().localPoints() + ); + } + // Check for bridge overlap if (!bridgeOverlap()) { @@ -864,7 +912,8 @@ void Foam::ggiPolyPatch::calcTransforms() ) { FatalErrorIn("label ggiPolyPatch::shadowIndex() const") - << "ggi patch " << name() << " has " + << "ggi patch " << name() << " with shadow " + << shadowName() << " has " << patchToPatch().uncoveredMasterFaces().size() << " uncovered master faces and " << patchToPatch().uncoveredSlaveFaces().size() From 99169ff574f8086e5c994da7ca60c16d99916a52 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 24 Aug 2011 11:30:00 +0100 Subject: [PATCH 4/6] Header update --- .../interFoam/laminar/damBreak/constant/polyMesh/boundary | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tutorials/multiphase/interFoam/laminar/damBreak/constant/polyMesh/boundary b/tutorials/multiphase/interFoam/laminar/damBreak/constant/polyMesh/boundary index 65fa6d3a7..8a00e1577 100644 --- a/tutorials/multiphase/interFoam/laminar/damBreak/constant/polyMesh/boundary +++ b/tutorials/multiphase/interFoam/laminar/damBreak/constant/polyMesh/boundary @@ -1,6 +1,6 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | -| \\ / F ield | OpenFOAM Extend Project: Open Source CFD | +| \\ / F ield | OpenFOAM Extend Project: Open source CFD | | \\ / O peration | Version: 1.6-ext | | \\ / A nd | Web: www.extend-project.de | | \\/ M anipulation | | From 105f6dd539411b5da3c126bb628675eb6c9ae0b0 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 24 Aug 2011 11:35:47 +0100 Subject: [PATCH 5/6] Block solver development --- .../blockCoupledScalarTransportFoam.C | 44 +-- .../blockMatrixTools.H | 130 --------- .../createFields.H | 36 +-- src/OpenFOAM/Make/files | 2 +- .../BlockLduMatrix/BlockLduMatrix.C | 3 +- .../BlockBiCGStab/BlockBiCGStabSolver.C | 2 +- .../BlockLduSolvers/BlockCG/BlockCGSolver.C | 2 +- .../BlockGMRES/BlockGMRESSolver.C | 2 +- .../BlockGaussSeidel/BlockGaussSeidelSolver.C | 6 +- .../BlockGaussSeidel/BlockGaussSeidelSolver.H | 2 +- .../BlockIterativeSolver.C | 3 +- .../BlockLduSolver/BlockLduSolver.C | 52 +++- .../BlockLduSolver/BlockLduSolver.H | 11 +- .../Segregated/SegregatedSolver.C | 43 +-- .../lduInterface/processorLduInterface.H | 3 +- .../GAMGAgglomerateLduAddressing.C | 2 +- .../GAMG/GAMGSolverAgglomerateMatrix.C | 13 +- .../primitiveMeshCheck/primitiveMeshCheck.C | 2 +- .../BlockLduMatrices/blockVectorNSolvers.C | 23 +- .../OpenFOAM/Fields/VectorNFieldTypes.H | 8 + .../blockMatrixTools}/blockMatrixTools.C | 270 ++++++++++++++++-- .../blockMatrixTools/blockMatrixTools.H | 179 ++++++++++++ .../processorFvPatchVectorNFields.C | 112 ++++---- .../fvsPatchFields/fvsPatchVectorNFieldsFwd.H | 4 + .../general/findRefCell/findRefCell.C | 5 +- .../fvMatrices/fvMatrix/fvMatrix.H | 14 + 26 files changed, 637 insertions(+), 336 deletions(-) delete mode 100644 applications/solvers/coupled/blockCoupledScalarTransportFoam/blockMatrixTools.H rename {applications/solvers/coupled/blockCoupledScalarTransportFoam => src/VectorN/finiteVolume/blockMatrixTools}/blockMatrixTools.C (53%) create mode 100644 src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.H diff --git a/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C b/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C index 6dc9fab52..76f92eb45 100644 --- a/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C +++ b/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C @@ -95,25 +95,28 @@ int main(int argc, char *argv[]) // Prepare block system BlockLduMatrix blockM(mesh); + //- Transfer the coupled interface list for processor/cyclic/etc. + // boundaries + blockM.interfaces() = blockT.boundaryField().blockInterfaces(); + // Grab block diagonal and set it to zero Field& d = blockM.diag().asSquare(); d = tensor2::zero; // Grab linear off-diagonal and set it to zero - Field& u = blockM.upper().asLinear(); Field& l = blockM.lower().asLinear(); + Field& u = blockM.upper().asLinear(); u = vector2::zero; l = vector2::zero; vector2Field& blockX = blockT.internalField(); - // vector2Field blockX(mesh.nCells(), vector2::zero); vector2Field blockB(mesh.nCells(), vector2::zero); //- Inset equations into block Matrix blockMatrixTools::insertEquation(0, TEqn, blockM, blockX, blockB); blockMatrixTools::insertEquation(1, TsEqn, blockM, blockX, blockB); - //- Add off-diagonal terms and remove from Block source + //- Add off-diagonal terms and remove from block source forAll(d, i) { d[i](0,1) = -alpha.value()*mesh.V()[i]; @@ -123,44 +126,13 @@ int main(int argc, char *argv[]) blockB[i][1] -= alpha.value()*blockX[i][0]*mesh.V()[i]; } - //- Transfer the coupled interface list for processor/cyclic/etc. - // boundaries - blockM.interfaces() = blockT.boundaryField().blockInterfaces(); - - //- Transfer the coupled interface coefficients - forAll(mesh.boundaryMesh(), patchI) - { - if (blockM.interfaces().set(patchI)) - { - Field& coupledLower = - blockM.coupleLower()[patchI].asLinear(); - - Field& coupledUpper = - blockM.coupleUpper()[patchI].asLinear(); - - const scalarField& TLower = TEqn.internalCoeffs()[patchI]; - const scalarField& TUpper = TEqn.boundaryCoeffs()[patchI]; - - const scalarField& TsLower = - TsEqn.internalCoeffs()[patchI]; - - const scalarField& TsUpper = - TsEqn.boundaryCoeffs()[patchI]; - - blockMatrixTools::blockInsert(0, TLower, coupledLower); - blockMatrixTools::blockInsert(1, TsLower, coupledLower); - blockMatrixTools::blockInsert(0, TUpper, coupledUpper); - blockMatrixTools::blockInsert(1, TsUpper, coupledUpper); - } - } - //- Block coupled solver call BlockSolverPerformance solverPerf = BlockLduSolver::New ( - word("blockVar"), + blockT.name(), blockM, - mesh.solutionDict().solver("blockVar") + mesh.solutionDict().solver(blockT.name()) )->solve(blockX, blockB); solverPerf.print(); diff --git a/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockMatrixTools.H b/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockMatrixTools.H deleted file mode 100644 index cb372a92a..000000000 --- a/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockMatrixTools.H +++ /dev/null @@ -1,130 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright held by original author - \\/ 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 2 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, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -InNamespace - Foam:: - -Description - Block matrix insertion and retrieval tools - -SourceFiles - blockMatrixTools.C - -\*---------------------------------------------------------------------------*/ - - -#ifndef blockMatrixTools_H -#define blockMatrixTools_H - -#include "blockLduMatrices.H" -#include "fvMatrices.H" -#include "surfaceFieldsFwd.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Namespace blockMatrixTools functions Declaration -\*---------------------------------------------------------------------------*/ - -namespace blockMatrixTools -{ - //- Insert field into block field - template - void blockInsert - ( - const direction dir, - const scalarField& x, - Field& blockX - ); - - //- Retrieve field from block field - template - void blockRetrieve - ( - const direction dir, - scalarField& x, - const Field& blockX - ); - - //- Insert matrix diagonal and source into the block system - template - void insertDiagSource - ( - const direction dir, - const fvScalarMatrix& m, - BlockLduMatrix& blockM, - Field& blockB - ); - - // Insert matrix into the block system - template - void insertUpperLower - ( - const direction dir, - const fvScalarMatrix& m, - BlockLduMatrix& blockM - ); - - // Insert matrix into the block system - template - void insertEquation - ( - const direction dir, - const fvScalarMatrix& m, - BlockLduMatrix& blockM, - Field& blockX, - Field& blockB - ); - - // 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 - template - void updateCoupling - ( - BlockLduMatrix& blockM, - Field& x, - Field& b - ); -} - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#ifdef NoRepository -# include "blockMatrixTools.C" -#endif - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H b/applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H index 69599c465..2f4eb6aaa 100644 --- a/applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H +++ b/applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H @@ -28,7 +28,7 @@ ); -// Working coupled solution field + // Working coupled solution field Info<< "Creating field blockT\n" << endl; volVector2Field blockT ( @@ -41,7 +41,7 @@ IOobject::NO_WRITE ), mesh, - dimensionedVector2(word(), dimless, vector2::zero) + dimensionedVector2("zero", dimless, vector2::zero) ); { @@ -53,38 +53,6 @@ blockT.write(); -// Removed comparison with the reference field. HJ, 17/Jun/2010 -#if (0) - Info<< "Reading field Tref\n" << endl; - volScalarField Tref - ( - IOobject - ( - "Tref", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - - Info<< "Reading field Tsref\n" << endl; - volScalarField Tsref - ( - IOobject - ( - "Tsref", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); -#endif - Info<< "Reading field U\n" << endl; volVectorField U diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index 435ad8c6b..2bd488e86 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -642,11 +642,11 @@ matrices/blockLduMatrix/BlockLduSmoothers/BlockILUSmoother/blockILUSmoothers.C matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/blockLduSolvers.C matrices/blockLduMatrix/BlockLduSolvers/BlockDiagonal/blockDiagonalSolvers.C -/*HJ matrices/blockLduMatrix/BlockLduSolvers/Segregated/segregatedSolvers.C */ matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/blockGaussSeidelSolvers.C matrices/blockLduMatrix/BlockLduSolvers/BlockCG/blockCGSolvers.C matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/blockBiCGStabSolvers.C matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/blockGMRESSolvers.C +matrices/blockLduMatrix/BlockLduSolvers/Segregated/segregatedSolvers.C matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterface/BlockLduInterfaceFields.C diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C index afb1dcf92..bdce6ccb4 100644 --- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C +++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C @@ -343,7 +343,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const BlockLduMatrix& ldum) { if (ldum.diagPtr_) { - os.writeKeyword("diagonal") << *ldum.diagPtr_ << token::END_STATEMENT << nl; + os.writeKeyword("diagonal") + << *ldum.diagPtr_ << token::END_STATEMENT << nl; } else { diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C index 81fe45c6b..993836dae 100644 --- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C +++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C @@ -93,7 +93,7 @@ Foam::BlockBiCGStabSolver::solve // Check convergence, solve if not converged - if (!solverPerf.checkConvergence(this->tolerance(), this->relTolerance())) + if (!stop(solverPerf)) { scalar rho = this->great_; scalar rhoOld = rho; diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C index 1c07bb7b0..ec8c2b295 100644 --- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C +++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C @@ -92,7 +92,7 @@ typename Foam::BlockSolverPerformance Foam::BlockCGSolver::solve // Check convergence, solve if not converged - if (!solverPerf.checkConvergence(this->tolerance(), this->relTolerance())) + if (!stop(solverPerf)) { scalar rho = this->great_; scalar rhoOld = rho; diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C index 6415cf7bd..1a651ab32 100644 --- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C +++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C @@ -128,7 +128,7 @@ Foam::BlockGMRESSolver::solve // Check convergence, solve if not converged - if (!solverPerf.checkConvergence(this->tolerance(), this->relTolerance())) + if (!stop(solverPerf)) { // Create the Hesenberg matrix scalarSquareMatrix H(nDirs_, 0); diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C index 0819dc9d3..e69fe326f 100644 --- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C +++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C @@ -49,7 +49,7 @@ Foam::BlockGaussSeidelSolver::BlockGaussSeidelSolver dict ), gs_(matrix), - nResSweeps_(readLabel(this->dict().lookup("nResSweeps"))) + nSweeps_(readLabel(this->dict().lookup("nSweeps"))) {} @@ -86,13 +86,13 @@ Foam::BlockGaussSeidelSolver::solve // Check convergence, solve if not converged - if (!solverPerf.checkConvergence(this->tolerance(), this->relTolerance())) + if (!stop(solverPerf)) { // Iteration loop do { - for (label i = 0; i < nResSweeps_; i++) + for (label i = 0; i < nSweeps_; i++) { gs_.precondition(x, b); diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.H b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.H index a29fbcc0c..4585b6a37 100644 --- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.H +++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.H @@ -64,7 +64,7 @@ class BlockGaussSeidelSolver BlockGaussSeidelPrecon gs_; //- Number of sweeps before evaluating residual - label nResSweeps_; + label nSweeps_; diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.C index d4e8dd0e4..bb7a8d690 100644 --- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.C +++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.C @@ -82,7 +82,8 @@ Foam::scalar Foam::BlockIterativeSolver::normFactor if (BlockLduMatrix::debug >= 2) { - Info<< "Iterative solver normalisation factor = " << normFactor << endl; + Info<< "Iterative solver normalisation factor = " + << normFactor << endl; } return normFactor; diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.C index bd3ec1be7..3e53dd209 100644 --- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.C +++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.C @@ -65,7 +65,7 @@ void Foam::BlockLduSolver::readControl ")", dict ) << "Cannot read control " << controlName - << abort(FatalIOError); + << exit(FatalIOError); } } @@ -112,6 +112,52 @@ Foam::BlockLduSolver::New { word solverName(dict.lookup("solver")); + if (matrix.diagonal()) + { + return autoPtr > + ( + new BlockDiagonalSolver(fieldName, matrix, dict) + ); + } + else if (matrix.symmetric() || matrix.asymmetric()) + { + return BlockLduSolver::New + ( + solverName, + fieldName, + matrix, + dict + ); + } + else + { + FatalErrorIn + ( + "BlockLduSolver::New\n" + "(\n" + " const word& fieldName,\n" + " BlockLduMatrix& matrix,\n" + " const dictionary& dict\n" + ")" + ) << "cannot solve incomplete matrix, " + "no diagonal or off-diagonal coefficient" + << exit(FatalError); + + return autoPtr >(NULL); + } +} + + +template +Foam::autoPtr > +Foam::BlockLduSolver::New +( + const word& solverName, + const word& fieldName, + BlockLduMatrix& matrix, + const dictionary& dict +) +{ if (matrix.diagonal()) { return autoPtr > @@ -143,6 +189,7 @@ Foam::BlockLduSolver::New ( "BlockLduSolver::New\n" "(\n" + " const word& solverName,\n" " const word& fieldName,\n" " BlockLduMatrix& matrix,\n" " const dictionary& dict\n" @@ -173,6 +220,7 @@ Foam::BlockLduSolver::New ( "BlockLduSolver::New\n" "(\n" + " const word& solverName,\n" " const word& fieldName,\n" " BlockLduMatrix& matrix,\n" " const dictionary& dict\n" @@ -189,6 +237,7 @@ Foam::BlockLduSolver::New ( "BlockLduSolver::New\n" "(\n" + " const word& solverName,\n" " const word& fieldName,\n" " BlockLduMatrix& matrix,\n" " const dictionary& dict\n" @@ -217,6 +266,7 @@ Foam::BlockLduSolver::New ( "BlockLduSolver::New\n" "(\n" + " const word& solverName,\n" " const word& fieldName,\n" " BlockLduMatrix& matrix,\n" " const dictionary& dict\n" diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H index 62aa78379..dc99b7e8b 100644 --- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H +++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H @@ -172,7 +172,7 @@ public: // Selectors - //- Return a new solver + //- Return a new solver from a dictionary static autoPtr > New ( const word& fieldName, @@ -180,6 +180,15 @@ public: const dictionary& dict ); + //- Return a new solver given type + static autoPtr > New + ( + const word& solverName, + const word& fieldName, + BlockLduMatrix& matrix, + const dictionary& dict + ); + // Destructor diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/Segregated/SegregatedSolver.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/Segregated/SegregatedSolver.C index 46aa30edf..a5cdeb804 100644 --- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/Segregated/SegregatedSolver.C +++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/Segregated/SegregatedSolver.C @@ -74,21 +74,6 @@ Foam::BlockSolverPerformance Foam::SegregatedSolver::solve switchingDiag = true; } - // Check switching upper - bool switchingUpper = false; - - if (blockMatrix.thereIsUpper()) - { - if (blockMatrix.upper().activeType() == blockCoeffBase::SCALAR) - { - scalarMatrix_.upper() = blockMatrix.upper().asScalar(); - } - else - { - switchingUpper = true; - } - } - // Check switching lower bool switchingLower = false; @@ -104,6 +89,21 @@ Foam::BlockSolverPerformance Foam::SegregatedSolver::solve } } + // Check switching upper + bool switchingUpper = false; + + if (blockMatrix.thereIsUpper()) + { + if (blockMatrix.upper().activeType() == blockCoeffBase::SCALAR) + { + scalarMatrix_.upper() = blockMatrix.upper().asScalar(); + } + else + { + switchingUpper = true; + } + } + // Decouple quadratic coupling by multiplying out the square coefficient // coupling Field* dBPtr = NULL; @@ -120,7 +120,7 @@ Foam::BlockSolverPerformance Foam::SegregatedSolver::solve } // Prepare solver performance - word segSolverName(this->dict().lookup("solver")); + word segSolverName(this->dict().lookup("segSolver")); BlockSolverPerformance solverPerf ( @@ -156,20 +156,21 @@ Foam::BlockSolverPerformance Foam::SegregatedSolver::solve scalarMatrix_.diag() = blockMatrix.diag().component(cmpt); } - if (switchingUpper) - { - scalarMatrix_.upper() = blockMatrix.upper().component(cmpt); - } - if (switchingLower) { scalarMatrix_.lower() = blockMatrix.lower().component(cmpt); } + if (switchingUpper) + { + scalarMatrix_.upper() = blockMatrix.upper().component(cmpt); + } + // Call the scalar solver BlockSolverPerformance scalarPerf = blockScalarSolver::New ( + segSolverName, this->fieldName(), scalarMatrix_, this->dict() diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/processorLduInterface.H b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/processorLduInterface.H index b5083cd56..252bc23ec 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/processorLduInterface.H +++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/processorLduInterface.H @@ -139,7 +139,8 @@ public: UList& ) const; - //- Raw field receive function with data compression returning field + //- Raw field receive function with data compression, + // returning a field template tmp > compressedReceive ( diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C index 3b8c25416..47967a673 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C @@ -290,7 +290,7 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing fineInterfaceAddr[inti] ).ptr() ); - + coarseInterfaceAddr[inti] = coarseInterfaces[inti].faceCells(); } } diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C index bcc721d29..0a4dc97e3 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C @@ -46,9 +46,14 @@ void Foam::GAMGSolver::agglomerateMatrix(const label fineLevelIndex) const labelList& faceRestrictAddr = agglomeration_.faceRestrictAddressing(fineLevelIndex); - // Coarse matrix diagonal initialised by restricting the finer mesh diagonal + // Coarse matrix diagonal initialised by restricting the fine mesh diagonal scalarField& coarseDiag = coarseMatrix.diag(); - agglomeration_.restrictField(coarseDiag, fineMatrix.diag(), fineLevelIndex); + agglomeration_.restrictField + ( + coarseDiag, + fineMatrix.diag(), + fineLevelIndex + ); // Get reference to fine-level interfaces const lduInterfaceFieldPtrsList& fineInterfaces = @@ -96,7 +101,7 @@ void Foam::GAMGSolver::agglomerateMatrix(const label fineLevelIndex) { if (fineInterfaces.set(inti)) { - const GAMGInterface& coarseInterface = + const GAMGInterface& coarseInterface = refCast ( agglomeration_.interfaceLevel(fineLevelIndex + 1)[inti] @@ -182,7 +187,7 @@ void Foam::GAMGSolver::agglomerateMatrix(const label fineLevelIndex) } } } - else // ... Otherwise it is symmetric so agglomerate just the upper + else // ... Otherwise it is symmetric so agglomerate just the upper { // Get off-diagonal matrix coefficients const scalarField& fineUpper = fineMatrix.upper(); diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C index 50b7d967b..f1c8d5fa5 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C @@ -537,7 +537,7 @@ bool Foam::primitiveMesh::checkFaceOrthogonality << ::acos(minDDotS)/mathematicalConstant::pi*180.0 << " average: " << ::acos(sumDDotS/neiSize)/mathematicalConstant::pi*180.0 - << " Threshold = " << nonOrthThreshold_ + << " Threshold = " << nonOrthThreshold_ << endl; } } diff --git a/src/VectorN/OpenFOAM/BlockLduMatrices/blockVectorNSolvers.C b/src/VectorN/OpenFOAM/BlockLduMatrices/blockVectorNSolvers.C index 44336037e..37e1cf78f 100644 --- a/src/VectorN/OpenFOAM/BlockLduMatrices/blockVectorNSolvers.C +++ b/src/VectorN/OpenFOAM/BlockLduMatrices/blockVectorNSolvers.C @@ -43,6 +43,7 @@ License #include "BlockCGSolver.H" #include "BlockGaussSeidelSolver.H" #include "BlockGMRESSolver.H" +#include "SegregatedSolver.H" #include "VectorTensorNFields.H" #include "ExpandTensorN.H" @@ -105,23 +106,29 @@ defineTemplateRunTimeSelectionTable \ typedef BlockDiagonalSolver block##Type##DiagonalSolver; \ defineNamedTemplateTypeNameAndDebug(block##Type##DiagonalSolver, 0); \ \ -typedef BlockBiCGStabSolver block##Type##BiCGStabSolver; \ -makeBlockSolverTypeName(block##Type##BiCGStabSolver); \ -addSolverToBlockMatrix(Type, block##Type##BiCGStabSolver, symMatrix); \ -addSolverToBlockMatrix(Type, block##Type##BiCGStabSolver, asymMatrix); \ +typedef BlockGaussSeidelSolver block##Type##GaussSeidelSolver; \ +makeBlockSolverTypeName(block##Type##GaussSeidelSolver); \ +addSolverToBlockMatrix(Type, block##Type##GaussSeidelSolver, symMatrix); \ +addSolverToBlockMatrix(Type, block##Type##GaussSeidelSolver, asymMatrix); \ \ typedef BlockCGSolver block##Type##CGSolver; \ makeBlockSolverTypeName(block##Type##CGSolver); \ addSolverToBlockMatrix(Type, block##Type##CGSolver, symMatrix); \ \ -typedef BlockGaussSeidelSolver block##Type##GaussSeidelSolver; \ -makeBlockSolverTypeName(block##Type##GaussSeidelSolver); \ -addSolverToBlockMatrix(Type, block##Type##GaussSeidelSolver, symMatrix); \ +typedef BlockBiCGStabSolver block##Type##BiCGStabSolver; \ +makeBlockSolverTypeName(block##Type##BiCGStabSolver); \ +addSolverToBlockMatrix(Type, block##Type##BiCGStabSolver, symMatrix); \ +addSolverToBlockMatrix(Type, block##Type##BiCGStabSolver, asymMatrix); \ \ typedef BlockGMRESSolver block##Type##GMRESSolver; \ makeBlockSolverTypeName(block##Type##GMRESSolver); \ addSolverToBlockMatrix(Type, block##Type##GMRESSolver, symMatrix); \ -addSolverToBlockMatrix(Type, block##Type##GMRESSolver, asymMatrix); +addSolverToBlockMatrix(Type, block##Type##GMRESSolver, asymMatrix); \ + \ +typedef SegregatedSolver block##Type##SegregatedSolver; \ +makeBlockSolverTypeName(block##Type##SegregatedSolver); \ +addSolverToBlockMatrix(Type, block##Type##SegregatedSolver, symMatrix); \ +addSolverToBlockMatrix(Type, block##Type##SegregatedSolver, asymMatrix); forAllVectorNTypes(makeSolver) diff --git a/src/VectorN/OpenFOAM/Fields/VectorNFieldTypes.H b/src/VectorN/OpenFOAM/Fields/VectorNFieldTypes.H index 971b6dd09..101db9542 100644 --- a/src/VectorN/OpenFOAM/Fields/VectorNFieldTypes.H +++ b/src/VectorN/OpenFOAM/Fields/VectorNFieldTypes.H @@ -32,46 +32,54 @@ Description #define VectorNFieldTypes_H #include "vector2.H" +#include "vector3.H" #include "vector4.H" #include "vector6.H" #include "vector8.H" #include "tensor2.H" +#include "tensor3.H" #include "tensor4.H" #include "tensor6.H" #include "tensor8.H" #include "diagTensor2.H" +#include "diagTensor3.H" #include "diagTensor4.H" #include "diagTensor6.H" #include "diagTensor8.H" #define forAllVectorNTypes(m, args...) \ m(vector2, Vector2, args) \ + m(vector3, Vector3, args) \ m(vector4, Vector4, args) \ m(vector6, Vector6, args) \ m(vector8, Vector8, args) #define forAllTensorNTypes(m, args...) \ m(tensor2, Tensor2, args) \ + m(tensor3, Tensor3, args) \ m(tensor4, Tensor4, args) \ m(tensor6, Tensor6, args) \ m(tensor8, Tensor8, args) #define forAllDiagTensorNTypes(m, args...) \ m(diagTensor2, DiagTensor2, args) \ + m(diagTensor3, DiagTensor3, args) \ m(diagTensor4, DiagTensor4, args) \ m(diagTensor6, DiagTensor6, args) \ m(diagTensor8, DiagTensor8, args) #define forAllSphericalTensorNTypes(m, args...) \ m(sphericalTensor2, SphericalTensor2, args) \ + m(sphericalTensor3, SphericalTensor3, args) \ m(sphericalTensor4, SphericalTensor4, args) \ m(sphericalTensor6, SphericalTensor6, args) \ m(sphericalTensor8, SphericalTensor8, args) #define forAllVectorTensorNTypes(m, args...) \ m(tensor2, diagTensor2, sphericalTensor2, vector2, scalar, args) \ + m(tensor3, diagTensor3, sphericalTensor3, vector3, scalar, args) \ m(tensor4, diagTensor4, sphericalTensor4, vector4, scalar, args) \ m(tensor6, diagTensor6, sphericalTensor6, vector6, scalar, args) \ m(tensor8, diagTensor8, sphericalTensor8, vector8, scalar, args) diff --git a/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockMatrixTools.C b/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.C similarity index 53% rename from applications/solvers/coupled/blockCoupledScalarTransportFoam/blockMatrixTools.C rename to src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.C index 601c8f640..e51fa657d 100644 --- a/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockMatrixTools.C +++ b/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.C @@ -51,6 +51,21 @@ void blockInsert } +template +void blockAdd +( + const direction dir, + const scalarField& x, + Field& blockX +) +{ + forAll (x, i) + { + blockX[i](dir) += x[i]; + } +} + + template void blockRetrieve ( @@ -131,6 +146,46 @@ void insertUpperLower return; } + if (m.symmetric() && blockM.symmetric()) + { + Info<< "Both m and blockM are symmetric: inserting only upper triangle" + << endl; + } + else + { + // Either scalar or block matrix is asymmetric: insert lower triangle + const scalarField& lower = m.lower(); + + if (blockM.lower().activeType() == blockCoeffBase::UNALLOCATED) + { + blockM.lower().asScalar() = lower; + } + else if + ( + blockM.lower().activeType() == blockCoeffBase::SCALAR + || blockM.lower().activeType() == blockCoeffBase::LINEAR + ) + { + typename CoeffField::linearTypeField& blockLower = + blockM.lower().asLinear(); + + forAll (lower, i) + { + blockLower[i](dir) = lower[i]; + } + } + else if (blockM.lower().activeType() == blockCoeffBase::SQUARE) + { + typename CoeffField::squareTypeField& blockLower = + blockM.lower().asSquare(); + + forAll (lower, i) + { + blockLower[i](dir, dir) = lower[i]; + } + } + } + if (m.hasUpper()) { const scalarField& upper = m.upper(); @@ -178,42 +233,61 @@ void insertUpperLower << abort(FatalError); } - if (m.symmetric() && blockM.symmetric()) + // Insert block interface fields + forAll(blockM.interfaces(), patchI) { - Info<< "Both m and blockM are symmetric: inserting only upper triangle" - << endl; - } - else - { - // Either scalar or block matrix is asymmetric: insert lower triangle - const scalarField& lower = m.lower(); - - if (blockM.lower().activeType() == blockCoeffBase::UNALLOCATED) + if (blockM.interfaces().set(patchI)) { - blockM.lower().asScalar() = lower; - } - else if - ( - blockM.lower().activeType() == blockCoeffBase::SCALAR - || blockM.lower().activeType() == blockCoeffBase::LINEAR - ) - { - typename CoeffField::linearTypeField& blockLower = - blockM.lower().asLinear(); + // Couple upper and lower + const scalarField& cUpper = m.boundaryCoeffs()[patchI]; + const scalarField& cLower = m.internalCoeffs()[patchI]; - forAll (lower, i) + if + ( + blockM.coupleUpper()[patchI].activeType() + == blockCoeffBase::UNALLOCATED + ) { - blockLower[i](dir) = lower[i]; + blockM.coupleUpper()[patchI].asScalar() = cUpper; + blockM.coupleLower()[patchI].asScalar() = cLower; } - } - else if (blockM.lower().activeType() == blockCoeffBase::SQUARE) - { - typename CoeffField::squareTypeField& blockLower = - blockM.lower().asSquare(); - - forAll (lower, i) + else if + ( + blockM.coupleUpper()[patchI].activeType() + == blockCoeffBase::SCALAR + || blockM.coupleUpper()[patchI].activeType() + == blockCoeffBase::LINEAR + ) { - blockLower[i](dir, dir) = lower[i]; + typename CoeffField::linearTypeField& blockUpper = + blockM.coupleUpper()[patchI].asLinear(); + + typename CoeffField::linearTypeField& blockLower = + blockM.coupleLower()[patchI].asLinear(); + + forAll (cUpper, i) + { + blockUpper[i](dir) = cUpper[i]; + blockLower[i](dir) = cLower[i]; + } + } + else if + ( + blockM.coupleUpper()[patchI].activeType() + == blockCoeffBase::SQUARE + ) + { + typename CoeffField::squareTypeField& blockUpper = + blockM.coupleUpper()[patchI].asSquare(); + + typename CoeffField::squareTypeField& blockLower = + blockM.coupleLower()[patchI].asSquare(); + + forAll (cUpper, i) + { + blockUpper[i](dir, dir) = cUpper[i]; + blockLower[i](dir, dir) = cLower[i]; + } } } } @@ -236,6 +310,142 @@ void insertEquation } +template +void insertCouplingDiagSource +( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& m, + BlockLduMatrix& blockM, + Field& blockB +) +{ + // Prepare the diagonal and source + + scalarField diag = m.diag(); + scalarField source = m.source(); + + // Add boundary source contribution + m.addBoundaryDiag(diag, 0); + m.addBoundarySource(source, false); + + + // Add off-diagonal block coefficients + typename CoeffField::squareTypeField& blockDiag = + blockM.diag().asSquare(); + + // Set off-diagonal coefficient + forAll (diag, i) + { + blockDiag[i](dirI, dirJ) = diag[i]; + } + + blockAdd(dirI, source, blockB); +} + + +template +void insertCouplingUpperLower +( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& m, + BlockLduMatrix& blockM +) +{ + if (m.diagonal()) + { + // Matrix for insertion is diagonal-only: nothing to do + return; + } + + if (m.symmetric() && blockM.symmetric()) + { + Info<< "Both m and blockM are symmetric: inserting only upper triangle" + << endl; + } + else + { + // Either scalar or block matrix is asymmetric: insert lower triangle + const scalarField& lower = m.lower(); + + typename CoeffField::squareTypeField& blockLower = + blockM.lower().asSquare(); + + forAll (lower, i) + { + blockLower[i](dirI, dirJ) = lower[i]; + } + } + + if (m.hasUpper()) + { + const scalarField& upper = m.upper(); + + typename CoeffField::squareTypeField& blockUpper = + blockM.upper().asSquare(); + + forAll (upper, i) + { + blockUpper[i](dirI, dirJ) = upper[i]; + } + } + else + { + FatalErrorIn + ( + "void insertCouplingUpperLower\n" + "(\n" + " const direction dirI,\n" + " const direction dirJ,\n" + " const fvScalarMatrix& m,\n" + " BlockLduMatrix& blockM\n" + ")" + ) << "Error in matrix insertion: problem with block structure" + << abort(FatalError); + } + + // Insert block interface fields + forAll(blockM.interfaces(), patchI) + { + if (blockM.interfaces().set(patchI)) + { + // Couple upper and lower + const scalarField& cUpper = m.boundaryCoeffs()[patchI]; + const scalarField& cLower = m.internalCoeffs()[patchI]; + + typename CoeffField::squareTypeField& blockUpper = + blockM.coupleUpper()[patchI].asSquare(); + + typename CoeffField::squareTypeField& blockLower = + blockM.coupleLower()[patchI].asSquare(); + + forAll (cUpper, i) + { + blockUpper[i](dirI, dirJ) = cUpper[i]; + blockLower[i](dirI, dirJ) = cLower[i]; + } + } + } +} + + +template +void insertCoupling +( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& m, + BlockLduMatrix& blockM, + Field& blockX, + Field& blockB +) +{ + insertCouplingDiagSource(dirI, dirJ, m, blockM, blockB); + insertCouplingUpperLower(dirI, dirJ, m, blockM); +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace blockMatrixTools diff --git a/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.H b/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.H new file mode 100644 index 000000000..9e275c781 --- /dev/null +++ b/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.H @@ -0,0 +1,179 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +InNamespace + Foam:: + +Description + Block matrix insertion and retrieval tools + +SourceFiles + blockMatrixTools.C + +\*---------------------------------------------------------------------------*/ + + +#ifndef blockMatrixTools_H +#define blockMatrixTools_H + +#include "blockLduMatrices.H" +#include "fvMatrices.H" +#include "surfaceFieldsFwd.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Namespace blockMatrixTools functions Declaration +\*---------------------------------------------------------------------------*/ + +namespace blockMatrixTools +{ + // Field operations + + //- Insert field into block field + template + void blockInsert + ( + const direction dir, + const scalarField& x, + Field& blockX + ); + + //- Add field into block field + template + void blockAdd + ( + const direction dir, + const scalarField& x, + Field& blockX + ); + + //- Retrieve field from block field + template + void blockRetrieve + ( + const direction dir, + scalarField& x, + const Field& blockX + ); + + + // Diagonal block operations + + //- Insert matrix diagonal and source into the block system + template + void insertDiagSource + ( + const direction dir, + const fvScalarMatrix& m, + BlockLduMatrix& blockM, + Field& blockB + ); + + // Insert upper and lower coefficients matrix into the block system + template + void insertUpperLower + ( + const direction dir, + const fvScalarMatrix& m, + BlockLduMatrix& blockM + ); + + // Insert matrix into the block system + template + void insertEquation + ( + const direction dir, + const fvScalarMatrix& m, + BlockLduMatrix& blockM, + Field& blockX, + Field& blockB + ); + + + // Coupling block operations + + //- Insert coupling matrix diagonal and source into the block system + template + void insertCouplingDiagSource + ( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& m, + BlockLduMatrix& blockM, + Field& blockB + ); + + // Insert coupling matrix into the block system + template + void insertCouplingUpperLower + ( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& m, + BlockLduMatrix& blockM + ); + + // Insert coupling matrix into the block system + template + void insertCoupling + ( + const direction dir, + const fvScalarMatrix& m, + BlockLduMatrix& blockM, + Field& blockX, + Field& blockB + ); + + // 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 + template + void updateCoupling + ( + BlockLduMatrix& blockM, + Field& x, + Field& b + ); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "blockMatrixTools.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/processorFvPatchVectorNFields.C b/src/VectorN/finiteVolume/fields/fvPatchFields/processorFvPatchVectorNFields.C index 18f4a9d92..9ef9a2246 100644 --- a/src/VectorN/finiteVolume/fields/fvPatchFields/processorFvPatchVectorNFields.C +++ b/src/VectorN/finiteVolume/fields/fvPatchFields/processorFvPatchVectorNFields.C @@ -34,65 +34,65 @@ namespace Foam // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -#define VectorNMatrixInterfaceFunc(Type) \ -template <> \ -void processorFvPatchField::initInterfaceMatrixUpdate \ -( \ - const Field& psiInternal, \ - Field&, \ - const BlockLduMatrix&, \ - const CoeffField&, \ - const Pstream::commsTypes commsType \ -) const \ -{ \ - procPatch_.compressedSend \ - ( \ - commsType, \ - this->patch().patchInternalField(psiInternal)() \ - ); \ -} \ - \ -template <> \ -void processorFvPatchField::updateInterfaceMatrix \ -( \ - const Field& psiInternal, \ - Field& result, \ - const BlockLduMatrix&, \ - const CoeffField& coeffs, \ - const Pstream::commsTypes commsType \ -) const \ -{ \ - Field pnf(this->size()); \ - \ - if (coeffs.activeType() == blockCoeffBase::SCALAR) \ - { \ - pnf = coeffs.asScalar() * \ - procPatch_.compressedReceive(commsType, this->size())(); \ - } \ - else if (coeffs.activeType() == blockCoeffBase::LINEAR) \ - { \ - pnf = cmptMultiply(coeffs.asLinear(), \ - procPatch_.compressedReceive(commsType, this->size())() \ - ); \ - } \ - else if (coeffs.activeType() == blockCoeffBase::SQUARE) \ - { \ - pnf = coeffs.asSquare() & \ - procPatch_.compressedReceive(commsType, this->size())(); \ - } \ - \ - const unallocLabelList& faceCells = this->patch().faceCells(); \ - \ - forAll(faceCells, facei) \ - { \ - result[faceCells[facei]] -= pnf[facei]; \ - } \ +#define VectorNMatrixInterfaceFunc(Type) \ +template <> \ +void processorFvPatchField::initInterfaceMatrixUpdate \ +( \ + const Field& psiInternal, \ + Field&, \ + const BlockLduMatrix&, \ + const CoeffField&, \ + const Pstream::commsTypes commsType \ +) const \ +{ \ + procPatch_.compressedSend \ + ( \ + commsType, \ + this->patch().patchInternalField(psiInternal)() \ + ); \ +} \ + \ +template <> \ +void processorFvPatchField::updateInterfaceMatrix \ +( \ + const Field& psiInternal, \ + Field& result, \ + const BlockLduMatrix&, \ + const CoeffField& coeffs, \ + const Pstream::commsTypes commsType \ +) const \ +{ \ + Field pnf(this->size()); \ + \ + if (coeffs.activeType() == blockCoeffBase::SCALAR) \ + { \ + pnf = coeffs.asScalar() * \ + procPatch_.compressedReceive(commsType, this->size())(); \ + } \ + else if (coeffs.activeType() == blockCoeffBase::LINEAR) \ + { \ + pnf = cmptMultiply(coeffs.asLinear(), \ + procPatch_.compressedReceive(commsType, this->size())() \ + ); \ + } \ + else if (coeffs.activeType() == blockCoeffBase::SQUARE) \ + { \ + pnf = coeffs.asSquare() & \ + procPatch_.compressedReceive(commsType, this->size())(); \ + } \ + \ + const unallocLabelList& faceCells = this->patch().faceCells(); \ + \ + forAll(faceCells, facei) \ + { \ + result[faceCells[facei]] -= pnf[facei]; \ + } \ } -#define doMakePatchTypeField(type, Type, args...) \ - VectorNMatrixInterfaceFunc(type) \ - \ +#define doMakePatchTypeField(type, Type, args...) \ + VectorNMatrixInterfaceFunc(type) \ + \ makePatchTypeField(fvPatch##Type##Field, processorFvPatch##Type##Field); forAllVectorNTypes(doMakePatchTypeField) diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFieldsFwd.H b/src/VectorN/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFieldsFwd.H index 557e8f51e..989733496 100644 --- a/src/VectorN/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFieldsFwd.H +++ b/src/VectorN/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFieldsFwd.H @@ -39,21 +39,25 @@ namespace Foam template class fvsPatchField; typedef fvsPatchField fvsPatchVector2Field; +typedef fvsPatchField fvsPatchVector3Field; typedef fvsPatchField fvsPatchVector4Field; typedef fvsPatchField fvsPatchVector6Field; typedef fvsPatchField fvsPatchVector8Field; typedef fvsPatchField fvsPatchTensor2Field; +typedef fvsPatchField fvsPatchTensor3Field; typedef fvsPatchField fvsPatchTensor4Field; typedef fvsPatchField fvsPatchTensor6Field; typedef fvsPatchField fvsPatchTensor8Field; typedef fvsPatchField fvsPatchDiagTensor2Field; +typedef fvsPatchField fvsPatchDiagTensor3Field; typedef fvsPatchField fvsPatchDiagTensor4Field; typedef fvsPatchField fvsPatchDiagTensor6Field; typedef fvsPatchField fvsPatchDiagTensor8Field; typedef fvsPatchField fvsPatchSphericalTensor2Field; +typedef fvsPatchField fvsPatchSphericalTensor3Field; typedef fvsPatchField fvsPatchSphericalTensor4Field; typedef fvsPatchField fvsPatchSphericalTensor6Field; typedef fvsPatchField fvsPatchSphericalTensor8Field; diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C index 96e95fd69..9d56b25ca 100644 --- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C +++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C @@ -90,8 +90,9 @@ void Foam::setRefCell " bool\n" ")", dict - ) << "Unable to set reference cell for field " << field.name() - << nl << " Reference point " << refPointName + ) << "Unable to set reference cell for field " + << field.name() << nl + << " Reference point " << refPointName << " " << refPointi << " found on " << sumHasRef << " domains (should be one)" << nl << exit(FatalIOError); diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H index 2953c98e4..4c3947b2d 100644 --- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H +++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H @@ -276,6 +276,13 @@ public: return source_; } + //- Access to fvBoundary scalar field containing + // pseudo-matrix coeffs for internal cells + const FieldField& internalCoeffs() const + { + return internalCoeffs_; + } + //- fvBoundary scalar field containing pseudo-matrix coeffs // for internal cells FieldField& internalCoeffs() @@ -283,6 +290,13 @@ public: return internalCoeffs_; } + //- Access to fvBoundary scalar field containing + // pseudo-matrix coeffs for boundary cells + const FieldField& boundaryCoeffs() const + { + return boundaryCoeffs_; + } + //- fvBoundary scalar field containing pseudo-matrix coeffs // for boundary cells FieldField& boundaryCoeffs() From c9523bb2f7dfeb7dd18f14238c82d7c536225bc1 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 24 Aug 2011 11:36:10 +0100 Subject: [PATCH 6/6] Block solver development --- .../blockCoupledSwirlTest/constant/polyMesh/boundary | 2 +- .../blockCoupledSwirlTest/system/fvSolution | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/constant/polyMesh/boundary b/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/constant/polyMesh/boundary index 664566406..aa9eaa611 100644 --- a/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/constant/polyMesh/boundary +++ b/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/constant/polyMesh/boundary @@ -1,6 +1,6 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | -| \\ / F ield | OpenFOAM Extend Project: Open Source CFD | +| \\ / F ield | OpenFOAM Extend Project: Open source CFD | | \\ / O peration | Version: 1.6-ext | | \\ / A nd | Web: www.extend-project.de | | \\/ M anipulation | | diff --git a/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/system/fvSolution b/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/system/fvSolution index a9243a0e6..db177fc76 100644 --- a/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/system/fvSolution +++ b/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/system/fvSolution @@ -16,7 +16,7 @@ FoamFile solvers { - blockVar + blockT { solver BiCGStab; preconditioner Cholesky; @@ -26,7 +26,7 @@ solvers minIter 1; maxIter 1000; - }; + } } SIMPLE