From 65cd7ee3f30bd7328c1ab6c5e4a57c709aa6dfe0 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 3 Oct 2014 15:28:15 +0100 Subject: [PATCH 01/45] Formatting --- src/solidModels/constitutiveModel/constitutiveModel.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solidModels/constitutiveModel/constitutiveModel.H b/src/solidModels/constitutiveModel/constitutiveModel.H index 3c1a10b07..bf9155b7b 100644 --- a/src/solidModels/constitutiveModel/constitutiveModel.H +++ b/src/solidModels/constitutiveModel/constitutiveModel.H @@ -116,8 +116,8 @@ public: //- Construct from components constitutiveModel ( - const volSymmTensorField& sigma, - const volVectorField& U + const volSymmTensorField& sigma, + const volVectorField& U ); From 7ba34de8e47e35e160103afd86851b4667991e3f Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 15 Oct 2014 11:22:29 +0100 Subject: [PATCH 02/45] Formatting --- .../postProcessing/stressField/solidStress/solidStress.C | 7 +++---- .../OutputFilterFunctionObject.C | 1 + src/sampling/probes/probes.C | 1 + 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/applications/utilities/postProcessing/stressField/solidStress/solidStress.C b/applications/utilities/postProcessing/stressField/solidStress/solidStress.C index 432cc283d..2bfa1f6d4 100644 --- a/applications/utilities/postProcessing/stressField/solidStress/solidStress.C +++ b/applications/utilities/postProcessing/stressField/solidStress/solidStress.C @@ -36,7 +36,6 @@ Description int main(int argc, char *argv[]) { - # include "addTimeOptions.H" # include "setRootCase.H" @@ -53,7 +52,7 @@ int main(int argc, char *argv[]) # include "createMesh.H" # include "readMechanicalProperties.H" - for (label i=startTime; i::destroyFilter() ptr_.reset(); } + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template diff --git a/src/sampling/probes/probes.C b/src/sampling/probes/probes.C index b6864cb6d..2596f5ac4 100644 --- a/src/sampling/probes/probes.C +++ b/src/sampling/probes/probes.C @@ -323,6 +323,7 @@ void Foam::probes::write() { // Check if the mesh is changing and if so, resample const fvMesh& mesh = refCast(obr_); + if (mesh.moving() || mesh.changing()) { cellList_.clear(); From ec21448e80511f5905aa4b311627dd9e23a04389 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 15 Oct 2014 15:11:57 +0100 Subject: [PATCH 03/45] Added turbulence modelling terms into fvSchemes --- tutorials/compressible/dbnsTurbFoam/naca0012/system/fvSchemes | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tutorials/compressible/dbnsTurbFoam/naca0012/system/fvSchemes b/tutorials/compressible/dbnsTurbFoam/naca0012/system/fvSchemes index 33d808a1f..9248dfd35 100644 --- a/tutorials/compressible/dbnsTurbFoam/naca0012/system/fvSchemes +++ b/tutorials/compressible/dbnsTurbFoam/naca0012/system/fvSchemes @@ -32,11 +32,14 @@ divSchemes div(phi,epsilon) Gauss upwind; div(devRhoReff) Gauss linear; div((devRhoReff&U)) Gauss linear; + + div((muEff*dev2(grad(U).T()))) Gauss linear; } laplacianSchemes { default none; + laplacian(muEff,U) Gauss linear limited 0.5; laplacian(DkEff,k) Gauss linear limited 0.5; laplacian(DepsilonEff,epsilon) Gauss linear limited 0.5; laplacian(alphaEff,e) Gauss linear limited 0.5; From 5305c8f4aed5c2012088ccbdbd97c045f6f93d8f Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 25 Nov 2014 13:34:24 +0000 Subject: [PATCH 04/45] Removed end Info --- .../utilities/postProcessing/velocityField/Mach/Mach.C | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/applications/utilities/postProcessing/velocityField/Mach/Mach.C b/applications/utilities/postProcessing/velocityField/Mach/Mach.C index a99140493..edef1eb36 100644 --- a/applications/utilities/postProcessing/velocityField/Mach/Mach.C +++ b/applications/utilities/postProcessing/velocityField/Mach/Mach.C @@ -36,6 +36,7 @@ Description #include "basicPsiThermo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh) { bool writeResults = !args.optionFound("noWrite"); @@ -134,8 +135,6 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh) { Info<< " Missing U or T" << endl; } - - Info<< "\nEnd\n" << endl; } From 953b9dea771648b638c45ee5771130e050281a29 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 27 Nov 2014 19:09:11 +0000 Subject: [PATCH 05/45] Updated for dynamic meshes; improved numerics --- .../compressible/sonicDyMFoam/Make/files | 3 + .../compressible/sonicDyMFoam/Make/options | 21 ++ .../solvers/compressible/sonicDyMFoam/UEqn.H | 24 +++ .../sonicDyMFoam/compressibleContinuityErrs.H | 20 ++ .../compressible/sonicDyMFoam/createFields.H | 64 ++++++ .../solvers/compressible/sonicDyMFoam/eEqn.H | 44 ++++ .../compressible/sonicDyMFoam/limitU.H | 17 ++ .../compressible/sonicDyMFoam/readControls.H | 14 ++ .../sonicDyMFoam/readFieldBounds.H | 20 ++ .../sonicDyMFoam/readTransportProperties.H | 18 ++ .../compressible/sonicDyMFoam/sonicDyMFoam.C | 194 ++++++++++++++++++ 11 files changed, 439 insertions(+) create mode 100644 applications/solvers/compressible/sonicDyMFoam/Make/files create mode 100644 applications/solvers/compressible/sonicDyMFoam/Make/options create mode 100644 applications/solvers/compressible/sonicDyMFoam/UEqn.H create mode 100644 applications/solvers/compressible/sonicDyMFoam/compressibleContinuityErrs.H create mode 100644 applications/solvers/compressible/sonicDyMFoam/createFields.H create mode 100644 applications/solvers/compressible/sonicDyMFoam/eEqn.H create mode 100644 applications/solvers/compressible/sonicDyMFoam/limitU.H create mode 100644 applications/solvers/compressible/sonicDyMFoam/readControls.H create mode 100644 applications/solvers/compressible/sonicDyMFoam/readFieldBounds.H create mode 100644 applications/solvers/compressible/sonicDyMFoam/readTransportProperties.H create mode 100644 applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C diff --git a/applications/solvers/compressible/sonicDyMFoam/Make/files b/applications/solvers/compressible/sonicDyMFoam/Make/files new file mode 100644 index 000000000..08fa87a2c --- /dev/null +++ b/applications/solvers/compressible/sonicDyMFoam/Make/files @@ -0,0 +1,3 @@ +sonicDyMFoam.C + +EXE = $(FOAM_APPBIN)/sonicDyMFoam diff --git a/applications/solvers/compressible/sonicDyMFoam/Make/options b/applications/solvers/compressible/sonicDyMFoam/Make/options new file mode 100644 index 000000000..5223e8e5b --- /dev/null +++ b/applications/solvers/compressible/sonicDyMFoam/Make/options @@ -0,0 +1,21 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ + -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -lbasicThermophysicalModels \ + -lspecie \ + -lcompressibleTurbulenceModel \ + -lcompressibleRASModels \ + -lcompressibleLESModels \ + -lfiniteVolume \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh \ + -ldynamicMesh \ + -lmeshTools \ + -llduSolvers diff --git a/applications/solvers/compressible/sonicDyMFoam/UEqn.H b/applications/solvers/compressible/sonicDyMFoam/UEqn.H new file mode 100644 index 000000000..377aa0a13 --- /dev/null +++ b/applications/solvers/compressible/sonicDyMFoam/UEqn.H @@ -0,0 +1,24 @@ + fvVectorMatrix UEqn + ( + fvm::ddt(rho, U) + + fvm::div(phi, U) + + turbulence->divDevRhoReff(U) + ); + + if (oCorr == nOuterCorr - 1) + { + if (mesh.solutionDict().relax("UFinal")) + { + UEqn.relax(mesh.solutionDict().relaxationFactor("UFinal")); + } + else + { + UEqn.relax(1); + } + } + else + { + UEqn.relax(); + } + + solve(UEqn == -fvc::grad(p)); diff --git a/applications/solvers/compressible/sonicDyMFoam/compressibleContinuityErrs.H b/applications/solvers/compressible/sonicDyMFoam/compressibleContinuityErrs.H new file mode 100644 index 000000000..eec848532 --- /dev/null +++ b/applications/solvers/compressible/sonicDyMFoam/compressibleContinuityErrs.H @@ -0,0 +1,20 @@ +{ +# include "rhoEqn.H" +} +{ + scalar sumLocalContErr = + sum + ( + mag(rho.internalField() - (psi*p)().internalField()) + )/sum(rho.internalField()); + + scalar globalContErr = + sum(rho.internalField() - (psi*p)().internalField()) + /sum(rho.internalField()); + + cumulativeContErr += globalContErr; + + Info<< "time step continuity errors : sum local = " << sumLocalContErr + << ", global = " << globalContErr + << ", cumulative = " << cumulativeContErr << endl; +} diff --git a/applications/solvers/compressible/sonicDyMFoam/createFields.H b/applications/solvers/compressible/sonicDyMFoam/createFields.H new file mode 100644 index 000000000..fac70b85e --- /dev/null +++ b/applications/solvers/compressible/sonicDyMFoam/createFields.H @@ -0,0 +1,64 @@ + Info<< "Reading thermophysical properties\n" << endl; + + autoPtr pThermo + ( + basicPsiThermo::New(mesh) + ); + basicPsiThermo& thermo = pThermo(); + + volScalarField& p = thermo.p(); + volScalarField& e = thermo.e(); + const volScalarField& psi = thermo.psi(); + + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh + ), + thermo.rho() + ); + + Info<< "Reading field U\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + +# include "compressibleCreatePhi.H" + + // Store energy source term for under-relaxation + volScalarField pDivU + ( + IOobject + ( + "pDivU", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + p*fvc::div(phi/fvc::interpolate(rho)) + ); + + Info<< "Creating turbulence model\n" << endl; + autoPtr turbulence + ( + compressible::turbulenceModel::New + ( + rho, + U, + phi, + thermo + ) + ); diff --git a/applications/solvers/compressible/sonicDyMFoam/eEqn.H b/applications/solvers/compressible/sonicDyMFoam/eEqn.H new file mode 100644 index 000000000..a1d0940f8 --- /dev/null +++ b/applications/solvers/compressible/sonicDyMFoam/eEqn.H @@ -0,0 +1,44 @@ +{ + fvScalarMatrix eEqn + ( + fvm::ddt(rho, e) + + fvm::div(phi, e) + - fvm::laplacian(turbulence->alphaEff(), e) + == +// - p*fvc::div(phi/fvc::interpolate(rho) + fvc::meshPhi(rho, U)) + - fvm::SuSp(pDivU/e, e) + // viscous heating? + ); + + if (oCorr == nOuterCorr - 1) + { + if (mesh.solutionDict().relax("eFinal")) + { + eEqn.relax(mesh.solutionDict().relaxationFactor("eFinal")); + } + else + { + eEqn.relax(1); + } + } + else + { + eEqn.relax(); + } + + eEqn.solve(); + + // Bound the energy using TMin and TMax + { + dimensionedScalar Tstd("Tstd", dimTemperature, specie::Tstd); + + volScalarField Cv = thermo.Cv(); + volScalarField R = thermo.Cp() - Cv; + + e = Foam::min(e, TMax*Cv + R*Tstd); + e = Foam::max(e, TMin*Cv + R*Tstd); + e.correctBoundaryConditions(); + } + + thermo.correct(); +} diff --git a/applications/solvers/compressible/sonicDyMFoam/limitU.H b/applications/solvers/compressible/sonicDyMFoam/limitU.H new file mode 100644 index 000000000..d7eaeb3d9 --- /dev/null +++ b/applications/solvers/compressible/sonicDyMFoam/limitU.H @@ -0,0 +1,17 @@ +{ + // Bound the velocity + volScalarField magU = mag(U); + + if (max(magU) > UMax) + { + 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/compressible/sonicDyMFoam/readControls.H b/applications/solvers/compressible/sonicDyMFoam/readControls.H new file mode 100644 index 000000000..3bd20c5c5 --- /dev/null +++ b/applications/solvers/compressible/sonicDyMFoam/readControls.H @@ -0,0 +1,14 @@ +# include "readTimeControls.H" +# include "readPIMPLEControls.H" + + bool correctPhi = false; + if (pimple.found("correctPhi")) + { + correctPhi = Switch(pimple.lookup("correctPhi")); + } + + bool checkMeshCourantNo = false; + if (pimple.found("checkMeshCourantNo")) + { + checkMeshCourantNo = Switch(pimple.lookup("checkMeshCourantNo")); + } diff --git a/applications/solvers/compressible/sonicDyMFoam/readFieldBounds.H b/applications/solvers/compressible/sonicDyMFoam/readFieldBounds.H new file mode 100644 index 000000000..9628254e3 --- /dev/null +++ b/applications/solvers/compressible/sonicDyMFoam/readFieldBounds.H @@ -0,0 +1,20 @@ + // Read field bounds + dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds"); + + // Pressure bounds + dimensionedScalar pMin("pMin", dimPressure, 0); + dimensionedScalar pMax("pMax", dimPressure, GREAT); + + fieldBounds.lookup("p") >> pMin.value() >> pMax.value(); + + // Temperature bounds + dimensionedScalar TMin("TMin", dimTemperature, 0); + dimensionedScalar TMax("TMax", dimTemperature, GREAT); + + fieldBounds.lookup("T") >> TMin.value() >> TMax.value(); + + // Velocity bound + dimensionedScalar UMax("UMax", dimVelocity, GREAT); + + fieldBounds.lookup(U.name()) >> UMax.value(); + dimensionedScalar smallU("smallU", dimVelocity, 1e-10); diff --git a/applications/solvers/compressible/sonicDyMFoam/readTransportProperties.H b/applications/solvers/compressible/sonicDyMFoam/readTransportProperties.H new file mode 100644 index 000000000..1502e2033 --- /dev/null +++ b/applications/solvers/compressible/sonicDyMFoam/readTransportProperties.H @@ -0,0 +1,18 @@ + Info<< "Reading transportProperties\n" << endl; + + IOdictionary transportProperties + ( + IOobject + ( + "transportProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + dimensionedScalar mu + ( + transportProperties.lookup("mu") + ); diff --git a/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C b/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C new file mode 100644 index 000000000..f14dd3367 --- /dev/null +++ b/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C @@ -0,0 +1,194 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Application + sonicDyMFoam + +Description + Transient solver for trans-sonic/supersonic, laminar flow of a + compressible gas with support for mesh motion and topological changes + + Updated from sonicFoamAutoMotion by Hrvoje Jasak + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "dynamicFvMesh.H" +#include "specie.H" +#include "basicPsiThermo.H" +#include "turbulenceModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ +# include "setRootCase.H" + +# include "createTime.H" +# include "createDynamicFvMesh.H" +# include "createFields.H" +# include "initContinuityErrs.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { +# include "readControls.H" +# include "readFieldBounds.H" +# include "compressibleCourantNo.H" +# include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + bool meshChanged = mesh.update(); + +# include "volContinuity.H" + + if (checkMeshCourantNo) + { +# include "meshCourantNo.H" + } + + // Mesh motion update +// if (correctPhi && meshChanged) +// { +// # include "correctPhi.H" +// } + + if (meshChanged) + { +# include "CourantNo.H" + } + + // --- PIMPLE loop + label oCorr = 0; + do + { + // Under-relax pDivU term + pDivU.storePrevIter(); + pDivU = p*fvc::div(phi/fvc::interpolate(rho)); + pDivU.relax(); + +# include "rhoEqn.H" +# include "eEqn.H" +# include "UEqn.H" + + // --- PISO loop + for (int corr = 0; corr < nCorr; corr++) + { + U = UEqn.H()/UEqn.A(); + +# include "limitU.H" + + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi)* + ( + (fvc::interpolate(U) & mesh.Sf()) + - fvc::meshPhi(rho, U) + ) + ); + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + { + // Store pressure for under-relaxation + p.storePrevIter(); + + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvm::div(phid, p) + - fvm::laplacian(rho/UEqn.A(), p) + ); + + if + ( +// oCorr == nOuterCorr - 1 + corr == nCorr - 1 + && nonOrth == nNonOrthCorr + ) + { + pEqn.solve + ( + mesh.solutionDict().solver(p.name() + "Final") + ); + } + else + { + pEqn.solve(mesh.solutionDict().solver(p.name())); + } + + if (nonOrth == nNonOrthCorr) + { + phi = pEqn.flux(); + } + + + // Bound the pressure + if (min(p) < pMin || max(p) > pMax) + { + p.max(pMin); + p.min(pMax); + p.correctBoundaryConditions(); + } + + // Relax the pressure + p.relax(); + } + +# include "compressibleContinuityErrs.H" + + U -= fvc::grad(p)/UEqn.A(); + U.correctBoundaryConditions(); + +# include "limitU.H" + } + + // Recalculate density + rho = thermo.rho(); + + turbulence->correct(); + } while (++oCorr < nOuterCorr); + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return(0); +} + + +// ************************************************************************* // From 03283d2edecc07e62d2de73e19fdce69740f31e3 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 5 Dec 2014 09:46:52 +0000 Subject: [PATCH 06/45] Formatting --- .../translatingWallVelocityFvPatchVectorField.C | 3 ++- .../translatingWallVelocityFvPatchVectorField.H | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.C index 81752df27..5d9c3b30b 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.C @@ -109,7 +109,8 @@ void translatingWallVelocityFvPatchVectorField::updateCoeffs() return; } - // Remove the component of U normal to the wall in case the wall is not flat + // Remove the component of U normal to the wall in case the wall + // is not flat vectorField n = patch().nf(); vectorField::operator=(U_ - n*(n & U_)); diff --git a/src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.H index 521933077..35d28d57c 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.H @@ -52,7 +52,7 @@ class translatingWallVelocityFvPatchVectorField { // Private data - //- Origin of the rotation + //- Translating velocity vector U_; From 135a28fa777e7f14deb9886ee302953d9e1cc5cc Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 8 Dec 2014 12:11:09 +0000 Subject: [PATCH 07/45] Formatting --- src/foam/primitives/ints/ulong/ulongIO.C | 1 + 1 file changed, 1 insertion(+) diff --git a/src/foam/primitives/ints/ulong/ulongIO.C b/src/foam/primitives/ints/ulong/ulongIO.C index f29f4593c..aba95599c 100644 --- a/src/foam/primitives/ints/ulong/ulongIO.C +++ b/src/foam/primitives/ints/ulong/ulongIO.C @@ -42,6 +42,7 @@ Foam::word Foam::name(const unsigned long val) return buf.str(); } + // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // Foam::Istream& Foam::operator>>(Istream& is, unsigned long& val) From abe380ca6a8832b3572ac0ab2c906dfd8644acd9 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 8 Dec 2014 12:11:39 +0000 Subject: [PATCH 08/45] LabelPairList --- .../polyMesh/globalMeshData/globalMeshData.H | 5 +- .../polyBoundaryMesh/polyBoundaryMesh.C | 6 +-- .../polyBoundaryMesh/polyBoundaryMesh.H | 16 +++--- src/foam/primitives/Lists/labelPairList.H | 51 +++++++++++++++++++ src/foam/primitives/Tuple2/Tuple2.H | 2 +- 5 files changed, 64 insertions(+), 16 deletions(-) create mode 100644 src/foam/primitives/Lists/labelPairList.H diff --git a/src/foam/meshes/polyMesh/globalMeshData/globalMeshData.H b/src/foam/meshes/polyMesh/globalMeshData/globalMeshData.H index e4881fa58..a2cd5464d 100644 --- a/src/foam/meshes/polyMesh/globalMeshData/globalMeshData.H +++ b/src/foam/meshes/polyMesh/globalMeshData/globalMeshData.H @@ -74,7 +74,7 @@ SourceFiles #include "polyMesh.H" #include "Switch.H" #include "processorTopology.H" -#include "labelPair.H" +#include "labelPairList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -119,9 +119,6 @@ class globalMeshData }; - typedef List labelPairList; - - // Private data //- Reference to mesh diff --git a/src/foam/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C b/src/foam/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C index 01edae0f4..b4a4e7365 100644 --- a/src/foam/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C +++ b/src/foam/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C @@ -168,7 +168,7 @@ void Foam::polyBoundaryMesh::calcGeometry() } -const Foam::List& +const Foam::labelPairListList& Foam::polyBoundaryMesh::neighbourEdges() const { if (Pstream::parRun()) @@ -180,8 +180,8 @@ Foam::polyBoundaryMesh::neighbourEdges() const if (!neighbourEdgesPtr_) { - neighbourEdgesPtr_ = new List(size()); - List& neighbourEdges = *neighbourEdgesPtr_; + neighbourEdgesPtr_ = new labelPairListList(size()); + labelPairListList& neighbourEdges = *neighbourEdgesPtr_; // Initialize. label nEdgePairs = 0; diff --git a/src/foam/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H b/src/foam/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H index e1a3cfc1e..0beb12d6c 100644 --- a/src/foam/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H +++ b/src/foam/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H @@ -37,14 +37,13 @@ SourceFiles #include "polyPatchList.H" #include "regIOobject.H" -#include "labelPair.H" +#include "labelPairList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { -typedef List labelPairList; class polyMesh; // Forward declaration of friend functions and operators @@ -69,7 +68,7 @@ class polyBoundaryMesh const polyMesh& mesh_; //- Edges of neighbouring patches - mutable List* neighbourEdgesPtr_; + mutable labelPairListList* neighbourEdgesPtr_; // Private Member Functions @@ -138,13 +137,14 @@ public: return mesh_; } - //- Per patch the edges on the neighbouring patch. Is for every external - // edge the neighbouring patch and neighbouring (external) patch edge - // label. Note that edge indices are offset by nInternalEdges to keep - // it as much as possible consistent with coupled patch addressing + //- Per patch the edges on the neighbouring patch + // For every external edge the neighbouring patch and + // neighbouring (external) patch edge label. Note that edge + // indices are offset by nInternalEdges to keep it as much as + // possible consistent with coupled patch addressing // (where coupling is by local patch face index). // Only valid for singly connected polyBoundaryMesh and not parallel - const List& neighbourEdges() const; + const labelPairListList& neighbourEdges() const; //- Return a list of patch names wordList names() const; diff --git a/src/foam/primitives/Lists/labelPairList.H b/src/foam/primitives/Lists/labelPairList.H new file mode 100644 index 000000000..8acd6a635 --- /dev/null +++ b/src/foam/primitives/Lists/labelPairList.H @@ -0,0 +1,51 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Typedef + Foam::labelPairList + +Description + List of labelPairs. + +\*---------------------------------------------------------------------------*/ + +#ifndef labelPairList_H +#define labelPairList_H + +#include "labelPair.H" +#include "List.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + typedef List labelPairList; + typedef List labelPairListList; + typedef List labelPairListListList; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/primitives/Tuple2/Tuple2.H b/src/foam/primitives/Tuple2/Tuple2.H index 88afa41b9..11e9755b9 100644 --- a/src/foam/primitives/Tuple2/Tuple2.H +++ b/src/foam/primitives/Tuple2/Tuple2.H @@ -25,7 +25,7 @@ Class Foam::Tuple2 Description - A 2-tuple. + A 2-tuple: storing two objects of different types. SeeAlso Foam::Pair for storing two objects of identical types. From 16d2a566f28517aaf952a3ab21265fe7d342c9b5 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 8 Dec 2014 12:22:58 +0000 Subject: [PATCH 09/45] Improved warning message --- .../polyPatches/constraint/ggi/ggiPolyPatch.C | 11 +++++++++++ .../polyPatches/constraint/ggi/ggiPolyPatch.H | 2 +- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C index 15a23a177..53f2a847e 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C @@ -962,6 +962,17 @@ void Foam::ggiPolyPatch::calcTransforms() const << abort(FatalError); } } + else + { + InfoIn("label ggiPolyPatch::shadowIndex() const") + << "ggi patch " << name() << " with shadow " + << shadowName() << " has " + << patchToPatch().uncoveredMasterFaces().size() + << " uncovered master faces and " + << patchToPatch().uncoveredSlaveFaces().size() + << " uncovered slave faces. Bridging is switched on. " + << abort(FatalError); + } } } diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H index 6a57f3f04..313971bae 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H @@ -151,7 +151,7 @@ class ggiPolyPatch // Memory management //- Clear geometry - void clearGeom(); + void clearGeom(); //- Clear out void clearOut(); From 5dcdf87b7921824f50fcbf407f0476c401986845 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 8 Dec 2014 12:23:33 +0000 Subject: [PATCH 10/45] LabelPairList --- src/foam/meshes/polyMesh/globalMeshData/globalMeshData.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/foam/meshes/polyMesh/globalMeshData/globalMeshData.H b/src/foam/meshes/polyMesh/globalMeshData/globalMeshData.H index a2cd5464d..33692aed7 100644 --- a/src/foam/meshes/polyMesh/globalMeshData/globalMeshData.H +++ b/src/foam/meshes/polyMesh/globalMeshData/globalMeshData.H @@ -88,7 +88,7 @@ Ostream& operator<<(Ostream&, const globalMeshData&); /*---------------------------------------------------------------------------*\ - Class globalMeshData Declaration + Class globalMeshData Declaration \*---------------------------------------------------------------------------*/ class globalMeshData From 00e8dce8dab56d004d17815f347bc90b34f3a00d Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 9 Dec 2014 17:57:24 +0000 Subject: [PATCH 11/45] Formatting --- .../slidingInterface/slidingInterfaceProjectPoints.C | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/slidingInterface/slidingInterfaceProjectPoints.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/slidingInterface/slidingInterfaceProjectPoints.C index a0c455ab1..4af61af0e 100644 --- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/slidingInterface/slidingInterfaceProjectPoints.C +++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/slidingInterface/slidingInterfaceProjectPoints.C @@ -297,7 +297,10 @@ bool Foam::slidingInterface::projectPoints() const } // Projected slave points are stored for mesh motion correction - if (projectedSlavePointsPtr_) delete projectedSlavePointsPtr_; + if (projectedSlavePointsPtr_) + { + delete projectedSlavePointsPtr_; + } projectedSlavePointsPtr_ = new pointField(slavePointFaceHits.size(), vector::zero); From 6af575dbe1c62773f2beaa274af976557d8f0efe Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 16 Dec 2014 22:37:18 +0000 Subject: [PATCH 12/45] Update for topo changes support --- .../compressible/sonicDyMFoam/createFields.H | 1 + .../solvers/compressible/sonicDyMFoam/pEqn.H | 69 +++++++++++++++ .../compressible/sonicDyMFoam/sonicDyMFoam.C | 83 +++---------------- 3 files changed, 83 insertions(+), 70 deletions(-) create mode 100644 applications/solvers/compressible/sonicDyMFoam/pEqn.H diff --git a/applications/solvers/compressible/sonicDyMFoam/createFields.H b/applications/solvers/compressible/sonicDyMFoam/createFields.H index fac70b85e..6d3406829 100644 --- a/applications/solvers/compressible/sonicDyMFoam/createFields.H +++ b/applications/solvers/compressible/sonicDyMFoam/createFields.H @@ -6,6 +6,7 @@ ); basicPsiThermo& thermo = pThermo(); + volScalarField& T = const_cast(thermo.T()); volScalarField& p = thermo.p(); volScalarField& e = thermo.e(); const volScalarField& psi = thermo.psi(); diff --git a/applications/solvers/compressible/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicDyMFoam/pEqn.H new file mode 100644 index 000000000..4be700755 --- /dev/null +++ b/applications/solvers/compressible/sonicDyMFoam/pEqn.H @@ -0,0 +1,69 @@ +{ + U = UEqn.H()/UEqn.A(); + +# include "limitU.H" + + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi)* + ( + (fvc::interpolate(U) & mesh.Sf()) + - fvc::meshPhi(rho, U) + ) + ); + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + { + // Store pressure for under-relaxation + p.storePrevIter(); + + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvm::div(phid, p) + - fvm::laplacian(rho/UEqn.A(), p) + ); + + if + ( +// oCorr == nOuterCorr - 1 + corr == nCorr - 1 + && nonOrth == nNonOrthCorr + ) + { + pEqn.solve + ( + mesh.solutionDict().solver(p.name() + "Final") + ); + } + else + { + pEqn.solve(mesh.solutionDict().solver(p.name())); + } + + if (nonOrth == nNonOrthCorr) + { + phi = pEqn.flux(); + } + + + // Bound the pressure + if (min(p) < pMin || max(p) > pMax) + { + p.max(pMin); + p.min(pMax); + p.correctBoundaryConditions(); + } + + // Relax the pressure + p.relax(); + } + +# include "compressibleContinuityErrs.H" + + U -= fvc::grad(p)/UEqn.A(); + U.correctBoundaryConditions(); + +# include "limitU.H" +} diff --git a/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C b/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C index f14dd3367..377f63da7 100644 --- a/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C +++ b/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C @@ -77,10 +77,19 @@ int main(int argc, char *argv[]) } // Mesh motion update -// if (correctPhi && meshChanged) -// { + if (meshChanged) + { + T.correctBoundaryConditions(); + p.correctBoundaryConditions(); + e.correctBoundaryConditions(); + thermo.correct(); + rho = thermo.rho(); + + if (correctPhi) + { // # include "correctPhi.H" -// } + } + } if (meshChanged) { @@ -103,73 +112,7 @@ int main(int argc, char *argv[]) // --- PISO loop for (int corr = 0; corr < nCorr; corr++) { - U = UEqn.H()/UEqn.A(); - -# include "limitU.H" - - surfaceScalarField phid - ( - "phid", - fvc::interpolate(psi)* - ( - (fvc::interpolate(U) & mesh.Sf()) - - fvc::meshPhi(rho, U) - ) - ); - - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) - { - // Store pressure for under-relaxation - p.storePrevIter(); - - fvScalarMatrix pEqn - ( - fvm::ddt(psi, p) - + fvm::div(phid, p) - - fvm::laplacian(rho/UEqn.A(), p) - ); - - if - ( -// oCorr == nOuterCorr - 1 - corr == nCorr - 1 - && nonOrth == nNonOrthCorr - ) - { - pEqn.solve - ( - mesh.solutionDict().solver(p.name() + "Final") - ); - } - else - { - pEqn.solve(mesh.solutionDict().solver(p.name())); - } - - if (nonOrth == nNonOrthCorr) - { - phi = pEqn.flux(); - } - - - // Bound the pressure - if (min(p) < pMin || max(p) > pMax) - { - p.max(pMin); - p.min(pMax); - p.correctBoundaryConditions(); - } - - // Relax the pressure - p.relax(); - } - -# include "compressibleContinuityErrs.H" - - U -= fvc::grad(p)/UEqn.A(); - U.correctBoundaryConditions(); - -# include "limitU.H" +# include "pEqn.H" } // Recalculate density From 8622f1fd19a8da595e1cf3b6d8084a30af3576c1 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 17 Dec 2014 09:47:48 +0000 Subject: [PATCH 13/45] Formatting and type --- ...rbulentMixingLengthDissipationRateInletFvPatchScalarField.C | 3 ++- ...rbulentMixingLengthDissipationRateInletFvPatchScalarField.H | 2 +- .../epsilonWallFunctionFvPatchScalarField.C | 2 +- .../omegaWallFunction/omegaWallFunctionFvPatchScalarField.C | 2 +- .../epsilonWallFunctionFvPatchScalarField.C | 2 +- .../omegaWallFunction/omegaWallFunctionFvPatchScalarField.C | 2 +- 6 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C index b996a598a..4dcb06522 100644 --- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C +++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C @@ -119,7 +119,8 @@ turbulentMixingLengthDissipationRateInletFvPatchScalarField // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void turbulentMixingLengthDissipationRateInletFvPatchScalarField::updateCoeffs() +void +turbulentMixingLengthDissipationRateInletFvPatchScalarField::updateCoeffs() { if (updated()) { diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.H index f50a8f291..b1c28259e 100644 --- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.H +++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.H @@ -32,7 +32,7 @@ Description \verbatim inlet { - type compressible::turbulentMixingLengthDissipationRateInlet; + type compressible::turbulentMixingLengthDissipationRateInlet; mixingLength 0.005; // 5 mm value uniform 200; // placeholder } diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C index f31c4b808..5ec811db5 100644 --- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C @@ -177,7 +177,7 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs() if (!db().foundObject(GName_)) { InfoIn("void epsilonWallFunctionFvPatchScalarField::updateCoeffs()") - << "Cannot access " << GName_ << " field. for patch " + << "Cannot access " << GName_ << " field for patch " << patch().name() << ". Evaluating as zeroGradient" << endl; diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C index 880901fce..4efbaa143 100644 --- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C @@ -182,7 +182,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs() if (!db().foundObject(GName_)) { InfoIn("void omegaWallFunctionFvPatchScalarField::updateCoeffs()") - << "Cannot access " << GName_ << " field. for patch " + << "Cannot access " << GName_ << " field for patch " << patch().name() << ". Evaluating as zeroGradient" << endl; diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C index 13571fbda..604310ebe 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C @@ -172,7 +172,7 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs() if (!db().foundObject(GName_)) { InfoIn("void epsilonWallFunctionFvPatchScalarField::updateCoeffs()") - << "Cannot access " << GName_ << " field. for patch " + << "Cannot access " << GName_ << " field for patch " << patch().name() << ". Evaluating as zeroGradient" << endl; diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C index 39d07d242..cf4bc6348 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C @@ -177,7 +177,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs() if (!db().foundObject(GName_)) { InfoIn("void omegaWallFunctionFvPatchScalarField::updateCoeffs()") - << "Cannot access " << GName_ << " field. for patch " + << "Cannot access " << GName_ << " field for patch " << patch().name() << ". Evaluating as zeroGradient" << endl; From d5b124d72f617975020c0296f55b77446bdecd8f Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 16 Jan 2015 11:08:02 +0000 Subject: [PATCH 14/45] Bug fix: out-of range access in implicit least squares Conflicts: src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C --- .../leastSquaresGrad/leastSquaresVectors.C | 135 +++++++----------- 1 file changed, 53 insertions(+), 82 deletions(-) diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C index ecbf38f43..8a083a01b 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C @@ -101,19 +101,15 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const // Set local references to mesh data const unallocLabelList& owner = mesh().owner(); const unallocLabelList& neighbour = mesh().neighbour(); - const volVectorField& C = mesh().C(); - const surfaceScalarField& w = mesh().weights(); -// const surfaceScalarField& magSf = mesh().magSf(); - // Set up temporary storage for the dd tensor (before inversion) symmTensorField dd(mesh().nCells(), symmTensor::zero); - forAll(owner, facei) + forAll(owner, faceI) { - label own = owner[facei]; - label nei = neighbour[facei]; + label own = owner[faceI]; + label nei = neighbour[faceI]; vector d = C[nei] - C[own]; symmTensor wdd = (1.0/magSqr(d))*sqr(d); @@ -122,113 +118,88 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const dd[nei] += wdd; } - - forAll(lsP.boundaryField(), patchi) + forAll(lsP.boundaryField(), patchI) { - const fvsPatchScalarField& pw = w.boundaryField()[patchi]; - // Note: least squares in 1.4.1 and other OpenCFD versions - // are wrong because of incorrect weighting. HJ, 23/Oct/2008 -// const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi]; - - const fvPatch& p = pw.patch(); + const fvPatch& p = mesh().boundary()[patchI]; const unallocLabelList& faceCells = p.patch().faceCells(); - // Build the d-vectors - - // Original version: closest distance to boundary - vectorField pd = - mesh().Sf().boundaryField()[patchi] - /( - mesh().magSf().boundaryField()[patchi] - *mesh().deltaCoeffs().boundaryField()[patchi] - ); - - if (!mesh().orthogonal()) - { - pd -= mesh().correctionVectors().boundaryField()[patchi] - /mesh().deltaCoeffs().boundaryField()[patchi]; - } - // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010 - // Experimental: review fixed gradient condition. HJ, 30/Sep/2010 -// vectorField pd = p.delta(); + const vectorField pd = p.delta(); - if (p.coupled()) + forAll(pd, patchFaceI) { - forAll(pd, patchFacei) - { - const vector& d = pd[patchFacei]; - - dd[faceCells[patchFacei]] += (1.0/magSqr(d))*sqr(d); - } - } - else - { - forAll(pd, patchFacei) - { - const vector& d = pd[patchFacei]; - - dd[faceCells[patchFacei]] += (1.0/magSqr(d))*sqr(d); - } + const vector& d = pd[patchFaceI]; + dd[faceCells[patchFaceI]] += (1.0/magSqr(d))*sqr(d); } } - - // Invert the dd tensor -// symmTensorField invDd = inv(dd); - // Fix: householder inverse. HJ, 3/Nov/2009 - symmTensorField invDd = hinv(dd); - + // For easy access of neighbour coupled patch field needed for + // lsN vectors on implicitly coupled boundaries. VV, 18/June/2014 + volSymmTensorField volInvDd + ( + IOobject + ( + "volInvDd", + mesh().pointsInstance(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh(), + dimensionedSymmTensor("zero", dimless, symmTensor::zero), + "zeroGradient" + ); + symmTensorField& invDd = volInvDd.internalField(); +// invDd = inv(dd); + invDd = hinv(dd); // Revisit all faces and calculate the lsP and lsN vectors - forAll(owner, facei) + forAll(owner, faceI) { - label own = owner[facei]; - label nei = neighbour[facei]; + label own = owner[faceI]; + label nei = neighbour[faceI]; vector d = C[nei] - C[own]; scalar magSfByMagSqrd = 1.0/magSqr(d); - lsP[facei] = magSfByMagSqrd*(invDd[own] & d); - lsN[facei] = -magSfByMagSqrd*(invDd[nei] & d); + lsP[faceI] = magSfByMagSqrd*(invDd[own] & d); + lsN[faceI] = -magSfByMagSqrd*(invDd[nei] & d); } - forAll(lsP.boundaryField(), patchi) + forAll(lsP.boundaryField(), patchI) { - fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchi]; - - const fvsPatchScalarField& pw = w.boundaryField()[patchi]; - // Note: least squares in 1.4.1 and other OpenCFD versions - // are wrong because of incorrect weighting. HJ, 23/Oct/2008 -// const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi]; - - const fvPatch& p = pw.patch(); + fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchI]; + fvsPatchVectorField& patchLsN = lsN.boundaryField()[patchI]; + const fvPatch& p = mesh().boundary()[patchI]; const unallocLabelList& faceCells = p.faceCells(); - // Build the d-vectors // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010 - vectorField pd = p.delta(); + const vectorField pd = p.delta(); if (p.coupled()) { - forAll(pd, patchFacei) - { - const vector& d = pd[patchFacei]; + const symmTensorField invDdNei = + volInvDd.boundaryField()[patchI].patchNeighbourField(); - patchLsP[patchFacei] = - (1.0/magSqr(d)) - *(invDd[faceCells[patchFacei]] & d); + forAll(pd, patchFaceI) + { + const vector& d = pd[patchFaceI]; + + patchLsP[patchFaceI] = (1.0/magSqr(d)) + *(invDd[faceCells[patchFaceI]] & d); + + patchLsN[patchFaceI] = - (1.0/magSqr(d)) + *(invDdNei[patchFaceI] & d); } } else { - forAll(pd, patchFacei) + forAll(pd, patchFaceI) { - const vector& d = pd[patchFacei]; + const vector& d = pd[patchFaceI]; - patchLsP[patchFacei] = - (1.0/magSqr(d)) - *(invDd[faceCells[patchFacei]] & d); + patchLsP[patchFaceI] = (1.0/magSqr(d)) + *(invDd[faceCells[patchFaceI]] & d); } } } From d8524310402c3c5368a595f9de2283158b506fc9 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 20 Jan 2015 10:58:44 +0000 Subject: [PATCH 15/45] Revert to Paraview-4.0.1 Conflicts: ThirdParty/AllMake.stage4 --- ThirdParty/AllMake.stage4 | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/ThirdParty/AllMake.stage4 b/ThirdParty/AllMake.stage4 index 6a3409bd7..2d0e77a37 100755 --- a/ThirdParty/AllMake.stage4 +++ b/ThirdParty/AllMake.stage4 @@ -56,10 +56,10 @@ echo Starting ThirdParty AllMake: Stage4 echo ======================================== echo -# qt-everywhere-opensource-src-4.8.5 +# qt-everywhere-opensource-src-4.8.6 if [ ! -z "$QT_THIRD_PARTY" ] then - ( rpm_make -p qt-everywhere-opensource-src-4.8.5 -s qt-everywhere-opensource-src-4.8.5.spec -u http://download.qt-project.org/official_releases/qt/4.8/4.8.5/qt-everywhere-opensource-src-4.8.5.tar.gz ) + ( rpm_make -p qt-everywhere-opensource-src-4.8.6 -s qt-everywhere-opensource-src-4.8.6.spec -u http://download.qt-project.org/official_releases/qt/4.8/4.8.6/qt-everywhere-opensource-src-4.8.6.tar.gz ) else echo "Using system installed QT" echo "" @@ -75,10 +75,14 @@ then then ( rpm_make -p ParaView-4.0.1 -s ParaView-4.0.1.spec -u http://downloads.sourceforge.net/project/openfoam-extend/foam-extend-3.1/ThirdParty/ParaView-v4.0.1-source.tgz \ -f --define='_qmakePath $QT_BIN_DIR/qmake' +# ( rpm_make -p ParaView-4.2.0 -s ParaView-4.2.0.spec -u http://www.paraview.org/files/v4.2/ParaView-v4.2.0-source.tar.gz \ +# -f --define='_qmakePath $QT_BIN_DIR/qmake' +# ( rpm_make -p ParaView-4.2.0 -s ParaView-4.2.0.spec -u http://downloads.sourceforge.net/project/openfoam-extend/foam-extend-3.1/ThirdParty/ParaView-v4.2.0-source.tgz \ +# -f --define='_qmakePath $QT_BIN_DIR/qmake' ) else echo "WARNING: " - echo "WARNING: Skipping the installation of ParaView-4.0.1." + echo "WARNING: Skipping the installation of ParaView-4.2.0." echo "WARNING: Please make sure the QT_BIN_DIR environment variable properly" echo "WARNING: initialized in the file prefs.sh or prefs.csh" echo "WARNING: The command \$QT_BIN_DIR/qmake needs to be valid" From 7ba950d12359823a90d8410b4bc8b477c1bd8982 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 21 Jan 2015 21:29:40 +0000 Subject: [PATCH 16/45] Formatting --- .../GGIInterpolationPolygonIntersection.C | 2 +- .../GGIInterpolationQuickRejectTests.C | 2 +- .../polyPatches/constraint/ggi/ggiPolyPatch.H | 5 ++- .../constraint/ggi/ggiPolyPatchTemplates.C | 33 +++++++++---------- 4 files changed, 20 insertions(+), 22 deletions(-) diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolationPolygonIntersection.C b/src/foam/interpolations/GGIInterpolation/GGIInterpolationPolygonIntersection.C index dff08b5d0..9ad0ab483 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolationPolygonIntersection.C +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolationPolygonIntersection.C @@ -192,7 +192,7 @@ GGIInterpolation::polygonIntersection ( "GGIInterpolation::" "polygonIntersection" - ) << "Intersection might be wrong wrong: clipping side " + ) << "Intersection might be wrong. Clipping side " << intersectionArea/clippingArea << " subject: " << intersectionArea/subjectArea << endl; } diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C b/src/foam/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C index b82fa0e99..82d50cc45 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C @@ -52,7 +52,7 @@ GGIInterpolation::faceBoundBoxExtendSpanFraction_ template const label - GGIInterpolation::octreeSearchMinNLevel_ +GGIInterpolation::octreeSearchMinNLevel_ ( debug::optimisationSwitch("GGIOctreeSearchMinNLevel", 3) ); diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H index 313971bae..0b3f73cb8 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H @@ -279,9 +279,8 @@ public: } - // Destructor - - virtual ~ggiPolyPatch(); + //- Destructor + virtual ~ggiPolyPatch(); // Member functions diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatchTemplates.C b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatchTemplates.C index 59e29599a..52b14ec62 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatchTemplates.C +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatchTemplates.C @@ -31,7 +31,6 @@ Author // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -// New. HJ, 12/Jun/2011 template Foam::tmp > Foam::ggiPolyPatch::fastExpand ( @@ -40,22 +39,6 @@ Foam::tmp > Foam::ggiPolyPatch::fastExpand { // Check and expand the field from patch size to zone size // with communication - - // Algorithm: - // 1) Master processor holds maps of all zone addressing (data provided) - // and all remote zone addressing (data required) - // 2) Each processor will send the locally active data to the master - // 3) Master assembles all the data - // 4) Master sends to all processors the data they need to receive - // - // Notes: - // A) If the size of zone addressing is zero, data is not sent - // B) Communicated data on each processor has the size of live faces - // C) Expanded data will be equal to actual data from other processors - // only for the faces marked in remote; for other faces, it will be - // equal to zero - // D) On processor zero, complete data is available - // HJ, 4/Jun/2011 if (ff.size() != size()) { FatalErrorIn @@ -81,6 +64,22 @@ Foam::tmp > Foam::ggiPolyPatch::fastExpand << abort(FatalError); } + // Algorithm: + // 1) Master processor holds maps of all zone addressing (data provided) + // and all remote zone addressing (data required) + // 2) Each processor will send the locally active data to the master + // 3) Master assembles all the data + // 4) Master sends to all processors the data they need to receive + // + // Notes: + // A) If the size of zone addressing is zero, data is not sent + // B) Communicated data on each processor has the size of live faces + // C) Expanded data will be equal to actual data from other processors + // only for the faces marked in remote; for other faces, it will be + // equal to zero + // D) On processor zero, complete data is available + // HJ, 4/Jun/2011 + // HR, 10/Jul/2013 // This function requires send-receive-addressing, but usage is not // symmetric across processors. Hence trigger re-calculate at this point From 2154e0dee1dd35ee278788bfcc9f8acff0214783 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 27 Jan 2015 12:12:36 +0000 Subject: [PATCH 17/45] Added mesh Co number check --- .../mesh/manipulation/moveDyMEngineMesh/moveDyMEngineMesh.C | 1 + .../mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C | 1 + 2 files changed, 2 insertions(+) diff --git a/applications/utilities/mesh/manipulation/moveDyMEngineMesh/moveDyMEngineMesh.C b/applications/utilities/mesh/manipulation/moveDyMEngineMesh/moveDyMEngineMesh.C index bfa198e7d..e0918b8a3 100644 --- a/applications/utilities/mesh/manipulation/moveDyMEngineMesh/moveDyMEngineMesh.C +++ b/applications/utilities/mesh/manipulation/moveDyMEngineMesh/moveDyMEngineMesh.C @@ -69,6 +69,7 @@ int main(int argc, char *argv[]) mesh.update(); # include "checkVolContinuity.H" +# include "meshCourantNo.H" if ( diff --git a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C index b5bd06d9c..842eec56f 100644 --- a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C +++ b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C @@ -62,6 +62,7 @@ int main(int argc, char *argv[]) mesh.update(); # include "checkVolContinuity.H" +# include "meshCourantNo.H" if (runTime.timeIndex() % checkFrequency == 0) { From 15a24604c7e118ffd7726d0cf5da9860ee381fd2 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 27 Jan 2015 12:12:53 +0000 Subject: [PATCH 18/45] Formatting Conflicts: applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C --- .../basic/potentialDyMFoam/potentialDyMFoam.C | 203 ++++++++++++++++++ .../compressible/sonicDyMFoam/sonicDyMFoam.C | 4 +- .../solvers/compressible/sonicFoam/UEqn.H | 4 +- .../sonicFoamAutoMotion/Make/files | 3 - .../sonicFoamAutoMotion/Make/options | 10 - .../compressibleContinuityErrs.H | 20 -- .../sonicFoamAutoMotion/createFields.H | 103 --------- .../readThermodynamicProperties.H | 23 -- .../readTransportProperties.H | 18 -- .../sonicFoamAutoMotion/sonicFoamAutoMotion.C | 138 ------------ .../steadyCompressibleFoam/readFieldBounds.H | 14 +- .../readFieldBounds.H | 14 +- .../readFieldBounds.H | 19 +- .../fluentMeshToFoam/fluentMeshToFoam.L | 2 +- .../manipulation/createPatch/createPatchDict | 2 +- 15 files changed, 233 insertions(+), 344 deletions(-) create mode 100644 applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C delete mode 100644 applications/solvers/compressible/sonicFoamAutoMotion/Make/files delete mode 100644 applications/solvers/compressible/sonicFoamAutoMotion/Make/options delete mode 100644 applications/solvers/compressible/sonicFoamAutoMotion/compressibleContinuityErrs.H delete mode 100644 applications/solvers/compressible/sonicFoamAutoMotion/createFields.H delete mode 100644 applications/solvers/compressible/sonicFoamAutoMotion/readThermodynamicProperties.H delete mode 100644 applications/solvers/compressible/sonicFoamAutoMotion/readTransportProperties.H delete mode 100644 applications/solvers/compressible/sonicFoamAutoMotion/sonicFoamAutoMotion.C diff --git a/applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C b/applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C new file mode 100644 index 000000000..c1707f5b5 --- /dev/null +++ b/applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C @@ -0,0 +1,203 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Application + potentialDyMFoam + +Description + Transient solver for potential flow with dynamic mesh. + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "dynamicFvMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + argList::validOptions.insert("resetU", ""); + argList::validOptions.insert("writep", ""); + +# include "setRootCase.H" +# include "createTime.H" +# include "createDynamicFvMesh.H" +# include "createFields.H" +# include "initTotalVolume.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.loop()) + { +# include "readPISOControls.H" +# include "checkTotalVolume.H" + + Info<< "Time = " << runTime.timeName() << nl << endl; + + bool meshChanged = mesh.update(); + reduce(meshChanged, orOp()); + + p.internalField() = 0; + + if (args.optionFound("resetU")) + { + U.internalField() = vector::zero; + } + +# include "volContinuity.H" +# include "meshCourantNo.H" + + // Solve potential flow equations + adjustPhi(phi, U, p); + + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) + { + p.storePrevIter(); + + fvScalarMatrix pEqn + ( + fvm::laplacian + ( + dimensionedScalar + ( + "1", + dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0), + 1 + ), + p + ) + == + fvc::div(phi) + ); + + pEqn.setReference(pRefCell, pRefValue); + pEqn.solve(); + + if (nonOrth == nNonOrthCorr) + { + phi -= pEqn.flux(); + } + else + { + p.relax(); + } + } + + Info<< "continuity error = " + << mag(fvc::div(phi))().weightedAverage(mesh.V()).value() + << endl; + + U = fvc::reconstruct(phi); + U.correctBoundaryConditions(); + + Info<< "Interpolated U error = " + << (sqrt(sum(sqr((fvc::interpolate(U) & mesh.Sf()) - phi))) + /sum(mesh.magSf())).value() + << endl; + + // Calculate velocity magnitude + { + volScalarField magU = mag(U); + + Info<< "mag(U): max: " << gMax(magU.internalField()) + << " min: " << gMin(magU.internalField()) << endl; + } + + if (args.optionFound("writep")) + { + // Find reference patch + label refPatch = -1; + scalar maxMagU = 0; + + // Go through all velocity patches and find the one that fixes + // velocity to the largest value + + forAll (U.boundaryField(), patchI) + { + const fvPatchVectorField& Upatch = U.boundaryField()[patchI]; + + if (Upatch.fixesValue()) + { + // Calculate mean velocity + scalar u = sum(mag(Upatch)); + label patchSize = Upatch.size(); + + reduce(u, sumOp()); + reduce(patchSize, sumOp