From a8c917fd4b69e0cc8e0295ce921cbacfa630365d Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 11 Dec 2013 17:26:34 +0000 Subject: [PATCH 1/3] DyM solvers: correct Uf using phi after construction --- applications/solvers/combustion/engineFoam/pEqn.H | 5 +---- .../compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H | 5 +---- .../solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H | 5 +---- .../incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H | 2 +- .../solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H | 5 +---- .../multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H | 2 +- .../compressibleInterFoam/compressibleInterDyMFoam/pEqn.H | 2 +- .../multiphase/interFoam/interDyMFoam/interDyMFoam.C | 2 +- .../solvers/multiphase/interFoam/interDyMFoam/pEqn.H | 2 +- .../interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C | 2 +- .../interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H | 2 +- src/finiteVolume/cfdTools/compressible/createRhoUf.H | 7 ++++++- src/finiteVolume/cfdTools/incompressible/createUf.H | 8 +++++++- 13 files changed, 24 insertions(+), 25 deletions(-) diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H index 2f6333cc97..bd12a82495 100644 --- a/applications/solvers/combustion/engineFoam/pEqn.H +++ b/applications/solvers/combustion/engineFoam/pEqn.H @@ -91,10 +91,7 @@ K = 0.5*magSqr(U); { rhoUf = fvc::interpolate(rho*U); surfaceVectorField n(mesh.Sf()/mesh.magSf()); - rhoUf += - mesh.Sf() - *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf)) - /sqr(mesh.magSf()); + rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf)); } if (thermo.dpdt()) diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H index fd892bb932..1c58a02c33 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H @@ -107,10 +107,7 @@ K = 0.5*magSqr(U); { rhoUf = fvc::interpolate(rho*U); surfaceVectorField n(mesh.Sf()/mesh.magSf()); - rhoUf += - mesh.Sf() - *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf)) - /sqr(mesh.magSf()); + rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf)); } if (thermo.dpdt()) diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H index 65afd1aeb8..4799f64567 100644 --- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H +++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H @@ -42,8 +42,5 @@ U.correctBoundaryConditions(); { rhoUf = fvc::interpolate(rho*U); surfaceVectorField n(mesh.Sf()/mesh.magSf()); - rhoUf += - mesh.Sf() - *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf)) - /sqr(mesh.magSf()); + rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf)); } diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H index 7ecfe1f662..e2938bdb94 100644 --- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H +++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H @@ -51,7 +51,7 @@ fvOptions.correct(U); { Uf = fvc::interpolate(U); surfaceVectorField n(mesh.Sf()/mesh.magSf()); - Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf()); + Uf += n*(phi/mesh.magSf() - (n & Uf)); } // Make the fluxes relative to the mesh motion diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H index b10c24aa0b..d23ec0f950 100644 --- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H +++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H @@ -93,10 +93,7 @@ K = 0.5*magSqr(U); { rhoUf = fvc::interpolate(rho*U); surfaceVectorField n(mesh.Sf()/mesh.magSf()); - rhoUf += - mesh.Sf() - *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf)) - /sqr(mesh.magSf()); + rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf)); } if (thermo.dpdt()) diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H index 0f0d576841..a0e0e6c74d 100644 --- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H +++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H @@ -85,6 +85,6 @@ { Uf = fvc::interpolate(U); surfaceVectorField n(mesh.Sf()/mesh.magSf()); - Uf += mesh.Sf()*(phiv - (mesh.Sf() & Uf))/sqr(mesh.magSf()); + Uf += n*(phiv/mesh.magSf() - (n & Uf)); } } diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H index c86921371f..a014a5ab4e 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H @@ -114,7 +114,7 @@ { Uf = fvc::interpolate(U); surfaceVectorField n(mesh.Sf()/mesh.magSf()); - Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf()); + Uf += n*(phi/mesh.magSf() - (n & Uf)); } // Make the fluxes relative to the mesh motion diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C index 0a3d7ed11e..ba7b4efa9a 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C @@ -55,10 +55,10 @@ int main(int argc, char *argv[]) pimpleControl pimple(mesh); #include "createFields.H" - #include "createUf.H" #include "readTimeControls.H" #include "createPrghCorrTypes.H" #include "correctPhi.H" + #include "createUf.H" #include "CourantNo.H" #include "setInitialDeltaT.H" diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H index c30b361e72..7644120791 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H @@ -67,7 +67,7 @@ { Uf = fvc::interpolate(U); surfaceVectorField n(mesh.Sf()/mesh.magSf()); - Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf()); + Uf += n*(phi/mesh.magSf() - (n & Uf)); } // Make the fluxes relative to the mesh motion diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C index 0bfd8ad56e..810edeff30 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C @@ -65,10 +65,10 @@ int main(int argc, char *argv[]) pimpleControl pimple(mesh); #include "createFields.H" - #include "createUf.H" #include "readTimeControls.H" #include "createPcorrTypes.H" #include "../interFoam/interDyMFoam/correctPhi.H" + #include "createUf.H" #include "CourantNo.H" #include "setInitialDeltaT.H" diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H index 7863f5d748..6fac1df95c 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H @@ -68,7 +68,7 @@ { Uf = fvc::interpolate(U); surfaceVectorField n(mesh.Sf()/mesh.magSf()); - Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf()); + Uf += n*(phi/mesh.magSf() - (n & Uf)); } // Make the fluxes relative to the mesh motion diff --git a/src/finiteVolume/cfdTools/compressible/createRhoUf.H b/src/finiteVolume/cfdTools/compressible/createRhoUf.H index e91115ae75..312aec60dc 100644 --- a/src/finiteVolume/cfdTools/compressible/createRhoUf.H +++ b/src/finiteVolume/cfdTools/compressible/createRhoUf.H @@ -46,9 +46,14 @@ surfaceVectorField rhoUf IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), - linearInterpolate(rho*U) + fvc::interpolate(rho*U) ); +{ + surfaceVectorField n(mesh.Sf()/mesh.magSf()); + rhoUf += n*(phi/mesh.magSf() - (n & rhoUf)); +} + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif diff --git a/src/finiteVolume/cfdTools/incompressible/createUf.H b/src/finiteVolume/cfdTools/incompressible/createUf.H index aab92ce0e2..0b65f2daed 100644 --- a/src/finiteVolume/cfdTools/incompressible/createUf.H +++ b/src/finiteVolume/cfdTools/incompressible/createUf.H @@ -46,9 +46,15 @@ surfaceVectorField Uf IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), - linearInterpolate(U) + fvc::interpolate(U) ); +{ + surfaceVectorField n(mesh.Sf()/mesh.magSf()); + Uf += n*(phi/mesh.magSf() - (n & Uf)); +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif From 79614840b3bcb74ed30aa688218c85890f50ee82 Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 11 Dec 2013 17:27:27 +0000 Subject: [PATCH 2/3] All ldu solvers: new optional parameter "minIter" to specify the minimum number of iterations the solver should perform. --- .../matrices/LduMatrix/LduMatrix/LduMatrix.H | 3 +++ .../matrices/LduMatrix/LduMatrix/LduMatrixSolver.C | 2 ++ .../matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C | 13 ++++++++++--- .../matrices/LduMatrix/Solvers/PCICG/PCICG.C | 13 ++++++++++--- .../LduMatrix/Solvers/SmoothSolver/SmoothSolver.C | 13 ++++++++++--- .../matrices/lduMatrix/lduMatrix/lduMatrix.C | 2 +- .../matrices/lduMatrix/lduMatrix/lduMatrix.H | 3 +++ .../matrices/lduMatrix/lduMatrix/lduMatrixSolver.C | 1 + .../lduMatrix/solvers/GAMG/GAMGSolverSolve.C | 13 ++++++++++--- .../matrices/lduMatrix/solvers/PBiCG/PBiCG.C | 13 ++++++++++--- src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C | 13 ++++++++++--- .../lduMatrix/solvers/smoothSolver/smoothSolver.C | 13 ++++++++++--- 12 files changed, 80 insertions(+), 22 deletions(-) diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H index c2414b0f8b..71aae68507 100644 --- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H +++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H @@ -125,6 +125,9 @@ public: //- Maximum number of iterations in the solver label maxIter_; + //- Minimum number of iterations in the solver + label minIter_; + //- Final convergence tolerance Type tolerance_; diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C index 50831a0ec9..ef3afff4ab 100644 --- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C +++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C @@ -135,6 +135,7 @@ Foam::LduMatrix::solver::solver controlDict_(solverDict), maxIter_(1000), + minIter_(0), tolerance_(1e-6*pTraits::one), relTol_(pTraits::zero) { @@ -148,6 +149,7 @@ template void Foam::LduMatrix::solver::readControls() { readControl(controlDict_, maxIter_, "maxIter"); + readControl(controlDict_, minIter_, "minIter"); readControl(controlDict_, tolerance_, "tolerance"); readControl(controlDict_, relTol_, "relTol"); } diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C index 7ee6da05dc..302785c1c4 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C @@ -104,7 +104,11 @@ Foam::PBiCCCG::solve solverPerf.finalResidual() = solverPerf.initialResidual(); // --- Check convergence, solve if not converged - if (!solverPerf.checkConvergence(this->tolerance_, this->relTol_)) + if + ( + this->minIter_ > 0 + || !solverPerf.checkConvergence(this->tolerance_, this->relTol_) + ) { // --- Select and construct the preconditioner autoPtr::preconditioner> @@ -182,8 +186,11 @@ Foam::PBiCCCG::solve } while ( - solverPerf.nIterations()++ < this->maxIter_ - && !(solverPerf.checkConvergence(this->tolerance_, this->relTol_)) + ( + solverPerf.nIterations()++ < this->maxIter_ + && !solverPerf.checkConvergence(this->tolerance_, this->relTol_) + ) + || solverPerf.nIterations() < this->minIter_ ); } diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C index cc8f1d1c72..87d087d199 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C @@ -92,7 +92,11 @@ Foam::PCICG::solve(Field& psi) const solverPerf.finalResidual() = solverPerf.initialResidual(); // --- Check convergence, solve if not converged - if (!solverPerf.checkConvergence(this->tolerance_, this->relTol_)) + if + ( + this->minIter_ > 0 + || !solverPerf.checkConvergence(this->tolerance_, this->relTol_) + ) { // --- Select and construct the preconditioner autoPtr::preconditioner> @@ -174,8 +178,11 @@ Foam::PCICG::solve(Field& psi) const } while ( - solverPerf.nIterations()++ < this->maxIter_ - && !(solverPerf.checkConvergence(this->tolerance_, this->relTol_)) + ( + solverPerf.nIterations()++ < this->maxIter_ + && !solverPerf.checkConvergence(this->tolerance_, this->relTol_) + ) + || solverPerf.nIterations() < this->minIter_ ); } diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C index 9cf5c81223..b6725506ed 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C @@ -113,7 +113,11 @@ Foam::SmoothSolver::solve(Field& psi) const // Check convergence, solve if not converged - if (!solverPerf.checkConvergence(this->tolerance_, this->relTol_)) + if + ( + this->minIter_ > 0 + || !solverPerf.checkConvergence(this->tolerance_, this->relTol_) + ) { autoPtr::smoother> smootherPtr = LduMatrix::smoother::New @@ -140,8 +144,11 @@ Foam::SmoothSolver::solve(Field& psi) const ); } while ( - (solverPerf.nIterations() += nSweeps_) < this->maxIter_ - && !(solverPerf.checkConvergence(this->tolerance_, this->relTol_)) + ( + (solverPerf.nIterations() += nSweeps_) < this->maxIter_ + && !solverPerf.checkConvergence(this->tolerance_, this->relTol_) + ) + || solverPerf.nIterations() < this->minIter_ ); } } diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C index fd9c43da08..b8bc5c930c 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C @@ -31,7 +31,7 @@ License namespace Foam { -defineTypeNameAndDebug(lduMatrix, 1); + defineTypeNameAndDebug(lduMatrix, 1); } diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H index a75668dfdd..fdaf834a74 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H @@ -107,6 +107,9 @@ public: //- Maximum number of iterations in the solver label maxIter_; + //- Minimum number of iterations in the solver + label minIter_; + //- Final convergence tolerance scalar tolerance_; diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C index cee5c1dc2f..aed21fe9f7 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C @@ -164,6 +164,7 @@ Foam::lduMatrix::solver::solver void Foam::lduMatrix::solver::readControls() { maxIter_ = controlDict_.lookupOrDefault