From 56e21d8f61bac64af2887f34a1227dd1036167c1 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 27 Apr 2015 15:16:23 +0100 Subject: [PATCH 01/21] Formatting --- src/finiteVolume/cfdTools/general/MRF/MRFZone.H | 4 ++-- src/finiteVolume/cfdTools/general/MRF/MRFZones.C | 3 +++ 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H index 56eea426c..1758a863e 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H @@ -96,7 +96,7 @@ class MRFZone //- Divide faces in frame according to patch void setMRFFaces(); - //- Make the given absolute mass/vol flux relative within the MRF region + //- Make the given absolute mass/vol flux relative within MRF region template void relativeRhoFlux ( @@ -104,7 +104,7 @@ class MRFZone surfaceScalarField& phi ) const; - //- Make the given relative mass/vol flux absolute within the MRF region + //- Make the given relative mass/vol flux absolute within MRF region template void absoluteRhoFlux ( diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZones.C b/src/finiteVolume/cfdTools/general/MRF/MRFZones.C index 2081eaaa3..b68601548 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZones.C +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZones.C @@ -85,6 +85,7 @@ Foam::tmp Foam::MRFZones::fluxCorrection() const return tMRFZonesPhiCorr; } + Foam::tmp Foam::MRFZones::meshPhi() const { tmp tMRFZonesPhiCorr @@ -113,6 +114,7 @@ Foam::tmp Foam::MRFZones::meshPhi() const return tMRFZonesPhiCorr; } + void Foam::MRFZones::addCoriolis(fvVectorMatrix& UEqn) const { forAll(*this, i) @@ -275,4 +277,5 @@ Foam::tmp Foam::MRFZones::Su return tPhiSource; } + // ************************************************************************* // From 76b941e7ff59526f1453305da2fb6cb144334e50 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 27 Apr 2015 15:34:02 +0100 Subject: [PATCH 02/21] Coupled MRF porous solver --- .../coupled/MRFPorousFoam/MRFPorousFoam.C | 120 ++++++++++++++++++ .../solvers/coupled/MRFPorousFoam/Make/files | 3 + .../coupled/MRFPorousFoam/Make/options | 12 ++ .../solvers/coupled/MRFPorousFoam/UEqn.H | 15 +++ .../solvers/coupled/MRFPorousFoam/boundPU.H | 32 +++++ .../coupled/MRFPorousFoam/convergenceCheck.H | 8 ++ .../coupled/MRFPorousFoam/couplingTerms.H | 15 +++ .../coupled/MRFPorousFoam/createFields.H | 52 ++++++++ .../solvers/coupled/MRFPorousFoam/createMRF.H | 2 + .../coupled/MRFPorousFoam/createPorous.H | 1 + .../MRFPorousFoam/initConvergenceCheck.H | 4 + .../solvers/coupled/MRFPorousFoam/pEqn.H | 23 ++++ .../MRFPorousFoam/readBlockSolverControls.H | 15 +++ .../coupled/MRFPorousFoam/readFieldBounds.H | 14 ++ 14 files changed, 316 insertions(+) create mode 100644 applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C create mode 100644 applications/solvers/coupled/MRFPorousFoam/Make/files create mode 100644 applications/solvers/coupled/MRFPorousFoam/Make/options create mode 100644 applications/solvers/coupled/MRFPorousFoam/UEqn.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/boundPU.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/convergenceCheck.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/couplingTerms.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/createFields.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/createMRF.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/createPorous.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/pEqn.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/readFieldBounds.H diff --git a/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C b/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C new file mode 100644 index 000000000..9ee0e7f41 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C @@ -0,0 +1,120 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Application + MRFPorousFoam + +Description + Steady-state solver for incompressible, turbulent flow, with implicit + coupling between pressure and velocity achieved by fvBlockMatrix. + Turbulence is in this version solved using the existing turbulence + structure. + + Added support for MRF and porous zones + +Authors + Hrvoje Jasak, Wikki Ltd. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "fvBlockMatrix.H" +#include "singlePhaseTransportModel.H" +#include "RASModel.H" +#include "MRFZones.H" +#include "porousZones.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" +# include "createFields.H" +# include "createMRF.H" +# include "createPorous.H" +# include "initContinuityErrs.H" +# include "initConvergenceCheck.H" + + Info<< "\nStarting time loop\n" << endl; + while (runTime.loop()) + { +# include "readBlockSolverControls.H" +# include "readFieldBounds.H" + + Info<< "Time = " << runTime.timeName() << nl << endl; + + p.storePrevIter(); + + // Initialize the Up block system (matrix, source and reference to Up) + fvBlockMatrix UpEqn(Up); + + // Assemble and insert momentum equation +# include "UEqn.H" + + // Assemble and insert pressure equation +# include "pEqn.H" + + // Assemble and insert coupling terms +# include "couplingTerms.H" + + // Solve the block matrix + maxResidual = cmptMax(UpEqn.solve().initialResidual()); + + // Retrieve solution + UpEqn.retrieveSolution(0, U.internalField()); + UpEqn.retrieveSolution(3, p.internalField()); + + U.correctBoundaryConditions(); + p.correctBoundaryConditions(); + + phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource; + +# include "continuityErrs.H" + + // Make flux relative in rotating zones + mrfZones.relativeFlux(phi); + +# include "continuityErrs.H" + +# include "boundPU.H" + + p.relax(); + + turbulence->correct(); + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + +# include "convergenceCheck.H" + } + + Info<< "End\n" << endl; + + return 0; +} diff --git a/applications/solvers/coupled/MRFPorousFoam/Make/files b/applications/solvers/coupled/MRFPorousFoam/Make/files new file mode 100644 index 000000000..2ae8a1e4f --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/Make/files @@ -0,0 +1,3 @@ +MRFPorousFoam.C + +EXE = $(FOAM_APPBIN)/MRFPorousFoam diff --git a/applications/solvers/coupled/MRFPorousFoam/Make/options b/applications/solvers/coupled/MRFPorousFoam/Make/options new file mode 100644 index 000000000..5d2b99e6d --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/Make/options @@ -0,0 +1,12 @@ +EXE_INC = \ + -I$(LIB_SRC)/turbulenceModels \ + -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +EXE_LIBS = \ + -lincompressibleRASModels \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -llduSolvers diff --git a/applications/solvers/coupled/MRFPorousFoam/UEqn.H b/applications/solvers/coupled/MRFPorousFoam/UEqn.H new file mode 100644 index 000000000..ed0e2bc44 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/UEqn.H @@ -0,0 +1,15 @@ + // Momentum equation + fvVectorMatrix UEqn + ( + fvm::ddt(U) + + fvm::div(phi, U) + + turbulence->divDevReff(U) + ); + + // Add MRF and porous sources + mrfZones.addCoriolis(UEqn); + pZones.addResistance(UEqn); + + UEqn.relax(); + + UpEqn.insertEquation(0, UEqn); diff --git a/applications/solvers/coupled/MRFPorousFoam/boundPU.H b/applications/solvers/coupled/MRFPorousFoam/boundPU.H new file mode 100644 index 000000000..eca655e8c --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/boundPU.H @@ -0,0 +1,32 @@ +{ + // Bound the pressure + dimensionedScalar p1 = min(p); + dimensionedScalar p2 = max(p); + + if (p1 < pMin || p2 > pMax) + { + Info<< "p: " << p1.value() << " " << p2.value() + << ". Bounding." << endl; + + p.max(pMin); + p.min(pMax); + p.correctBoundaryConditions(); + } + + // Bound the velocity + volScalarField magU = mag(U); + dimensionedScalar U1 = max(magU); + + if (U1 > UMax) + { + Info<< "U: " << U1.value() << ". Bounding." << endl; + + volScalarField Ulimiter = pos(magU - UMax)*UMax/(magU + smallU) + + neg(magU - UMax); + Ulimiter.max(scalar(0)); + Ulimiter.min(scalar(1)); + + U *= Ulimiter; + U.correctBoundaryConditions(); + } +} diff --git a/applications/solvers/coupled/MRFPorousFoam/convergenceCheck.H b/applications/solvers/coupled/MRFPorousFoam/convergenceCheck.H new file mode 100644 index 000000000..9c7734131 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/convergenceCheck.H @@ -0,0 +1,8 @@ +// Check convergence +if (maxResidual < convergenceCriterion) +{ + Info<< "reached convergence criterion: " << convergenceCriterion << endl; + runTime.writeAndEnd(); + Info<< "latestTime = " << runTime.timeName() << endl; +} + diff --git a/applications/solvers/coupled/MRFPorousFoam/couplingTerms.H b/applications/solvers/coupled/MRFPorousFoam/couplingTerms.H new file mode 100644 index 000000000..9aa6f9510 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/couplingTerms.H @@ -0,0 +1,15 @@ +{ + // Calculate grad p coupling matrix. Needs to be here if one uses + // gradient schemes with limiters. VV, 9/June/2014 + BlockLduSystem pInU(fvm::grad(p)); + + // Calculate div U coupling. Could be calculated only once since + // it is only geometry dependent. VV, 9/June/2014 + BlockLduSystem UInp(fvm::UDiv(U)); + + // Last argument in insertBlockCoupling says if the column direction + // should be incremented. This is needed for arbitrary positioning + // of U and p in the system. This could be better. VV, 30/April/2014 + UpEqn.insertBlockCoupling(0, 3, pInU, true); + UpEqn.insertBlockCoupling(3, 0, UInp, false); +} diff --git a/applications/solvers/coupled/MRFPorousFoam/createFields.H b/applications/solvers/coupled/MRFPorousFoam/createFields.H new file mode 100644 index 000000000..51403c4a7 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/createFields.H @@ -0,0 +1,52 @@ +Info << "Reading field p\n" << endl; +volScalarField p +( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +Info << "Reading field U\n" << endl; +volVectorField U +( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +#include "createPhi.H" + +singlePhaseTransportModel laminarTransport(U, phi); +autoPtr turbulence +( + incompressible::RASModel::New(U, phi, laminarTransport) +); + +// Block vector field for velocity (first entry) and pressure (second +// entry). +Info << "Creating field Up\n" << endl; +volVector4Field Up +( + IOobject + ( + "Up", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedVector4("zero", dimless, vector4::zero) +); diff --git a/applications/solvers/coupled/MRFPorousFoam/createMRF.H b/applications/solvers/coupled/MRFPorousFoam/createMRF.H new file mode 100644 index 000000000..161446a8e --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/createMRF.H @@ -0,0 +1,2 @@ + MRFZones mrfZones(mesh); + mrfZones.correctBoundaryVelocity(U); diff --git a/applications/solvers/coupled/MRFPorousFoam/createPorous.H b/applications/solvers/coupled/MRFPorousFoam/createPorous.H new file mode 100644 index 000000000..430b466aa --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/createPorous.H @@ -0,0 +1 @@ + porousZones pZones(mesh); diff --git a/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H b/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H new file mode 100644 index 000000000..482bc421d --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H @@ -0,0 +1,4 @@ +// initialize values for convergence checks + + scalar maxResidual = 0; + scalar convergenceCriterion = 0; diff --git a/applications/solvers/coupled/MRFPorousFoam/pEqn.H b/applications/solvers/coupled/MRFPorousFoam/pEqn.H new file mode 100644 index 000000000..9029b873f --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/pEqn.H @@ -0,0 +1,23 @@ +// Pressure parts of the continuity equation +surfaceScalarField rUAf +( + "rUAf", + fvc::interpolate(1.0/UEqn.A()) +); + +surfaceScalarField presSource +( + "presSource", + rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf()) +); + +fvScalarMatrix pEqn +( + - fvm::laplacian(rUAf, p) + == + - fvc::div(presSource) +); + +pEqn.setReference(pRefCell, pRefValue); + +UpEqn.insertEquation(3, pEqn); diff --git a/applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H b/applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H new file mode 100644 index 000000000..93c2bda47 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H @@ -0,0 +1,15 @@ + label pRefCell = 0; + scalar pRefValue = 0; + setRefCell + ( + p, + mesh.solutionDict().subDict("blockSolver"), + pRefCell, + pRefValue + ); + + mesh.solutionDict().subDict("blockSolver").readIfPresent + ( + "convergence", + convergenceCriterion + ); diff --git a/applications/solvers/coupled/MRFPorousFoam/readFieldBounds.H b/applications/solvers/coupled/MRFPorousFoam/readFieldBounds.H new file mode 100644 index 000000000..4a6c6f787 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/readFieldBounds.H @@ -0,0 +1,14 @@ + // Read field bounds + dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds"); + + // Pressure bounds + dimensionedScalar pMin("pMin", p.dimensions(), 0); + dimensionedScalar pMax("pMax", p.dimensions(), 0); + + fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value(); + + // Velocity bound + dimensionedScalar UMax("UMax", U.dimensions(), 0); + + fieldBounds.lookup(U.name()) >> UMax.value(); + dimensionedScalar smallU("smallU", dimVelocity, 1e-10); From 95fb99d2a6629be7e15ac4601afa84834ab685e1 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 27 Apr 2015 18:09:57 +0100 Subject: [PATCH 03/21] Feature: block coupled updates --- .../fvMatrices/fvBlockMatrix/fvBlockMatrix.C | 273 ++---------------- .../fvMatrices/fvBlockMatrix/fvBlockMatrix.H | 36 +-- 2 files changed, 44 insertions(+), 265 deletions(-) diff --git a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C index f2697f212..18dbee6ff 100644 --- a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C +++ b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C @@ -30,7 +30,7 @@ License template template -void fvBlockMatrix::insertSolutionVector +void Foam::fvBlockMatrix::insertSolutionVector ( const direction dir, const Field& xSingle @@ -59,7 +59,7 @@ void fvBlockMatrix::insertSolutionVector template template -void fvBlockMatrix::insertDiagSource +void Foam::fvBlockMatrix::insertDiagSource ( const direction dir, fvMatrix& matrix @@ -162,7 +162,7 @@ void fvBlockMatrix::insertDiagSource template template -void fvBlockMatrix::insertUpperLower +void Foam::fvBlockMatrix::insertUpperLower ( const direction dir, const fvMatrix& matrix @@ -290,7 +290,7 @@ void fvBlockMatrix::insertUpperLower template template -void fvBlockMatrix::updateCouplingCoeffs +void Foam::fvBlockMatrix::updateCouplingCoeffs ( const direction dir, const fvMatrix& matrix @@ -299,7 +299,9 @@ void fvBlockMatrix::updateCouplingCoeffs const direction nCmpts = pTraits::nComponents; direction localDir = dir; - const GeometricField& psi = matrix.psi(); + const GeometricField& psi = + matrix.psi(); + forAll(psi.boundaryField(), patchI) { const fvPatchField& pf = psi.boundaryField()[patchI]; @@ -373,7 +375,7 @@ void fvBlockMatrix::updateCouplingCoeffs template template -void fvBlockMatrix::insertBlock +void Foam::fvBlockMatrix::insertBlock ( const direction dirI, const direction dirJ, @@ -471,7 +473,7 @@ void fvBlockMatrix::insertBlock template template -void fvBlockMatrix::insertBoundaryContributions +void Foam::fvBlockMatrix::insertBoundaryContributions ( const direction dirI, const direction dirJ, @@ -561,39 +563,30 @@ void fvBlockMatrix::insertBoundaryContributions template -void fvBlockMatrix::insertCouplingDiagSource +void Foam::fvBlockMatrix::insertCouplingDiag ( const direction dirI, const direction dirJ, - const fvScalarMatrix& matrix + const scalarField& coeffIJ ) { - // Prepare the diagonal and source - scalarField diag = matrix.diag(); - scalarField source = matrix.source(); - - // Add boundary source contribution - matrix.addBoundaryDiag(diag, 0); - matrix.addBoundarySource(source, false); - // Get reference to block diagonal of the block system typename CoeffField::squareTypeField& blockDiag = this->diag().asSquare(); - // Get reference to this source field of the block system - Field& b = this->source(); - // Set off-diagonal coefficient - forAll(diag, cellI) + forAll(coeffIJ, cellI) { - blockDiag[cellI](dirI, dirJ) += diag[cellI]; - b[cellI](dirI) += source[cellI]; + blockDiag[cellI](dirI, dirJ) += coeffIJ[cellI]; } + + // Source compensation is done in function updateSourceCoupling() + // after all coupling terms are added. HJ, 27/Apr/2015 } template -void fvBlockMatrix::insertCouplingUpperLower +void Foam::fvBlockMatrix::insertCouplingUpperLower ( const direction dirI, const direction dirJ, @@ -707,7 +700,7 @@ Foam::fvBlockMatrix::fvBlockMatrix template template -void fvBlockMatrix::retrieveSolution +void Foam::fvBlockMatrix::retrieveSolution ( const direction dir, Field& xSingle @@ -736,7 +729,7 @@ void fvBlockMatrix::retrieveSolution template template -void fvBlockMatrix::insertEquation +void Foam::fvBlockMatrix::insertEquation ( const direction dir, fvMatrix& matrix @@ -751,7 +744,7 @@ void fvBlockMatrix::insertEquation template template -void fvBlockMatrix::insertBlockCoupling +void Foam::fvBlockMatrix::insertBlockCoupling ( const direction dirI, const direction dirJ, @@ -765,20 +758,22 @@ void fvBlockMatrix::insertBlockCoupling template -void fvBlockMatrix::insertEquationCoupling +void Foam::fvBlockMatrix::insertEquationCoupling ( const direction dirI, const direction dirJ, - const fvScalarMatrix& matrix + const scalarField& coeffIJ ) { - insertCouplingDiagSource(dirI, dirJ, matrix); - insertCouplingUpperLower(dirI, dirJ, matrix); + // Multiply coefficients by volume + scalarField coeffIJVol = coeffIJ*psi_.mesh().V(); + + insertCouplingDiag(dirI, dirJ, coeffIJVol); } template -void fvBlockMatrix::blockAdd +void Foam::fvBlockMatrix::blockAdd ( const direction dir, const scalarField& xSingle, @@ -793,7 +788,7 @@ void fvBlockMatrix::blockAdd template -void fvBlockMatrix::updateSourceCoupling() +void Foam::fvBlockMatrix::updateSourceCoupling() { // Eliminate off-diagonal block coefficients from the square diagonal // With this change, the segregated matrix can be assembled with complete @@ -827,215 +822,7 @@ void fvBlockMatrix::updateSourceCoupling() template -void fvBlockMatrix::insertPicardTensor -( - const direction UEqnDir, - const volVectorField& U, - const surfaceScalarField& phi -) -{ - const direction nCmpts = pTraits::nComponents; - direction localDirI = UEqnDir; - direction localDirJ = UEqnDir; - - // Sanity check - { - const direction blockMatrixSize = pTraits::nComponents; - - if (nCmpts > blockMatrixSize) - { - FatalErrorIn - ( - "void fvBlockMatrix::insertPicardTensor\n" - "(\n" - " const direction UEqnDir,\n" - " const volVectorField& U,\n" - " const surfaceScalarField& phi\n" - ")" - ) << "Trying to insert Picard tensor term into smaller " - << "fvBlockMatrix. Do you have momentum equation?" - << abort(FatalError); - } - } - - // Get weights for U which makes the implicit flux part - const fvMesh& mesh = U.mesh(); - - tmp > tinterpScheme = - fvc::scheme(mesh, U.name()); - - tmp tweights = tinterpScheme().weights(U); - const scalarField& wIn = tweights().internalField(); - - // Calculate the Pi tensor. Consider hard coding the interpolation scheme to - // correspond to the div(U, phi) interpolation scheme for consistency. - // VV, 21/July/2014. - const surfaceTensorField pi(mesh.Sf()*fvc::interpolate(U, phi, "(phi,U)")); -// const surfaceTensorField pi(fvc::interpolate(U, phi, "(phi,U)")*mesh.Sf()); - const tensorField& piIn = pi.internalField(); - - BlockLduSystem piSystem(mesh); - - // Get references to ldu fields of pi block system always as square - typename CoeffField::squareTypeField& piDiag = - piSystem.diag().asSquare(); - typename CoeffField::squareTypeField& piUpper = - piSystem.upper().asSquare(); - typename CoeffField::squareTypeField& piLower = - piSystem.lower().asSquare(); - - vectorField& piSource = piSystem.source(); - - piLower = -wIn*piIn; - piUpper = piLower + piIn; - piSystem.negSumDiag(); - - Info << "Diag coeff: " << piDiag[125] << nl - << "Lower coeff: " << piLower[125] << nl - << "Upper coeff: " << piUpper[125] << endl; - - // Boundary contributions - treated explicitly because of the problem with - // inconsistent return type of internalCoeffs. VV, 21/July/2014. - forAll(U.boundaryField(), patchI) - { - const fvPatchVectorField& Ub = U.boundaryField()[patchI]; - const fvPatch& patch = Ub.patch(); - const tensorField& pib = pi.boundaryField()[patchI]; - const fvsPatchScalarField& wb = tweights().boundaryField()[patchI]; - const unallocLabelList& fc = patch.faceCells(); - -// const vectorField internalCoeffs(Ub.valueInternalCoeffs(wb)); - - // Explicit diag boundary contribution -// forAll(Ub, faceI) -// { -// piSource[fc[faceI]] -= pib[faceI] & Ub[faceI]; -// piSource[fc[faceI]] -= Ub[faceI] & pib[faceI]; -// } - - // Hard coded zeroGradient if patch does not fix value - if (!U.boundaryField()[patchI].fixesValue()) - { - forAll(Ub, faceI) - { - piDiag[fc[faceI]] += pib[faceI]; - } - } - // Coupled patches treated implicitly - else if (patch.coupled()) - { - typename CoeffField::squareTypeField& pipCoupleUpper = - piSystem.coupleUpper()[patchI].asSquare(); - typename CoeffField::squareTypeField& pipCoupleLower = - piSystem.coupleLower()[patchI].asSquare(); - - const tensorField pcl = -wb*pib; - const tensorField pcu = pcl + pib; - - // Coupling contributions - pipCoupleLower -= pcl; - pipCoupleUpper -= pcu; - } - else - { - const vectorField boundaryCoeffs(Ub.valueBoundaryCoeffs(wb)); - - // Boundary contribution - forAll(Ub, faceI) - { - piSource[fc[faceI]] -= pib[faceI] & boundaryCoeffs[faceI]; - } - } - } - - // Consider chucking the above code into fvm::picardTensor operator. - // VV, 21/July/2014. - - // Finally add Picard piSystem into this fvBlockMatrix - typename CoeffField::squareTypeField& blockDiag = - this->diag().asSquare(); - typename CoeffField::squareTypeField& blockUpper = - this->upper().asSquare(); - typename CoeffField::squareTypeField& blockLower = - this->lower().asSquare(); - - Field& blockSource = this->source(); - - // Add diagonal, source, lower and upper - for (direction cmptI = 0; cmptI < nCmpts; cmptI++) - { - for (direction cmptJ = 0; cmptJ < nCmpts; cmptJ++) - { - forAll(piDiag, cellI) - { - blockDiag[cellI](localDirI, localDirJ) += - piDiag[cellI](cmptI, cmptJ); -// blockSource[cellI](localDirI, localDirJ) += -// piSource[cellI](cmptI, cmptJ); - } - - forAll(piUpper, faceI) - { - blockUpper[faceI](localDirI, localDirJ) += - piUpper[faceI](cmptI, cmptJ); - blockLower[faceI](localDirI, localDirJ) += - piLower[faceI](cmptI, cmptJ); - } - - localDirJ++; - } - - forAll(piSource, cellI) - { - blockSource[cellI](localDirI) += piSource[cellI](cmptI); - } - - localDirI++; - } - - // Reset local direction for coupling contributions - localDirI = UEqnDir; - localDirJ = UEqnDir; - - // Add coupling contributions - forAll(U.boundaryField(), patchI) - { - if (U.boundaryField()[patchI].patch().coupled()) - { - typename CoeffField::squareTypeField& pcoupleUpper = - this->coupleUpper()[patchI].asSquare(); - typename CoeffField::squareTypeField& pcoupleLower = - this->coupleLower()[patchI].asSquare(); - - const typename CoeffField::squareTypeField& pipcu = - piSystem.coupleUpper()[patchI].asSquare(); - const typename CoeffField::squareTypeField& pipcl = - piSystem.coupleLower()[patchI].asSquare(); - - for (direction cmptI = 0; cmptI < nCmpts; cmptI++) - { - for (direction cmptJ = 0; cmptJ < nCmpts; cmptJ++) - { - forAll(pipcu, faceI) - { - pcoupleUpper[faceI](localDirI, localDirJ) += - pipcu[faceI](cmptI, cmptJ); - pcoupleLower[faceI](localDirI, localDirJ) += - pipcl[faceI](cmptI, cmptJ); - } - } - } - - // Reset local directions for other patches - localDirI = UEqnDir; - localDirJ = UEqnDir; - } - } -} - - -template -Foam::BlockSolverPerformance fvBlockMatrix::solve +Foam::BlockSolverPerformance Foam::fvBlockMatrix::solve ( const dictionary& solverControls ) @@ -1057,7 +844,7 @@ Foam::BlockSolverPerformance fvBlockMatrix::solve template -Foam::BlockSolverPerformance fvBlockMatrix::solve() +Foam::BlockSolverPerformance Foam::fvBlockMatrix::solve() { return solve(psi_.mesh().solutionDict().solverDict(psi_.name())); } diff --git a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H index d212a7253..972db44d8 100644 --- a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H +++ b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H @@ -26,11 +26,12 @@ Class Description fvBlockMatrix is an extension of fvMatrix for block coupled types. It holds - general insertion and retrieval tools for block coupling along with specific + general insertion and retrieval tools for block coupling and specific functions for p-U coupling. Author Vuko Vukcevic, FMENA Zagreb. + Update by Hrvoje Jasak SourceFiles fvBlockMatrix.C @@ -112,14 +113,13 @@ class fvBlockMatrix // Insertion functions for fvScalarMatrix into off-diagonal positions // (coupling matrices) - // These two functions are dubious. Reconsider. HJ, 17/July/2014. - //- Insert coupling matrix diag and source into this fvBlockMatrix - void insertCouplingDiagSource + //- Insert coupling matrix diag element into this fvBlockMatrix + void insertCouplingDiag ( const direction dirI, const direction dirJ, - const fvScalarMatrix& matrix + const scalarField& coeffIJ ); //- Insert coupling matrix lower and upper into this fvBlockMatrix @@ -144,8 +144,8 @@ class fvBlockMatrix const bool incrementColumnDir ); - //- Insert source and coupling coeffs of BlockLduSystem (obtained by - // implicit grad/div operator) into this fvBlockMatrix + //- Insert source and coupling coeffs of BlockLduSystem + // (eg. obtained by implicit grad/div operator) template void insertBoundaryContributions ( @@ -192,7 +192,7 @@ public: // Insertion and retrieval public tools - //- Retrieve part of internal (solved) field from this fvBlockMatrix + //- Retrieve part of internal field from this fvBlockMatrix template void retrieveSolution ( @@ -220,12 +220,13 @@ public: ); //- Insert scalar equation coupling into this fvBlockMatrix - // Not tested: VV, 10/July/2014 + // Source compensation is done in function updateSourceCoupling() + // after all coupling terms are added. HJ, 27/Apr/2015 void insertEquationCoupling ( const direction dirI, const direction dirJ, - const fvScalarMatrix& matrix + const scalarField& coeffIJ ); //- Add field into block field @@ -237,20 +238,11 @@ public: ); //- Update coupling of block system. - // Subtracts the block-coefficient coupling as specified by the user - // from the source, leaving the implicit update given by - // linearisation + // Subtracts the block-coefficient coupling as specified by the + // user from the source, leaving the implicit update given by + // linearisation void updateSourceCoupling(); - //- Insert Picard tensor term that comes from Picard linearisation - // of convection term in momentum equation. VV, 21/July/2014. - void insertPicardTensor - ( - const direction UEqnDir, - const volVectorField& U, - const surfaceScalarField& phi - ); - // Solver calls for fvBlockMatrix From acb22cc5ccb3b7b32866204e624401c1f8aff0ed Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 27 Apr 2015 18:10:18 +0100 Subject: [PATCH 04/21] Formatting --- src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C | 1 + 1 file changed, 1 insertion(+) diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C index 857adc7ee..48761845a 100644 --- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C +++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C @@ -69,6 +69,7 @@ tmp kOmegaSST::F1(const volScalarField& CDkOmega) const return tanh(pow4(arg1)); } + tmp kOmegaSST::F2() const { volScalarField arg2 = min From 2e0286b2c1ab9e8d2560dc0980429a2b273d05e7 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 27 Apr 2015 18:11:12 +0100 Subject: [PATCH 05/21] Feature: coupled implicit k-epsilon model --- .../incompressible/RAS/Make/files | 2 + .../RAS/coupledKEpsilon/coupledKEpsilon.C | 323 ++++++++++++++++++ .../RAS/coupledKEpsilon/coupledKEpsilon.H | 175 ++++++++++ 3 files changed, 500 insertions(+) create mode 100644 src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.C create mode 100644 src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.H diff --git a/src/turbulenceModels/incompressible/RAS/Make/files b/src/turbulenceModels/incompressible/RAS/Make/files index 51a3bbfa1..fede0056b 100644 --- a/src/turbulenceModels/incompressible/RAS/Make/files +++ b/src/turbulenceModels/incompressible/RAS/Make/files @@ -18,6 +18,8 @@ NonlinearKEShih/NonlinearKEShih.C LienLeschzinerLowRe/LienLeschzinerLowRe.C LamBremhorstKE/LamBremhorstKE.C +coupledKEpsilon/coupledKEpsilon.C + /* Wall functions */ wallFunctions = derivedFvPatchFields/wallFunctions diff --git a/src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.C b/src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.C new file mode 100644 index 000000000..4050938dc --- /dev/null +++ b/src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.C @@ -0,0 +1,323 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +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 "coupledKEpsilon.H" +#include "fvBlockMatrix.H" +#include "addToRunTimeSelectionTable.H" + +#include "backwardsCompatibilityWallFunctions.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace incompressible +{ +namespace RASModels +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(coupledKEpsilon, 0); +addToRunTimeSelectionTable(RASModel, coupledKEpsilon, dictionary); + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +coupledKEpsilon::coupledKEpsilon +( + const volVectorField& U, + const surfaceScalarField& phi, + transportModel& lamTransportModel +) +: + RASModel(typeName, U, phi, lamTransportModel), + + Cmu_ + ( + dimensioned::lookupOrAddToDict + ( + "Cmu", + coeffDict_, + 0.09 + ) + ), + C1_ + ( + dimensioned::lookupOrAddToDict + ( + "C1", + coeffDict_, + 1.44 + ) + ), + C2_ + ( + dimensioned::lookupOrAddToDict + ( + "C2", + coeffDict_, + 1.92 + ) + ), + sigmaEps_ + ( + dimensioned::lookupOrAddToDict + ( + "sigmaEps", + coeffDict_, + 1.3 + ) + ), + + k_ + ( + IOobject + ( + "k", + runTime_.timeName(), + U_.db(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + autoCreateK("k", mesh_) + ), + epsilon_ + ( + IOobject + ( + "epsilon", + runTime_.timeName(), + U_.db(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + autoCreateEpsilon("epsilon", mesh_) + ), + nut_ + ( + IOobject + ( + "nut", + runTime_.timeName(), + U_.db(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + autoCreateNut("nut", mesh_) + ), + ke_ + ( + IOobject + ( + "ke", + runTime_.timeName(), + U_.db(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedVector2("zero", dimless, vector2::zero) + ) +{ + nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); + nut_ = min(nut_, nuRatio()*nu()); + nut_.correctBoundaryConditions(); + + printCoeffs(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +tmp coupledKEpsilon::R() const +{ + return tmp + ( + new volSymmTensorField + ( + IOobject + ( + "R", + runTime_.timeName(), + U_.db(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)), + k_.boundaryField().types() + ) + ); +} + + +tmp coupledKEpsilon::devReff() const +{ + return tmp + ( + new volSymmTensorField + ( + IOobject + ( + "devRhoReff", + runTime_.timeName(), + U_.db(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + -nuEff()*dev(twoSymm(fvc::grad(U_))) + ) + ); +} + + +tmp coupledKEpsilon::divDevReff(volVectorField& U) const +{ + return + ( + - fvm::laplacian(nuEff(), U) + - fvc::div(nuEff()*dev(fvc::grad(U)().T())) + ); +} + + +bool coupledKEpsilon::read() +{ + if (RASModel::read()) + { + Cmu_.readIfPresent(coeffDict()); + C1_.readIfPresent(coeffDict()); + C2_.readIfPresent(coeffDict()); + sigmaEps_.readIfPresent(coeffDict()); + + return true; + } + else + { + return false; + } +} + + +void coupledKEpsilon::correct() +{ + // Bound in case of topological change + // HJ, 22/Aug/2007 + if (mesh_.changing()) + { + bound(k_, k0_); + bound(epsilon_, epsilon0_); + } + + RASModel::correct(); + + if (!turbulence_) + { + return; + } + + + fvBlockMatrix keEqn(ke_); + + volScalarField G("RASModel::G", nut_*2*magSqr(symm(fvc::grad(U_)))); + + // Update epsilon and G at the wall + epsilon_.boundaryField().updateCoeffs(); + + // Dissipation equation + { + fvScalarMatrix epsEqn + ( + fvm::ddt(epsilon_) + + fvm::div(phi_, epsilon_) + + fvm::SuSp(-fvc::div(phi_), epsilon_) + - fvm::laplacian(DepsilonEff(), epsilon_) + == + C1_*G*epsilon_/k_ + - fvm::Sp(C2_*epsilon_/k_, epsilon_) + ); + + epsEqn.relax(); + epsEqn.completeAssembly(); + + keEqn.insertEquation(1, epsEqn); + } + + // Turbulent kinetic energy equation + { + fvScalarMatrix kEqn + ( + fvm::ddt(k_) + + fvm::div(phi_, k_) + + fvm::SuSp(-fvc::div(phi_), k_) + - fvm::laplacian(DkEff(), k_) + == + G + - fvm::Sp(epsilon_/k_, k_) + ); + + kEqn.relax(); + + keEqn.insertEquation(0, kEqn); + } + + // Add coupling term: C1*Cmu*(symm(grad(U))) k but with wall function + // corrections: must be calculated from G. HJ, 27/Apr/2015 + + // Add coupling term: epsilon source depends on k + // k, e sink terms cannot be changed because of boundedness + keEqn.insertEquationCoupling + ( + 1, 0, -C1_*G*epsilon_/sqr(k_) + ); + + // Update source coupling: coupling terms eliminated from source + keEqn.updateSourceCoupling(); + + keEqn.solve(); + + // Retrieve solution + keEqn.retrieveSolution(0, k_.internalField()); + keEqn.retrieveSolution(1, epsilon_.internalField()); + + k_.correctBoundaryConditions(); + epsilon_.correctBoundaryConditions(); + + bound(epsilon_, epsilon0_); + bound(k_, k0_); + + // Re-calculate viscosity + nut_ = Cmu_*sqr(k_)/epsilon_; + nut_ = min(nut_, nuRatio()*nu()); + nut_.correctBoundaryConditions(); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace RASModels +} // End namespace incompressible +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.H b/src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.H new file mode 100644 index 000000000..755e9a0b7 --- /dev/null +++ b/src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.H @@ -0,0 +1,175 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +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 + Foam::incompressible::RASModels::coupledKEpsilon + +Description + Standard k-epsilon turbulence model for incompressible flows using + a block-coupled solution. + + The default model coefficients correspond to the following: + @verbatim + coupledKEpsilonCoeffs + { + Cmu 0.09; + C1 1.44; + C2 1.92; + sigmaEps 1.3; + } + @endverbatim + +SourceFiles + coupledKEpsilon.C + +\*---------------------------------------------------------------------------*/ + +#ifndef coupledKEpsilon_H +#define coupledKEpsilon_H + +#include "RASModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace incompressible +{ +namespace RASModels +{ + +/*---------------------------------------------------------------------------*\ + Class coupledKEpsilon Declaration +\*---------------------------------------------------------------------------*/ + +class coupledKEpsilon +: + public RASModel +{ + // Private data + + // Model coefficients + + dimensionedScalar Cmu_; + dimensionedScalar C1_; + dimensionedScalar C2_; + dimensionedScalar sigmaEps_; + + + // Fields + + volScalarField k_; + volScalarField epsilon_; + volScalarField nut_; + + // Block vector field for k-epsilon + volVector2Field ke_; + + +public: + + //- Runtime type information + TypeName("coupledKEpsilon"); + + // Constructors + + //- Construct from components + coupledKEpsilon + ( + const volVectorField& U, + const surfaceScalarField& phi, + transportModel& transport + ); + + + //- Destructor + virtual ~coupledKEpsilon() + {} + + + // Member Functions + + //- Return the turbulence viscosity + virtual tmp nut() const + { + return nut_; + } + + //- Return the effective diffusivity for k + tmp DkEff() const + { + return tmp + ( + new volScalarField("DkEff", nut_ + nu()) + ); + } + + //- Return the effective diffusivity for epsilon + tmp DepsilonEff() const + { + return tmp + ( + new volScalarField("DepsilonEff", nut_/sigmaEps_ + nu()) + ); + } + + //- Return the turbulence kinetic energy + virtual tmp k() const + { + return k_; + } + + //- Return the turbulence kinetic energy dissipation rate + virtual tmp epsilon() const + { + return epsilon_; + } + + //- Return the Reynolds stress tensor + virtual tmp R() const; + + //- Return the effective stress tensor including the laminar stress + virtual tmp devReff() const; + + //- Return the source term for the momentum equation + virtual tmp divDevReff(volVectorField& U) const; + + //- Solve the turbulence equations and correct the turbulence viscosity + virtual void correct(); + + //- Read RASProperties dictionary + virtual bool read(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace RASModels +} // End namespace incompressible +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // From 75259a047bf3c0397d74811bf6bd74e0f7b06179 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 27 Apr 2015 18:12:34 +0100 Subject: [PATCH 06/21] Feature: block coupling rewrite --- .../blockCoupledScalarTransportFoam.C | 59 +++++-------------- .../createFields.H | 9 --- 2 files changed, 14 insertions(+), 54 deletions(-) diff --git a/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C b/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C index cc2b1f7e0..25b389b87 100644 --- a/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C +++ b/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C @@ -38,13 +38,8 @@ Description #include "fieldTypes.H" #include "Time.H" #include "fvMesh.H" +#include "fvBlockMatrix.H" -#include "blockLduSolvers.H" -#include "VectorNFieldTypes.H" -#include "volVectorNFields.H" -#include "blockVectorNMatrices.H" - -#include "blockMatrixTools.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -92,53 +87,27 @@ int main(int argc, char *argv[]) TsEqn.relax(); // Prepare block system - BlockLduMatrix blockM(mesh); - - //- Transfer the coupled interface list for processor/cyclic/etc. - // boundaries - blockM.interfaces() = blockT.boundaryField().blockInterfaces(); - - // Grab block diagonal and set it to zero - Field& d = blockM.diag().asSquare(); - d = tensor2::zero; - - // Grab linear off-diagonal and set it to zero - Field& l = blockM.lower().asLinear(); - Field& u = blockM.upper().asLinear(); - u = vector2::zero; - l = vector2::zero; - - vector2Field& blockX = blockT.internalField(); - vector2Field blockB(mesh.nCells(), vector2::zero); + fvBlockMatrix blockM(blockT); //- Inset equations into block Matrix - blockMatrixTools::insertEquation(0, TEqn, blockM, blockX, blockB); - blockMatrixTools::insertEquation(1, TsEqn, blockM, blockX, blockB); + blockM.insertEquation(0, TEqn); + blockM.insertEquation(1, TsEqn); - //- Add off-diagonal terms and remove from block source - forAll(d, i) - { - d[i](0, 1) = -alpha.value()*mesh.V()[i]; - d[i](1, 0) = -alpha.value()*mesh.V()[i]; + //- Add off-diagonal coupling terms + scalarField coupling(mesh.nCells(), -alpha.value()); - blockB[i][0] -= alpha.value()*blockX[i][1]*mesh.V()[i]; - blockB[i][1] -= alpha.value()*blockX[i][0]*mesh.V()[i]; - } + blockM.insertEquationCoupling(0, 1, coupling); + blockM.insertEquationCoupling(1, 0, coupling); + + // Update source coupling: coupling terms eliminated from source + blockM.updateSourceCoupling(); //- Block coupled solver call - BlockSolverPerformance solverPerf = - BlockLduSolver::New - ( - blockT.name(), - blockM, - mesh.solutionDict().solver(blockT.name()) - )->solve(blockX, blockB); - - solverPerf.print(); + blockM.solve(); // Retrieve solution - blockMatrixTools::blockRetrieve(0, T.internalField(), blockX); - blockMatrixTools::blockRetrieve(1, Ts.internalField(), blockX); + blockM.retrieveSolution(0, T.internalField()); + blockM.retrieveSolution(1, Ts.internalField()); T.correctBoundaryConditions(); Ts.correctBoundaryConditions(); diff --git a/applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H b/applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H index 2f4eb6aaa..d3d4b42a8 100644 --- a/applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H +++ b/applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H @@ -44,15 +44,6 @@ dimensionedVector2("zero", dimless, vector2::zero) ); - { - vector2Field& blockX = blockT.internalField(); - - blockMatrixTools::blockInsert(0, T.internalField(), blockX); - blockMatrixTools::blockInsert(1, Ts.internalField(), blockX); - } - - blockT.write(); - Info<< "Reading field U\n" << endl; volVectorField U From de605a3b04869983b7906c8893032a9686cea855 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 27 Apr 2015 18:13:14 +0100 Subject: [PATCH 07/21] Feature: remove double continuity check --- applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C | 2 -- 1 file changed, 2 deletions(-) diff --git a/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C b/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C index 9ee0e7f41..f9fb25c1d 100644 --- a/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C +++ b/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C @@ -93,8 +93,6 @@ int main(int argc, char *argv[]) phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource; -# include "continuityErrs.H" - // Make flux relative in rotating zones mrfZones.relativeFlux(phi); From fd1921f95543024c6b51e0d9022cf08bb7d15fb4 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 27 Apr 2015 18:13:30 +0100 Subject: [PATCH 08/21] Feature: Update build script --- applications/solvers/coupled/Allwmake | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/applications/solvers/coupled/Allwmake b/applications/solvers/coupled/Allwmake index 9ce03ddae..e6320de43 100755 --- a/applications/solvers/coupled/Allwmake +++ b/applications/solvers/coupled/Allwmake @@ -8,6 +8,6 @@ wmake blockCoupledScalarTransportFoam wmake conjugateHeatFoam wmake conjugateHeatSimpleFoam wmake pUCoupledFoam -wmake pUCoupledFullPicardFoam -wmake pUCoupledSemiPicardFoam +wmake MRFPorousFoam + # ----------------------------------------------------------------- end-of-file From 8a8afda6ac3b94f4127d6e72db0ec0240f599ad2 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 28 Apr 2015 10:19:40 +0100 Subject: [PATCH 09/21] Bugfix: compiler error typename outside of template --- .../basic/basicSymmetry/basicSymmetryFvPatchScalarField.C | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/finiteVolume/fields/fvPatchFields/basic/basicSymmetry/basicSymmetryFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/basic/basicSymmetry/basicSymmetryFvPatchScalarField.C index 97bb0a32f..fc0670bd8 100644 --- a/src/finiteVolume/fields/fvPatchFields/basic/basicSymmetry/basicSymmetryFvPatchScalarField.C +++ b/src/finiteVolume/fields/fvPatchFields/basic/basicSymmetry/basicSymmetryFvPatchScalarField.C @@ -50,7 +50,7 @@ void basicSymmetryFvPatchField::evaluate(const Pstream::commsTypes) { // Local typedefs typedef scalar Type; - typedef typename outerProduct::type gradType; + typedef outerProduct::type gradType; typedef GeometricField gradFieldType; if (!updated()) @@ -140,7 +140,7 @@ tmp basicSymmetryFvPatchField::snGrad() const { // Local typedefs typedef vector Type; - typedef typename outerProduct::type gradType; + typedef outerProduct::type gradType; typedef GeometricField gradFieldType; vectorField nHat = this->patch().nf(); @@ -213,7 +213,7 @@ void basicSymmetryFvPatchField::evaluate(const Pstream::commsTypes) { // Local typedefs typedef vector Type; - typedef typename outerProduct::type gradType; + typedef outerProduct::type gradType; typedef GeometricField gradFieldType; if (!updated()) From a863a5fc19d0aa3bce4a845c92a9ba84c4b9febf Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 29 Apr 2015 10:17:25 +0100 Subject: [PATCH 10/21] Feature: coupled solver residual scripts --- bin/plotForces.py | 62 +++++++++++++++++++++++++ bin/plotResCoupled.py | 79 ++++++++++++++++++++++++++++++++ bin/plotResidual.py | 103 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 244 insertions(+) create mode 100755 bin/plotForces.py create mode 100755 bin/plotResCoupled.py create mode 100755 bin/plotResidual.py diff --git a/bin/plotForces.py b/bin/plotForces.py new file mode 100755 index 000000000..910f367c1 --- /dev/null +++ b/bin/plotForces.py @@ -0,0 +1,62 @@ +#!/usr/bin/python + +import re +forceRegex=r"([0-9.Ee\-+]+)\s+\(+([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\)\s\(([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\s([0-9 .Ee\-+]+)\)+\s\(+([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\)\s\(([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\)+" +t = [] +fpx = []; fpy = []; fpz = [] +fvx = []; fvy = []; fvz = [] +mpx = []; mpy = []; mpz = [] +mvx = []; mvy = []; mvz = [] +pipefile=open('forces/0/forces.dat','r') +lines = pipefile.readlines() +for line in lines: + match=re.search(forceRegex,line) + if match: + t.append(float(match.group(1))) + fpx.append(float(match.group(2))) + fpy.append(float(match.group(3))) + fpz.append(float(match.group(4))) + fvx.append(float(match.group(5))) + fvy.append(float(match.group(6))) + fvz.append(float(match.group(7))) + mpx.append(float(match.group(8))) + mpy.append(float(match.group(9))) + mpz.append(float(match.group(10))) + mvx.append(float(match.group(11))) + mvy.append(float(match.group(12))) + mvz.append(float(match.group(13))) + + +# combine pressure and viscous forces and moments +fx = fpx +fy = fpy +fz = fpz + +mx = mpx +my = mpy +mz = mpz + +for i in range(1,len(t)): + fx[i] += fvx[i] + fy[i] += fvy[i] + fz[i] += fvz[i] + mx[i] += mvx[i] + my[i] += mvy[i] + mz[i] += mvz[i] + +# plot forces +import pylab +pylab.xlabel('iteration') +pylab.ylabel('force (N)') +pylab.grid(True) +# +pylab.plot(t,fx,'-',label="fx") +pylab.plot(t,fy,'-',label="fy") +pylab.plot(t,fz,'-',label="fz") + +#pylab.plot(t,mx,'-',label="mx") +#pylab.plot(t,my,'-',label="my") +#pylab.plot(t,mz,'-',label="mz") + +pylab.legend() +pylab.show() diff --git a/bin/plotResCoupled.py b/bin/plotResCoupled.py new file mode 100755 index 000000000..7050c7e51 --- /dev/null +++ b/bin/plotResCoupled.py @@ -0,0 +1,79 @@ +#!/usr/bin/python + +import sys +if len(sys.argv) != 2: + print 'script requires name of log file' + sys.exit() + +logfilename = sys.argv[1] +print 'Reading file', logfilename + +import re +UpRegex=r"([A-Z,a-z]*):*.*Solving for Up, Initial residual = \(([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\), Final residual = \(([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\), No Iterations ([0-9]*)" +komegaRegex=r"([A-Z,a-z]*):*.*Solving for kOmega, Initial residual = \(([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\), Final residual = \(([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\), No Iterations ([0-9]*)" + +tUp = [] +Ux = [] +Uy = [] +Uz = [] +p = [] +iUp = 0 + +tkomega = [] +k = [] +omega = [] +ikomega = 0 + +#HJ take name of log file as script argument +pipefile=open(logfilename,'r') +lines = pipefile.readlines() + +for line in lines: + matchUp=re.search(UpRegex,line) + if matchUp: + iUp = iUp + 1 + tUp.append(iUp) + Ux.append(float(matchUp.group(2))) + Uy.append(float(matchUp.group(3))) + Uz.append(float(matchUp.group(4))) + p.append(float(matchUp.group(5))) + matchkomega=re.search(komegaRegex,line) + if matchkomega: + ikomega = ikomega + 1 + tkomega.append(ikomega) + k.append(float(matchkomega.group(2))) + omega.append(float(matchkomega.group(3))) + +outfile=open('residual.dat','w') + +print 'hits = ', ikomega + +#HJ need better way of combining lists +if ikomega > 0: + for index in range(0,ikomega): + outfile.write(str(tUp[index])+' '+str(Ux[index])+' '+str(Uy[index])+' '+str(Uz[index])+' '+str(p[index])+' '+str(k[index])+' '+str(omega[index])+'\n') +elif iUp > 0: + for index in range(0,iUp): + outfile.write(str(tUp[index])+' '+str(Ux[index])+' '+str(Uy[index])+' '+str(Uz[index])+' '+str(p[index])+'\n') + +outfile.close() + +# prepare plot +import pylab +pylab.xlabel('iteration') +pylab.ylabel('residual') +pylab.grid(True) + +# plot in log scale +if iUp > 0: + pylab.semilogy(tUp,Ux,'-',label="Ux") + pylab.semilogy(tUp,Uy,'-',label="Uy") + pylab.semilogy(tUp,Uz,'-',label="Uz") + pylab.semilogy(tUp,p,'-',label="p") + +if ikomega > 0: + pylab.semilogy(tkomega,k,'-',label="k") + pylab.semilogy(tkomega,omega,'-',label="omega") + +pylab.legend() +pylab.show() diff --git a/bin/plotResidual.py b/bin/plotResidual.py new file mode 100755 index 000000000..53d650415 --- /dev/null +++ b/bin/plotResidual.py @@ -0,0 +1,103 @@ +#!/usr/bin/python + +import sys +if len(sys.argv) != 2: + print 'script requires name of log file' + sys.exit() + +logfilename = sys.argv[1] +print 'Reading file', logfilename + +import re +UpRegex=r"([A-Z,a-z]*):*.*Solving for Up, Initial residual = \(([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\), Final residual = \(([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\s([0-9.Ee\-+]*)\), No Iterations ([0-9]*)" +kRegex=r"([A-Z,a-z]*):*.*Solving for k, Initial residual = ([0-9.Ee\-+]*), Final residual = ([0-9.Ee\-+]*), No Iterations ([0-9]*)" +omegaRegex=r"([A-Z,a-z]*):*.*Solving for omega, Initial residual = ([0-9.Ee\-+]*), Final residual = ([0-9.Ee\-+]*), No Iterations ([0-9]*)" +epsilonRegex=r"([A-Z,a-z]*):*.*Solving for epsilon, Initial residual = ([0-9.Ee\-+]*), Final residual = ([0-9.Ee\-+]*), No Iterations ([0-9]*)" + +tUp = [] +Ux = [] +Uy = [] +Uz = [] +p = [] +iUp = 0 + +tk = [] +k = [] +ik = 0 + +tomega = [] +omega = [] +iomega = 0 + +tepsilon = [] +epsilon = [] +iepsilon = 0 + +#HJ take name of log file as script argument +pipefile=open(logfilename,'r') +lines = pipefile.readlines() + +for line in lines: + matchUp=re.search(UpRegex,line) + if matchUp: + iUp = iUp + 1 + tUp.append(iUp) + Ux.append(float(matchUp.group(2))) + Uy.append(float(matchUp.group(3))) + Uz.append(float(matchUp.group(4))) + p.append(float(matchUp.group(5))) + matchk=re.search(kRegex,line) + if matchk: + ik = ik + 1 + tk.append(ik) + k.append(float(matchk.group(2))) + matchomega=re.search(omegaRegex,line) + if matchomega: + iomega = iomega + 1 + tomega.append(iomega) + omega.append(float(matchomega.group(2))) + matchepsilon=re.search(epsilonRegex,line) + if matchepsilon: + iepsilon = iepsilon + 1 + tepsilon.append(iepsilon) + epsilon.append(float(matchepsilon.group(2))) + +outfile=open('residual.dat','w') + +#HJ need better way of combining lists +if iomega > 0: + for index in range(0,iomega): + outfile.write(str(tUp[index])+' '+str(Ux[index])+' '+str(Uy[index])+' '+str(Uz[index])+' '+str(p[index])+' '+str(k[index])+' '+str(omega[index])+'\n') +elif iepsilon > 0: + for index in range(0,iepsilon): + outfile.write(str(tUp[index])+' '+str(Ux[index])+' '+str(Uy[index])+' '+str(Uz[index])+' '+str(p[index])+' '+str(k[index])+' '+str(epsilon[index])+'\n') +elif iUp > 0: + for index in range(0,iUp): + outfile.write(str(tUp[index])+' '+str(Ux[index])+' '+str(Uy[index])+' '+str(Uz[index])+' '+str(p[index])+'\n') + +outfile.close() + +# prepare plot +import pylab +pylab.xlabel('iteration') +pylab.ylabel('residual') +pylab.grid(True) + +# plot in log scale +if iUp > 0: + pylab.semilogy(tUp,Ux,'-',label="Ux") + pylab.semilogy(tUp,Uy,'-',label="Uy") + pylab.semilogy(tUp,Uz,'-',label="Uz") + pylab.semilogy(tUp,p,'-',label="p") + +if ik > 0: + pylab.semilogy(tk,k,'-',label="k") + +if iomega > 0: + pylab.semilogy(tomega,omega,'-',label="omega") + +if iepsilon > 0: + pylab.semilogy(tepsilon,epsilon,'-',label="epsilon") + +pylab.legend() +pylab.show() From c080289664410ee6dc86f49b58de7097fbfeca96 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 29 Apr 2015 13:58:36 +0100 Subject: [PATCH 11/21] Bugfix: formatting --- .../constant/viscoelasticProperties | 91 +++++++++---------- 1 file changed, 45 insertions(+), 46 deletions(-) diff --git a/tutorials/viscoelastic/viscoelasticFluidFoam/PTT-Exponential/constant/viscoelasticProperties b/tutorials/viscoelastic/viscoelasticFluidFoam/PTT-Exponential/constant/viscoelasticProperties index ac13aa608..0dacfc5a5 100644 --- a/tutorials/viscoelastic/viscoelasticFluidFoam/PTT-Exponential/constant/viscoelasticProperties +++ b/tutorials/viscoelastic/viscoelasticFluidFoam/PTT-Exponential/constant/viscoelasticProperties @@ -33,56 +33,55 @@ Tr = 170 C; activation energy: E0 = 48.2 kJ/mol; WLF-shift parameters: C1 = 14.3 rheology { + type multiMode; - type multiMode; + models + ( + first + { + type PTT-Exponential; + rho rho [1 -3 0 0 0 0 0] 850; + etaS etaS [1 -1 -1 0 0 0 0] 0.0; + etaP etaP [1 -1 -1 0 0 0 0] 280.43457; + lambda lambda [0 0 1 0 0 0 0] 3.8946e-3; + epsilon epsilon [0 0 0 0 0 0 0] 0.3; + zeta zeta [0 0 0 0 0 0 0] 0.08; + } - models - ( - first - { - type PTT-Exponential; - rho rho [1 -3 0 0 0 0 0] 850; - etaS etaS [1 -1 -1 0 0 0 0] 0.0; - etaP etaP [1 -1 -1 0 0 0 0] 280.43457; - lambda lambda [0 0 1 0 0 0 0] 3.8946e-3; - epsilon epsilon [0 0 0 0 0 0 0] 0.3; - zeta zeta [0 0 0 0 0 0 0] 0.08; - } + second + { + type PTT-Exponential; + rho rho [1 -3 0 0 0 0 0] 850; + etaS etaS [1 -1 -1 0 0 0 0] 0.0; + etaP etaP [1 -1 -1 0 0 0 0] 810.4203; + lambda lambda [0 0 1 0 0 0 0] 5.1390e-2; + epsilon epsilon [0 0 0 0 0 0 0] 0.2; + zeta zeta [0 0 0 0 0 0 0] 0.08; + } - second - { - type PTT-Exponential; - rho rho [1 -3 0 0 0 0 0] 850; - etaS etaS [1 -1 -1 0 0 0 0] 0.0; - etaP etaP [1 -1 -1 0 0 0 0] 810.4203; - lambda lambda [0 0 1 0 0 0 0] 5.1390e-2; - epsilon epsilon [0 0 0 0 0 0 0] 0.2; - zeta zeta [0 0 0 0 0 0 0] 0.08; - } - - third - { - type PTT-Exponential; - rho rho [1 -3 0 0 0 0 0] 850; - etaS etaS [1 -1 -1 0 0 0 0] 0.0; - etaP etaP [1 -1 -1 0 0 0 0] 1678.6357; - lambda lambda [0 0 1 0 0 0 0] 5.0349e-1; - epsilon epsilon [0 0 0 0 0 0 0] 0.02; - zeta zeta [0 0 0 0 0 0 0] 0.08; - } - - fourth - { - type PTT-Exponential; - rho rho [1 -3 0 0 0 0 0] 850; - etaS etaS [1 -1 -1 0 0 0 0] 0.0; - etaP etaP [1 -1 -1 0 0 0 0] 1381.0029; - lambda lambda [0 0 1 0 0 0 0] 4.5911e0; - epsilon epsilon [0 0 0 0 0 0 0] 0.02; - zeta zeta [0 0 0 0 0 0 0] 0.08; - } - ); + third + { + type PTT-Exponential; + rho rho [1 -3 0 0 0 0 0] 850; + etaS etaS [1 -1 -1 0 0 0 0] 0.0; + etaP etaP [1 -1 -1 0 0 0 0] 1678.6357; + lambda lambda [0 0 1 0 0 0 0] 5.0349e-1; + epsilon epsilon [0 0 0 0 0 0 0] 0.02; + zeta zeta [0 0 0 0 0 0 0] 0.08; + } + fourth + { + type PTT-Exponential; + rho rho [1 -3 0 0 0 0 0] 850; + etaS etaS [1 -1 -1 0 0 0 0] 0.0; + etaP etaP [1 -1 -1 0 0 0 0] 1381.0029; + lambda lambda [0 0 1 0 0 0 0] 4.5911e0; + epsilon epsilon [0 0 0 0 0 0 0] 0.02; + zeta zeta [0 0 0 0 0 0 0] 0.08; + } + ); } + // ************************************************************************* // From a6d1c4021fff47324f4943fea508aa18a6198c4e Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 29 Apr 2015 13:59:22 +0100 Subject: [PATCH 12/21] Bug fix: Maxwell transport. M. Nobrega --- .../viscoelastic/viscoelasticLaws/Maxwell/Maxwell.C | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/transportModels/viscoelastic/viscoelasticLaws/Maxwell/Maxwell.C b/src/transportModels/viscoelastic/viscoelasticLaws/Maxwell/Maxwell.C index c53a75326..5beb5dbdd 100644 --- a/src/transportModels/viscoelastic/viscoelasticLaws/Maxwell/Maxwell.C +++ b/src/transportModels/viscoelastic/viscoelasticLaws/Maxwell/Maxwell.C @@ -77,7 +77,6 @@ Foam::tmp Foam::Maxwell::divTau(volVectorField& U) const - fvc::laplacian(etaPEff/rho_, U, "laplacian(etaPEff,U)") + fvm::laplacian( (etaPEff + etaS_)/rho_, U, "laplacian(etaPEff+etaS,U)") ); - } @@ -86,6 +85,9 @@ void Foam::Maxwell::correct() // Velocity gradient tensor volTensorField L = fvc::grad(U()); + // Convected derivate term + volTensorField C = tau_ & L; + // Twice the rate of deformation tensor volSymmTensorField twoD = twoSymm(L); @@ -93,9 +95,11 @@ void Foam::Maxwell::correct() fvSymmTensorMatrix tauEqn ( fvm::ddt(tau_) + + fvm::div(phi(), tau_) == etaP_/lambda_*twoD - - fvm::Sp( 1/lambda_, tau_ ) + + twoSymm(C) + - fvm::Sp( 1/lambda_, tau_) ); tauEqn.relax(); From f19018541ee762207c156c96217ce685742a5477 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 30 Apr 2015 14:49:40 +0100 Subject: [PATCH 13/21] Formatting --- .../blockVectorNLduInterfaceFields.C | 10 +++++----- .../processorBlockLduInterfaceField.H | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField/blockVectorNLduInterfaceFields.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField/blockVectorNLduInterfaceFields.C index 5b340677e..a3c8fa30a 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField/blockVectorNLduInterfaceFields.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField/blockVectorNLduInterfaceFields.C @@ -28,7 +28,7 @@ Description Author -\*----------------------------------------------------------------------------*/ +\*---------------------------------------------------------------------------*/ #include "BlockLduInterfaceField.H" #include "processorBlockLduInterfaceField.H" @@ -41,10 +41,10 @@ namespace Foam // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -#define makeTemplateTypeNameAndDebug(type, Type, args...) \ - defineTemplateTypeNameAndDebug(BlockLduInterfaceField, 0); \ - \ - defineTemplateTypeNameAndDebug(processorBlockLduInterfaceField, 0); \ +#define makeTemplateTypeNameAndDebug(type, Type, args...) \ + defineTemplateTypeNameAndDebug(BlockLduInterfaceField, 0); \ + \ + defineTemplateTypeNameAndDebug(processorBlockLduInterfaceField, 0); forAllVectorNTypes(makeTemplateTypeNameAndDebug); diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField/processorBlockLduInterfaceField/processorBlockLduInterfaceField.H b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField/processorBlockLduInterfaceField/processorBlockLduInterfaceField.H index e3487dc08..7af2ef84a 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField/processorBlockLduInterfaceField/processorBlockLduInterfaceField.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField/processorBlockLduInterfaceField/processorBlockLduInterfaceField.H @@ -48,7 +48,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class processorBlockLduInterfaceField Declaration + Class processorBlockLduInterfaceField Declaration \*---------------------------------------------------------------------------*/ template From 60f060da52c8c3108c025ed61b3de5d225da7e2f Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 1 May 2015 13:16:50 +0100 Subject: [PATCH 14/21] Update tolerances for SP and LDP --- .../iglooWithFridges/system/controlDict | 2 +- .../wingMotion/wingMotion_snappyHexMesh/system/controlDict | 2 +- .../incompressible/simpleFoam/motorBike/system/controlDict | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/controlDict b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/controlDict index 5bee4900a..c38b5c728 100644 --- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/controlDict +++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/controlDict @@ -37,7 +37,7 @@ purgeWrite 0; writeFormat ascii; -writePrecision 6; +writePrecision 10; writeCompression uncompressed; diff --git a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/system/controlDict b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/system/controlDict index 95f30f5c7..0c9467f78 100644 --- a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/system/controlDict +++ b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/system/controlDict @@ -35,7 +35,7 @@ purgeWrite 0; writeFormat ascii; -writePrecision 6; +writePrecision 10; writeCompression uncompressed; diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict index a3bbc84aa..2ed2dc070 100644 --- a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict +++ b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict @@ -36,7 +36,7 @@ purgeWrite 0; writeFormat ascii; -writePrecision 6; +writePrecision 10; writeCompression compressed; From 7eefef1a1062b55894557763c252f21816246ceb Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 1 May 2015 13:17:11 +0100 Subject: [PATCH 15/21] Bugfix: Removed --- .../mesh/cfMesh/cartesianMesh/ship5415Octree/ship5415Octree.foam | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 tutorials/mesh/cfMesh/cartesianMesh/ship5415Octree/ship5415Octree.foam diff --git a/tutorials/mesh/cfMesh/cartesianMesh/ship5415Octree/ship5415Octree.foam b/tutorials/mesh/cfMesh/cartesianMesh/ship5415Octree/ship5415Octree.foam deleted file mode 100644 index e69de29bb..000000000 From 51cd34cc0e14ecee3f6dda64411d59f5ca5a9201 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 1 May 2015 13:18:09 +0100 Subject: [PATCH 16/21] Formatting --- .../fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.C | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.C index 015f56640..42a9b4a53 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.C @@ -315,8 +315,7 @@ void cyclicGgiFvPatchField::updateInterfaceMatrix const CoeffField&, const Pstream::commsTypes commsType ) const -{ -} +{} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // From 6b022758d1b15a8d08718a78d3f68879e95bcf90 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 1 May 2015 13:21:04 +0100 Subject: [PATCH 17/21] Feature: Single precision and long double precision port --- .../ersPointSource/ersPointSource.C | 2 +- .../contactStressFoam/contactPatchPair.C | 12 +- .../newContactStressFoam/contactPatchPair.C | 2 +- .../elasticAcpSolidFoam/updateCrack.H | 50 ++- .../elasticIncrAcpSolidFoam/updateCrack.H | 130 ++++-- .../elasticOrthoAcpSolidFoam/updateCrack.H | 409 ++++++++++-------- .../testFadField/testFadField.C | 2 +- etc/bashrc | 2 +- src/Pstream/mpi/Pstream.C | 2 + src/cudaSolvers/cudaSolver/cudaTypes.H | 14 +- src/dbns/multigrid/mgMeshLevel/mgMeshLevel.C | 43 ++ src/foam/Make/files | 1 + src/foam/db/IOstreams/Fstreams/IFstream.H | 4 +- src/foam/db/IOstreams/IOstreams/Istream.H | 7 +- src/foam/db/IOstreams/IOstreams/Ostream.H | 7 +- src/foam/db/IOstreams/Pstreams/IPstream.C | 22 + src/foam/db/IOstreams/Pstreams/IPstream.H | 5 +- src/foam/db/IOstreams/Pstreams/OPstream.C | 8 + src/foam/db/IOstreams/Pstreams/OPstream.H | 5 +- src/foam/db/IOstreams/Sstreams/ISstream.C | 8 + src/foam/db/IOstreams/Sstreams/ISstream.H | 9 +- src/foam/db/IOstreams/Sstreams/OSstream.C | 8 + src/foam/db/IOstreams/Sstreams/OSstream.H | 7 +- .../db/IOstreams/Sstreams/prefixOSstream.C | 7 + .../db/IOstreams/Sstreams/prefixOSstream.H | 7 +- .../IOstreams/StringStreams/IStringStream.H | 8 +- .../IOstreams/StringStreams/OStringStream.H | 4 +- src/foam/db/IOstreams/Tstreams/ITstream.C | 7 + src/foam/db/IOstreams/Tstreams/ITstream.H | 7 +- src/foam/db/IOstreams/hashes/OSHA1stream.H | 4 +- src/foam/db/IOstreams/token/token.H | 11 + src/foam/db/IOstreams/token/tokenI.H | 70 ++- src/foam/db/IOstreams/token/tokenIO.C | 12 + src/foam/db/Time/Time.C | 2 +- .../processor/ProcessorPointPatchField.C | 4 +- src/foam/meshes/meshShapes/face/face.C | 4 +- .../primitiveShapes/triangle/triangleI.H | 10 +- src/foam/primitives/Scalar/Scalar.H | 4 +- src/foam/primitives/Scalar/doubleFloat.H | 7 + .../longDoubleScalar/longDoubleScalar.C | 43 ++ .../longDoubleScalar/longDoubleScalar.H | 115 +++++ src/foam/primitives/Scalar/scalar/scalar.H | 19 + src/foam/primitives/Tensor/tensor/tensor.C | 4 +- src/foam/primitives/transform/symmTransform.H | 291 +++++++------ src/foam/primitives/transform/transform.H | 267 ++++++------ .../MGridGenGAMGAgglomerate.C | 42 ++ .../solidCohesiveFvPatchVectorField.C | 2 +- .../dirichletNeumann/dirichletNeumann.C | 19 +- 48 files changed, 1129 insertions(+), 600 deletions(-) create mode 100644 src/foam/primitives/Scalar/longDoubleScalar/longDoubleScalar.C create mode 100644 src/foam/primitives/Scalar/longDoubleScalar/longDoubleScalar.H diff --git a/applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C b/applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C index f8b206c57..b66a430df 100644 --- a/applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C +++ b/applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C @@ -46,7 +46,7 @@ ersPointSource::ersPointSource alpha_(readScalar(dict.lookup("alpha"))), direction_(dict.lookup("direction")) { - q() = -alpha_*qmax_*min(direction_ & p.Sf()/p.magSf(), 0.0); + q() = -alpha_*qmax_*min(direction_ & p.Sf()/p.magSf(), scalar(0)); } diff --git a/applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C b/applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C index 4804383a1..2af8a4675 100644 --- a/applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C +++ b/applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C @@ -30,18 +30,10 @@ Description #include "contactPatchPair.H" #include "volFields.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // Construct from components -contactPatchPair::contactPatchPair + Foam::contactPatchPair::contactPatchPair ( const fvMesh& m, const label master, @@ -75,6 +67,4 @@ contactPatchPair::contactPatchPair {} -} // End namespace Foam - // ************************************************************************* // diff --git a/applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C b/applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C index 9e5cf7992..59a70f0e4 100644 --- a/applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C +++ b/applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C @@ -302,7 +302,7 @@ void Foam::contactPatchPair::correct ( slavePressure ), - 0.0 + scalar(0) ); // Calculate master traction, using the positive part of diff --git a/applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H b/applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H index 50bb49ddd..a9f4c0b83 100644 --- a/applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H +++ b/applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H @@ -6,13 +6,15 @@ nCoupledFacesToBreak = 0; // scalarField effTraction = // cohesiveZone.internalField() * // mag(traction.internalField()); - scalarField normalTraction = - cohesiveZone.internalField() * - ( n.internalField() & traction.internalField() ); - normalTraction = max(normalTraction, 0.0); // only consider tensile tractions - scalarField shearTraction = - cohesiveZone.internalField() * - mag( (I - Foam::sqr(n.internalField())) & traction.internalField() ); + scalarField normalTraction = + cohesiveZone.internalField()* + ( n.internalField() & traction.internalField() ); + + // only consider tensile tractions + normalTraction = max(normalTraction, scalar(0)); + scalarField shearTraction = + cohesiveZone.internalField() * + mag( (I - Foam::sqr(n.internalField())) & traction.internalField() ); // the traction fraction is monitored to decide which faces to break: // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face @@ -27,16 +29,18 @@ nCoupledFacesToBreak = 0; scalarField effTractionFraction(normalTraction.size(), 0.0); if(cohesivePatchUPtr) - { - effTractionFraction = - (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/tauMaxI)*(shearTraction/tauMaxI); - } + { + effTractionFraction = + (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + + (shearTraction/tauMaxI)*(shearTraction/tauMaxI); + } else - { - // solidCohesiveFixedModeMix only uses sigmaMax - effTractionFraction = - (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI); - } + { + // solidCohesiveFixedModeMix only uses sigmaMax + effTractionFraction = + (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI); + } maxEffTractionFraction = gMax(effTractionFraction); SLList