From 052a91e32aafae3ea3127a58db79074b707feb4b Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 13 Feb 2017 16:10:00 +0000 Subject: [PATCH 001/333] Bugfix: Named release --- ReleaseNotes | 1 + 1 file changed, 1 insertion(+) diff --git a/ReleaseNotes b/ReleaseNotes index 2293b6209..e0c886326 100644 --- a/ReleaseNotes +++ b/ReleaseNotes @@ -1,6 +1,7 @@ # -*- mode: org; -*- # #+TITLE: *Release notes for foam-extend-4.0* +#+TITLE: *Version 4.0 - Guimaraes* #+AUTHOR: foam-extend administrators: #+AUTHOR: Hrvoje Jasak #+AUTHOR: Håkan Nilsson From 38bc5f681f07f8bda9a8a5bd8c67d8188220deea Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 13 Feb 2017 16:10:43 +0000 Subject: [PATCH 002/333] Bugfix: Remove debug message --- src/finiteVolume/cfdTools/general/MRF/MRFZone.C | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C index 9633fff7a..ca008576d 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C @@ -84,7 +84,7 @@ void Foam::MRFZone::setMRFFaces() labelHashSet excludedPatches(excludedPatchLabels_); - Info<< "Excluded patches: " << excludedPatchLabels_ << endl; + forAll (patches, patchI) { const polyPatch& pp = patches[patchI]; @@ -246,12 +246,7 @@ Foam::vector Foam::MRFZone::Omega() const const scalar t = mesh_.time().value(); const scalar ramp = sin(2*pi/(4*rampTime_)*Foam::min(rampTime_, t)); - Info<< "MRF Ramp: " << ramp << endl; - return ramp*omega_.value()*axis_.value(); - - // return Foam::min(mesh_.time().value()/rampTime_, 1.0)* - // omega_.value()*axis_.value(); } } From 195f5b6094895bace0ec0a33c73dd84b476e5037 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 13 Feb 2017 16:10:53 +0000 Subject: [PATCH 003/333] Bugfix: Formatting --- .../solvers/incompressible/MRFSimpleFoam/createFields.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/solvers/incompressible/MRFSimpleFoam/createFields.H b/applications/solvers/incompressible/MRFSimpleFoam/createFields.H index e43f40963..c7021ecb7 100644 --- a/applications/solvers/incompressible/MRFSimpleFoam/createFields.H +++ b/applications/solvers/incompressible/MRFSimpleFoam/createFields.H @@ -12,7 +12,7 @@ mesh ); - Info << "Reading field U\n" << endl; + Info<< "Reading field U\n" << endl; volVectorField U ( IOobject From 1b05cb13dc318011ab1675caf3a3f3cdf843fee5 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 13 Feb 2017 16:11:19 +0000 Subject: [PATCH 004/333] Bugfix: Mac port: check for array size --- src/sampling/sampledSurface/sampledPatch/sampledPatch.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sampling/sampledSurface/sampledPatch/sampledPatch.C b/src/sampling/sampledSurface/sampledPatch/sampledPatch.C index c50315b2a..217d191ff 100644 --- a/src/sampling/sampledSurface/sampledPatch/sampledPatch.C +++ b/src/sampling/sampledSurface/sampledPatch/sampledPatch.C @@ -154,7 +154,7 @@ void Foam::sampledPatch::remapFaces ) { // recalculate the cells cut - if (&faceMap && faceMap.size()) + if (!faceMap.empty()) { MeshStorage::remapFaces(faceMap); } From 9fb69e5fd8e7a7d2fd3337a9eb3f57770e132cd2 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 7 Mar 2017 14:43:25 +0000 Subject: [PATCH 005/333] Bugfix: Typo in comment --- .../muSgsWallFunction/muSgsWallFunctionFvPatchScalarField.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/muSgsWallFunctions/muSgsWallFunction/muSgsWallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/muSgsWallFunctions/muSgsWallFunction/muSgsWallFunctionFvPatchScalarField.H index 90801690d..d08491a34 100644 --- a/src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/muSgsWallFunctions/muSgsWallFunction/muSgsWallFunctionFvPatchScalarField.H +++ b/src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/muSgsWallFunctions/muSgsWallFunction/muSgsWallFunctionFvPatchScalarField.H @@ -26,7 +26,7 @@ Class muSgsWallFunctionFvPatchScalarField Description - Spalart Allmaas wall function boundary condition for compressible flows + Spalart Allmaras wall function boundary condition for compressible flows SourceFiles muSgsWallFunctionFvPatchScalarField.C From cb59b93ba6354ec4cb27eb9d60fb0304def27f16 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 21 Mar 2017 19:51:00 +0000 Subject: [PATCH 006/333] Formatting --- .../solvers/compressible/steadyCompressibleMRFFoam/iEqn.H | 1 - 1 file changed, 1 deletion(-) diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H index 2fa37cfcc..910176930 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H @@ -19,7 +19,6 @@ ); iEqn.relax(); - iEqn.solve(); // Calculate enthalpy out of rothalpy From 92a208f6c295a1e8779437b5c7a9bad0f993bb4b Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:18:30 +0100 Subject: [PATCH 007/333] Clean-up of fam solver --- src/finiteArea/faMatrices/faMatrix/faMatrix.H | 8 -- .../faMatrices/faMatrix/faMatrixSolve.C | 6 -- .../faScalarMatrix/faScalarMatrix.C | 77 ------------------- .../faScalarMatrix/faScalarMatrix.H | 12 --- 4 files changed, 103 deletions(-) diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrix.H b/src/finiteArea/faMatrices/faMatrix/faMatrix.H index 2ebd26a62..27f6983fa 100644 --- a/src/finiteArea/faMatrices/faMatrix/faMatrix.H +++ b/src/finiteArea/faMatrices/faMatrix/faMatrix.H @@ -296,14 +296,6 @@ public: // alpha is read from controlDict void relax(); - //- Construct and return the solver - // Solver controls read from Istream - autoPtr solver(const dictionary&); - - //- Construct and return the solver - // Solver controls read from faSolution - autoPtr solver(); - //- Solve returning the solution statistics. // Solver controls read Istream lduSolverPerformance solve(const dictionary&); diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C b/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C index 948c04b51..5c0a3da25 100644 --- a/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C +++ b/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C @@ -153,12 +153,6 @@ lduSolverPerformance faMatrix::solve(const dictionary& solverControls) } -template -autoPtr::faSolver> faMatrix::solver() -{ - return solver(psi_.mesh().solutionDict().solverDict(psi_.name())); -} - template lduSolverPerformance faMatrix::faSolver::solve() { diff --git a/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C b/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C index 5b89e8903..184ddbcdb 100644 --- a/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C +++ b/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C @@ -57,83 +57,6 @@ void faMatrix::setComponentReference } -template<> -autoPtr::faSolver> faMatrix::solver -( - const dictionary& solverControls -) -{ - if (debug) - { - Info<< "faMatrix::solver(const dictionary&) : " - "solver for faMatrix" - << endl; - } - - scalarField saveDiag = diag(); - addBoundaryDiag(diag(), 0); - - // Make a copy of interfaces: no longer a reference - // HJ, 20/Nov/2007 - lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces(); - - autoPtr::faSolver> solverPtr - ( - new faMatrix::faSolver - ( - *this, - lduSolver::New - ( - psi_.name(), - *this, - boundaryCoeffs_, - internalCoeffs_, - interfaces, - solverControls - ) - ) - ); - - diag() = saveDiag; - - return solverPtr; -} - - -template<> -lduSolverPerformance faMatrix::faSolver::solve -( - const dictionary& solverControls -) -{ - scalarField saveDiag = faMat_.diag(); - faMat_.addBoundaryDiag(faMat_.diag(), 0); - - scalarField totalSource = faMat_.source(); - faMat_.addBoundarySource(totalSource, false); - - solver_->read(solverControls); - - // Cast into a non-const to solve. HJ, 6/May/2016 - GeometricField& psi = - const_cast&> - ( - faMat_.psi() - ); - - lduSolverPerformance solverPerf = - solver_->solve(psi.internalField(), totalSource); - - solverPerf.print(); - - faMat_.diag() = saveDiag; - - psi.correctBoundaryConditions(); - - return solverPerf; -} - - template<> lduSolverPerformance faMatrix::solve ( diff --git a/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.H b/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.H index d4294f1bf..e7e4c2593 100644 --- a/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.H +++ b/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.H @@ -54,18 +54,6 @@ void faMatrix::setComponentReference const scalar value ); -template<> -autoPtr::faSolver> faMatrix::solver -( - const dictionary& -); - -template<> -lduSolverPerformance faMatrix::faSolver::solve -( - const dictionary& -); - template<> lduSolverPerformance faMatrix::solve(const dictionary&); From a1309b87f500330b346b07f1d34dc51ef05147f2 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:18:44 +0100 Subject: [PATCH 008/333] Formatting --- .../basic/zeroGradient/zeroGradientFaPatchField.H | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/finiteArea/fields/faPatchFields/basic/zeroGradient/zeroGradientFaPatchField.H b/src/finiteArea/fields/faPatchFields/basic/zeroGradient/zeroGradientFaPatchField.H index 125f78a9b..34999eb24 100644 --- a/src/finiteArea/fields/faPatchFields/basic/zeroGradient/zeroGradientFaPatchField.H +++ b/src/finiteArea/fields/faPatchFields/basic/zeroGradient/zeroGradientFaPatchField.H @@ -42,7 +42,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class zeroGradientFaPatch Declaration + Class zeroGradientFaPatchField Declaration \*---------------------------------------------------------------------------*/ template @@ -119,6 +119,11 @@ public: } + //- Destructor + virtual ~zeroGradientFaPatchField() + {} + + // Member functions // Evaluation functions From f67e2e0c4fe79bd3e351b83f1efacf503ea22de0 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:24:25 +0100 Subject: [PATCH 009/333] Improved formulation of isentropic boundary conditions --- ...sentropicTotalPressureFvPatchScalarField.C | 57 +++++++++++++++++-- ...sentropicTotalPressureFvPatchScalarField.H | 3 + ...tropicTotalTemperatureFvPatchScalarField.C | 46 +++++---------- ...tropicTotalTemperatureFvPatchScalarField.H | 12 ---- 4 files changed, 70 insertions(+), 48 deletions(-) diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalPressure/isentropicTotalPressureFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalPressure/isentropicTotalPressureFvPatchScalarField.C index 4287ba7cc..41c501189 100644 --- a/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalPressure/isentropicTotalPressureFvPatchScalarField.C +++ b/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalPressure/isentropicTotalPressureFvPatchScalarField.C @@ -39,6 +39,7 @@ isentropicTotalPressureFvPatchScalarField ) : fixedValueFvPatchScalarField(p, iF), + phiName_("phi"), UName_("U"), TName_("T"), p0_(p.size(), 0.0) @@ -54,6 +55,7 @@ isentropicTotalPressureFvPatchScalarField ) : fixedValueFvPatchScalarField(p, iF), + phiName_(dict.lookupOrDefault("phi", "phi")), UName_(dict.lookupOrDefault("U", "U")), TName_(dict.lookupOrDefault("T", "T")), p0_("p0", dict, p.size()) @@ -82,6 +84,7 @@ isentropicTotalPressureFvPatchScalarField ) : fixedValueFvPatchScalarField(ptf, p, iF, mapper), + phiName_(ptf.phiName_), UName_(ptf.UName_), TName_(ptf.TName_), p0_(ptf.p0_, mapper) @@ -95,6 +98,7 @@ isentropicTotalPressureFvPatchScalarField ) : fixedValueFvPatchScalarField(ptf), + phiName_(ptf.phiName_), UName_(ptf.UName_), TName_(ptf.TName_), p0_(ptf.p0_) @@ -109,6 +113,7 @@ isentropicTotalPressureFvPatchScalarField ) : fixedValueFvPatchScalarField(ptf, iF), + phiName_(ptf.phiName_), UName_(ptf.UName_), TName_(ptf.TName_), p0_(ptf.p0_) @@ -149,9 +154,38 @@ void Foam::isentropicTotalPressureFvPatchScalarField::updateCoeffs() return; } + if + ( + !this->db().objectRegistry::found(phiName_) + || !this->db().objectRegistry::found(UName_) + || !this->db().objectRegistry::found(TName_) + ) + { + // Flux not available, do not update + InfoIn + ( + "void isentropicTotalPressureFvPatchScalarField::updateCoeffs()" + ) << "Flux, pressure or velocity field not found. " + << "Performing fixed value update" << endl; + + fixedValueFvPatchScalarField::updateCoeffs(); + + return; + } + + // Get flux + const surfaceScalarField& phi = + db().lookupObject(phiName_); + + const scalarField& phip = + patch().patchField(phi); + // Get velocity - const fvPatchVectorField& U = - patch().lookupPatchField(UName_); + // const fvPatchVectorField& U = + // patch().lookupPatchField(UName_); + + // Get pressure + const scalarField& p = *this; // Get temperature const fvPatchScalarField& T = @@ -166,15 +200,28 @@ void Foam::isentropicTotalPressureFvPatchScalarField::updateCoeffs() scalarField gamma = Cp/Cv; scalarField R = Cp - Cv; - scalarField Ma = -(patch().nf() & U)/sqrt(gamma*R*T); + scalarField rho = p/(R*T); + + // Velocity + // scalarField Un = -(patch().nf() & U); + + // Velocity from incoming flux + scalarField Un = max(-phip/(patch().magSf()*rho), scalar(0)); + + scalarField Ma = Un/sqrt(gamma*R*T); scalarField a = 1 + 0.5*(gamma - 1)*sqr(Ma); - operator==(p0_*pow(a, -gamma/(gamma - 1))); + operator== + ( + pos(Ma - 1)*p0_ + + neg(Ma - 1)*p0_*pow(a, -gamma/(gamma - 1)) + ); fixedValueFvPatchScalarField::updateCoeffs(); } + void Foam::isentropicTotalPressureFvPatchScalarField::updateCoeffs ( const vectorField& Up @@ -190,6 +237,8 @@ void Foam::isentropicTotalPressureFvPatchScalarField::write ) const { fixedValueFvPatchScalarField::write(os); + writeEntryIfDifferent(os, "phi", "phi", phiName_); + writeEntryIfDifferent(os, "U", "U", UName_); writeEntryIfDifferent(os, "T", "T", TName_); p0_.writeEntry("p0", os); } diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalPressure/isentropicTotalPressureFvPatchScalarField.H b/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalPressure/isentropicTotalPressureFvPatchScalarField.H index 78134c4d9..288e44ae6 100644 --- a/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalPressure/isentropicTotalPressureFvPatchScalarField.H +++ b/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalPressure/isentropicTotalPressureFvPatchScalarField.H @@ -52,6 +52,9 @@ class isentropicTotalPressureFvPatchScalarField { // Private data + //- Name of the flux field + word phiName_; + //- Name of the velocity field word UName_; diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalTemperature/isentropicTotalTemperatureFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalTemperature/isentropicTotalTemperatureFvPatchScalarField.C index 517a37df7..cdf7bfe8e 100644 --- a/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalTemperature/isentropicTotalTemperatureFvPatchScalarField.C +++ b/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalTemperature/isentropicTotalTemperatureFvPatchScalarField.C @@ -144,6 +144,20 @@ void Foam::isentropicTotalTemperatureFvPatchScalarField::updateCoeffs() return; } + if (!this->db().objectRegistry::found(pName_)) + { + // Flux not available, do not update + InfoIn + ( + "void isentropicTotalPressureFvPatchScalarField::updateCoeffs()" + ) << "Pressure field " << pName_ << " not found. " + << "Performing fixed value update" << endl; + + fixedValueFvPatchScalarField::updateCoeffs(); + + return; + } + // Get pressure and temperature const scalarField& T = *this; @@ -174,38 +188,6 @@ void Foam::isentropicTotalTemperatureFvPatchScalarField::updateCoeffs } -Foam::tmp -Foam::isentropicTotalTemperatureFvPatchScalarField::snGrad() const -{ - return tmp - ( - new scalarField(this->size(), 0.0) - ); -} - - -Foam::tmp -Foam::isentropicTotalTemperatureFvPatchScalarField:: -gradientInternalCoeffs() const -{ - return tmp - ( - new scalarField(this->size(), 0.0) - ); -} - - -Foam::tmp -Foam::isentropicTotalTemperatureFvPatchScalarField:: -gradientBoundaryCoeffs() const -{ - return tmp - ( - new scalarField(this->size(), 0.0) - ); -} - - void Foam::isentropicTotalTemperatureFvPatchScalarField::write ( Ostream& os diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalTemperature/isentropicTotalTemperatureFvPatchScalarField.H b/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalTemperature/isentropicTotalTemperatureFvPatchScalarField.H index e622d8187..15ca93016 100644 --- a/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalTemperature/isentropicTotalTemperatureFvPatchScalarField.H +++ b/src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalTemperature/isentropicTotalTemperatureFvPatchScalarField.H @@ -170,18 +170,6 @@ public: //- Update the coefficients associated with the patch field virtual void updateCoeffs(); - //- Return patch-normal gradient: set to zero - virtual tmp snGrad() const; - - //- Return the matrix diagonal coefficients corresponding to the - // evaluation of the gradient of this patchField - // Removed gradient contribution - virtual tmp gradientInternalCoeffs() const; - - //- Return the matrix source coefficients corresponding to the - // evaluation of the gradient of this patchField - virtual tmp gradientBoundaryCoeffs() const; - //- Write virtual void write(Ostream&) const; From 7a89d0627897afc5e85ce11b8f86f557b6a2da52 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:26:05 +0100 Subject: [PATCH 010/333] Moved crMatrix into foam library --- src/{lduSolvers => foam/matrices}/crMatrix/crAddressing.C | 0 src/{lduSolvers => foam/matrices}/crMatrix/crAddressing.H | 0 src/{lduSolvers => foam/matrices}/crMatrix/crMatrix.C | 0 src/{lduSolvers => foam/matrices}/crMatrix/crMatrix.H | 0 4 files changed, 0 insertions(+), 0 deletions(-) rename src/{lduSolvers => foam/matrices}/crMatrix/crAddressing.C (100%) rename src/{lduSolvers => foam/matrices}/crMatrix/crAddressing.H (100%) rename src/{lduSolvers => foam/matrices}/crMatrix/crMatrix.C (100%) rename src/{lduSolvers => foam/matrices}/crMatrix/crMatrix.H (100%) diff --git a/src/lduSolvers/crMatrix/crAddressing.C b/src/foam/matrices/crMatrix/crAddressing.C similarity index 100% rename from src/lduSolvers/crMatrix/crAddressing.C rename to src/foam/matrices/crMatrix/crAddressing.C diff --git a/src/lduSolvers/crMatrix/crAddressing.H b/src/foam/matrices/crMatrix/crAddressing.H similarity index 100% rename from src/lduSolvers/crMatrix/crAddressing.H rename to src/foam/matrices/crMatrix/crAddressing.H diff --git a/src/lduSolvers/crMatrix/crMatrix.C b/src/foam/matrices/crMatrix/crMatrix.C similarity index 100% rename from src/lduSolvers/crMatrix/crMatrix.C rename to src/foam/matrices/crMatrix/crMatrix.C diff --git a/src/lduSolvers/crMatrix/crMatrix.H b/src/foam/matrices/crMatrix/crMatrix.H similarity index 100% rename from src/lduSolvers/crMatrix/crMatrix.H rename to src/foam/matrices/crMatrix/crMatrix.H From f91bdbd2e8b57e7849d330636a0d044a6453ea7d Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:28:06 +0100 Subject: [PATCH 011/333] Updates in AMG and block AMG coarsening --- .../BlockMatrixAgglomeration.C | 76 +- .../scalarBlockMatrixAgglomeration.H | 5 +- .../tensorBlockMatrixAgglomeration.H | 3 +- .../BlockMatrixClustering.C | 1415 +++++++++++++++++ .../BlockMatrixClustering.H | 218 +++ .../blockMatrixClusterings.C | 43 + .../blockMatrixClusterings.H | 65 + .../scalarBlockMatrixClustering.H | 86 + .../tensorBlockMatrixClustering.H | 83 + .../BlockAMG/coarseBlockAMGLevel.H | 2 +- .../BlockLduSolver/BlockLduSolver.H | 7 +- .../BlockLduSolvers/blockVectorNSolvers.C | 4 + .../BlockCoeffComponentNorm.C | 10 +- .../BlockCoeffComponentNorm.H | 13 +- .../BlockCoeffMaxNorm/BlockCoeffMaxNorm.C | 19 +- .../BlockCoeffMaxNorm/BlockCoeffMaxNorm.H | 13 +- .../scalarBlockCoeffMaxNorm.H | 6 +- .../tensorBlockCoeffMaxNorm.H | 12 +- .../BlockCoeffNorm/BlockCoeffNorm.H | 26 +- .../BlockCoeffTwoNorm/BlockCoeffTwoNorm.C | 15 +- .../BlockCoeffTwoNorm/BlockCoeffTwoNorm.H | 12 +- .../tensorBlockCoeffTwoNorm.H | 10 +- src/lduSolvers/amg/amgPolicy/amgPolicy.H | 7 +- .../amg/amgPolicy/clusterAmgPolicy.C | 1021 ++++++++++++ .../amg/amgPolicy/clusterAmgPolicy.H | 165 ++ src/lduSolvers/amg/amgPolicy/pamgPolicy.C | 20 +- src/lduSolvers/amg/amgPolicy/pamgPolicy.H | 4 +- 27 files changed, 3246 insertions(+), 114 deletions(-) create mode 100644 src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/blockMatrixClusterings.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/blockMatrixClusterings.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/scalarBlockMatrixClustering.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/tensorBlockMatrixClustering.H create mode 100644 src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.C create mode 100644 src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.H diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.C index 25c22b8fa..bf0c934aa 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.C +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.C @@ -36,7 +36,6 @@ Author #include "boolList.H" #include "tolerancesSwitch.H" #include "coeffFields.H" -#include "addToRunTimeSelectionTable.H" #include "BlockAMGInterfaceField.H" #include "coarseBlockAMGLevel.H" @@ -153,21 +152,32 @@ void Foam::BlockMatrixAgglomeration::calcAgglomeration() scalarField magDiag(diag.size()); scalarField magOffDiag(upperAddr.size()); - normPtr_->coeffMag(diag, magDiag); + // Calculate norm of block coefficient. Note: the norm + // may be signed, ie. it will take the sign of the coefficient + normPtr_->normalize(magDiag, diag); + + // Take magnitude of the norm + mag(magDiag, magDiag); if (matrix_.asymmetric()) { scalarField magUpper(upperAddr.size()); scalarField magLower(upperAddr.size()); - normPtr_->coeffMag(matrix_.upper(), magUpper); - normPtr_->coeffMag(matrix_.lower(), magLower); + // Take max of norm magnitude in upper and lower triangle + normPtr_->normalize(magUpper, matrix_.upper()); + mag(magUpper, magUpper); - magOffDiag = Foam::max(magUpper, magLower); + normPtr_->normalize(magLower, matrix_.lower()); + mag(magLower, magLower); + + Foam::max(magOffDiag, magUpper, magLower); } else if (matrix_.symmetric()) { - normPtr_->coeffMag(matrix_.upper(), magOffDiag); + // Take max of norm magnitude + normPtr_->normalize(magOffDiag, matrix_.upper()); + mag(magOffDiag, magOffDiag); } else { @@ -223,13 +233,9 @@ void Foam::BlockMatrixAgglomeration::calcAgglomeration() if (nSolo_ > 0) { // Found solo equations - nCoarseEqns_++; + sizeOfGroups[nCoarseEqns_] = nSolo_; - if (BlockLduMatrix::debug >= 3) - { - Pout<< "Found " << nSolo_ << " weakly connected equations." - << endl; - } + nCoarseEqns_++; } } @@ -340,11 +346,14 @@ void Foam::BlockMatrixAgglomeration::calcAgglomeration() } } + // Resize group size count + sizeOfGroups.setSize(nCoarseEqns_); + // The decision on parallel agglomeration needs to be made for the // whole gang of processes; otherwise I may end up with a different // number of agglomeration levels on different processors. - // If the number of coarse equations is les than minimum and + // If the number of coarse equations is less than minimum and // if the matrix has reduced in size by at least 1/3, coarsen if ( @@ -357,9 +366,32 @@ void Foam::BlockMatrixAgglomeration::calcAgglomeration() reduce(coarsen_, andOp()); - if (BlockLduMatrix::debug >= 3) + // if (BlockLduMatrix::debug >= 3) { - Pout << "Coarse level size: " << nCoarseEqns_; + // Count singleton clusters + label nSingleClusters = 0; + + // Adjust start based on solo cells cluster status + label start = 0; + + if (nSolo_ > 0) + { + start = 1; + } + + for (label i = start; i < sizeOfGroups.size(); i++) + { + if (sizeOfGroups[i] == 1) + { + nSingleClusters++; + } + } + + Pout<< "Coarse level size: " << nCoarseEqns_ + << " nSolo = " << nSolo_ + << " cluster (" << min(sizeOfGroups) << " " + << max(sizeOfGroups) + << "). N singleton clusters = " << nSingleClusters; if (coarsen_) { @@ -369,20 +401,6 @@ void Foam::BlockMatrixAgglomeration::calcAgglomeration() { Pout << ". Rejected" << endl; } - - // Count cluster size - labelList clusterSize(nCoarseEqns_, 0); - - forAll (agglomIndex_, eqnI) - { - clusterSize[agglomIndex_[eqnI]]++; - } - - label minClusterSize = gMin(clusterSize); - label maxClusterSize = gMax(clusterSize); - - Info<< "Cluster size: min = " << minClusterSize - << " max = " << maxClusterSize << endl; } } diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/scalarBlockMatrixAgglomeration.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/scalarBlockMatrixAgglomeration.H index 42226f6e5..05477172e 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/scalarBlockMatrixAgglomeration.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/scalarBlockMatrixAgglomeration.H @@ -52,11 +52,10 @@ namespace Foam /*---------------------------------------------------------------------------*\ - Class BlockMatrixAgglomeration Declaration + Class BlockMatrixAgglomeration Declaration \*---------------------------------------------------------------------------*/ -//- Restrict matrix - +// Disable restrict matrix: no square type template<> inline autoPtr > BlockMatrixAgglomeration::restrictMatrix() const diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/tensorBlockMatrixAgglomeration.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/tensorBlockMatrixAgglomeration.H index 648465be3..6a19fdeae 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/tensorBlockMatrixAgglomeration.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/tensorBlockMatrixAgglomeration.H @@ -49,8 +49,7 @@ namespace Foam Class BlockMatrixAgglomeration Declaration \*---------------------------------------------------------------------------*/ -//- Restrict matrix - +// Disable restrict matrix: no square type template<> inline autoPtr > BlockMatrixAgglomeration::restrictMatrix diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C new file mode 100644 index 000000000..11f027efc --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C @@ -0,0 +1,1415 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend 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 foam-extend. If not, see . + +Class + BlockMatrixClustering + +Description + Block matrix AMG coarsening by Jasak clustering algorithm + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#include "BlockMatrixClustering.H" +#include "boolList.H" +#include "tolerancesSwitch.H" +#include "coeffFields.H" +#include "BlockAMGInterfaceField.H" +#include "coarseBlockAMGLevel.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +template +const Foam::debug::tolerancesSwitch +Foam::BlockMatrixClustering::weightFactor_ +( + "aamgWeightFactor", + 0.65 +); + + +template +const Foam::debug::tolerancesSwitch +Foam::BlockMatrixClustering::diagFactor_ +( + "aamgDiagFactor", + 1e-8 +); + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +void Foam::BlockMatrixClustering::calcClustering() +{ + if (matrix_.diagonal()) + { + // Diag only matrix. Reset and return + agglomIndex_ = 0; + nCoarseEqns_ = 1; + + return; + } + + // Algorithm: + // 1) Calculate appropriate connection strength for symmetric and asymmetric + // matrix + // 2) Collect solo (disconnected) cells and remove them from agglomeration + // by placing them into cluster zero (solo group) + // 3) Loop through all equations and for each equation find the best fit + // neighbour. If all neighbours are grouped, add equation + // to best group + + // Initialise child array + agglomIndex_ = -1; + + const label nRows = matrix_.lduAddr().size(); + + // Get matrix addressing + const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); + const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); + const unallocLabelList& losortAddr = matrix_.lduAddr().losortAddr(); + + const unallocLabelList& ownerStartAddr = + matrix_.lduAddr().ownerStartAddr(); + const unallocLabelList& losortStartAddr = + matrix_.lduAddr().losortStartAddr(); + + + // Calculate clustering + + // Get matrix coefficients and norms Note: the norm + // may be signed, ie. it will take the sign of the coefficient + + const CoeffField& diag = matrix_.diag(); + scalarField normDiag(diag.size()); + normPtr_->normalize(normDiag, diag); + + scalarField normUpper(upperAddr.size(), 0); + scalarField normLower(upperAddr.size(), 0); + + // Note: + // Matrix properties are no longer assumed, eg. the diag and off-diag + // sign is checked + // HJ, 29/Mar/2017 + + // Note: negative connections are eliminated in max(...) below + // HJ, 30/Mar/2017 + + if (matrix_.thereIsUpper()) + { + normPtr_->normalize(normUpper, matrix_.upper()); + + // Owner: upper triangle + forAll (lowerAddr, coeffI) + { + // Sign of strong positive upper is opposite of the sign of + // its diagonal coefficient + normUpper[coeffI] = Foam::max + ( + -1*sign(normDiag[lowerAddr[coeffI]])*normUpper[coeffI], + 0 + ); + } + } + + if (matrix_.thereIsLower()) + { + normPtr_->normalize(normLower, matrix_.lower()); + + // Neighbour: lower triangle + forAll (lowerAddr, coeffI) + { + // Sign of strong positive upper is opposite of the sign of + // its diagonal coefficient + normLower[coeffI] = Foam::max + ( + -1*sign(normDiag[upperAddr[coeffI]])*normLower[coeffI], + 0 + ); + } + } + else + { + normLower = normUpper; + } + + // Take magnitude of the diagonal norm after the normalisation + // of upper and lower is complete + // Note: upper and lower are already normalised + normDiag = mag(normDiag); + + labelList sizeOfGroups(nRows, 0); + + nCoarseEqns_ = 0; + + // Gather disconnected and weakly connected equations into cluster zero + // Weak connection is assumed to be the one where the off-diagonal + // coefficient is smaller than diagFactor_()*diag + { + // Algorithm + // Mark all cells to belong to zero cluster + // Go through all upper and lower coefficients and for the ones + // larger than threshold mark the equations out of cluster zero + + scalarField scaledNormDiag = diagFactor_()*normDiag; + + boolList zeroCluster(normDiag.size(), true); + + if (matrix_.symmetric()) + { + // Owner: upper triangle + forAll (lowerAddr, coeffI) + { + if (normUpper[coeffI] > scaledNormDiag[lowerAddr[coeffI]]) + { + zeroCluster[lowerAddr[coeffI]] = false; + } + } + + // Neighbour: lower triangle with symm coefficients + forAll (upperAddr, coeffI) + { + if (normUpper[coeffI] > scaledNormDiag[upperAddr[coeffI]]) + { + zeroCluster[upperAddr[coeffI]] = false; + } + } + } + else if (matrix_.asymmetric()) + { + // Owner: upper triangle + forAll (lowerAddr, coeffI) + { + if (normUpper[coeffI] > scaledNormDiag[lowerAddr[coeffI]]) + { + zeroCluster[lowerAddr[coeffI]] = false; + } + } + + // Neighbour: lower triangle with lower coeffs + forAll (upperAddr, coeffI) + { + if (normLower[coeffI] > scaledNormDiag[upperAddr[coeffI]]) + { + zeroCluster[upperAddr[coeffI]] = false; + } + } + } + + // Collect solo equations + forAll (zeroCluster, eqnI) + { + if (zeroCluster[eqnI]) + { + // Found solo equation + nSolo_++; + + agglomIndex_[eqnI] = nCoarseEqns_; + } + } + + if (nSolo_ > 0) + { + // Found solo equations + sizeOfGroups[nCoarseEqns_] = nSolo_; + + nCoarseEqns_++; + } + } + + + // Loop through cells + // - if the cell is not already grouped, open a new group (seed) + // - find stroungest grouped and ungrouped connection: + // - grouped connection is towards an already grouped cell + // - ungrouped connection is towards an ungrouped cell + // - decide if a grouped or ungrouped connection is better + + label curEqn, indexUngrouped, indexGrouped, colI, groupPassI, + nextUngrouped, nextGrouped; + + scalar curWeight, weightUngrouped, weightGrouped; + + for (label rowI = 0; rowI < nRows; rowI++) + { + if (agglomIndex_[rowI] == -1) + { + // Found new ungrouped equation + curEqn = rowI; + + // Reset grouped and upgrouped index + indexUngrouped = -1; + indexGrouped = -1; + + // Make next group (coarse equation) and search for neighbours + agglomIndex_[curEqn] = nCoarseEqns_; + + // Work on the group until the min group size is satisfied + // As each new element of the group is found, group search starts + // from this element until group cluster size is satisfied + for + ( + groupPassI = 1; + groupPassI < minClusterSize_; + groupPassI++ + ) + { + weightUngrouped = 0; + weightGrouped = 0; + + indexUngrouped = -1; + indexGrouped = -1; + + nextUngrouped = -1; + nextGrouped = -1; + + // Visit all neighbour equations + + // Visit upper triangle coefficients from the equation + for + ( + label rowCoeffI = ownerStartAddr[curEqn]; + rowCoeffI < ownerStartAddr[curEqn + 1]; + rowCoeffI++ + ) + { + // Get column index, upper triangle + colI = upperAddr[rowCoeffI]; + + curWeight = normUpper[rowCoeffI]/normDiag[curEqn]; + + if (agglomIndex_[colI] == -1) + { + // Found ungrouped neighbour + if + ( + indexUngrouped == -1 + || curWeight > weightUngrouped + ) + { + // Found or updated ungrouped neighbour + indexUngrouped = rowCoeffI; + weightUngrouped = curWeight; + + // Record possible next equation for search + nextUngrouped = colI; + } + } + else if (agglomIndex_[curEqn] != agglomIndex_[colI]) + { + // Check for neighbour in solo group + if (nSolo_ == 0 || agglomIndex_[colI] != 0) + { + // Found neighbour belonging to other group + if (indexGrouped == -1 || curWeight > weightGrouped) + { + // Found or updated grouped neighbour + indexGrouped = rowCoeffI; + weightGrouped = curWeight; + + // Record possible next equation for search + nextGrouped = colI; + } + } + } + } + + // Visit lower triangle coefficients from the equation + for + ( + label rowCoeffI = losortStartAddr[curEqn]; + rowCoeffI < losortStartAddr[curEqn + 1]; + rowCoeffI++ + ) + { + // Get column index, lower triangle + colI = lowerAddr[losortAddr[rowCoeffI]]; + + curWeight = normLower[losortAddr[rowCoeffI]]/ + normDiag[curEqn]; + + if (agglomIndex_[colI] == -1) + { + // Found first or better ungrouped neighbour + if + ( + indexUngrouped == -1 + || curWeight > weightUngrouped + ) + { + // Found or updated ungrouped neighbour + indexUngrouped = rowCoeffI; + weightUngrouped = curWeight; + + // Record possible next equation for search + nextUngrouped = colI; + } + } + else if (agglomIndex_[curEqn] != agglomIndex_[colI]) + { + // Check for neighbour in solo group + if (nSolo_ == 0 || agglomIndex_[colI] != 0) + { + // Found first of better neighbour belonging to + // other group + if + ( + indexGrouped == -1 + || curWeight > weightGrouped + ) + { + // Found or updated grouped neighbour + indexGrouped = rowCoeffI; + weightGrouped = curWeight; + + // Record possible next equation for search + nextGrouped = colI; + } + } + } + } + + // Decide what to do with current equation + + // If ungrouped candidate exists and it is stronger + // than best weighted grouped candidate, use it + if + ( + indexUngrouped != -1 + && ( + indexGrouped == -1 + || weightUngrouped >= weightFactor_()*weightGrouped + ) + ) + { + // Found new element of group. Add it and use as + // start of next search + + agglomIndex_[nextUngrouped] = agglomIndex_[curEqn]; + sizeOfGroups[agglomIndex_[curEqn]]++; + + // Search from nextEqn + curEqn = nextUngrouped; + } + else + { + // Group full or cannot be extended with new + // ungrouped candidates + break; + } + } + + + // Finished group passes + if + ( + groupPassI > 1 + || indexGrouped == -1 + || sizeOfGroups[agglomIndex_[nextGrouped]] > maxClusterSize_ + ) + { + // There is no group to put this equation into + sizeOfGroups[agglomIndex_[rowI]]++; + nCoarseEqns_++; + } + else + { + // Dump current cell into the best group available + agglomIndex_[rowI] = agglomIndex_[nextGrouped]; + sizeOfGroups[agglomIndex_[nextGrouped]]++; + } + } + } + + // Resize group size count + sizeOfGroups.setSize(nCoarseEqns_); + + // The decision on parallel agglomeration needs to be made for the + // whole gang of processes; otherwise I may end up with a different + // number of agglomeration levels on different processors. + + // If the number of coarse equations is less than minimum and + // if the matrix has reduced in size by at least 1/3, coarsen + if + ( + nCoarseEqns_ > BlockMatrixCoarsening::minCoarseEqns() + && 3*nCoarseEqns_ <= 2*nRows + ) + { + coarsen_ = true; + } + + reduce(coarsen_, andOp()); + + // if (BlockLduMatrix::debug >= 3) + { + // Count singleton clusters + label nSingleClusters = 0; + + // Adjust start based on solo cells cluster status + label start = 0; + + if (nSolo_ > 0) + { + start = 1; + } + + for (label i = start; i < sizeOfGroups.size(); i++) + { + if (sizeOfGroups[i] == 1) + { + nSingleClusters++; + } + } + + Pout<< "Coarse level size: " << nCoarseEqns_ + << " nSolo = " << nSolo_ + << " cluster (" << min(sizeOfGroups) << " " + << max(sizeOfGroups) + << "). N singleton clusters = " << nSingleClusters; + + if (coarsen_) + { + Pout << ". Accepted" << endl; + } + else + { + Pout << ". Rejected" << endl; + } + } +} + + +template +void Foam::BlockMatrixClustering::restrictDiag +( + const CoeffField& Coeff, + CoeffField& coarseCoeff +) const +{ + typedef CoeffField TypeCoeffField; + + if + ( + Coeff.activeType() == blockCoeffBase::SQUARE + && coarseCoeff.activeType() == blockCoeffBase::SQUARE + ) + { + typedef typename TypeCoeffField::squareType squareType; + typedef typename TypeCoeffField::squareTypeField squareTypeField; + + squareTypeField& activeCoarseCoeff = coarseCoeff.asSquare(); + const squareTypeField& activeCoeff = Coeff.asSquare(); + + forAll (coarseCoeff, i) + { + activeCoarseCoeff[i] = pTraits::zero; + } + + forAll (Coeff, i) + { + activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; + } + } + else if + ( + Coeff.activeType() == blockCoeffBase::LINEAR + && coarseCoeff.activeType() == blockCoeffBase::LINEAR + ) + { + typedef typename TypeCoeffField::linearType linearType; + typedef typename TypeCoeffField::linearTypeField linearTypeField; + + linearTypeField& activeCoarseCoeff = coarseCoeff.asLinear(); + const linearTypeField& activeCoeff = Coeff.asLinear(); + + forAll (coarseCoeff, i) + { + activeCoarseCoeff[i] = pTraits::zero; + } + + forAll (Coeff, i) + { + activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; + } + } + else if + ( + Coeff.activeType() == blockCoeffBase::SCALAR + && coarseCoeff.activeType() == blockCoeffBase::SCALAR + ) + { + typedef typename TypeCoeffField::scalarType scalarType; + typedef typename TypeCoeffField::scalarTypeField scalarTypeField; + + scalarTypeField& activeCoarseCoeff = coarseCoeff.asScalar(); + const scalarTypeField& activeCoeff = Coeff.asScalar(); + + forAll (coarseCoeff, i) + { + activeCoarseCoeff[i] = pTraits::zero; + } + + forAll (Coeff, i) + { + activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; + } + } + else + { + FatalErrorIn + ( + "void BlockMatrixClustering::restrictDiag() const" + ) << "Problem in coeff type morphing" + << abort(FatalError); + } +} + + +template +template +void Foam::BlockMatrixClustering::agglomerateCoeffs +( + const labelList& coeffRestrictAddr, + Field& activeCoarseDiag, + Field& activeCoarseUpper, + const Field& activeFineUpper, + const Field& activeFineUpperTranspose +) const +{ + // Does the matrix have solo equations + bool soloEqns = nSolo_ > 0; + + // Get addressing + const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); + + forAll(coeffRestrictAddr, fineCoeffI) + { + label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffI]]; + label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffI]]; + + // If the coefficient touches block zero and + // solo equations are present, skip it + if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) + { + continue; + } + + label cCoeff = coeffRestrictAddr[fineCoeffI]; + + if (cCoeff >= 0) + { + activeCoarseUpper[cCoeff] += activeFineUpper[fineCoeffI]; + } + else + { + // Add the fine face coefficient into the diagonal + // Note: upper and lower coeffs are transpose of + // each other. HJ, 28/May/2014 + activeCoarseDiag[-1 - cCoeff] += + activeFineUpper[fineCoeffI] + + activeFineUpperTranspose[fineCoeffI]; + } + } +} + + +template +template +void Foam::BlockMatrixClustering::agglomerateCoeffs +( + const labelList& coeffRestrictAddr, + Field& activeCoarseDiag, + Field& activeCoarseUpper, + const Field& activeFineUpper, + Field& activeCoarseLower, + const Field& activeFineLower +) const +{ + // Does the matrix have solo equations + bool soloEqns = nSolo_ > 0; + + // Get addressing + const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); + + forAll(coeffRestrictAddr, fineCoeffI) + { + label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffI]]; + label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffI]]; + + // If the coefficient touches block zero and + // solo equations are present, skip it + if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) + { + continue; + } + + label cCoeff = coeffRestrictAddr[fineCoeffI]; + + if (cCoeff >= 0) + { + activeCoarseUpper[cCoeff] += activeFineUpper[fineCoeffI]; + activeCoarseLower[cCoeff] += activeFineLower[fineCoeffI]; + } + else + { + // Add the fine face coefficients into the diagonal. + activeCoarseDiag[-1 - cCoeff] += + activeFineUpper[fineCoeffI] + + activeFineLower[fineCoeffI]; + } + } +} + + +template +void Foam::BlockMatrixClustering::restrictDiagDecoupled +( + const CoeffField& Coeff, + CoeffField& coarseCoeff +) const +{ + typedef CoeffField TypeCoeffField; + + if + ( + Coeff.activeType() == blockCoeffBase::LINEAR + && coarseCoeff.activeType() == blockCoeffBase::LINEAR + ) + { + typedef typename TypeCoeffField::linearType linearType; + typedef typename TypeCoeffField::linearTypeField linearTypeField; + + linearTypeField& activeCoarseCoeff = coarseCoeff.asLinear(); + const linearTypeField& activeCoeff = Coeff.asLinear(); + + forAll (coarseCoeff, i) + { + activeCoarseCoeff[i] = pTraits::zero; + } + + forAll (Coeff, i) + { + activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; + } + } + else if + ( + Coeff.activeType() == blockCoeffBase::SCALAR + && coarseCoeff.activeType() == blockCoeffBase::SCALAR + ) + { + typedef typename TypeCoeffField::scalarType scalarType; + typedef typename TypeCoeffField::scalarTypeField scalarTypeField; + + scalarTypeField& activeCoarseCoeff = coarseCoeff.asScalar(); + const scalarTypeField& activeCoeff = Coeff.asScalar(); + + forAll (coarseCoeff, i) + { + activeCoarseCoeff[i] = pTraits::zero; + } + + forAll (Coeff, i) + { + activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; + } + } + else + { + FatalErrorIn + ( + "void BlockMatrixClustering::restrictDiagDecoupled()" + " const" + ) << "Problem in coeff type morphing" + << abort(FatalError); + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::BlockMatrixClustering::BlockMatrixClustering +( + const BlockLduMatrix& matrix, + const dictionary& dict, + const label groupSize, + const label minCoarseEqns +) +: + BlockMatrixCoarsening(matrix, dict, groupSize, minCoarseEqns), + matrix_(matrix), + minClusterSize_(readLabel(dict.lookup("minClusterSize"))), + maxClusterSize_(readLabel(dict.lookup("maxClusterSize"))), + normPtr_(BlockCoeffNorm::New(dict)), + agglomIndex_(matrix_.lduAddr().size()), + groupSize_(groupSize), + nSolo_(0), + nCoarseEqns_(0), + coarsen_(false), + lTime_() +{ + calcClustering(); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::BlockMatrixClustering::~BlockMatrixClustering() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +Foam::autoPtr > +Foam::BlockMatrixClustering::restrictMatrix() const +{ + if (!coarsen_) + { + FatalErrorIn + ( + "autoPtr > " + "BlockMatrixClustering::restrictMatrix() const" + ) << "Requesting coarse matrix when it cannot be created" + << abort(FatalError); + } + + // Construct the coarse matrix and ldu addressing for the next level + // Algorithm: + // 1) Loop through all fine coeffs. If the agglom labels on two sides are + // different, this creates a coarse coeff. Define owner and neighbour + // for this coeff based on cluster IDs. + // 2) Check if the coeff has been seen before. If yes, add the coefficient + // to the appropriate field (stored with the equation). If no, create + // a new coeff with neighbour ID and add the coefficient + // 3) Once all the coeffs have been created, loop through all clusters and + // insert the coeffs in the upper order. At the same time, collect the + // owner and neighbour addressing. + // 4) Agglomerate the diagonal by summing up the fine diagonal + + // Get addressing + const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); + + const label nFineCoeffs = upperAddr.size(); + +# ifdef FULLDEBUG + if (agglomIndex_.size() != matrix_.lduAddr().size()) + { + FatalErrorIn + ( + "autoPtr >" + "BlockMatrixClustering::restrictMatrix() const" + ) << "agglomIndex array does not correspond to fine level. " << endl + << " Size: " << agglomIndex_.size() + << " number of equations: " << matrix_.lduAddr().size() + << abort(FatalError); + } +# endif + + // Does the matrix have solo equations + bool soloEqns = nSolo_ > 0; + + // Storage for block neighbours and coefficients + + // Guess initial maximum number of neighbours in block + label maxNnbrs = 10; + + // Number of neighbours per block + labelList blockNnbrs(nCoarseEqns_, 0); + + // Setup initial packed storage for neighbours and coefficients + labelList blockNbrsData(maxNnbrs*nCoarseEqns_); + + // Create face-restriction addressing + labelList coeffRestrictAddr(nFineCoeffs); + + // Initial neighbour array (not in upper-triangle order) + labelList initCoarseNeighb(nFineCoeffs); + + // Counter for coarse coeffs + label nCoarseCoeffs = 0; + + // Loop through all fine coeffs + forAll (upperAddr, fineCoeffi) + { + label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]]; + label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]]; + + // If the coefficient touches block zero and solo equations are + // present, skip it + if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) + { + continue; + } + + if (rmUpperAddr == rmLowerAddr) + { + // For each fine coeff inside of a coarse cluster keep the address + // of the cluster corresponding to the coeff in the + // coeffRestrictAddr as a negative index + coeffRestrictAddr[fineCoeffi] = -(rmUpperAddr + 1); + } + else + { + // This coeff is a part of a coarse coeff + + label cOwn = rmUpperAddr; + label cNei = rmLowerAddr; + + // Get coarse owner and neighbour + if (rmUpperAddr > rmLowerAddr) + { + cOwn = rmLowerAddr; + cNei = rmUpperAddr; + } + + // Check the neighbour to see if this coeff has already been found + bool nbrFound = false; + label& ccnCoeffs = blockNnbrs[cOwn]; + + for (int i = 0; i < ccnCoeffs; i++) + { + if (initCoarseNeighb[blockNbrsData[maxNnbrs*cOwn + i]] == cNei) + { + nbrFound = true; + coeffRestrictAddr[fineCoeffi] = + blockNbrsData[maxNnbrs*cOwn + i]; + break; + } + } + + if (!nbrFound) + { + if (ccnCoeffs >= maxNnbrs) + { + // Double the size of list and copy data + label oldMaxNnbrs = maxNnbrs; + maxNnbrs *= 2; + + // Resize and copy list + const labelList oldBlockNbrsData = blockNbrsData; + blockNbrsData.setSize(maxNnbrs*nCoarseEqns_); + + forAll (blockNnbrs, i) + { + for (int j = 0; j < blockNnbrs[i]; j++) + { + blockNbrsData[maxNnbrs*i + j] = + oldBlockNbrsData[oldMaxNnbrs*i + j]; + } + } + } + + blockNbrsData[maxNnbrs*cOwn + ccnCoeffs] = nCoarseCoeffs; + initCoarseNeighb[nCoarseCoeffs] = cNei; + coeffRestrictAddr[fineCoeffi] = nCoarseCoeffs; + ccnCoeffs++; + + // New coarse coeff created + nCoarseCoeffs++; + } + } + } // End for all fine coeffs + + + // Renumber into upper-triangular order + + // All coarse owner-neighbour storage + labelList coarseOwner(nCoarseCoeffs); + labelList coarseNeighbour(nCoarseCoeffs); + labelList coarseCoeffMap(nCoarseCoeffs); + + label coarseCoeffi = 0; + + forAll (blockNnbrs, cci) + { + label* cCoeffs = &blockNbrsData[maxNnbrs*cci]; + label ccnCoeffs = blockNnbrs[cci]; + + for (int i = 0; i < ccnCoeffs; i++) + { + coarseOwner[coarseCoeffi] = cci; + coarseNeighbour[coarseCoeffi] = initCoarseNeighb[cCoeffs[i]]; + coarseCoeffMap[cCoeffs[i]] = coarseCoeffi; + coarseCoeffi++; + } + } + + forAll(coeffRestrictAddr, fineCoeffi) + { + label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]]; + label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]]; + + // If the coefficient touches block zero and solo equations are + // present, skip it + if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) + { + continue; + } + + if (coeffRestrictAddr[fineCoeffi] >= 0) + { + coeffRestrictAddr[fineCoeffi] = + coarseCoeffMap[coeffRestrictAddr[fineCoeffi]]; + } + } + + // Clear the temporary storage for the coarse matrix data + blockNnbrs.setSize(0); + blockNbrsData.setSize(0); + initCoarseNeighb.setSize(0); + coarseCoeffMap.setSize(0); + + + // Create coarse-level coupled interfaces + + // Create coarse interfaces, addressing and coefficients + const label interfaceSize = + const_cast& >(matrix_).interfaces().size(); + + const typename BlockLduInterfaceFieldPtrsList::Type& + interfaceFields = + const_cast&>(matrix_).interfaces(); + + // Set the coarse interfaces and coefficients + lduInterfacePtrsList coarseInterfaces(interfaceSize); + + labelListList coarseInterfaceAddr(interfaceSize); + + // Add the coarse level + + // Set the coarse ldu addressing onto the list + autoPtr coarseAddrPtr + ( + new lduPrimitiveMesh + ( + nCoarseEqns_, + coarseOwner, + coarseNeighbour, + Pstream::worldComm, //HJ, AMG Comm fineMesh.comm(), + true + ) + ); + + // Initialise transfer of restrict addressing on the interface + // HJ, consider blocking comms. HJ, 9/Jun/2016 + forAll (interfaceFields, intI) + { + if (interfaceFields.set(intI)) + { + interfaceFields[intI].coupledInterface().initInternalFieldTransfer + ( + Pstream::blocking, + agglomIndex_ + ); + } + } + + // Store coefficients to avoid tangled communications + // HJ, 1/Apr/2009 + FieldField fineInterfaceAddr(interfaceFields.size()); + + forAll (interfaceFields, intI) + { + if (interfaceFields.set(intI)) + { + const lduInterface& fineInterface = + interfaceFields[intI].coupledInterface(); + + fineInterfaceAddr.set + ( + intI, + new labelField + ( + fineInterface.internalFieldTransfer + ( + Pstream::blocking, + agglomIndex_ + ) + ) + ); + } + } + + // Create AMG interfaces + forAll (interfaceFields, intI) + { + if (interfaceFields.set(intI)) + { + const lduInterface& fineInterface = + interfaceFields[intI].coupledInterface(); + + coarseInterfaces.set + ( + intI, + AMGInterface::New + ( + coarseAddrPtr(), + coarseInterfaces, + fineInterface, + fineInterface.interfaceInternalField(agglomIndex_), + fineInterfaceAddr[intI] + ).ptr() + ); + } + } + + forAll (interfaceFields, intI) + { + if (interfaceFields.set(intI)) + { + const AMGInterface& coarseInterface = + refCast(coarseInterfaces[intI]); + + coarseInterfaceAddr[intI] = coarseInterface.faceCells(); + } + } + + // Add interfaces + coarseAddrPtr->addInterfaces + ( + coarseInterfaces, + coarseInterfaceAddr, + matrix_.patchSchedule() + ); + + // Set the coarse level matrix + autoPtr > coarseMatrixPtr + ( + new BlockLduMatrix(coarseAddrPtr()) + ); + BlockLduMatrix& coarseMatrix = coarseMatrixPtr(); + + // Get interfaces from coarse matrix + typename BlockLduInterfaceFieldPtrsList::Type& + coarseInterfaceFieldsTransfer = coarseMatrix.interfaces(); + + // Aggolmerate the upper and lower coupled coefficients + forAll (interfaceFields, intI) + { + if (interfaceFields.set(intI)) + { + const AMGInterface& coarseInterface = + refCast(coarseInterfaces[intI]); + + coarseInterfaceFieldsTransfer.set + ( + intI, + BlockAMGInterfaceField::New + ( + coarseInterface, + interfaceFields[intI] + ).ptr() + ); + + // Since the type of clustering is now templated, clustering + // of block coefficients must be done by a FIELD (not interface) + // via a new set of virtual functions + // HJ, 16/Mar/2016 + + // Note: in the scalar AMG, clustering is done by the interface + // (always scalar) but in the block matrix it is done by a + // templated block interface field + // HJ, 16/Mar/2016 + + // Cast the interface into AMG type + const BlockAMGInterfaceField& coarseField = + refCast > + ( + coarseInterfaceFieldsTransfer[intI] + ); + + coarseMatrix.coupleUpper().set + ( + intI, + coarseField.agglomerateBlockCoeffs + ( + matrix_.coupleUpper()[intI] + ) + ); + + coarseMatrix.coupleLower().set + ( + intI, + coarseField.agglomerateBlockCoeffs + ( + matrix_.coupleLower()[intI] + ) + ); + } + } + + // Matrix restriction done! + + typedef CoeffField TypeCoeffField; + + typedef typename TypeCoeffField::squareTypeField squareTypeField; + typedef typename TypeCoeffField::linearTypeField linearTypeField; + typedef typename TypeCoeffField::scalarTypeField scalarTypeField; + + TypeCoeffField& coarseUpper = coarseMatrix.upper(); + TypeCoeffField& coarseDiag = coarseMatrix.diag(); + const TypeCoeffField& fineUpper = matrix_.upper(); + const TypeCoeffField& fineDiag = matrix_.diag(); + + // KRJ: 2013-01-31: Many cases needed as there are different combinations + + // Note: + // In coarsening, the off-diagonal coefficient type should be preserved + // and the case of taking out of off-diag from diag may need to be handled + // separately (expand the off-diag coeff to diag type before removing it + // from the diag coefficient). Since this has not been encountered yet + // only matching diag/off-diag types are handled. + // HJ, 15/Feb/2016 + if (matrix_.symmetric()) + { + if + ( + fineDiag.activeType() == blockCoeffBase::SQUARE + || fineUpper.activeType() == blockCoeffBase::SQUARE + ) + { + squareTypeField& activeCoarseDiag = coarseDiag.asSquare(); + + squareTypeField& activeCoarseUpper = coarseUpper.asSquare(); + const squareTypeField& activeFineUpper = fineUpper.asSquare(); + + // Use lower as transpose of upper + squareTypeField activeFineUpperTranspose = + activeFineUpper.T(); + + restrictDiag(fineDiag, coarseDiag); + + agglomerateCoeffs + ( + coeffRestrictAddr, + activeCoarseDiag, + activeCoarseUpper, + activeFineUpper, + activeFineUpperTranspose + ); + } + else if + ( + fineDiag.activeType() == blockCoeffBase::LINEAR + || fineUpper.activeType() == blockCoeffBase::LINEAR + ) + { + linearTypeField& activeCoarseDiag = coarseDiag.asLinear(); + + linearTypeField& activeCoarseUpper = coarseUpper.asLinear(); + const linearTypeField& activeFineUpper = fineUpper.asLinear(); + + restrictDiag(fineDiag, coarseDiag); + + agglomerateCoeffs + ( + coeffRestrictAddr, + activeCoarseDiag, + activeCoarseUpper, + activeFineUpper, + activeFineUpper + ); + } + else if + ( + fineDiag.activeType() == blockCoeffBase::SCALAR + || fineUpper.activeType() == blockCoeffBase::SCALAR + ) + { + scalarTypeField& activeCoarseDiag = coarseDiag.asScalar(); + + scalarTypeField& activeCoarseUpper = coarseUpper.asScalar(); + const scalarTypeField& activeFineUpper = fineUpper.asScalar(); + + restrictDiag(fineDiag, coarseDiag); + + agglomerateCoeffs + ( + coeffRestrictAddr, + activeCoarseDiag, + activeCoarseUpper, + activeFineUpper, + activeFineUpper + ); + } + else + { + FatalErrorIn + ( + "autoPtr >" + "BlockMatrixClustering::restrictMatrix() const" + ) << "Matrix coeff type morphing error, symmetric matrix" + << abort(FatalError); + } + } + else // asymmetric matrix + { + TypeCoeffField& coarseLower = coarseMatrix.lower(); + const TypeCoeffField& fineLower = matrix_.lower(); + + if + ( + fineDiag.activeType() == blockCoeffBase::SQUARE + || fineUpper.activeType() == blockCoeffBase::SQUARE + ) + { + squareTypeField& activeCoarseDiag = coarseDiag.asSquare(); + + squareTypeField& activeCoarseUpper = coarseUpper.asSquare(); + const squareTypeField& activeFineUpper = fineUpper.asSquare(); + + squareTypeField& activeCoarseLower = coarseLower.asSquare(); + const squareTypeField& activeFineLower = fineLower.asSquare(); + + restrictDiag(fineDiag, coarseDiag); + + agglomerateCoeffs + ( + coeffRestrictAddr, + activeCoarseDiag, + activeCoarseUpper, + activeFineUpper, + activeCoarseLower, + activeFineLower + ); + } + else if + ( + fineDiag.activeType() == blockCoeffBase::LINEAR + || fineUpper.activeType() == blockCoeffBase::LINEAR + ) + { + linearTypeField& activeCoarseDiag = coarseDiag.asLinear(); + + linearTypeField& activeCoarseUpper = coarseUpper.asLinear(); + const linearTypeField& activeFineUpper = fineUpper.asLinear(); + + linearTypeField& activeCoarseLower = coarseLower.asLinear(); + const linearTypeField& activeFineLower = fineLower.asLinear(); + + restrictDiag(fineDiag, coarseDiag); + + agglomerateCoeffs + ( + coeffRestrictAddr, + activeCoarseDiag, + activeCoarseUpper, + activeFineUpper, + activeCoarseLower, + activeFineLower + ); + } + else if + ( + fineDiag.activeType() == blockCoeffBase::SCALAR + || fineUpper.activeType() == blockCoeffBase::SCALAR + ) + { + scalarTypeField& activeCoarseDiag = coarseDiag.asScalar(); + + scalarTypeField& activeCoarseUpper = coarseUpper.asScalar(); + const scalarTypeField& activeFineUpper = fineUpper.asScalar(); + + scalarTypeField& activeCoarseLower = coarseLower.asScalar(); + const scalarTypeField& activeFineLower = fineLower.asScalar(); + + restrictDiag(fineDiag, coarseDiag); + + agglomerateCoeffs + ( + coeffRestrictAddr, + activeCoarseDiag, + activeCoarseUpper, + activeFineUpper, + activeCoarseLower, + activeFineLower + ); + } + else + { + FatalErrorIn + ( + "autoPtr >" + "BlockMatrixClustering::restrictMatrix() const" + ) << "Matrix coeff type morphing error, asymmetric matrix" + << abort(FatalError); + } + } + + // Create and return BlockAMGLevel + return autoPtr > + ( + new coarseBlockAMGLevel + ( + coarseAddrPtr, + coarseMatrixPtr, + this->dict(), + this->type(), + this->groupSize(), + this->minCoarseEqns() + ) + ); +} + + +template +void Foam::BlockMatrixClustering::restrictResidual +( + const Field& res, + Field& coarseRes +) const +{ + coarseRes = pTraits::zero; + + forAll (res, i) + { + coarseRes[agglomIndex_[i]] += res[i]; + } +} + + +template +void Foam::BlockMatrixClustering::prolongateCorrection +( + Field& x, + const Field& coarseX +) const +{ + forAll (x, i) + { + x[i] += coarseX[agglomIndex_[i]]; + } +} + + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.H new file mode 100644 index 000000000..6e56ec25b --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.H @@ -0,0 +1,218 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend 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 foam-extend. If not, see . + +Class + BlockMatrixClustering + +Description + Block matrix AMG coarsening by Jasak clustering algorithm + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved. + +SourceFiles + BlockMatrixClustering.C + +\*---------------------------------------------------------------------------*/ + +#ifndef BlockMatrixClustering_H +#define BlockMatrixClustering_H + +#include "BlockMatrixCoarsening.H" +#include "BlockLduMatrix.H" +#include "BlockCoeffNorm.H" +#include "BlockCoeff.H" +#include "tolerancesSwitch.H" +#include "cpuTime.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +/*---------------------------------------------------------------------------*\ + Class BlockMatrixClustering Declaration +\*---------------------------------------------------------------------------*/ + +template +class BlockMatrixClustering +: + public BlockMatrixCoarsening +{ + // Private Data + + //- Reference to matrix + const BlockLduMatrix& matrix_; + + //- Min cluster size + const label minClusterSize_; + + //- Max cluster size + const label maxClusterSize_; + + //- Reference to a templated norm calculator + autoPtr > normPtr_; + + //- Child array: for each fine equation give clustering cluster + // index + labelField agglomIndex_; + + //- Group size + label groupSize_; + + //- Number of solo cells + label nSolo_; + + //- Number of coarse equations + label nCoarseEqns_; + + //- Can a coarse level be constructed? + bool coarsen_; + + cpuTime lTime_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + BlockMatrixClustering(const BlockMatrixClustering&); + + // Disallow default bitwise assignment + void operator=(const BlockMatrixClustering&); + + + //- Calculate clustering index (child) + void calcClustering(); + + //- Restrict CoeffField. Used for diag coefficient + void restrictDiag + ( + const CoeffField& Coeff, + CoeffField& coarseCoeff + ) const; + + //- Agglomerate coeffs, symmetric matrix + template + void agglomerateCoeffs + ( + const labelList& coeffRestrictAddr, + Field& activeCoarseDiag, + Field& activeCoarseUpper, + const Field& activeFineUpper, + const Field& activeFineUpperTranspose + ) const; + + //- Agglomerate coeffs, assymmetric matrix + template + void agglomerateCoeffs + ( + const labelList& coeffRestrictAddr, + Field& activeCoarseDiag, + Field& activeCoarseUpper, + const Field& activeFineUpper, + Field& activeCoarseLower, + const Field& activeFineLower + ) const; + + //- Restrict CoeffField, decoupled version. Used for diag coefficient + void restrictDiagDecoupled + ( + const CoeffField& Coeff, + CoeffField& coarseCoeff + ) const; + + // Private Static Data + + //- Weighting factor + static const debug::tolerancesSwitch weightFactor_; + + //- Diagonal scaling factor + static const debug::tolerancesSwitch diagFactor_; + + +public: + + //- Runtime type information + TypeName("cluster"); + + + // Constructors + + //- Construct from matrix and group size + BlockMatrixClustering + ( + const BlockLduMatrix& matrix, + const dictionary& dict, + const label groupSize, + const label minCoarseEqns + ); + + + //- Destructor + virtual ~BlockMatrixClustering(); + + + // Member Functions + + //- Can a coarse level be constructed? + virtual bool coarsen() const + { + return coarsen_; + } + + //- Restrict matrix + virtual autoPtr > restrictMatrix() const; + + //- Restrict residual + virtual void restrictResidual + ( + const Field& res, + Field& coarseRes + ) const; + + //- Prolongate correction + virtual void prolongateCorrection + ( + Field& x, + const Field& coarseX + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "BlockMatrixClustering.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/blockMatrixClusterings.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/blockMatrixClusterings.C new file mode 100644 index 000000000..b919fb77b --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/blockMatrixClusterings.C @@ -0,0 +1,43 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend 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 foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "blockMatrixClusterings.H" +#include "blockMatrixCoarsenings.H" +#include "coarseBlockAMGLevel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +makeBlockMatrixCoarsenings(blockMatrixClustering); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/blockMatrixClusterings.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/blockMatrixClusterings.H new file mode 100644 index 000000000..2eeb1f390 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/blockMatrixClusterings.H @@ -0,0 +1,65 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend 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 foam-extend. If not, see . + +Class + BlockMatrixClustering + +Description + Typedefs for block matrix AMG coarsening by Jasak clustering algorithm. + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved. + +SourceFiles + blockMatrixClustering.C + +\*---------------------------------------------------------------------------*/ + +#ifndef blockMatrixClusterings_H +#define blockMatrixClusterings_H + +#include "scalarBlockMatrixClustering.H" +#include "tensorBlockMatrixClustering.H" +#include "BlockMatrixClustering.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +typedef BlockMatrixClustering blockMatrixClusteringScalar; +typedef BlockMatrixClustering blockMatrixClusteringVector; +typedef BlockMatrixClustering blockMatrixClusteringTensor; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/scalarBlockMatrixClustering.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/scalarBlockMatrixClustering.H new file mode 100644 index 000000000..a274599fe --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/scalarBlockMatrixClustering.H @@ -0,0 +1,86 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend 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 foam-extend. If not, see . + +Class + BlockMatrixClustering + +Description + Specialisation of the BlockMatrixClustering for scalars. + +Author + Klas Jareteg, 2013-01-31 + +\*---------------------------------------------------------------------------*/ + +#ifndef scalarBlockMatrixClustering_H +#define scalarBlockMatrixClustering_H + +#include "blockMatrixCoarsenings.H" +#include "blockMatrixClusterings.H" +#include "BlockMatrixClustering.H" +#include "BlockMatrixCoarsening.H" +#include "runTimeSelectionTables.H" +#include "scalarBlockLduMatrix.H" +#include "scalarBlockConstraint.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * Forward declaration of template friend fuctions * * * * * * * // + + +/*---------------------------------------------------------------------------*\ + Class BlockMatrixClustering Declaration +\*---------------------------------------------------------------------------*/ + +// Disable restrict matrix: no square type +template<> +inline autoPtr > +BlockMatrixClustering::restrictMatrix() const +{ + FatalErrorIn + ( + "autoPtr > " + "BlockMatrixClustering::restrictMatrix() const" + ) << "Function not implemented for Type=scalar. " << endl + << abort(FatalError); + + // Dummy return to keep compiler happy + return autoPtr >(NULL); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/tensorBlockMatrixClustering.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/tensorBlockMatrixClustering.H new file mode 100644 index 000000000..5d3f31946 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/tensorBlockMatrixClustering.H @@ -0,0 +1,83 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend 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 foam-extend. If not, see . + +Class + BlockMatrixClustering + +Description + Specialisation of the BlockMatrixClustering for tensors. + +Author + Klas Jareteg, 2013-01-31 + +\*---------------------------------------------------------------------------*/ + +#ifndef tensorBlockMatrixClustering_H +#define tensorBlockMatrixClustering_H + +#include "blockMatrixCoarsenings.H" +#include "blockMatrixClusterings.H" +#include "BlockMatrixClustering.H" +#include "runTimeSelectionTables.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class BlockMatrixClustering Declaration +\*---------------------------------------------------------------------------*/ + +// Disable restrict matrix: no square type +template<> +inline autoPtr > +BlockMatrixClustering::restrictMatrix +() const +{ + FatalErrorIn + ( + "autoPtr > " + "BlockMatrixClustering::" + "restrictMatrix() const" + ) << "Function not implemented for Type=tensor. " << endl + << abort(FatalError); + + // Dummy return to keep compiler happy + return autoPtr >(NULL); +} + + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.H b/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.H index a2d6f4711..fbeb70d75 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.H @@ -49,7 +49,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class coarseBlockAMGLevel Declaration + Class coarseBlockAMGLevel Declaration \*---------------------------------------------------------------------------*/ template diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H index 3ff3de320..b748adb23 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H @@ -189,10 +189,9 @@ public: ); - // Destructor - - virtual ~BlockLduSolver() - {} + //- Destructor + virtual ~BlockLduSolver() + {} // Member functions diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C index a1318306b..b8ccb4f40 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C @@ -50,6 +50,7 @@ License #include "blockAMGPrecons.H" #include "blockMatrixCoarsenings.H" #include "blockMatrixAgglomerations.H" +#include "blockMatrixClusterings.H" #include "blockCoeffNorms.H" #include "blockCoeffTwoNorms.H" #include "blockCoeffMaxNorms.H" @@ -147,6 +148,9 @@ defineTemplateRunTimeSelectionTable(block##Type##MatrixCoarsening, matrix); \ typedef BlockMatrixAgglomeration block##Type##MatrixAgglomeration; \ makeBlockMatrixCoarsening(block##Type##MatrixCoarsening, block##Type##MatrixAgglomeration); \ \ +typedef BlockMatrixClustering block##Type##MatrixClustering; \ +makeBlockMatrixCoarsening(block##Type##MatrixCoarsening, block##Type##MatrixClustering); \ + \ typedef BlockCoeffNorm block##Type##CoeffNorm; \ defineNamedTemplateTypeNameAndDebug(block##Type##CoeffNorm, 0); \ defineTemplateRunTimeSelectionTable(block##Type##CoeffNorm, dictionary); \ diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffComponentNorm/BlockCoeffComponentNorm.C b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffComponentNorm/BlockCoeffComponentNorm.C index d4e71f5ea..f9fc86ba8 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffComponentNorm/BlockCoeffComponentNorm.C +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffComponentNorm/BlockCoeffComponentNorm.C @@ -50,18 +50,18 @@ Foam::scalar Foam::BlockCoeffComponentNorm::normalize const Foam::BlockCoeff& a ) { - return mag(a.component(cmpt_)); + return a.component(cmpt_); } template -void Foam::BlockCoeffComponentNorm::coeffMag +void Foam::BlockCoeffComponentNorm::normalize ( - const Foam::CoeffField& a, - Foam::Field& b + Foam::Field& b, + const Foam::CoeffField& a ) { - b = mag(a.component(cmpt_)); + b = a.component(cmpt_); } diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffComponentNorm/BlockCoeffComponentNorm.H b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffComponentNorm/BlockCoeffComponentNorm.H index 518fc9089..bb0c8689b 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffComponentNorm/BlockCoeffComponentNorm.H +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffComponentNorm/BlockCoeffComponentNorm.H @@ -95,10 +95,9 @@ public: const dictionary& dict ); - // Destructor - - virtual ~BlockCoeffComponentNorm() - {} + //- Destructor + virtual ~BlockCoeffComponentNorm() + {} // Member functions @@ -109,10 +108,10 @@ public: const BlockCoeff& a ); - virtual void coeffMag + virtual void normalize ( - const CoeffField& a, - Field& b + Field& b, + const CoeffField& a ); }; diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/BlockCoeffMaxNorm.C b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/BlockCoeffMaxNorm.C index 51bc5f275..7b10fc89d 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/BlockCoeffMaxNorm.C +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/BlockCoeffMaxNorm.C @@ -46,9 +46,14 @@ Foam::BlockCoeffMaxNorm::BlockCoeffMaxNorm template Foam::scalar Foam::BlockCoeffMaxNorm::normalize ( - const Foam::BlockCoeff& a + const BlockCoeff& a ) { + // Note. This does not properly account for the sign of the + // off-diagonal coefficient. If the off-diag is negative + // the function should look for cmptMin and vice-versa + // HJ, 28/Feb/2017 + if (a.activeType() == BlockCoeff::SCALAR) { return mag(a.asScalar()); @@ -74,23 +79,23 @@ Foam::scalar Foam::BlockCoeffMaxNorm::normalize template -void Foam::BlockCoeffMaxNorm::coeffMag +void Foam::BlockCoeffMaxNorm::normalize ( - const Foam::CoeffField& a, - Foam::Field& b + Field& b, + const CoeffField& a ) { if (a.activeType() == BlockCoeff::SCALAR) { - b = mag(a.asScalar()); + b = a.asScalar(); } else if (a.activeType() == BlockCoeff::LINEAR) { - b = cmptMax(cmptMag(a.asLinear())); + b = cmptMax(a.asLinear()); } else if (a.activeType() == BlockCoeff::SQUARE) { - b = cmptMax(cmptMag(a.asSquare())); + b = cmptMax(a.asSquare()); } else { diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/BlockCoeffMaxNorm.H b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/BlockCoeffMaxNorm.H index 97404d4e5..acc7ef137 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/BlockCoeffMaxNorm.H +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/BlockCoeffMaxNorm.H @@ -92,10 +92,9 @@ public: const dictionary& dict ); - // Destructor - - virtual ~BlockCoeffMaxNorm() - {} + //- Destructor + virtual ~BlockCoeffMaxNorm() + {} // Member functions @@ -106,10 +105,10 @@ public: const BlockCoeff& a ); - virtual void coeffMag + virtual void normalize ( - const CoeffField& a, - Field& b + Field& b, + const CoeffField& a ); }; diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/scalarBlockCoeffMaxNorm.H b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/scalarBlockCoeffMaxNorm.H index 398acef39..4ca39f88c 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/scalarBlockCoeffMaxNorm.H +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/scalarBlockCoeffMaxNorm.H @@ -68,10 +68,10 @@ inline scalar BlockCoeffMaxNorm::normalize template<> -inline void BlockCoeffMaxNorm::coeffMag +inline void BlockCoeffMaxNorm::normalize ( - const CoeffField& a, - Field& b + Field& b, + const CoeffField& a ) { b = mag(a.asScalar()); diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/tensorBlockCoeffMaxNorm.H b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/tensorBlockCoeffMaxNorm.H index 7a0c3a51e..1e99f7632 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/tensorBlockCoeffMaxNorm.H +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffMaxNorm/tensorBlockCoeffMaxNorm.H @@ -82,10 +82,10 @@ inline scalar BlockCoeffMaxNorm::normalize template<> -inline void BlockCoeffMaxNorm::coeffMag +inline void BlockCoeffMaxNorm::normalize ( - const CoeffField& a, - Field& b + Field& b, + const CoeffField& a ) { if (a.activeType() == BlockCoeff::SCALAR) @@ -94,6 +94,12 @@ inline void BlockCoeffMaxNorm::coeffMag } else if (a.activeType() == BlockCoeff::LINEAR) { + // Note. This does not properly account for the sign of the + // off-diagonal coefficient. If the off-diag is negative, the return + // should look for cmptMin and vice-versa + // Consider this unreliable + // HJ, 28/Feb/2017 + b = cmptMax(cmptMag(a.asLinear())); } else diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffNorm/BlockCoeffNorm.H b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffNorm/BlockCoeffNorm.H index 453efc4a2..e5aeea3fd 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffNorm/BlockCoeffNorm.H +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffNorm/BlockCoeffNorm.H @@ -86,22 +86,15 @@ public: autoPtr, BlockCoeffNorm, dictionary, - ( - const dictionary& dict - ), - ( - dict - ) + (const dictionary& dict), + (dict) ); // Constructors //- Construct from dictionary - BlockCoeffNorm - ( - const dictionary& dict - ); + BlockCoeffNorm(const dictionary& dict); // Selectors @@ -113,10 +106,9 @@ public: ); - // Destructor - - virtual ~BlockCoeffNorm() - {} + //- Destructor + virtual ~BlockCoeffNorm() + {} // Member functions @@ -127,10 +119,10 @@ public: const BlockCoeff& a ) = 0; - virtual void coeffMag + virtual void normalize ( - const CoeffField& a, - Field& b + Field& b, + const CoeffField& a ) = 0; }; diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/BlockCoeffTwoNorm.C b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/BlockCoeffTwoNorm.C index 1ff743b76..dfdb1f025 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/BlockCoeffTwoNorm.C +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/BlockCoeffTwoNorm.C @@ -51,21 +51,25 @@ Foam::scalar Foam::BlockCoeffTwoNorm::normalize { if (a.activeType() == Foam::BlockCoeff::SCALAR) { + // Note: for two-norm, use mag return mag(a.asScalar()); } else if (a.activeType() == Foam::BlockCoeff::LINEAR) { + // Note: for two-norm, use mag return mag(a.asLinear()); } else if (a.activeType() == Foam::BlockCoeff::SQUARE) { + // Note: for two-norm, use mag return mag(a.asSquare()); } else { FatalErrorIn ( - "scalar BlockCoeffTwoNorm(const BlockCoeff& a)" + "scalar BlockCoeffTwoNorm::normalize" + "(const BlockCoeff& a)" ) << "Unknown type" << abort(FatalError); return 0; @@ -77,22 +81,25 @@ Foam::scalar Foam::BlockCoeffTwoNorm::normalize template -void Foam::BlockCoeffTwoNorm::coeffMag +void Foam::BlockCoeffTwoNorm::normalize ( - const Foam::CoeffField& a, - Foam::Field& b + Field& b, + const CoeffField& a ) { if (a.activeType() == Foam::BlockCoeff::SCALAR) { + // Note: for two-norm, use mag b = mag(a.asScalar()); } else if (a.activeType() == Foam::BlockCoeff::LINEAR) { + // Note: for two-norm, use mag b = mag(a.asLinear()); } else if (a.activeType() == Foam::BlockCoeff::SQUARE) { + // Note: for two-norm, use mag b = mag(a.asSquare()); } else diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/BlockCoeffTwoNorm.H b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/BlockCoeffTwoNorm.H index cad9966e5..0bb5ae150 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/BlockCoeffTwoNorm.H +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/BlockCoeffTwoNorm.H @@ -92,10 +92,10 @@ public: const dictionary& dict ); - // Destructor - virtual ~BlockCoeffTwoNorm() - {} + //- Destructor + virtual ~BlockCoeffTwoNorm() + {} // Member functions @@ -106,10 +106,10 @@ public: const BlockCoeff& a ); - virtual void coeffMag + virtual void normalize ( - const CoeffField& a, - Field& b + Field& b, + const CoeffField& a ); }; diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/tensorBlockCoeffTwoNorm.H b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/tensorBlockCoeffTwoNorm.H index f564a930c..db2a2bc7d 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/tensorBlockCoeffTwoNorm.H +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/tensorBlockCoeffTwoNorm.H @@ -65,10 +65,12 @@ inline scalar BlockCoeffTwoNorm::normalize { if (a.activeType() == BlockCoeff::SCALAR) { + // Note: for two-norm, use mag return mag(a.asScalar()); } else if (a.activeType() == BlockCoeff::LINEAR) { + // Note: for two-norm, use mag return mag(a.asLinear()); } else @@ -84,18 +86,20 @@ inline scalar BlockCoeffTwoNorm::normalize template<> -inline void BlockCoeffTwoNorm::coeffMag +inline void BlockCoeffTwoNorm::normalize ( - const CoeffField& a, - Field& b + Field& b, + const CoeffField& a ) { if (a.activeType() == BlockCoeff::SCALAR) { + // Note: for two-norm, use mag b = mag(a.asScalar()); } else if (a.activeType() == BlockCoeff::LINEAR) { + // Note: for two-norm, use mag b = mag(a.asLinear()); } else diff --git a/src/lduSolvers/amg/amgPolicy/amgPolicy.H b/src/lduSolvers/amg/amgPolicy/amgPolicy.H index 2ddaab5da..d92e0d4a7 100644 --- a/src/lduSolvers/amg/amgPolicy/amgPolicy.H +++ b/src/lduSolvers/amg/amgPolicy/amgPolicy.H @@ -123,10 +123,9 @@ public: {} - // Destructor - - virtual ~amgPolicy() - {} + //- Destructor + virtual ~amgPolicy() + {} // Member Functions diff --git a/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.C b/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.C new file mode 100644 index 000000000..8c3822966 --- /dev/null +++ b/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.C @@ -0,0 +1,1021 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend 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 foam-extend. If not, see . + +Class + clusterAmgPolicy + +Description + Block matrix AMG coarsening by Jasak clustering algorithm + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#include "clusterAmgPolicy.H" +#include "amgMatrix.H" +#include "boolList.H" +#include "addToRunTimeSelectionTable.H" +#include "AMGInterfaceField.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(clusterAmgPolicy, 0); + + addToRunTimeSelectionTable(amgPolicy, clusterAmgPolicy, matrix); + +} // End namespace Foam + + +const Foam::debug::tolerancesSwitch Foam::clusterAmgPolicy::weightFactor_ +( + "aamgWeightFactor", + 0.65 +); + + +const Foam::debug::tolerancesSwitch Foam::clusterAmgPolicy::diagFactor_ +( + "aamgDiagFactor", + 1e-8 +); + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::clusterAmgPolicy::calcChild() +{ + if (matrix_.diagonal()) + { + // Diag only matrix. Reset and return + child_ = 0; + nCoarseEqns_ = 1; + + return; + } + + // Algorithm: + // 1) Calculate appropriate connection strength for symmetric and asymmetric + // matrix + // 2) Collect solo (disconnected) cells and remove them from agglomeration + // by placing them into cluster zero (solo group) + // 3) Loop through all equations and for each equation find the best fit + // neighbour. If all neighbours are grouped, add equation + // to best group + + // Initialise child array + child_ = -1; + + const label nRows = matrix_.lduAddr().size(); + + // Get matrix addressing + const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); + const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); + const unallocLabelList& losortAddr = matrix_.lduAddr().losortAddr(); + + const unallocLabelList& ownerStartAddr = + matrix_.lduAddr().ownerStartAddr(); + const unallocLabelList& losortStartAddr = + matrix_.lduAddr().losortStartAddr(); + + + // Calculate agglomeration + + // Get magnitudes of matrix coefficients + + const scalarField& diag = matrix_.diag(); + const scalarField magDiag = mag(matrix_.diag()); + + // Calculate magnitude of strong positive connections + scalarField magUpper(upperAddr.size(), 0); + scalarField magLower(upperAddr.size(), 0); + + // Note: + // Matrix properties are no longer assumed, eg. the diag and off-diag + // sign is checked + // HJ, 29/Mar/2017 + + // Note: negative connections are eliminated in max(...) below + // HJ, 30/Mar/2017 + + if (matrix_.hasUpper()) + { + const scalarField& upper = matrix_.upper(); + + // Owner: upper triangle + forAll (lowerAddr, coeffI) + { + // Sign of strong positive upper is opposite of the sign of + // its diagonal coefficient + magUpper[coeffI] = + Foam::max(-1*sign(diag[lowerAddr[coeffI]])*upper[coeffI], 0); + } + } + + if (matrix_.hasLower()) + { + const scalarField& lower = matrix_.lower(); + + // Neighbour: lower triangle + forAll (lowerAddr, coeffI) + { + // Sign of strong positive upper is opposite of the sign of + // its diagonal coefficient + magLower[coeffI] = + Foam::max(-1*sign(diag[upperAddr[coeffI]])*lower[coeffI], 0); + } + } + else + { + magLower = magUpper; + } + + labelList sizeOfGroups(nRows, 0); + + nCoarseEqns_ = 0; + + // Gather disconnected and weakly connected equations into cluster zero + // Weak connection is assumed to be the one where the off-diagonal + // coefficient is smaller than diagFactor_()*diag + { + // Algorithm + // Mark all cells to belong to zero cluster + // Go through all upper and lower coefficients and for the ones + // larger than threshold mark the equations out of cluster zero + + scalarField magScaledDiag = diagFactor_()*magDiag; + + boolList zeroCluster(magDiag.size(), true); + + if (matrix_.symmetric()) + { + // Owner: upper triangle + forAll (lowerAddr, coeffI) + { + if (magUpper[coeffI] > magScaledDiag[lowerAddr[coeffI]]) + { + zeroCluster[lowerAddr[coeffI]] = false; + } + } + + // Neighbour: lower triangle with symm coefficients + forAll (upperAddr, coeffI) + { + if (magUpper[coeffI] > magScaledDiag[upperAddr[coeffI]]) + { + zeroCluster[upperAddr[coeffI]] = false; + } + } + } + else if (matrix_.asymmetric()) + { + // Owner: upper triangle + forAll (lowerAddr, coeffI) + { + if (magUpper[coeffI] > magScaledDiag[lowerAddr[coeffI]]) + { + zeroCluster[lowerAddr[coeffI]] = false; + } + } + + // Neighbour: lower triangle with lower coeffs + forAll (upperAddr, coeffI) + { + if (magLower[coeffI] > magScaledDiag[upperAddr[coeffI]]) + { + zeroCluster[upperAddr[coeffI]] = false; + } + } + } + + // Collect solo equations + forAll (zeroCluster, eqnI) + { + if (zeroCluster[eqnI]) + { + // Found solo equation + nSolo_++; + + child_[eqnI] = nCoarseEqns_; + } + } + + if (nSolo_ > 0) + { + // Found solo equations + sizeOfGroups[nCoarseEqns_] = nSolo_; + + nCoarseEqns_++; + } + } + + + // Loop through cells + // - if the cell is not already grouped, open a new group (seed) + // - find stroungest grouped and ungrouped connection: + // - grouped connection is towards an already grouped cell + // - ungrouped connection is towards an ungrouped cell + // - decide if a grouped or ungrouped connection is better + + label curEqn, indexUngrouped, indexGrouped, colI, groupPassI, + nextUngrouped, nextGrouped; + + scalar curWeight, weightUngrouped, weightGrouped; + + for (label rowI = 0; rowI < nRows; rowI++) + { + if (child_[rowI] == -1) + { + // Found new ungrouped equation + curEqn = rowI; + + // Reset grouped and upgrouped index + indexUngrouped = -1; + indexGrouped = -1; + + // Make next group (coarse equation) and search for neighbours + child_[curEqn] = nCoarseEqns_; + + // Work on the group until the min group size is satisfied + // As each new element of the group is found, group search starts + // from this element until group cluster size is satisfied + for + ( + groupPassI = 1; + groupPassI < minClusterSize_; + groupPassI++ + ) + { + weightUngrouped = 0; + weightGrouped = 0; + + indexUngrouped = -1; + indexGrouped = -1; + + nextUngrouped = -1; + nextGrouped = -1; + + // Visit all neighbour equations + + // Visit upper triangle coefficients from the equation + for + ( + label rowCoeffI = ownerStartAddr[curEqn]; + rowCoeffI < ownerStartAddr[curEqn + 1]; + rowCoeffI++ + ) + { + // Get column index, upper triangle + colI = upperAddr[rowCoeffI]; + + curWeight = magUpper[rowCoeffI]/magDiag[curEqn]; + + if (child_[colI] == -1) + { + // Found ungrouped neighbour + if + ( + indexUngrouped == -1 + || curWeight > weightUngrouped + ) + { + // Found or updated ungrouped neighbour + indexUngrouped = rowCoeffI; + weightUngrouped = curWeight; + + // Record possible next equation for search + nextUngrouped = colI; + } + } + else if (child_[curEqn] != child_[colI]) + { + // Check for neighbour in solo group + if (nSolo_ == 0 || child_[colI] != 0) + { + // Found neighbour belonging to other group + if (indexGrouped == -1 || curWeight > weightGrouped) + { + // Found or updated grouped neighbour + indexGrouped = rowCoeffI; + weightGrouped = curWeight; + + // Record possible next equation for search + nextGrouped = colI; + } + } + } + } + + // Visit lower triangle coefficients from the equation + for + ( + label rowCoeffI = losortStartAddr[curEqn]; + rowCoeffI < losortStartAddr[curEqn + 1]; + rowCoeffI++ + ) + { + // Get column index, lower triangle + colI = lowerAddr[losortAddr[rowCoeffI]]; + + curWeight = magLower[losortAddr[rowCoeffI]]/ + magDiag[curEqn]; + + if (child_[colI] == -1) + { + // Found first or better ungrouped neighbour + if + ( + indexUngrouped == -1 + || curWeight > weightUngrouped + ) + { + // Found or updated ungrouped neighbour + indexUngrouped = rowCoeffI; + weightUngrouped = curWeight; + + // Record possible next equation for search + nextUngrouped = colI; + } + } + else if (child_[curEqn] != child_[colI]) + { + // Check for neighbour in solo group + if (nSolo_ == 0 || child_[colI] != 0) + { + // Found first of better neighbour belonging to + // other group + if + ( + indexGrouped == -1 + || curWeight > weightGrouped + ) + { + // Found or updated grouped neighbour + indexGrouped = rowCoeffI; + weightGrouped = curWeight; + + // Record possible next equation for search + nextGrouped = colI; + } + } + } + } + + // Decide what to do with current equation + + // If ungrouped candidate exists and it is stronger + // than best weighted grouped candidate, use it + if + ( + indexUngrouped != -1 + && ( + indexGrouped == -1 + || weightUngrouped >= weightFactor_()*weightGrouped + ) + ) + { + // Found new element of group. Add it and use as + // start of next search + + child_[nextUngrouped] = child_[curEqn]; + sizeOfGroups[child_[curEqn]]++; + + // Search from nextEqn + curEqn = nextUngrouped; + } + else + { + // Group full or cannot be extended with new + // ungrouped candidates + break; + } + } + + + // Finished group passes + if + ( + groupPassI > 1 + || indexGrouped == -1 + || sizeOfGroups[child_[nextGrouped]] > maxClusterSize_ + ) + { + // There is no group to put this equation into + sizeOfGroups[child_[rowI]]++; + nCoarseEqns_++; + } + else + { + // Dump current cell into the best group available + child_[rowI] = child_[nextGrouped]; + sizeOfGroups[child_[nextGrouped]]++; + } + } + } + + // Resize group size count + sizeOfGroups.setSize(nCoarseEqns_); + + // The decision on parallel agglomeration needs to be made for the + // whole gang of processes; otherwise I may end up with a different + // number of agglomeration levels on different processors. + + // If the number of coarse equations is less than minimum and + // if the matrix has reduced in size by at least 1/3, coarsen + if (nCoarseEqns_ > minCoarseEqns() && 3*nCoarseEqns_ <= 2*nRows) + { + coarsen_ = true; + } + + reduce(coarsen_, andOp()); + + // if (lduMatrix::debug >= 2) + { + // Count solo cells + label nSingleClusters = 0; + + // Adjust start based on solo cells cluster status + label start = 0; + + if (nSolo_ > 0) + { + start = 1; + } + + for (label i = start; i < sizeOfGroups.size(); i++) + { + if (sizeOfGroups[i] == 1) + { + nSingleClusters++; + } + } + + Pout<< "Coarse level size: " << nCoarseEqns_ + << " nSolo = " << nSolo_ + << " cluster (" << min(sizeOfGroups) << " " + << max(sizeOfGroups) + << "). N singleton clusters = " << nSingleClusters; + + if (coarsen_) + { + Pout << ". Accepted" << endl; + } + else + { + Pout << ". Rejected" << endl; + } + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::clusterAmgPolicy::clusterAmgPolicy +( + const lduMatrix& matrix, + const label groupSize, + const label minCoarseEqns +) +: + amgPolicy(groupSize, minCoarseEqns), + matrix_(matrix), + minClusterSize_(groupSize), + maxClusterSize_(2*groupSize), + child_(matrix_.lduAddr().size()), + nSolo_(0), + nCoarseEqns_(0), + coarsen_(false) +{ + if (groupSize < 2) + { + FatalErrorIn + ( + "clusterAmgPolicy::clusterAmgPolicy\n" + "(\n" + " const lduMatrix& matrix,\n" + " const label groupSize,\n" + " const label minCoarseEqns\n" + ")" + ) << "Group size smaller than 2 is not allowed" + << abort(FatalError); + } + + calcChild(); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::clusterAmgPolicy::~clusterAmgPolicy() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::autoPtr Foam::clusterAmgPolicy::restrictMatrix +( + const FieldField& bouCoeffs, + const FieldField& intCoeffs, + const lduInterfaceFieldPtrsList& interfaceFields +) const +{ + if (!coarsen_) + { + FatalErrorIn + ( + "autoPtr clusterAmgPolicy::restrictMatrix() const" + ) << "Requesting coarse matrix when it cannot be created" + << abort(FatalError); + } + + // Construct the coarse matrix and ldu addressing for the next level + // Algorithm: + // 1) Loop through all fine coeffs. If the child labels on two sides are + // different, this creates a coarse coeff. Define owner and neighbour + // for this coeff based on cluster IDs. + // 2) Check if the coeff has been seen before. If yes, add the coefficient + // to the appropriate field (stored with the equation). If no, create + // a new coeff with neighbour ID and add the coefficient + // 3) Once all the coeffs have been created, loop through all clusters and + // insert the coeffs in the upper order. At the same time, collect the + // owner and neighbour addressing. + // 4) Agglomerate the diagonal by summing up the fine diagonal + + // Get addressing + const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); + + const label nFineCoeffs = upperAddr.size(); + +# ifdef FULLDEBUG + if (child_.size() != matrix_.lduAddr().size()) + { + FatalErrorIn + ( + "autoPtr clusterAmgPolicy::restrictMatrix() const" + ) << "Child array does not correspond to fine level. " << endl + << " Child size: " << child_.size() + << " number of equations: " << matrix_.lduAddr().size() + << abort(FatalError); + } +# endif + + + // Does the matrix have solo equations + bool soloEqns = nSolo_ > 0; + + // Storage for block neighbours and coefficients + + // Guess initial maximum number of neighbours in block + label maxNnbrs = 10; + + // Number of neighbours per block + labelList blockNnbrs(nCoarseEqns_, 0); + + // Setup initial packed storage for neighbours and coefficients + labelList blockNbrsData(maxNnbrs*nCoarseEqns_); + + // Create face-restriction addressing + // Note: value of coeffRestrictAddr for off-diagonal coefficients + // touching solo cells will be invalid + // HJ, 7/Apr/2015 + labelList coeffRestrictAddr(nFineCoeffs); + + // Initial neighbour array (not in upper-triangle order) + labelList initCoarseNeighb(nFineCoeffs); + + // Counter for coarse coeffs + label nCoarseCoeffs = 0; + + // Note on zero cluster coarsening + // If the matrix contains a solo group, it will be in index zero. + // Since solo equations are disconnected from the rest of the matrix + // they do not create new off-diagonal coefficients + + // Loop through all fine coeffs + forAll (upperAddr, fineCoeffI) + { + label rmUpperAddr = child_[upperAddr[fineCoeffI]]; + label rmLowerAddr = child_[lowerAddr[fineCoeffI]]; + + // If the coefficient touches block zero and solo equations are + // present, skip it + if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) + { + continue; + } + + if (rmUpperAddr == rmLowerAddr) + { + // For each fine coeff inside of a coarse cluster keep the address + // of the cluster corresponding to the coeff in the + // coeffRestrictAddr as a negative index + coeffRestrictAddr[fineCoeffI] = -(rmUpperAddr + 1); + } + else + { + // This coeff is a part of a coarse coeff + + label cOwn = rmUpperAddr; + label cNei = rmLowerAddr; + + // Get coarse owner and neighbour + if (rmUpperAddr > rmLowerAddr) + { + cOwn = rmLowerAddr; + cNei = rmUpperAddr; + } + + // Check the neighbour to see if this coeff has already been found + bool nbrFound = false; + label& ccnCoeffs = blockNnbrs[cOwn]; + + for (int i = 0; i < ccnCoeffs; i++) + { + if (initCoarseNeighb[blockNbrsData[maxNnbrs*cOwn + i]] == cNei) + { + nbrFound = true; + coeffRestrictAddr[fineCoeffI] = + blockNbrsData[maxNnbrs*cOwn + i]; + break; + } + } + + if (!nbrFound) + { + if (ccnCoeffs >= maxNnbrs) + { + // Double the size of list and copy data + label oldMaxNnbrs = maxNnbrs; + maxNnbrs *= 2; + + // Resize and copy list + const labelList oldBlockNbrsData = blockNbrsData; + blockNbrsData.setSize(maxNnbrs*nCoarseEqns_); + + forAll (blockNnbrs, i) + { + for (int j = 0; j < blockNnbrs[i]; j++) + { + blockNbrsData[maxNnbrs*i + j] = + oldBlockNbrsData[oldMaxNnbrs*i + j]; + } + } + } + + blockNbrsData[maxNnbrs*cOwn + ccnCoeffs] = nCoarseCoeffs; + initCoarseNeighb[nCoarseCoeffs] = cNei; + coeffRestrictAddr[fineCoeffI] = nCoarseCoeffs; + ccnCoeffs++; + + // New coarse coeff created + nCoarseCoeffs++; + } + } + } // End for all fine coeffs + + + // Renumber into upper-triangular order + + // All coarse owner-neighbour storage + labelList coarseOwner(nCoarseCoeffs); + labelList coarseNeighbour(nCoarseCoeffs); + labelList coarseCoeffMap(nCoarseCoeffs); + + label coarseCoeffi = 0; + + forAll (blockNnbrs, cci) + { + label* cCoeffs = &blockNbrsData[maxNnbrs*cci]; + label ccnCoeffs = blockNnbrs[cci]; + + for (int i = 0; i < ccnCoeffs; i++) + { + coarseOwner[coarseCoeffi] = cci; + coarseNeighbour[coarseCoeffi] = initCoarseNeighb[cCoeffs[i]]; + coarseCoeffMap[cCoeffs[i]] = coarseCoeffi; + coarseCoeffi++; + } + } + + forAll (coeffRestrictAddr, fineCoeffI) + { + label rmUpperAddr = child_[upperAddr[fineCoeffI]]; + label rmLowerAddr = child_[lowerAddr[fineCoeffI]]; + + // If the coefficient touches block zero and solo equations are + // present, skip it + if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) + { + continue; + } + + if (coeffRestrictAddr[fineCoeffI] >= 0) + { + coeffRestrictAddr[fineCoeffI] = + coarseCoeffMap[coeffRestrictAddr[fineCoeffI]]; + } + } + + // Clear the temporary storage for the coarse matrix data + blockNnbrs.setSize(0); + blockNbrsData.setSize(0); + initCoarseNeighb.setSize(0); + coarseCoeffMap.setSize(0); + + + // Create coarse-level coupled interfaces + + // Create coarse interfaces, addressing and coefficients + + // Set the coarse interfaces and coefficients + lduInterfacePtrsList* coarseInterfacesPtr = + new lduInterfacePtrsList(interfaceFields.size()); + lduInterfacePtrsList& coarseInterfaces = *coarseInterfacesPtr; + + // Set the coarse interfaceFields and coefficients + lduInterfaceFieldPtrsList* coarseInterfaceFieldsPtr = + new lduInterfaceFieldPtrsList(interfaceFields.size()); + lduInterfaceFieldPtrsList& coarseInterfaceFields = + *coarseInterfaceFieldsPtr; + + FieldField* coarseBouCoeffsPtr = + new FieldField(interfaceFields.size()); + FieldField& coarseBouCoeffs = *coarseBouCoeffsPtr; + + FieldField* coarseIntCoeffsPtr = + new FieldField(interfaceFields.size()); + FieldField& coarseIntCoeffs = *coarseIntCoeffsPtr; + + labelListList coarseInterfaceAddr(interfaceFields.size()); + + // Add the coarse level + + // Set the coarse ldu addressing onto the list + lduPrimitiveMesh* coarseAddrPtr = + new lduPrimitiveMesh + ( + nCoarseEqns_, + coarseOwner, + coarseNeighbour, + true + ); + + // Initialise transfer of restrict addressing on the interface + forAll (interfaceFields, intI) + { + if (interfaceFields.set(intI)) + { + interfaceFields[intI].coupledInterface().initInternalFieldTransfer + ( + Pstream::blocking, + child_ + ); + } + } + + // Store coefficients to avoid tangled communications + // HJ, 1/Apr/2009 + FieldField fineInterfaceAddr(interfaceFields.size()); + + forAll (interfaceFields, intI) + { + if (interfaceFields.set(intI)) + { + const lduInterface& fineInterface = + interfaceFields[intI].coupledInterface(); + + fineInterfaceAddr.set + ( + intI, + new labelField + ( + fineInterface.internalFieldTransfer + ( + Pstream::blocking, + child_ + ) + ) + ); + } + } + + // Create AMG interfaces + forAll (interfaceFields, intI) + { + if (interfaceFields.set(intI)) + { + const lduInterface& fineInterface = + interfaceFields[intI].coupledInterface(); + + coarseInterfaces.set + ( + intI, + AMGInterface::New + ( + *coarseAddrPtr, + coarseInterfaces, + fineInterface, + fineInterface.interfaceInternalField(child_), + fineInterfaceAddr[intI] + ).ptr() + ); + } + } + + forAll (interfaceFields, intI) + { + if (interfaceFields.set(intI)) + { + const AMGInterface& coarseInterface = + refCast(coarseInterfaces[intI]); + + coarseInterfaceFields.set + ( + intI, + AMGInterfaceField::New + ( + coarseInterface, + interfaceFields[intI] + ).ptr() + ); + + // Note: scalar agglomeration is done by the interface + // (always scalar) but in the block matrix it is done by a + // templated block interface field + // HJ, 16/Mar/2016 + coarseBouCoeffs.set + ( + intI, + coarseInterface.agglomerateCoeffs(bouCoeffs[intI]) + ); + + coarseIntCoeffs.set + ( + intI, + coarseInterface.agglomerateCoeffs(intCoeffs[intI]) + ); + + coarseInterfaceAddr[intI] = coarseInterface.faceCells(); + } + } + + // Add interfaces + coarseAddrPtr->addInterfaces + ( + *coarseInterfacesPtr, + coarseInterfaceAddr, + matrix_.patchSchedule() + ); + + // Matrix restriction done! + + // Set the coarse level matrix + lduMatrix* coarseMatrixPtr = new lduMatrix(*coarseAddrPtr); + lduMatrix& coarseMatrix = *coarseMatrixPtr; + + // Coarse matrix diagonal initialised by restricting the + // finer mesh diagonal + scalarField& coarseDiag = coarseMatrix.diag(); + restrictResidual(matrix_.diag(), coarseDiag); + + // Check if matrix is assymetric and if so agglomerate both upper and lower + // coefficients ... + if (matrix_.hasLower()) + { + // Get off-diagonal matrix coefficients + const scalarField& fineUpper = matrix_.upper(); + const scalarField& fineLower = matrix_.lower(); + + // Coarse matrix upper coefficients + scalarField& coarseUpper = coarseMatrix.upper(); + scalarField& coarseLower = coarseMatrix.lower(); + + forAll (coeffRestrictAddr, fineCoeffI) + { + label rmUpperAddr = child_[upperAddr[fineCoeffI]]; + label rmLowerAddr = child_[lowerAddr[fineCoeffI]]; + + // If the coefficient touches block zero and solo equations are + // present, skip it + if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) + { + continue; + } + + label cCoeff = coeffRestrictAddr[fineCoeffI]; + + if (cCoeff >= 0) + { + coarseUpper[cCoeff] += fineUpper[fineCoeffI]; + coarseLower[cCoeff] += fineLower[fineCoeffI]; + } + else + { + // Add the fine face coefficients into the diagonal. + coarseDiag[-1 - cCoeff] += + fineUpper[fineCoeffI] + fineLower[fineCoeffI]; + } + } + } + else // ... Otherwise it is symmetric so agglomerate just the upper + { + // Get off-diagonal matrix coefficients + const scalarField& fineUpper = matrix_.upper(); + + // Coarse matrix upper coefficients + scalarField& coarseUpper = coarseMatrix.upper(); + + forAll (coeffRestrictAddr, fineCoeffI) + { + label rmUpperAddr = child_[upperAddr[fineCoeffI]]; + label rmLowerAddr = child_[lowerAddr[fineCoeffI]]; + + // If the coefficient touches block zero and solo equations are + // present, skip it + if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) + { + continue; + } + + label cCoeff = coeffRestrictAddr[fineCoeffI]; + + if (cCoeff >= 0) + { + coarseUpper[cCoeff] += fineUpper[fineCoeffI]; + } + else + { + // Add the fine face coefficient into the diagonal. + coarseDiag[-1 - cCoeff] += 2*fineUpper[fineCoeffI]; + } + } + } + + // Create and return amgMatrix + return autoPtr + ( + new amgMatrix + ( + coarseAddrPtr, + coarseInterfacesPtr, + coarseMatrixPtr, + coarseBouCoeffsPtr, + coarseIntCoeffsPtr, + coarseInterfaceFieldsPtr + ) + ); +} + + +void Foam::clusterAmgPolicy::restrictResidual +( + const scalarField& res, + scalarField& coarseRes +) const +{ + coarseRes = 0; + + forAll (res, i) + { + coarseRes[child_[i]] += res[i]; + } +} + + +void Foam::clusterAmgPolicy::prolongateCorrection +( + scalarField& x, + const scalarField& coarseX +) const +{ + forAll (x, i) + { + x[i] += coarseX[child_[i]]; + } +} + + +// ************************************************************************* // diff --git a/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.H b/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.H new file mode 100644 index 000000000..3ba6dd905 --- /dev/null +++ b/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.H @@ -0,0 +1,165 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend 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 foam-extend. If not, see . + +Class + clusterAmgPolicy + +Description + Clustering AMG policy + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + +SourceFiles + clusterAmgPolicy.C + +\*---------------------------------------------------------------------------*/ + +#ifndef clusterAmgPolicy_H +#define clusterAmgPolicy_H + +#include "amgPolicy.H" +#include "lduMatrix.H" +#include "tolerancesSwitch.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class clusterAmgPolicy Declaration +\*---------------------------------------------------------------------------*/ + +class clusterAmgPolicy +: + public amgPolicy +{ + // Private Data + + //- Reference to matrix + const lduMatrix& matrix_; + + //- Min cluster size + const label minClusterSize_; + + //- Max cluster size + const label maxClusterSize_; + + //- Child array: for each fine equation give coarse cluster index + labelField child_; + + //- Number of solo cells + label nSolo_; + + //- Number of coarse equations + label nCoarseEqns_; + + //- Can a coarse level be constructed? + bool coarsen_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + clusterAmgPolicy(const clusterAmgPolicy&); + + //- Disallow default bitwise assignment + void operator=(const clusterAmgPolicy&); + + //- Calculate child + void calcChild(); + + + // Private Static Data + + //- Weighting factor + static const debug::tolerancesSwitch weightFactor_; + + //- Diagonal scaling factor + static const debug::tolerancesSwitch diagFactor_; + + +public: + + //- Runtime type information + TypeName("cluster"); + + + // Constructors + + //- Construct from matrix and group size + clusterAmgPolicy + ( + const lduMatrix& matrix, + const label groupSize, + const label minCoarseEqns + ); + + + //- Destructor + virtual ~clusterAmgPolicy(); + + + // Member Functions + + //- Can a coarse level be constructed? + virtual bool coarsen() const + { + return coarsen_; + } + + //- Restrict matrix + virtual autoPtr restrictMatrix + ( + const FieldField& bouCoeffs, + const FieldField& intCoeffs, + const lduInterfaceFieldPtrsList& interfaces + ) const; + + //- Restrict residual + virtual void restrictResidual + ( + const scalarField& res, + scalarField& coarseRes + ) const; + + //- Prolongate correction + virtual void prolongateCorrection + ( + scalarField& x, + const scalarField& coarseX + ) const; + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lduSolvers/amg/amgPolicy/pamgPolicy.C b/src/lduSolvers/amg/amgPolicy/pamgPolicy.C index a4d302e32..8617fc302 100644 --- a/src/lduSolvers/amg/amgPolicy/pamgPolicy.C +++ b/src/lduSolvers/amg/amgPolicy/pamgPolicy.C @@ -61,6 +61,15 @@ Foam::pamgPolicy::diagFactor_ void Foam::pamgPolicy::calcChild() { + if (matrix_.diagonal()) + { + // Diag only matrix. Reset and return + child_ = 0; + nCoarseEqns_ = 1; + + return; + } + // Algorithm: // 1) Create temporary equation addressing using a double-pass algorithm. // to create the offset table. @@ -196,12 +205,6 @@ void Foam::pamgPolicy::calcChild() { // Found solo equations nCoarseEqns_++; - - if (lduMatrix::debug >= 2) - { - Pout<< "Found " << nSolo_ << " weakly connected equations." - << endl; - } } } @@ -281,6 +284,8 @@ void Foam::pamgPolicy::calcChild() { // This is a Dirichlet point: no neighbours // Put it into its own cluster + // This should have been already handled by the zero + // cluster. HJ, 29/Mar/2017 child_[eqnI] = nCoarseEqns_; nCoarseEqns_++; } @@ -308,7 +313,8 @@ void Foam::pamgPolicy::calcChild() if (lduMatrix::debug >= 2) { - Pout << "Coarse level size: " << nCoarseEqns_; + Pout<< "Coarse level size: " << nCoarseEqns_ + << " nSolo = " << nSolo_; if (coarsen_) { diff --git a/src/lduSolvers/amg/amgPolicy/pamgPolicy.H b/src/lduSolvers/amg/amgPolicy/pamgPolicy.H index 91e3cbeae..67d19c69e 100644 --- a/src/lduSolvers/amg/amgPolicy/pamgPolicy.H +++ b/src/lduSolvers/amg/amgPolicy/pamgPolicy.H @@ -107,9 +107,9 @@ public: const label minCoarseEqns ); - // Destructor - virtual ~pamgPolicy(); + //- Destructor + virtual ~pamgPolicy(); // Member Functions From 775661dfdfeca9df68b1641430baed06ee266afe Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:28:48 +0100 Subject: [PATCH 012/333] Formatting --- .../FieldFields/oneFieldField/oneFieldField.H | 2 +- .../vectorNFieldFields/VectorNFieldFields.H | 30 +++++++++---------- .../GeometricField/GeometricField.H | 8 ++--- 3 files changed, 18 insertions(+), 22 deletions(-) diff --git a/src/foam/fields/FieldFields/oneFieldField/oneFieldField.H b/src/foam/fields/FieldFields/oneFieldField/oneFieldField.H index 5d380d9b8..4a2fb88ff 100644 --- a/src/foam/fields/FieldFields/oneFieldField/oneFieldField.H +++ b/src/foam/fields/FieldFields/oneFieldField/oneFieldField.H @@ -45,7 +45,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class oneField Declaration + Class oneFieldField Declaration \*---------------------------------------------------------------------------*/ class oneFieldField diff --git a/src/foam/fields/FieldFields/vectorNFieldFields/VectorNFieldFields.H b/src/foam/fields/FieldFields/vectorNFieldFields/VectorNFieldFields.H index 3e4533a31..9bfac1286 100644 --- a/src/foam/fields/FieldFields/vectorNFieldFields/VectorNFieldFields.H +++ b/src/foam/fields/FieldFields/vectorNFieldFields/VectorNFieldFields.H @@ -43,21 +43,21 @@ SourceFiles // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -#define VectorN_FieldFunctions(tensorType, diagTensorType, sphericalTensorType, \ - vectorType, CmptType, args...) \ - \ -UNARY_FUNCTION(CmptType, vectorType, cmptSum) \ - \ -BINARY_FUNCTION(vectorType, vectorType, vectorType, cmptMultiply) \ -BINARY_TYPE_FUNCTION(vectorType, vectorType, vectorType, cmptMultiply) \ - \ -BINARY_OPERATOR(vectorType, vectorType, CmptType, /,divide) \ -BINARY_TYPE_OPERATOR(vectorType, vectorType, CmptType, /,divide) \ - \ -BINARY_OPERATOR(vectorType, vectorType, vectorType, +, add) \ -BINARY_OPERATOR(vectorType, vectorType, vectorType, -, subtract) \ - \ -BINARY_TYPE_OPERATOR(vectorType, vectorType, vectorType, +, add) \ +#define VectorN_FieldFunctions(tensorType, diagTensorType, sphericalTensorType, \ + vectorType, CmptType, args...) \ + \ +UNARY_FUNCTION(CmptType, vectorType, cmptSum) \ + \ +BINARY_FUNCTION(vectorType, vectorType, vectorType, cmptMultiply) \ +BINARY_TYPE_FUNCTION(vectorType, vectorType, vectorType, cmptMultiply) \ + \ +BINARY_OPERATOR(vectorType, vectorType, CmptType, /,divide) \ +BINARY_TYPE_OPERATOR(vectorType, vectorType, CmptType, /,divide) \ + \ +BINARY_OPERATOR(vectorType, vectorType, vectorType, +, add) \ +BINARY_OPERATOR(vectorType, vectorType, vectorType, -, subtract) \ + \ +BINARY_TYPE_OPERATOR(vectorType, vectorType, vectorType, +, add) \ BINARY_TYPE_OPERATOR(vectorType, vectorType, vectorType, -, subtract) diff --git a/src/foam/fields/GeometricFields/GeometricField/GeometricField.H b/src/foam/fields/GeometricFields/GeometricField/GeometricField.H index 10cd0776a..c2e7059a6 100644 --- a/src/foam/fields/GeometricFields/GeometricField/GeometricField.H +++ b/src/foam/fields/GeometricFields/GeometricField/GeometricField.H @@ -304,8 +304,7 @@ public: const word& patchFieldType=PatchField::calculatedType() ); - //- Constructor given IOobject, mesh, dimensioned - // and patch types. + //- Construct given IOobject, mesh, dimensioned and patch types GeometricField ( const IOobject&, @@ -314,9 +313,7 @@ public: const wordList& patchFieldTypes ); - //- Constructor from components - -#if ( !defined(SWIG) || (SWIG_VERSION > 0x010340) ) + //- Construct from components GeometricField ( const IOobject&, @@ -325,7 +322,6 @@ public: const Field&, const PtrList >& ); -#endif //- Construct and read given IOobject GeometricField From b21c18342269939a26b29c6dd06a0ffeb2735e24 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:29:29 +0100 Subject: [PATCH 013/333] Clean-up --- src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H | 3 +++ src/finiteVolume/fvMesh/fvMeshLduAddressing.H | 7 +++---- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H b/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H index b022cb433..4fbf36543 100644 --- a/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H +++ b/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H @@ -101,6 +101,9 @@ public: ); + // Destructor - default + + // Member functions // Access diff --git a/src/finiteVolume/fvMesh/fvMeshLduAddressing.H b/src/finiteVolume/fvMesh/fvMeshLduAddressing.H index 62fadaa37..2a2f21aae 100644 --- a/src/finiteVolume/fvMesh/fvMeshLduAddressing.H +++ b/src/finiteVolume/fvMesh/fvMeshLduAddressing.H @@ -102,10 +102,9 @@ public: } - // Destructor - - virtual ~fvMeshLduAddressing() - {} + //- Destructor + virtual ~fvMeshLduAddressing() + {} // Member Functions From 287e02f674195a65a56883cdb851f56c4a91619b Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:29:48 +0100 Subject: [PATCH 014/333] Make comm() a virtual function in fvMesh --- src/finiteVolume/fvMesh/fvMesh.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/finiteVolume/fvMesh/fvMesh.H b/src/finiteVolume/fvMesh/fvMesh.H index 86d083a80..1ea79365c 100644 --- a/src/finiteVolume/fvMesh/fvMesh.H +++ b/src/finiteVolume/fvMesh/fvMesh.H @@ -288,7 +288,7 @@ public: // Communication support //- Return communicator used for parallel communication - label comm() const + virtual int comm() const { return polyMesh::comm(); } From 113e5189069692d0203fd0b3a4b741977b670771 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:31:39 +0100 Subject: [PATCH 015/333] Clean-up --- .../pressureInletVelocityFvPatchVectorField.C | 2 +- src/finiteVolume/finiteVolume/fvSolution/fvSolution.H | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/finiteVolume/fields/fvPatchFields/derived/pressureInletVelocity/pressureInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/pressureInletVelocity/pressureInletVelocityFvPatchVectorField.C index 5d0a70820..1a0270fef 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/pressureInletVelocity/pressureInletVelocityFvPatchVectorField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/pressureInletVelocity/pressureInletVelocityFvPatchVectorField.C @@ -125,7 +125,7 @@ void Foam::pressureInletVelocityFvPatchVectorField::updateCoeffs() const surfaceScalarField& phi = db().lookupObject(phiName_); - const fvsPatchField& phip = + const fvsPatchScalarField& phip = patch().patchField(phi); vectorField n = patch().nf(); diff --git a/src/finiteVolume/finiteVolume/fvSolution/fvSolution.H b/src/finiteVolume/finiteVolume/fvSolution/fvSolution.H index e854fb28f..bfd08a810 100644 --- a/src/finiteVolume/finiteVolume/fvSolution/fvSolution.H +++ b/src/finiteVolume/finiteVolume/fvSolution/fvSolution.H @@ -51,8 +51,10 @@ class fvSolution { // Private Member Functions - //- Disallow default bitwise copy construct and assignment + //- Disallow default bitwise assignment fvSolution(const fvSolution&); + + //- Disallow default bitwise copy construct void operator=(const fvSolution&); From 6b24ec03166067f3e23e5e336fa6bfa0851f7224 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:32:19 +0100 Subject: [PATCH 016/333] crMatrix move --- src/foam/Make/files | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/foam/Make/files b/src/foam/Make/files index f0f928e7d..c00becfe5 100644 --- a/src/foam/Make/files +++ b/src/foam/Make/files @@ -238,6 +238,10 @@ $(LUscalarMatrix)/LUscalarMatrix.C $(LUscalarMatrix)/procLduMatrix.C $(LUscalarMatrix)/procLduInterface.C +crMatrix = matrices/crMatrix +$(crMatrix)/crAddressing.C +$(crMatrix)/crMatrix.C + lduMatrix = matrices/lduMatrix $(lduMatrix)/lduMatrix/lduMatrix.C $(lduMatrix)/lduMatrix/lduMatrixOperations.C @@ -697,6 +701,8 @@ $(BlockAMGInterfaceFields)/GGIBlockAMGInterfaceField/GGIBlockAMGInterfaceFields. BlockMatrixCoarsening = $(BlockAMG)/BlockMatrixCoarsening $(BlockMatrixCoarsening)/BlockMatrixCoarsening/blockMatrixCoarsenings.C $(BlockMatrixCoarsening)/BlockMatrixAgglomeration/blockMatrixAgglomerations.C +$(BlockMatrixCoarsening)/BlockMatrixClustering/blockMatrixClusterings.C +$(BlockMatrixCoarsening)/BlockSelectiveAMG/blockSelectiveAMGs.C BlockLduPrecons = matrices/blockLduMatrix/BlockLduPrecons $(BlockLduPrecons)/BlockLduPrecon/blockLduPrecons.C From 395ecdbe9dad5fd2ce95bb86d408ffd6dbd7d972 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:32:37 +0100 Subject: [PATCH 017/333] Formatting --- src/foam/meshes/MeshObject/meshObjectBase.H | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/foam/meshes/MeshObject/meshObjectBase.H b/src/foam/meshes/MeshObject/meshObjectBase.H index d4269cc9a..8bbd9c7de 100644 --- a/src/foam/meshes/MeshObject/meshObjectBase.H +++ b/src/foam/meshes/MeshObject/meshObjectBase.H @@ -175,10 +175,9 @@ public: } - // Destructor - - virtual ~meshObjectBase() - {} + //- Destructor + virtual ~meshObjectBase() + {} // Member Functions From 5aaf1f0bece9b7241d6464e08860684263cc3c67 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:33:02 +0100 Subject: [PATCH 018/333] Report internal field min/max separately --- .../field/minMaxField/minMaxField.C | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/postProcessing/functionObjects/field/minMaxField/minMaxField.C b/src/postProcessing/functionObjects/field/minMaxField/minMaxField.C index 70cafbc73..0b9d21853 100644 --- a/src/postProcessing/functionObjects/field/minMaxField/minMaxField.C +++ b/src/postProcessing/functionObjects/field/minMaxField/minMaxField.C @@ -92,8 +92,12 @@ bool Foam::minMaxField::execute() fieldName_ ); - Info<< "Field " << fieldName_ << " min = " << Foam::min(f).value() - << " max = " << Foam::max(f).value() << endl; + Info<< "Field " << fieldName_ + << " min = " << Foam::min(f).value() + << " (" << gMin(f.internalField()) << ")" + << " max = " << max(f).value() + << " (" << gMax(f.internalField()) << ")" + << endl; return true; } @@ -103,9 +107,12 @@ bool Foam::minMaxField::execute() volScalarField magF = mag(f); - Info<< "Field " << fieldName_ << " magnitude min = " - << Foam::min(magF).value() - << " max = " << Foam::max(magF).value() << endl; + Info<< "Field " << fieldName_ + << " magnitude min = " << Foam::min(magF).value() + << " (" << gMin(magF.internalField()) << ")" + << " max = " << max(magF).value() + << " (" << gMax(magF.internalField()) << ")" + << endl; return true; } From 29b8708ceafa78908a7017888d8a9c67631c2651 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:34:10 +0100 Subject: [PATCH 019/333] Bugfix: Update time loop syntax --- .../surfaceTracking/interTrackFoam/interTrackFoam.C | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/applications/solvers/surfaceTracking/interTrackFoam/interTrackFoam.C b/applications/solvers/surfaceTracking/interTrackFoam/interTrackFoam.C index 8dbdb73fc..0b9e62069 100644 --- a/applications/solvers/surfaceTracking/interTrackFoam/interTrackFoam.C +++ b/applications/solvers/surfaceTracking/interTrackFoam/interTrackFoam.C @@ -56,9 +56,13 @@ int main(int argc, char *argv[]) Info << "\nStarting time loop\n" << endl; - for (runTime++; !runTime.end(); runTime++) + while (runTime.run()) { - Info << "Time = " << runTime.value() << endl << endl; +# include "readTimeControls.H" +# include "CourantNo.H" +# include "setSurfaceStabilityDeltaT.H" + +# include "readPISOControls.H" interface.moveMeshPointsForOldFreeSurfDisplacement(); @@ -116,6 +120,8 @@ int main(int argc, char *argv[]) { phi -= pEqn.flux(); } + + p.relax(); } # include "continuityErrs.H" From ba5893ba4a88b56697f880747e5a84be4ce8510b Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:34:18 +0100 Subject: [PATCH 020/333] Formatting --- .../solvers/incompressible/pimpleDyMFoam/correctPhi.H | 2 +- .../solvers/surfaceTracking/freeSurface/freeSurface.H | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/applications/solvers/incompressible/pimpleDyMFoam/correctPhi.H b/applications/solvers/incompressible/pimpleDyMFoam/correctPhi.H index 23a957bb5..faf2ab8b8 100644 --- a/applications/solvers/incompressible/pimpleDyMFoam/correctPhi.H +++ b/applications/solvers/incompressible/pimpleDyMFoam/correctPhi.H @@ -26,7 +26,7 @@ IOobject::NO_WRITE ), mesh, - dimensionedScalar("pcorr", p.dimensions(), 0.0), + dimensionedScalar("pcorr", p.dimensions(), 0), pcorrTypes ); diff --git a/applications/solvers/surfaceTracking/freeSurface/freeSurface.H b/applications/solvers/surfaceTracking/freeSurface/freeSurface.H index e977e56a6..7ca6ba835 100644 --- a/applications/solvers/surfaceTracking/freeSurface/freeSurface.H +++ b/applications/solvers/surfaceTracking/freeSurface/freeSurface.H @@ -135,6 +135,7 @@ class freeSurface //- Interface smoothing at the begining of time step Switch smoothing_; + // Demand-driven data //- Patch to patch interpolation object which deals with @@ -229,6 +230,7 @@ class freeSurface const vector& axis ) const; + public: // Declare name of the class and it's debug switch @@ -247,9 +249,8 @@ public: ); - // Destructor - - ~freeSurface(); + //- Destructor + virtual ~freeSurface(); // Member Functions From c66ebe7b9150df759cc464a611ab1386773f08b0 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:35:36 +0100 Subject: [PATCH 021/333] Bugfix: surface stability deltaT control. Dirk Grunding --- .../setSurfaceStabilityDeltaT.H | 55 +++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 applications/solvers/surfaceTracking/interTrackFoam/setSurfaceStabilityDeltaT.H diff --git a/applications/solvers/surfaceTracking/interTrackFoam/setSurfaceStabilityDeltaT.H b/applications/solvers/surfaceTracking/interTrackFoam/setSurfaceStabilityDeltaT.H new file mode 100644 index 000000000..818a0cdea --- /dev/null +++ b/applications/solvers/surfaceTracking/interTrackFoam/setSurfaceStabilityDeltaT.H @@ -0,0 +1,55 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend 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 foam-extend. If not, see . + +Global + setSurfaceStabilityDeltaT + +Description + Reset the timestep to maintain a constant maximum courant Number. + Reduction of time-step is immediate, but increase is damped to avoid + unstable oscillations. + +\*---------------------------------------------------------------------------*/ + +if (adjustTimeStep) +{ + scalar maxDeltaTFact = maxCo/(CoNum + SMALL); + scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2); + + runTime.setDeltaT + ( + Foam::min + ( + interface.surfStabCrit().value(), + Foam::min + ( + deltaTFact*runTime.deltaT().value(), + maxDeltaT + ) + ) + ); + + Info<< "deltaT = " << runTime.deltaT().value() << endl; +} + +// ************************************************************************* // From 3f31a7d4c9c40117076b579773ab9f2e72e4e9c3 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 10:36:09 +0100 Subject: [PATCH 022/333] Improved formulation of steady compressible solvers. Jasak, Cvijetic, Rusche, Zuzul --- .../steadyCompressibleFoam/UEqn.H | 4 ++-- .../steadyCompressibleFoam/hEqn.H | 8 +++---- .../steadyCompressibleFoam/pEqn.H | 6 +----- .../steadyCompressibleMRFFoam/UEqn.H | 7 ++----- .../steadyCompressibleMRFFoam/iEqn.H | 21 +++++++++++++++++-- .../steadyCompressibleMRFFoam/pEqn.H | 14 ++----------- .../steadyCompressibleSRFFoam/UEqn.H | 4 ++-- .../steadyCompressibleSRFFoam/iEqn.H | 6 +++--- .../steadyCompressibleSRFFoam/pEqn.H | 14 ++----------- 9 files changed, 37 insertions(+), 47 deletions(-) diff --git a/applications/solvers/compressible/steadyCompressibleFoam/UEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/UEqn.H index 10f11a421..2833a63e9 100644 --- a/applications/solvers/compressible/steadyCompressibleFoam/UEqn.H +++ b/applications/solvers/compressible/steadyCompressibleFoam/UEqn.H @@ -3,8 +3,8 @@ fvVectorMatrix UEqn ( - fvm::ddt(rho, U) - + fvm::div(phi, U) + fvm::div(phi, U) + + fvm::SuSp(-fvc::div(phi), U) + turbulence->divDevRhoReff() ); diff --git a/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H index c93a3b888..d9378a743 100644 --- a/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H +++ b/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H @@ -1,13 +1,14 @@ { // Solve the enthalpy equation in total enthalpy formulation (see K) + T.storePrevIter(); + K = 0.5*(magSqr(U)); fvScalarMatrix hEqn ( - fvm::ddt(rho, h) - + fvc::ddt(rho, K) - + fvm::div(phi, h) + fvm::div(phi, h) + + fvm::SuSp(-fvc::div(phi), h) + fvc::div(phi, K) - fvm::laplacian(turbulence->alphaEff(), h) == @@ -16,7 +17,6 @@ ); hEqn.relax(); - hEqn.solve(); // Bounding of enthalpy taken out diff --git a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H index e954c7b13..f8fa0d40b 100644 --- a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H @@ -31,11 +31,7 @@ { fvScalarMatrix pEqn ( - fvm::ddt(psis, p) - + fvm::div(phid, p) - // Convective flux relaxation terms - + fvm::SuSp(-divPhid, p) - + divPhid*p + fvm::div(phid, p) + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H index 7f769de78..3e2f20ea0 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H @@ -1,10 +1,7 @@ - // Solve the momentum equation - U.storePrevIter(); - fvVectorMatrix UEqn ( - fvm::ddt(rho, U) - + fvm::div(phi, U) + fvm::div(phi, U) + + fvm::SuSp(-fvc::div(phi), U) + turbulence->divDevRhoReff() ); diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H index 910176930..720d3928b 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H @@ -5,13 +5,30 @@ Urel == U; mrfZones.relativeVelocity(Urel); + // Bound the relative velocity to preserve the i to h conversion bound + // HJ, 22/Mar/2017 + volScalarField magUrel = mag(Urel); + + if (max(magUrel) > UMax) + { + volScalarField UrelLimiter = pos(magUrel - UMax)*UMax/(magUrel + smallU) + + neg(magUrel - UMax); + UrelLimiter.max(scalar(0)); + UrelLimiter.min(scalar(1)); + + Urel *= UrelLimiter; + Urel.correctBoundaryConditions(); + } + // Create rotational velocity (= omega x r) Urot == U - Urel; + T.storePrevIter(); + fvScalarMatrix iEqn ( - fvm::ddt(rho, i) - + fvm::div(phi, i) + fvm::div(phi, i) + + fvm::SuSp(-fvc::div(phi), i) - fvm::laplacian(turbulence->alphaEff(), i) == // Viscous heating: note sign (devRhoReff has a minus in it) diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H index 61fed64e8..bedf65b19 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H @@ -16,7 +16,7 @@ mrfZones.correctBoundaryVelocity(U); // Calculate phi for boundary conditions - phi = rhof*fvc::interpolate(U) & mesh.Sf(); + phi = rhof*(fvc::interpolate(U) & mesh.Sf()); surfaceScalarField phid2 = rhoReff/rhof*phi; @@ -29,21 +29,11 @@ p.storePrevIter(); - volScalarField divPhid - ( - "divPhid", - fvc::div(phid) - ); - while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( - fvm::ddt(psis, p) - + fvm::div(phid, p) - // Convective flux relaxation terms - + fvm::SuSp(-divPhid, p) - + divPhid*p + fvm::div(phid, p) + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H index 727548f3b..9e09b5c42 100644 --- a/applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H +++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H @@ -3,8 +3,8 @@ fvVectorMatrix UrelEqn ( - fvm::ddt(rho, Urel) - + fvm::div(phi, Urel) + fvm::div(phi, Urel) + + fvm::SuSp(-fvc::div(phi), Urel) + turbulence->divDevRhoReff() + rho*SRF->Su() ); diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H index c37796f06..36f61d359 100644 --- a/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H +++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H @@ -1,10 +1,11 @@ { // Solve the rothalpy equation + T.storePrevIter(); fvScalarMatrix iEqn ( - fvm::ddt(rho, i) - + fvm::div(phi, i) + fvm::div(phi, i) + + fvm::SuSp(-fvc::div(phi), i) - fvm::laplacian(turbulence->alphaEff(), i) == // Viscous heating: note sign (devRhoReff has a minus in it) @@ -12,7 +13,6 @@ ); iEqn.relax(); - iEqn.solve(); // Calculate enthalpy out of rothalpy diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H index 59e03e2c4..762b129dc 100644 --- a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H @@ -20,23 +20,13 @@ p.storePrevIter(); - volScalarField divPhid - ( - "divPhid", - fvc::div(phid) - ); - while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( - fvm::ddt(psis, p) - + fvm::div(phid, p) - // Convective flux relaxation terms - + fvm::SuSp(-divPhid, p) - + divPhid*p - + fvc::div(phid2) + fvm::div(phid, p) - fvm::laplacian(rho*rUrelA, p) + + fvc::div(phid2) ); pEqn.solve(); From a8654cedcd786589fbf2f8110088e660443432ea Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 13:28:52 +0100 Subject: [PATCH 023/333] Clean-up of files files --- src/foam/Make/files | 1 - src/lduSolvers/Make/files | 3 --- 2 files changed, 4 deletions(-) diff --git a/src/foam/Make/files b/src/foam/Make/files index c00becfe5..3f3b90bfb 100644 --- a/src/foam/Make/files +++ b/src/foam/Make/files @@ -702,7 +702,6 @@ BlockMatrixCoarsening = $(BlockAMG)/BlockMatrixCoarsening $(BlockMatrixCoarsening)/BlockMatrixCoarsening/blockMatrixCoarsenings.C $(BlockMatrixCoarsening)/BlockMatrixAgglomeration/blockMatrixAgglomerations.C $(BlockMatrixCoarsening)/BlockMatrixClustering/blockMatrixClusterings.C -$(BlockMatrixCoarsening)/BlockSelectiveAMG/blockSelectiveAMGs.C BlockLduPrecons = matrices/blockLduMatrix/BlockLduPrecons $(BlockLduPrecons)/BlockLduPrecon/blockLduPrecons.C diff --git a/src/lduSolvers/Make/files b/src/lduSolvers/Make/files index 1a94853f8..2267afc05 100644 --- a/src/lduSolvers/Make/files +++ b/src/lduSolvers/Make/files @@ -1,6 +1,3 @@ -crMatrix/crAddressing.C -crMatrix/crMatrix.C - lduPrecon = lduPrecon $(lduPrecon)/CholeskyPrecon/CholeskyPrecon.C $(lduPrecon)/ILU0/ILU0.C From df50c6e4a9c842787414f77d13633d4e6d0a66f6 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 13:46:41 +0100 Subject: [PATCH 024/333] Rename group size in AMG agglomeration --- .../BlockMatrixClustering/BlockMatrixClustering.C | 8 ++++---- .../BlockMatrixClustering/BlockMatrixClustering.H | 11 +++++------ src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.C | 8 ++++---- src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.H | 4 ++-- 4 files changed, 15 insertions(+), 16 deletions(-) diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C index 11f027efc..587904a74 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C @@ -273,7 +273,7 @@ void Foam::BlockMatrixClustering::calcClustering() for ( groupPassI = 1; - groupPassI < minClusterSize_; + groupPassI < minGroupSize_; groupPassI++ ) { @@ -428,7 +428,7 @@ void Foam::BlockMatrixClustering::calcClustering() ( groupPassI > 1 || indexGrouped == -1 - || sizeOfGroups[agglomIndex_[nextGrouped]] > maxClusterSize_ + || sizeOfGroups[agglomIndex_[nextGrouped]] > maxGroupSize_ ) { // There is no group to put this equation into @@ -765,8 +765,8 @@ Foam::BlockMatrixClustering::BlockMatrixClustering : BlockMatrixCoarsening(matrix, dict, groupSize, minCoarseEqns), matrix_(matrix), - minClusterSize_(readLabel(dict.lookup("minClusterSize"))), - maxClusterSize_(readLabel(dict.lookup("maxClusterSize"))), + minGroupSize_(readLabel(dict.lookup("minGroupSize"))), + maxGroupSize_(readLabel(dict.lookup("maxGroupSize"))), normPtr_(BlockCoeffNorm::New(dict)), agglomIndex_(matrix_.lduAddr().size()), groupSize_(groupSize), diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.H index 6e56ec25b..804d323c2 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.H @@ -66,17 +66,16 @@ class BlockMatrixClustering //- Reference to matrix const BlockLduMatrix& matrix_; - //- Min cluster size - const label minClusterSize_; + //- Min group size + const label minGroupSize_; - //- Max cluster size - const label maxClusterSize_; + //- Max group size + const label maxGroupSize_; //- Reference to a templated norm calculator autoPtr > normPtr_; - //- Child array: for each fine equation give clustering cluster - // index + //- Child array: for each fine equation give a clustering index labelField agglomIndex_; //- Group size diff --git a/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.C b/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.C index 8c3822966..73a2c7fea 100644 --- a/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.C +++ b/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.C @@ -264,7 +264,7 @@ void Foam::clusterAmgPolicy::calcChild() for ( groupPassI = 1; - groupPassI < minClusterSize_; + groupPassI < minGroupSize_; groupPassI++ ) { @@ -419,7 +419,7 @@ void Foam::clusterAmgPolicy::calcChild() ( groupPassI > 1 || indexGrouped == -1 - || sizeOfGroups[child_[nextGrouped]] > maxClusterSize_ + || sizeOfGroups[child_[nextGrouped]] > maxGroupSize_ ) { // There is no group to put this equation into @@ -501,8 +501,8 @@ Foam::clusterAmgPolicy::clusterAmgPolicy : amgPolicy(groupSize, minCoarseEqns), matrix_(matrix), - minClusterSize_(groupSize), - maxClusterSize_(2*groupSize), + minGroupSize_(groupSize), + maxGroupSize_(2*groupSize), child_(matrix_.lduAddr().size()), nSolo_(0), nCoarseEqns_(0), diff --git a/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.H b/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.H index 3ba6dd905..75dfcc2f6 100644 --- a/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.H +++ b/src/lduSolvers/amg/amgPolicy/clusterAmgPolicy.H @@ -61,10 +61,10 @@ class clusterAmgPolicy const lduMatrix& matrix_; //- Min cluster size - const label minClusterSize_; + const label minGroupSize_; //- Max cluster size - const label maxClusterSize_; + const label maxGroupSize_; //- Child array: for each fine equation give coarse cluster index labelField child_; From 7f020db73e03e6c86ea1a774a01f307f813acf66 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Mar 2017 16:09:47 +0100 Subject: [PATCH 025/333] Improved handling of solo cells in clustering --- .../BlockMatrixClustering.C | 53 ++++++++++++++++--- 1 file changed, 46 insertions(+), 7 deletions(-) diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C index 587904a74..681affb55 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C @@ -299,7 +299,15 @@ void Foam::BlockMatrixClustering::calcClustering() // Get column index, upper triangle colI = upperAddr[rowCoeffI]; - curWeight = normUpper[rowCoeffI]/normDiag[curEqn]; + // curWeight = normUpper[rowCoeffI]/ + // // normDiag[curEqn]; + // max(normDiag[curEqn], normDiag[colI]); + + curWeight = Foam::min + ( + normUpper[rowCoeffI]/normDiag[curEqn], + normLower[rowCoeffI]/normDiag[colI] + ); if (agglomIndex_[colI] == -1) { @@ -348,8 +356,15 @@ void Foam::BlockMatrixClustering::calcClustering() // Get column index, lower triangle colI = lowerAddr[losortAddr[rowCoeffI]]; - curWeight = normLower[losortAddr[rowCoeffI]]/ - normDiag[curEqn]; + // curWeight = normLower[losortAddr[rowCoeffI]]/ + // // normDiag[curEqn]; + // max(normDiag[curEqn], normDiag[colI]); + + curWeight = Foam::min + ( + normLower[losortAddr[rowCoeffI]]/normDiag[curEqn], + normUpper[losortAddr[rowCoeffI]]/normDiag[colI] + ); if (agglomIndex_[colI] == -1) { @@ -428,12 +443,36 @@ void Foam::BlockMatrixClustering::calcClustering() ( groupPassI > 1 || indexGrouped == -1 - || sizeOfGroups[agglomIndex_[nextGrouped]] > maxGroupSize_ + || sizeOfGroups[agglomIndex_[nextGrouped]] >= maxGroupSize_ ) { - // There is no group to put this equation into - sizeOfGroups[agglomIndex_[rowI]]++; - nCoarseEqns_++; + // If this is a solo cell and a group is available + // force it into the group irrespective of size + if + ( + sizeOfGroups[agglomIndex_[rowI]] == 0 + && indexGrouped != -1 + + ) + { + // Group exists, but it's too big. Add it anyway + agglomIndex_[rowI] = agglomIndex_[nextGrouped]; + sizeOfGroups[agglomIndex_[nextGrouped]]++; + } + else if (nSolo_ > 0) + { + // There is no group, but there is a solo group. + // Add it there + agglomIndex_[rowI] = 0; + sizeOfGroups[0]++; + nSolo_ ++; + } + else + { + // No group and no solo group. Make its own group + sizeOfGroups[agglomIndex_[rowI]]++; + nCoarseEqns_++; + } } else { From ef86f503f737b401df46ad07c324ab75e5d4b332 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 4 Apr 2017 10:36:16 +0100 Subject: [PATCH 026/333] Improved formulation of steady compressible solver suite --- .../steadyCompressibleFoam/UEqn.H | 1 - .../steadyCompressibleFoam/pEqn.H | 9 +--- .../steadyCompressibleMRFFoam/UEqn.H | 1 - .../steadyCompressibleMRFFoam/hEqn.H | 1 + .../steadyCompressibleMRFFoam/pEqn.H | 3 +- .../steadyCompressibleSRFFoam/UEqn.H | 1 - .../steadyCompressibleSRFFoam/pEqn.H | 5 +- .../compressible/steadyUniversalFoam/UEqn.H | 3 +- .../compressible/steadyUniversalFoam/hEqn.H | 5 +- .../compressible/steadyUniversalFoam/pEqn.H | 7 ++- .../steadyUniversalFoam/rhoFromP.H | 1 + .../universalContinuityErrs.H | 49 ------------------- .../steadyUniversalMRFFoam/UEqn.H | 3 +- .../steadyUniversalMRFFoam/iEqn.H | 23 +++++++-- .../steadyUniversalMRFFoam/pEqn.H | 7 ++- .../steadyUniversalMRFFoam/rhoFromP.H | 1 + .../universalContinuityErrs.H | 49 ------------------- 17 files changed, 39 insertions(+), 130 deletions(-) delete mode 100644 applications/solvers/compressible/steadyUniversalFoam/universalContinuityErrs.H delete mode 100644 applications/solvers/compressible/steadyUniversalMRFFoam/universalContinuityErrs.H diff --git a/applications/solvers/compressible/steadyCompressibleFoam/UEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/UEqn.H index 2833a63e9..3d8b68d58 100644 --- a/applications/solvers/compressible/steadyCompressibleFoam/UEqn.H +++ b/applications/solvers/compressible/steadyCompressibleFoam/UEqn.H @@ -4,7 +4,6 @@ fvVectorMatrix UEqn ( fvm::div(phi, U) - + fvm::SuSp(-fvc::div(phi), U) + turbulence->divDevRhoReff() ); diff --git a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H index f8fa0d40b..3262ebc32 100644 --- a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H @@ -21,12 +21,6 @@ p.storePrevIter(); - volScalarField divPhid - ( - "divPhid", - fvc::div(phid) - ); - while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn @@ -45,7 +39,8 @@ } } -# include "compressibleContinuityErrs.H" + // Use incompressible continuity error check: div(rho U) = 0 +# include "continuityErrs.H" // Relax the pressure p.relax(); diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H index 3e2f20ea0..f8408bf59 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H @@ -1,7 +1,6 @@ fvVectorMatrix UEqn ( fvm::div(phi, U) - + fvm::SuSp(-fvc::div(phi), U) + turbulence->divDevRhoReff() ); diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H index 3aebca26b..6cf655b45 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H @@ -24,6 +24,7 @@ ); hEqn.relax(); + hEqn.solve(); // Bounding of enthalpy taken out thermo.correct(); diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H index bedf65b19..21d7d1df3 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H @@ -47,7 +47,8 @@ } } -# include "compressibleContinuityErrs.H" + // Use incompressible continuity error check: div(rho U) = 0 +# include "continuityErrs.H" // Relax the pressure p.relax(); diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H index 9e09b5c42..6f9413765 100644 --- a/applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H +++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H @@ -4,7 +4,6 @@ fvVectorMatrix UrelEqn ( fvm::div(phi, Urel) - + fvm::SuSp(-fvc::div(phi), Urel) + turbulence->divDevRhoReff() + rho*SRF->Su() ); diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H index 762b129dc..eed8e6966 100644 --- a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H @@ -25,8 +25,8 @@ fvScalarMatrix pEqn ( fvm::div(phid, p) - - fvm::laplacian(rho*rUrelA, p) + fvc::div(phid2) + - fvm::laplacian(rho*rUrelA, p) ); pEqn.solve(); @@ -38,7 +38,8 @@ } } -# include "compressibleContinuityErrs.H" + // Use incompressible continuity error check: div(rho U) = 0 +# include "continuityErrs.H" // Relax the pressure p.relax(); diff --git a/applications/solvers/compressible/steadyUniversalFoam/UEqn.H b/applications/solvers/compressible/steadyUniversalFoam/UEqn.H index 10f11a421..3d8b68d58 100644 --- a/applications/solvers/compressible/steadyUniversalFoam/UEqn.H +++ b/applications/solvers/compressible/steadyUniversalFoam/UEqn.H @@ -3,8 +3,7 @@ fvVectorMatrix UEqn ( - fvm::ddt(rho, U) - + fvm::div(phi, U) + fvm::div(phi, U) + turbulence->divDevRhoReff() ); diff --git a/applications/solvers/compressible/steadyUniversalFoam/hEqn.H b/applications/solvers/compressible/steadyUniversalFoam/hEqn.H index 0a5c1e1db..d024e7df1 100644 --- a/applications/solvers/compressible/steadyUniversalFoam/hEqn.H +++ b/applications/solvers/compressible/steadyUniversalFoam/hEqn.H @@ -10,8 +10,8 @@ fvScalarMatrix hEqn ( - fvm::ddt(rho, h) - + fvm::div(phi, h) + fvm::div(phi, h) + + fvm::SuSp(-fvc::div(phi), h) - fvm::laplacian(turbulence->alphaEff(), h) == fvc::div(faceU, p, "div(U,p)") @@ -22,7 +22,6 @@ ); hEqn.relax(); - hEqn.solve(); // Bounding of enthalpy taken out diff --git a/applications/solvers/compressible/steadyUniversalFoam/pEqn.H b/applications/solvers/compressible/steadyUniversalFoam/pEqn.H index 8eec3318e..667bc89bb 100644 --- a/applications/solvers/compressible/steadyUniversalFoam/pEqn.H +++ b/applications/solvers/compressible/steadyUniversalFoam/pEqn.H @@ -24,8 +24,7 @@ { fvScalarMatrix pEqn ( - fvm::ddt(psis, p) - + fvm::div(phid, p) + fvm::div(phid, p) + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); @@ -39,8 +38,8 @@ } } - // Use custom continuity error check -# include "universalContinuityErrs.H" + // Use incompressible continuity error check: div(rho U) = 0 +# include "continuityErrs.H" // Relax the pressure p.relax(); diff --git a/applications/solvers/compressible/steadyUniversalFoam/rhoFromP.H b/applications/solvers/compressible/steadyUniversalFoam/rhoFromP.H index 31278008e..5dcfe3bf7 100644 --- a/applications/solvers/compressible/steadyUniversalFoam/rhoFromP.H +++ b/applications/solvers/compressible/steadyUniversalFoam/rhoFromP.H @@ -12,4 +12,5 @@ rho = Foam::min(rho, rhoMax); rho = Foam::max(rho, rhoMin); rho.relax(); + rho.correctBoundaryConditions(); } diff --git a/applications/solvers/compressible/steadyUniversalFoam/universalContinuityErrs.H b/applications/solvers/compressible/steadyUniversalFoam/universalContinuityErrs.H deleted file mode 100644 index 723f773e1..000000000 --- a/applications/solvers/compressible/steadyUniversalFoam/universalContinuityErrs.H +++ /dev/null @@ -1,49 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend 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 foam-extend. If not, see . - -Global - continuityErrs - -Description - Calculates and prints the continuity errors. - -\*---------------------------------------------------------------------------*/ - -{ - volScalarField contErr = fvc::ddt(rho) + fvc::div(phi); - - sumLocalContErr = runTime.deltaT().value()* - mag(contErr)().weightedAverage(mesh.V()).value(); - - globalContErr = runTime.deltaT().value()* - contErr.weightedAverage(mesh.V()).value(); - - cumulativeContErr += globalContErr; - - Info<< "time step continuity errors : sum local = " << sumLocalContErr - << ", global = " << globalContErr - << ", cumulative = " << cumulativeContErr - << endl; -} - -// ************************************************************************* // diff --git a/applications/solvers/compressible/steadyUniversalMRFFoam/UEqn.H b/applications/solvers/compressible/steadyUniversalMRFFoam/UEqn.H index 7f769de78..919631569 100644 --- a/applications/solvers/compressible/steadyUniversalMRFFoam/UEqn.H +++ b/applications/solvers/compressible/steadyUniversalMRFFoam/UEqn.H @@ -3,8 +3,7 @@ fvVectorMatrix UEqn ( - fvm::ddt(rho, U) - + fvm::div(phi, U) + fvm::div(phi, U) + turbulence->divDevRhoReff() ); diff --git a/applications/solvers/compressible/steadyUniversalMRFFoam/iEqn.H b/applications/solvers/compressible/steadyUniversalMRFFoam/iEqn.H index 915787b4b..1fc70ba27 100644 --- a/applications/solvers/compressible/steadyUniversalMRFFoam/iEqn.H +++ b/applications/solvers/compressible/steadyUniversalMRFFoam/iEqn.H @@ -3,14 +3,29 @@ Urel == U; mrfZones.relativeVelocity(Urel); + // Bound the relative velocity to preserve the i to h conversion bound + // HJ, 22/Mar/2017 + volScalarField magUrel = mag(Urel); + + if (max(magUrel) > UMax) + { + volScalarField UrelLimiter = pos(magUrel - UMax)*UMax/(magUrel + smallU) + + neg(magUrel - UMax); + UrelLimiter.max(scalar(0)); + UrelLimiter.min(scalar(1)); + + Urel *= UrelLimiter; + Urel.correctBoundaryConditions(); + } + // Create rotational velocity (= omega x r) Urot == U - Urel; fvScalarMatrix iEqn ( - fvm::ddt(rho, i) - + fvm::div(phi, i) - - fvm::laplacian(turbulence->alphaEff(), i) + fvm::div(phi, i) + + fvm::SuSp(-fvc::div(phi), i) + - fvm::laplacian(turbulence->alphaEff(), i) == // Viscous heating: note sign (devRhoReff has a minus in it) - (turbulence->devRhoReff() && fvc::grad(U)) @@ -23,8 +38,8 @@ // From rothalpy, calculate enthalpy after solution of rothalpy equation h = i + 0.5*(magSqr(Urot) - magSqr(Urel)); h.correctBoundaryConditions(); + // Update thermo for new h thermo.correct(); psis = thermo.psi()/thermo.Cp()*thermo.Cv(); } - diff --git a/applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H b/applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H index bba75ac7f..23973366c 100644 --- a/applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H @@ -33,8 +33,7 @@ { fvScalarMatrix pEqn ( - fvm::ddt(psis, p) - + fvm::div(phid, p) + fvm::div(phid, p) + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); @@ -48,8 +47,8 @@ } } - // Use custom continuity error check -# include "universalContinuityErrs.H" + // Use incompressible continuity error check: div(rho U) = 0 +# include "continuityErrs.H" // Relax the pressure p.relax(); diff --git a/applications/solvers/compressible/steadyUniversalMRFFoam/rhoFromP.H b/applications/solvers/compressible/steadyUniversalMRFFoam/rhoFromP.H index 31278008e..5dcfe3bf7 100644 --- a/applications/solvers/compressible/steadyUniversalMRFFoam/rhoFromP.H +++ b/applications/solvers/compressible/steadyUniversalMRFFoam/rhoFromP.H @@ -12,4 +12,5 @@ rho = Foam::min(rho, rhoMax); rho = Foam::max(rho, rhoMin); rho.relax(); + rho.correctBoundaryConditions(); } diff --git a/applications/solvers/compressible/steadyUniversalMRFFoam/universalContinuityErrs.H b/applications/solvers/compressible/steadyUniversalMRFFoam/universalContinuityErrs.H deleted file mode 100644 index 723f773e1..000000000 --- a/applications/solvers/compressible/steadyUniversalMRFFoam/universalContinuityErrs.H +++ /dev/null @@ -1,49 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend 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 foam-extend. If not, see . - -Global - continuityErrs - -Description - Calculates and prints the continuity errors. - -\*---------------------------------------------------------------------------*/ - -{ - volScalarField contErr = fvc::ddt(rho) + fvc::div(phi); - - sumLocalContErr = runTime.deltaT().value()* - mag(contErr)().weightedAverage(mesh.V()).value(); - - globalContErr = runTime.deltaT().value()* - contErr.weightedAverage(mesh.V()).value(); - - cumulativeContErr += globalContErr; - - Info<< "time step continuity errors : sum local = " << sumLocalContErr - << ", global = " << globalContErr - << ", cumulative = " << cumulativeContErr - << endl; -} - -// ************************************************************************* // From e74a1c666ec948fb872e3b8858ce07ad08706a73 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 4 Apr 2017 16:34:43 +0100 Subject: [PATCH 027/333] Temporary compilation fix --- .../surfaceTracking/interTrackFoam/interTrackFoam.C | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/applications/solvers/surfaceTracking/interTrackFoam/interTrackFoam.C b/applications/solvers/surfaceTracking/interTrackFoam/interTrackFoam.C index 0b9e62069..d91ac9242 100644 --- a/applications/solvers/surfaceTracking/interTrackFoam/interTrackFoam.C +++ b/applications/solvers/surfaceTracking/interTrackFoam/interTrackFoam.C @@ -41,16 +41,16 @@ Description int main(int argc, char *argv[]) { # include "setRootCase.H" - # include "createTime.H" - # include "createDynamicFvMesh.H" pimpleControl pimple(mesh); # include "createFields.H" - +# include "createTimeControls.H" # include "initContinuityErrs.H" +# include "CourantNo.H" +# include "setInitialDeltaT.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -60,7 +60,7 @@ int main(int argc, char *argv[]) { # include "readTimeControls.H" # include "CourantNo.H" -# include "setSurfaceStabilityDeltaT.H" + //# include "setSurfaceStabilityDeltaT.H" # include "readPISOControls.H" From 3968a5c9005bc29299efe03c7c2b3e95765a69df Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 4 Apr 2017 16:35:14 +0100 Subject: [PATCH 028/333] Added support for block-coupled coefficients in fvPatchFields --- .../fvPatchFields/fvPatchField/fvPatchField.H | 163 +++++++++++++----- 1 file changed, 123 insertions(+), 40 deletions(-) diff --git a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H index daed7adcb..5d4a7b2c8 100644 --- a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H @@ -46,6 +46,7 @@ SourceFiles #include "fvPatch.H" #include "DimensionedField.H" +#include "coeffFields.H" #include "debugSwitch.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -329,6 +330,14 @@ public: return false; } + //- Return true if this patch field is block-coupled in the + // boundary condition, ie. if the coupling coefficient is a + // rank x rank implicit block + virtual bool blockCoupled() const + { + return false; + } + //- Return true if the boundary condition has already been updated bool updated() const { @@ -388,51 +397,125 @@ public: ); - //- Return the matrix diagonal coefficients corresponding to the - // evaluation of the value of this patchField with given weights - virtual tmp > valueInternalCoeffs - ( - const tmp >& - ) const - { - notImplemented + // Decoupled matrix coefficients + + //- Return the matrix diagonal coefficients corresponding to the + // evaluation of the value of this patchField with + // given weights + virtual tmp > valueInternalCoeffs ( - type() - + "::valueInternalCoeffs(const tmp >&)" - ); - return *this; - } + const tmp >& + ) const + { + notImplemented + ( + type() + + "::valueInternalCoeffs(const tmp >&)" + ); - //- Return the matrix source coefficients corresponding to the - // evaluation of the value of this patchField with given weights - virtual tmp > valueBoundaryCoeffs - ( - const tmp >& - ) const - { - notImplemented + return *this; + } + + //- Return the matrix source coefficients corresponding to the + // evaluation of the value of this patchField with given + // weights + virtual tmp > valueBoundaryCoeffs ( - type() - + "::valueBoundaryCoeffs(const tmp >&)" - ); - return *this; - } + const tmp >& + ) const + { + notImplemented + ( + type() + + "::valueBoundaryCoeffs(const tmp >&)" + ); + return *this; + } - //- Return the matrix diagonal coefficients corresponding to the - // evaluation of the gradient of this patchField - virtual tmp > gradientInternalCoeffs() const - { - notImplemented(type() + "::gradientInternalCoeffs()"); - return *this; - } + //- Return the matrix diagonal coefficients corresponding to the + // evaluation of the gradient of this patchField + virtual tmp > gradientInternalCoeffs() const + { + notImplemented(type() + "::gradientInternalCoeffs()"); + return *this; + } - //- Return the matrix source coefficients corresponding to the - // evaluation of the gradient of this patchField - virtual tmp > gradientBoundaryCoeffs() const - { - notImplemented(type() + "::gradientBoundaryCoeffs()"); - return *this; - } + //- Return the matrix source coefficients corresponding to the + // evaluation of the gradient of this patchField + virtual tmp > gradientBoundaryCoeffs() const + { + notImplemented(type() + "::gradientBoundaryCoeffs()"); + return *this; + } + + + // Block-coupled matrix coefficients + + //- Return the matrix diagonal coefficients corresponding to the + // evaluation of the value of this patchField with given + // weights + virtual tmp > blockValueInternalCoeffs + ( + const tmp& + ) const + { + notImplemented + ( + type() + + "::blockValueInternalCoeffs(const tmp >&)" + ); + + return tmp > + ( + new CoeffField(this->size()) + ); + } + + //- Return the matrix source coefficients corresponding to the + // evaluation of the value of this patchField with given + // weights + virtual tmp > blockValueBoundaryCoeffs + ( + const tmp& + ) const + { + notImplemented + ( + type() + + "::blockValueBoundaryCoeffs(const tmp >&)" + ); + + return *this; + } + + //- Return the matrix diagonal coefficients corresponding to the + // evaluation of the gradient of this patchField + virtual tmp > + blockGradientInternalCoeffs() const + { + notImplemented + ( + type() + "blockGradientInternalCoeffs() const" + ); + + return tmp > + ( + new CoeffField(this->size()) + ); + } + + //- Return the matrix source coefficients corresponding to the + // evaluation of the gradient of this patchField + virtual tmp > + blockGradientBoundaryCoeffs() const + { + notImplemented + ( + type() + "::blockGradientBoundaryCoeffs()" + ); + + return *this; + } // Patch interpolation and patch flux From 4b63cb70698911527d9097f53c24232c2b7fb1e3 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 5 Apr 2017 12:06:44 +0200 Subject: [PATCH 029/333] Switch off debug messages --- .../BlockMatrixAgglomeration/BlockMatrixAgglomeration.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.C index bf0c934aa..ba897f61d 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.C +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.C @@ -366,7 +366,7 @@ void Foam::BlockMatrixAgglomeration::calcAgglomeration() reduce(coarsen_, andOp()); - // if (BlockLduMatrix::debug >= 3) + if (BlockLduMatrix::debug >= 3) { // Count singleton clusters label nSingleClusters = 0; From 9fca371a485dab042655f0a0d276dffb4c443dbd Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 5 Apr 2017 12:06:59 +0200 Subject: [PATCH 030/333] Change handling of zero cluster --- .../BlockMatrixClustering/BlockMatrixClustering.C | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C index 681affb55..c07cf753c 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C @@ -459,14 +459,6 @@ void Foam::BlockMatrixClustering::calcClustering() agglomIndex_[rowI] = agglomIndex_[nextGrouped]; sizeOfGroups[agglomIndex_[nextGrouped]]++; } - else if (nSolo_ > 0) - { - // There is no group, but there is a solo group. - // Add it there - agglomIndex_[rowI] = 0; - sizeOfGroups[0]++; - nSolo_ ++; - } else { // No group and no solo group. Make its own group @@ -503,7 +495,7 @@ void Foam::BlockMatrixClustering::calcClustering() reduce(coarsen_, andOp()); - // if (BlockLduMatrix::debug >= 3) + if (BlockLduMatrix::debug >= 3) { // Count singleton clusters label nSingleClusters = 0; From 0e04d88628a9c899cc33e9bdc45446bb4dad99f3 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 5 Apr 2017 11:08:42 +0100 Subject: [PATCH 031/333] Added zero CoeffField --- src/foam/fields/CoeffField/CoeffField.C | 4 ++++ src/foam/fields/CoeffField/CoeffField.H | 24 ++++++++++++++++++------ 2 files changed, 22 insertions(+), 6 deletions(-) diff --git a/src/foam/fields/CoeffField/CoeffField.C b/src/foam/fields/CoeffField/CoeffField.C index 6fb7a6201..ea3233901 100644 --- a/src/foam/fields/CoeffField/CoeffField.C +++ b/src/foam/fields/CoeffField/CoeffField.C @@ -31,6 +31,10 @@ License template const char* const Foam::CoeffField::typeName("CoeffField"); +template +const Foam::CoeffField Foam::CoeffField::zero(0); + + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // diff --git a/src/foam/fields/CoeffField/CoeffField.H b/src/foam/fields/CoeffField/CoeffField.H index 159e73c38..a46ee682d 100644 --- a/src/foam/fields/CoeffField/CoeffField.H +++ b/src/foam/fields/CoeffField/CoeffField.H @@ -122,6 +122,19 @@ public: static const char* const typeName; + //- Empty field + static const CoeffField zero; + + + // Static Member Functions + + //- Return a null field + inline static const CoeffField& null() + { + return zero; + } + + // Constructors @@ -138,12 +151,8 @@ public: tmp > clone() const; - // Destructor - - ~CoeffField(); - - //- Clear data - void clear(); + //- Destructor + ~CoeffField(); // Member functions @@ -163,6 +172,9 @@ public: //- Return the field transpose tmp > transpose() const; + //- Clear data + void clear(); + // Return as typed. Fails when asked for the incorrect type From ed73abd931b91128f971a11f4d4e5bd8639fcfd3 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 10 Apr 2017 11:19:39 +0100 Subject: [PATCH 032/333] Improved form of BlockCoeffTwoNorm --- .../BlockCoeffTwoNorm/BlockCoeffTwoNorm.C | 20 +++-- .../BlockCoeffTwoNorm/blockCoeffTwoNorms.H | 1 + .../scalarBlockCoeffTwoNorm.H | 80 +++++++++++++++++++ .../tensorBlockCoeffTwoNorm.H | 11 --- 4 files changed, 96 insertions(+), 16 deletions(-) create mode 100644 src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/scalarBlockCoeffTwoNorm.H diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/BlockCoeffTwoNorm.C b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/BlockCoeffTwoNorm.C index dfdb1f025..1dc865113 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/BlockCoeffTwoNorm.C +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/BlockCoeffTwoNorm.C @@ -61,8 +61,13 @@ Foam::scalar Foam::BlockCoeffTwoNorm::normalize } else if (a.activeType() == Foam::BlockCoeff::SQUARE) { - // Note: for two-norm, use mag - return mag(a.asSquare()); + // Note: for two-norm, use mag of diag + typedef typename BlockCoeff::linearType linearType; + + linearType diag; + contractLinear(diag, a.asSquare()); + + return mag(diag); } else { @@ -94,13 +99,18 @@ void Foam::BlockCoeffTwoNorm::normalize } else if (a.activeType() == Foam::BlockCoeff::LINEAR) { - // Note: for two-norm, use mag + // Note: for two-norm, use mag = sqrt(sum magSqr(cmpt)) b = mag(a.asLinear()); } else if (a.activeType() == Foam::BlockCoeff::SQUARE) { - // Note: for two-norm, use mag - b = mag(a.asSquare()); + // Note: for two-norm, use mag of diag + typedef typename BlockCoeff::linearTypeField linearTypeField; + + linearTypeField diag(a.size()); + contractLinear(diag, a.asSquare()); + + b = mag(diag); } else { diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/blockCoeffTwoNorms.H b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/blockCoeffTwoNorms.H index fd70487b3..524766733 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/blockCoeffTwoNorms.H +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/blockCoeffTwoNorms.H @@ -38,6 +38,7 @@ SourceFiles #ifndef blockCoeffTwoNorms_H #define blockCoeffTwoNorms_H +#include "scalarBlockCoeffTwoNorm.H" #include "tensorBlockCoeffTwoNorm.H" #include "BlockCoeffTwoNorm.H" diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/scalarBlockCoeffTwoNorm.H b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/scalarBlockCoeffTwoNorm.H new file mode 100644 index 000000000..3a10881b6 --- /dev/null +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/scalarBlockCoeffTwoNorm.H @@ -0,0 +1,80 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend 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 foam-extend. If not, see . + +Class + BlockCoeffTwoNorm + +Description + Class for two norm of block coeffs. Specilization for scalar. Implemented + to avoid issues with asScalar, asSquare etc. + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#ifndef scalarBlockCoeffTwoNorm_H +#define scalarBlockCoeffTwoNorm_H + +#include "blockCoeffs.H" +#include "blockCoeffNorms.H" +#include "BlockCoeffNorm.H" +#include "BlockCoeffTwoNorm.H" +#include "runTimeSelectionTables.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +inline scalar BlockCoeffTwoNorm::normalize +( + const BlockCoeff& a +) +{ + // Note: for two-norm, use mag + return mag(a.asScalar()); +} + + +template<> +inline void BlockCoeffTwoNorm::normalize +( + Field& b, + const CoeffField& a +) +{ + b = mag(a.asScalar()); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/tensorBlockCoeffTwoNorm.H b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/tensorBlockCoeffTwoNorm.H index db2a2bc7d..e765a5a65 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/tensorBlockCoeffTwoNorm.H +++ b/src/foam/primitives/BlockCoeff/BlockCoeffNorm/BlockCoeffTwoNorm/tensorBlockCoeffTwoNorm.H @@ -33,9 +33,6 @@ Description Author Klas Jareteg, 2013-01-30 -SourceFiles - BlockCoeffTwoNorm.C - \*---------------------------------------------------------------------------*/ #ifndef tensorBlockCoeffTwoNorm_H @@ -52,11 +49,6 @@ SourceFiles namespace Foam { - -/*---------------------------------------------------------------------------*\ - Class BlockCoeffTwoNorm Declaration -\*---------------------------------------------------------------------------*/ - template<> inline scalar BlockCoeffTwoNorm::normalize ( @@ -116,9 +108,6 @@ inline void BlockCoeffTwoNorm::normalize } // End namespace Foam -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif From 568a927388f1f1d7d571e5cf4fddba0d8776807f Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 10 Apr 2017 11:42:59 +0100 Subject: [PATCH 033/333] Better error messaging for coupled patches --- .../fields/fvPatchFields/constraint/ggi/ggiFvPatchField.H | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.H index 315e2d1e0..85c9079e6 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.H @@ -198,11 +198,13 @@ public: Field& f ) const { - notImplemented + FatalErrorIn ( "ggiFvPatchField::transformCoupleField" "(Field& f) const" - ); + ) << "transformCoupleField not implemented for patch " + << this->patch().name() + << abort(FatalError); } //- Initialise neighbour matrix update From 1d1bc4dcf39d109612d03efda66ba0de5098c0c4 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 10 Apr 2017 11:43:38 +0100 Subject: [PATCH 034/333] Changes in master/slave virtual mamber function for coupled patches --- .../manipulation/createPatch/createPatch.C | 4 +- .../processorMeshesRebuild.C | 4 +- .../basic/coupled/coupledPolyPatch.H | 14 +- .../constraint/cyclic/cyclicPolyPatch.H | 11 +- .../polyPatches/constraint/ggi/ggiPolyPatch.H | 8 +- .../mixingPlane/mixingPlanePolyPatch.H | 8 +- .../overlapGgi/overlapGgiPolyPatch.C | 13 -- .../overlapGgi/overlapGgiPolyPatch.H | 158 ++++++++++-------- .../constraint/processor/processorPolyPatch.C | 6 +- .../constraint/processor/processorPolyPatch.H | 10 +- .../regionCouple/regionCouplePolyPatch.H | 8 +- .../meshes/polyMesh/syncTools/syncTools.C | 6 +- .../field/fieldValues/faceSource/faceSource.C | 2 +- 13 files changed, 119 insertions(+), 133 deletions(-) diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatch.C b/applications/utilities/mesh/manipulation/createPatch/createPatch.C index b3d567b7e..246893b53 100644 --- a/applications/utilities/mesh/manipulation/createPatch/createPatch.C +++ b/applications/utilities/mesh/manipulation/createPatch/createPatch.C @@ -343,7 +343,7 @@ void syncPoints ( isA(pp) && pp.nPoints() > 0 - && refCast(pp).owner() + && refCast(pp).master() ) { const processorPolyPatch& procPatch = @@ -380,7 +380,7 @@ void syncPoints ( isA(pp) && pp.nPoints() > 0 - && !refCast(pp).owner() + && !refCast(pp).master() ) { const processorPolyPatch& procPatch = diff --git a/applications/utilities/parallelProcessing/reconstructParMesh/processorMeshesRebuild.C b/applications/utilities/parallelProcessing/reconstructParMesh/processorMeshesRebuild.C index 2be72b483..8d66ceb4b 100644 --- a/applications/utilities/parallelProcessing/reconstructParMesh/processorMeshesRebuild.C +++ b/applications/utilities/parallelProcessing/reconstructParMesh/processorMeshesRebuild.C @@ -735,7 +735,7 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db) // If patch is a neighbour, its master has already inserted // the points - if (procPatch.neighbour()) + if (procPatch.slave()) { const label masterProcID = procPatch.neighbProcNo(); @@ -870,7 +870,7 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db) // If patch is a master, drop the faces and fill the // owner side addressing - if (procPatch.owner()) + if (procPatch.master()) { for ( diff --git a/src/foam/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H b/src/foam/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H index 4b5ec4f07..867d26196 100644 --- a/src/foam/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H +++ b/src/foam/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H @@ -146,9 +146,8 @@ public: ); - // Destructor - - virtual ~coupledPolyPatch(); + //- Destructor + virtual ~coupledPolyPatch(); // Member Functions @@ -161,6 +160,15 @@ public: return true; } + //- Does this side own the patch ? + virtual bool master() const = 0; + + //- Does the coupled side own the patch ? + virtual bool slave() const + { + return !master(); + } + //- Are the coupled planes separated bool separated() const diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H b/src/foam/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H index c1792453a..fb37812c8 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H @@ -268,13 +268,18 @@ public: } - // Destructor - - virtual ~cyclicPolyPatch(); + //- Destructor + virtual ~cyclicPolyPatch(); // Member Functions + //- Is this the master side? Yes: it contains both sets of faces + virtual bool master() const + { + return true; + } + //- Return connected points (in patch local point indexing). Demand // driven calculation. Does primitivePatch::clearOut // after calculation! diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H index 8e585b6dc..065a7507a 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H @@ -322,17 +322,11 @@ public: // Interpolation data //- Is this the master side? - bool master() const + virtual bool master() const { return index() < shadowIndex(); } - //- Is this the slave side? - bool slave() const - { - return !master(); - } - //- Use bridging to fix overlap error in interpolation bool bridgeOverlap() const { diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/mixingPlane/mixingPlanePolyPatch.H b/src/foam/meshes/polyMesh/polyPatches/constraint/mixingPlane/mixingPlanePolyPatch.H index 58173e15d..535f04e61 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/mixingPlane/mixingPlanePolyPatch.H +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/mixingPlane/mixingPlanePolyPatch.H @@ -314,17 +314,11 @@ public: // Interpolation data //- Is this the master side? - bool master() const + virtual bool master() const { return index() < shadowIndex(); } - //- Is this the slave side? - bool slave() const - { - return !master(); - } - //- Return coordinate system const coordinateSystem& cs() const; diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatch.C b/src/foam/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatch.C index a174ffa7b..e70be8859 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatch.C +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatch.C @@ -360,19 +360,6 @@ const Foam::label& Foam::overlapGgiPolyPatch::nCopies() const } -bool Foam::overlapGgiPolyPatch::master() const -{ - // The first overlapggi interface is master,second one is slave - if (angle() == shadow().angle()) - { - return index() < shadowIndex(); - } - - // Master is the one with the larger angle - return angle() > shadow().angle(); -} - - // Write void Foam::overlapGgiPolyPatch::write(Ostream& os) const { diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatch.H b/src/foam/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatch.H index 183a204cd..142eb9a64 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatch.H +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatch.H @@ -265,103 +265,113 @@ public: } - // Destructor - - virtual ~overlapGgiPolyPatch(); + //- Destructor + virtual ~overlapGgiPolyPatch(); // Member functions - //- Return shadow patch name - const word& shadowName() const - { - return shadowName_; - } + // Basic info - //- Return name of interpolation face zone - const word& zoneName() const - { - return zoneName_; - } + //- Return shadow patch name + const word& shadowName() const + { + return shadowName_; + } - //- Return shadow patch index - label shadowIndex() const; + //- Return name of interpolation face zone + const word& zoneName() const + { + return zoneName_; + } - //- Return zone patch index - label zoneIndex() const; + //- Return shadow patch index + label shadowIndex() const; - //- Return shadow patch - const overlapGgiPolyPatch& shadow() const; + //- Return zone patch index + label zoneIndex() const; - //- Return interpolation face zone - const faceZone& zone() const; + //- Return shadow patch + const overlapGgiPolyPatch& shadow() const; - //- Return rotation axis - const vector& rotationAxis() const - { - return rotationAxis_; - } + //- Return interpolation face zone + const faceZone& zone() const; - //- Return wedge angle - scalar angle() const - { - return 360.0/scalar(nCopies()); - } + //- Return rotation axis + const vector& rotationAxis() const + { + return rotationAxis_; + } - //- Return number of slave copies - const label& nCopies() const; + //- Return wedge angle + scalar angle() const + { + return 360.0/scalar(nCopies()); + } - //- Is this the master side? - bool master() const; + //- Return number of slave copies + const label& nCopies() const; - //- Is this the slave side? - bool slave() const - { - return !master(); - } - //- Is the patch localised on a single processor - bool localParallel() const; + // Interpolation data - //- Expand face field to full for 360 degrees coverage - template - tmp > expandData(const Field& spf) const; + //- Is this the master side? + virtual bool master() const + { + return index() < shadowIndex(); + } - //- Interpolate face field: given field on the shadow side, - // create an interpolated field on this side - template - tmp > interpolate(const Field& pf) const; + //- Is the patch localised on a single processor + bool localParallel() const; - template - tmp > interpolate(const tmp >& tpf) const; - //- Filter zone field to patch size - template - tmp > filter(const Field& zf) const; + // Interpolation functions - //- Return reconstructed cell centres - const vectorField& reconFaceCellCentres() const; + //- Expand face field to full for 360 degrees coverage + template + tmp > expandData(const Field& spf) const; - //- Initialize ordering for primitivePatch. Does not - // refer to *this (except for name() and type() etc.) - virtual void initOrder(const primitivePatch&) const; + //- Interpolate face field: given field on the shadow side, + // create an interpolated field on this side + template + tmp > interpolate(const Field& pf) const; - //- Return new ordering for primitivePatch. - // Ordering is -faceMap: for every face - // index of the new face -rotation: for every new face the clockwise - // shift of the original face. Return false if nothing changes - // (faceMap is identity, rotation is 0), true otherwise. - virtual bool order - ( - const primitivePatch&, - labelList& faceMap, - labelList& rotation - ) const; + template + tmp > interpolate(const tmp >& tpf) const; - //- Synchronise communications of ordering for primitivePatch - // Used in cases when no topological change happens locally, - // but is happening on other processors - virtual void syncOrder() const; + //- Filter zone field to patch size + template + tmp > filter(const Field& zf) const; + + + // Geometric data + + //- Return reconstructed cell centres + const vectorField& reconFaceCellCentres() const; + + + // Patch ordering + + //- Initialize ordering for primitivePatch. Does not + // refer to *this (except for name() and type() etc.) + virtual void initOrder(const primitivePatch&) const; + + //- Return new ordering for primitivePatch. + // Ordering is -faceMap: for every face + // index of the new face -rotation: for every new face the + // clockwise shift of the original face. Return false if nothing + // changes (faceMap is identity, rotation is 0), true otherwise. + virtual bool order + ( + const primitivePatch&, + labelList& faceMap, + labelList& rotation + ) const; + + //- Synchronise communications of ordering for primitivePatch + // Used in cases when no topological change happens locally, + // but is happening on other processors + virtual void syncOrder() const; //- Write diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C b/src/foam/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C index 48b1e9397..b353039c5 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C @@ -471,7 +471,7 @@ void Foam::processorPolyPatch::initOrder(const primitivePatch& pp) const // Master side sends the data and slave side rotates the faces // HJ, 10/Mar/2011 - if (owner()) + if (master()) { pointField ctrs(calcFaceCentres(pp, pp.points())); @@ -538,7 +538,7 @@ bool Foam::processorPolyPatch::order rotation.setSize(pp.size()); rotation = 0; - if (owner()) + if (master()) { // Do nothing (i.e. identical mapping, zero rotation). // See comment at top. @@ -807,7 +807,7 @@ void Foam::processorPolyPatch::syncOrder() const } // Read and discard info from owner side from the neighbour - if (neighbour()) + if (slave()) { vectorField masterCtrs; vectorField masterAnchors; diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H b/src/foam/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H index f90a286d1..79f0c9d79 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H @@ -217,18 +217,12 @@ public: virtual int tag() const; - //- Does the processor own the patch ? - bool owner() const + //- Is this the master side? + virtual bool master() const { return (myProcNo_ < neighbProcNo_); } - //- Is the processor the patch neighbour ? - bool neighbour() const - { - return !owner(); - } - //- Return processor-neighbbour patch face centres const vectorField& neighbFaceCentres() const { diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H b/src/foam/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H index b6cb05589..f0f6fa420 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H @@ -363,17 +363,11 @@ public: // Interpolation data //- Is this the master side? - bool master() const + virtual bool master() const { return master_; } - //- Is this the slave side? - bool slave() const - { - return !master(); - } - //- Use bridging to fix overlap error in interpolation bool bridgeOverlap() const { diff --git a/src/foam/meshes/polyMesh/syncTools/syncTools.C b/src/foam/meshes/polyMesh/syncTools/syncTools.C index 646e28e9b..d3e0c9598 100644 --- a/src/foam/meshes/polyMesh/syncTools/syncTools.C +++ b/src/foam/meshes/polyMesh/syncTools/syncTools.C @@ -132,7 +132,7 @@ Foam::PackedBoolList Foam::syncTools::getMasterPoints(const polyMesh& mesh) { donePoint.set(pointI, 1u); - if (pp.owner()) + if (pp.master()) { isMasterPoint.set(pointI, 1u); } @@ -251,7 +251,7 @@ Foam::PackedBoolList Foam::syncTools::getMasterEdges(const polyMesh& mesh) { doneEdge.set(edgeI, 1u); - if (pp.owner()) + if (pp.master()) { isMasterEdge.set(edgeI, 1u); } @@ -325,7 +325,7 @@ Foam::PackedBoolList Foam::syncTools::getMasterFaces(const polyMesh& mesh) const processorPolyPatch& pp = refCast(patches[patchI]); - if (!pp.owner()) + if (!pp.master()) { forAll(pp, i) { diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C index f614617d8..237ae1420 100644 --- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C +++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C @@ -102,7 +102,7 @@ void Foam::fieldValues::faceSource::setFaceZoneFaces() const polyPatch& pp = mesh().boundaryMesh()[facePatchId]; if (isA(pp)) { - if (refCast(pp).owner()) + if (refCast(pp).master()) { faceId = pp.whichFace(faceI); } From 21d8588cc88ba0a12efa1512fdbc2c183903875e Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 10 Apr 2017 11:44:58 +0100 Subject: [PATCH 035/333] Added view factor radiation model and tools --- .../preProcessing/faceAgglomerate/Make/files | 3 + .../faceAgglomerate/Make/options | 7 + .../faceAgglomerate/faceAgglomerate.C | 188 ++++ .../preProcessing/viewFactorsGen/Make/files | 3 + .../preProcessing/viewFactorsGen/Make/options | 7 + .../viewFactorsGen/searchingEngine.H | 54 ++ .../preProcessing/viewFactorsGen/shootRays.H | 83 ++ .../viewFactorsGen/viewFactorsGen.C | 851 ++++++++++++++++++ .../polyMesh/globalMeshData/globalIndex.H | 7 +- .../mapPolyMesh/mapDistribute/mapDistribute.C | 458 +++++++++- .../mapPolyMesh/mapDistribute/mapDistribute.H | 74 ++ src/fvAgglomerationMethods/Allwmake | 1 + .../pairPatchAgglomeration/Make/files | 3 + .../pairPatchAgglomeration/Make/options | 0 .../pairPatchAgglomeration.C | 533 +++++++++++ .../pairPatchAgglomeration.H | 217 +++++ .../pairPatchAgglomerationTemplates.C | 79 ++ 17 files changed, 2545 insertions(+), 23 deletions(-) create mode 100644 applications/utilities/preProcessing/faceAgglomerate/Make/files create mode 100644 applications/utilities/preProcessing/faceAgglomerate/Make/options create mode 100644 applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C create mode 100644 applications/utilities/preProcessing/viewFactorsGen/Make/files create mode 100644 applications/utilities/preProcessing/viewFactorsGen/Make/options create mode 100644 applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H create mode 100644 applications/utilities/preProcessing/viewFactorsGen/shootRays.H create mode 100644 applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C create mode 100644 src/fvAgglomerationMethods/pairPatchAgglomeration/Make/files create mode 100644 src/fvAgglomerationMethods/pairPatchAgglomeration/Make/options create mode 100644 src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C create mode 100644 src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.H create mode 100644 src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomerationTemplates.C diff --git a/applications/utilities/preProcessing/faceAgglomerate/Make/files b/applications/utilities/preProcessing/faceAgglomerate/Make/files new file mode 100644 index 000000000..34bd95a84 --- /dev/null +++ b/applications/utilities/preProcessing/faceAgglomerate/Make/files @@ -0,0 +1,3 @@ +faceAgglomerate.C + +EXE = $(FOAM_APPBIN)/faceAgglomerate diff --git a/applications/utilities/preProcessing/faceAgglomerate/Make/options b/applications/utilities/preProcessing/faceAgglomerate/Make/options new file mode 100644 index 000000000..3d9b744f2 --- /dev/null +++ b/applications/utilities/preProcessing/faceAgglomerate/Make/options @@ -0,0 +1,7 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/fvAgglomerationMethods/pairPatchAgglomeration/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lpairPatchAgglomeration diff --git a/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C b/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C new file mode 100644 index 000000000..c59e7a924 --- /dev/null +++ b/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C @@ -0,0 +1,188 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend 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 foam-extend. If not, see . + +Application + faceAgglomerate + +Description + Agglomerate boundary faces using the pairPatchAgglomeration algorithm. + It writes a map from the fine to coarse grid. + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "foamTime.H" +#include "fvMesh.H" +#include "volFields.H" +#include "unitConversion.H" +#include "pairPatchAgglomeration.H" +#include "labelIOList.H" +#include "syncTools.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ +# include "addRegionOption.H" + +# include "setRootCase.H" +# include "createTime.H" +# include "createNamedMesh.H" + + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + + labelListIOList finalAgglom + ( + IOobject + ( + "finalAgglom", + mesh.facesInstance(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + patches.size() + ); + + + // Read control dictionary + IOdictionary agglomDict + ( + IOobject + ( + "faceAgglomerateDict", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + bool writeAgglom = readBool(agglomDict.lookup("writeFacesAgglomeration")); + + const polyBoundaryMesh& boundary = mesh.boundaryMesh(); + + forAll(boundary, patchId) + { + const polyPatch& pp = boundary[patchId]; + + label patchI = pp.index(); + finalAgglom[patchI].setSize(pp.size(), 0); + + if (!pp.coupled()) + { + if (agglomDict.found(pp.name())) + { + Info << "\nAgglomerating patch : " << pp.name() << endl; + pairPatchAgglomeration agglomObject + ( + pp, + agglomDict.subDict(pp.name()) + ); + agglomObject.agglomerate(); + finalAgglom[patchI] = + agglomObject.restrictTopBottomAddressing(); + } + else + { + FatalErrorIn(args.executable()) + << "Patch " << pp.name() << " not found in dictionary: " + << agglomDict.name() << exit(FatalError); + } + } + } + + // Sync agglomeration across coupled patches + labelList nbrAgglom(mesh.nFaces() - mesh.nInternalFaces(), -1); + + forAll(boundary, patchId) + { + const polyPatch& pp = boundary[patchId]; + if (pp.coupled()) + { + finalAgglom[patchId] = identity(pp.size()); + forAll(pp, i) + { + nbrAgglom[pp.start() - mesh.nInternalFaces() + i] = + finalAgglom[patchId][i]; + } + } + } + + syncTools::swapBoundaryFaceList(mesh, nbrAgglom, false); + forAll(boundary, patchId) + { + const polyPatch& pp = boundary[patchId]; + if (pp.coupled() && !refCast(pp).master()) + { + forAll(pp, i) + { + finalAgglom[patchId][i] = + nbrAgglom[pp.start() - mesh.nInternalFaces() + i]; + } + } + } + + finalAgglom.write(); + + if (writeAgglom) + { + volScalarField facesAgglomeration + ( + IOobject + ( + "facesAgglomeration", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("facesAgglomeration", dimless, 0) + ); + + forAll(boundary, patchId) + { + + fvPatchScalarField& bFacesAgglomeration = + facesAgglomeration.boundaryField()[patchId]; + + forAll(bFacesAgglomeration, j) + { + bFacesAgglomeration[j] = finalAgglom[patchId][j]; + } + } + + Info << "\nWriting facesAgglomeration" << endl; + facesAgglomeration.write(); + } + + Info<< "End\n" << endl; + return 0; +} + + +// ************************************************************************* // diff --git a/applications/utilities/preProcessing/viewFactorsGen/Make/files b/applications/utilities/preProcessing/viewFactorsGen/Make/files new file mode 100644 index 000000000..aa0bc4481 --- /dev/null +++ b/applications/utilities/preProcessing/viewFactorsGen/Make/files @@ -0,0 +1,3 @@ +viewFactorsGen.C + +EXE = $(FOAM_APPBIN)/viewFactorsGen diff --git a/applications/utilities/preProcessing/viewFactorsGen/Make/options b/applications/utilities/preProcessing/viewFactorsGen/Make/options new file mode 100644 index 000000000..d27c95d03 --- /dev/null +++ b/applications/utilities/preProcessing/viewFactorsGen/Make/options @@ -0,0 +1,7 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools diff --git a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H new file mode 100644 index 000000000..fad039b3b --- /dev/null +++ b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H @@ -0,0 +1,54 @@ +Random rndGen(653213); + +// Determine mesh bounding boxes: +List meshBb +( + 1, + treeBoundBox + ( + boundBox(coarseMesh.points(), false) + ).extend(rndGen, 1E-3) +); + +// Dummy bounds dictionary +dictionary dict; +dict.add("bounds", meshBb); +dict.add +( + "distributionType", + distributedTriSurfaceMesh::distributionTypeNames_ + [ + distributedTriSurfaceMesh::FROZEN + ] +); +dict.add("mergeDistance", SMALL); + +labelHashSet includePatches; +forAll(patches, patchI) +{ + if (!isA(patches[patchI])) + { + includePatches.insert(patchI); + } +} + +distributedTriSurfaceMesh surfacesMesh +( + IOobject + ( + "wallSurface.stl", + runTime.constant(), // directory + "triSurface", // instance + runTime, // registry + IOobject::NO_READ, + IOobject::NO_WRITE + ), + triSurfaceTools::triangulate + ( + patches, + includePatches + ), + dict +); + +//surfacesMesh.searchableSurface::write(); diff --git a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H new file mode 100644 index 000000000..3cdca36d6 --- /dev/null +++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H @@ -0,0 +1,83 @@ +// All rays expressed as start face (local) index end end face (global) +// Pre-size by assuming a certain percentage is visible. + +// Maximum lenght for dynamicList +const label maxDynListLength = 10000; + +for (label procI = 0; procI < Pstream::nProcs(); procI++) +{ + // Shoot rays from me to procI. Note that even if processor has + // 0 faces we still need to call findLine to keep calls synced. + + DynamicField start(coarseMesh.nFaces()); + DynamicField end(start.size()); + DynamicList