From 04033b6968c5931107b131de29e9085c92d3937b Mon Sep 17 00:00:00 2001 From: andy Date: Fri, 2 Mar 2012 17:13:46 +0000 Subject: [PATCH 1/4] ENH: Re-inistated weighted-average fieldAverage usage --- .../field/fieldValues/cellSource/cellSource.C | 5 +++-- .../field/fieldValues/cellSource/cellSource.H | 9 +++++++-- .../fieldValues/cellSource/cellSourceTemplates.C | 10 ++++++++-- .../field/fieldValues/faceSource/faceSource.C | 5 +++-- .../field/fieldValues/faceSource/faceSource.H | 7 +++++-- .../fieldValues/faceSource/faceSourceTemplates.C | 11 +++++++++-- 6 files changed, 35 insertions(+), 12 deletions(-) diff --git a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C index b42f15e6b3..9caa7a85b9 100644 --- a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C +++ b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C @@ -50,12 +50,13 @@ namespace Foam const char* Foam::NamedEnum < Foam::fieldValues::cellSource::operationType, - 8 + 9 >::names[] = { "none", "sum", "average", + "weightedAverage", "volAverage", "volIntegrate", "min", @@ -68,7 +69,7 @@ namespace Foam const Foam::NamedEnum Foam::fieldValues::cellSource::sourceTypeNames_; -const Foam::NamedEnum +const Foam::NamedEnum Foam::fieldValues::cellSource::operationTypeNames_; diff --git a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.H b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.H index 573b53f522..5e857e02f2 100644 --- a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.H +++ b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.H @@ -51,9 +51,12 @@ Description - none - sum - average + - weightedAverage - volAverage - volIntegrate - CoV (Coefficient of variation: standard deviation/mean) + - min + - max SourceFiles cellSource.C @@ -105,6 +108,7 @@ public: opNone, opSum, opAverage, + opWeightedAverage, opVolAverage, opVolIntegrate, opMin, @@ -113,7 +117,7 @@ public: }; //- Operation type names - static const NamedEnum operationTypeNames_; + static const NamedEnum operationTypeNames_; private: @@ -169,7 +173,8 @@ protected: Type processValues ( const Field& values, - const scalarField& V + const scalarField& V, + const scalarField& weightField ) const; //- Output file header information diff --git a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceTemplates.C b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceTemplates.C index 61d39a43c9..942ae6a0a9 100644 --- a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceTemplates.C +++ b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceTemplates.C @@ -78,7 +78,8 @@ template Type Foam::fieldValues::cellSource::processValues ( const Field& values, - const scalarField& V + const scalarField& V, + const scalarField& weightField ) const { Type result = pTraits::zero; @@ -94,6 +95,11 @@ Type Foam::fieldValues::cellSource::processValues result = sum(values)/values.size(); break; } + case opWeightedAverage: + { + result = sum(values)/sum(weightField); + break; + } case opVolAverage: { result = sum(values*V)/sum(V); @@ -169,7 +175,7 @@ bool Foam::fieldValues::cellSource::writeValues(const word& fieldName) if (Pstream::master()) { - Type result = processValues(values, V); + Type result = processValues(values, V, weightField); if (valueOutput_) { diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C index 68befe8fd9..97f244c60b 100644 --- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C +++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C @@ -53,12 +53,13 @@ namespace Foam const char* Foam::NamedEnum < Foam::fieldValues::faceSource::operationType, - 8 + 9 >::names[] = { "none", "sum", "average", + "weightedAverage", "areaAverage", "areaIntegrate", "min", @@ -72,7 +73,7 @@ namespace Foam const Foam::NamedEnum Foam::fieldValues::faceSource::sourceTypeNames_; -const Foam::NamedEnum +const Foam::NamedEnum Foam::fieldValues::faceSource::operationTypeNames_; diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H index 3c5ab78624..b7a2924035 100644 --- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H +++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H @@ -59,6 +59,7 @@ Description - none - sum - average + - weightedAverage - areaAverage - areaIntegrate - min @@ -132,6 +133,7 @@ public: opNone, opSum, opAverage, + opWeightedAverage, opAreaAverage, opAreaIntegrate, opMin, @@ -140,7 +142,7 @@ public: }; //- Operation type names - static const NamedEnum operationTypeNames_; + static const NamedEnum operationTypeNames_; private: @@ -215,7 +217,8 @@ protected: Type processValues ( const Field& values, - const scalarField& magSf + const scalarField& magSf, + const scalarField& weightField ) const; //- Output file header information diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C index bf83e811bc..79c031818f 100644 --- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C +++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C @@ -97,7 +97,9 @@ template Type Foam::fieldValues::faceSource::processValues ( const Field& values, - const scalarField& magSf + const scalarField& magSf, + const scalarField& weightField + ) const { Type result = pTraits::zero; @@ -113,6 +115,11 @@ Type Foam::fieldValues::faceSource::processValues result = sum(values)/values.size(); break; } + case opWeightedAverage: + { + result = sum(values)/sum(weightField); + break; + } case opAreaAverage: { result = sum(values*magSf)/sum(magSf); @@ -203,7 +210,7 @@ bool Foam::fieldValues::faceSource::writeValues(const word& fieldName) if (Pstream::master()) { - Type result = processValues(values, magSf); + Type result = processValues(values, magSf, weightField); if (valueOutput_) { From 658f0a0680219999007a772fc5ab067f16ce7478 Mon Sep 17 00:00:00 2001 From: andy Date: Fri, 2 Mar 2012 17:14:30 +0000 Subject: [PATCH 2/4] ENH: Corrected/cleaned smallPoolFire2D tutorial input files --- tutorials/combustion/fireFoam/les/smallPoolFire2D/0/CH4 | 6 +++++- tutorials/combustion/fireFoam/les/smallPoolFire2D/0/G | 5 +++++ .../combustion/fireFoam/les/smallPoolFire2D/0/IDefault | 5 +++++ tutorials/combustion/fireFoam/les/smallPoolFire2D/0/N2 | 6 +++++- tutorials/combustion/fireFoam/les/smallPoolFire2D/0/O2 | 6 +++++- tutorials/combustion/fireFoam/les/smallPoolFire2D/0/T | 6 +++++- tutorials/combustion/fireFoam/les/smallPoolFire2D/0/U | 6 +++++- .../combustion/fireFoam/les/smallPoolFire2D/0/Ydefault | 3 ++- .../combustion/fireFoam/les/smallPoolFire2D/0/alphaSgs | 6 +++++- tutorials/combustion/fireFoam/les/smallPoolFire2D/0/k | 6 +++++- tutorials/combustion/fireFoam/les/smallPoolFire2D/0/muSgs | 6 +++++- tutorials/combustion/fireFoam/les/smallPoolFire2D/0/omega | 7 +++++-- tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p | 5 ++++- tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p_rgh | 6 +++++- 14 files changed, 66 insertions(+), 13 deletions(-) diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/CH4 b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/CH4 index dd13fa9a57..0c17615d4c 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/CH4 +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/CH4 @@ -27,22 +27,26 @@ boundaryField inletValue uniform 0; value uniform 0; } + sides { type inletOutlet; inletValue uniform 0; value uniform 0; } + base { type zeroGradient; } + inlet { type fixedValue; value uniform 1.0; } - frontBack + + frontAndBack { type empty; } diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/G b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/G index fe75242d76..5dfaf2e402 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/G +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/G @@ -28,6 +28,11 @@ boundaryField emissivity uniform 1.0; value uniform 0; } + + frontAndBack + { + type empty; + } } // ************************************************************************* // diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/IDefault b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/IDefault index e70bf3dec0..75186697cd 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/IDefault +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/IDefault @@ -28,6 +28,11 @@ boundaryField emissivity uniform 1.0; value uniform 0; } + + frontAndBack + { + type empty; + } } // ************************************************************************* // diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/N2 b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/N2 index ff1a3b1c42..248e3fa02d 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/N2 +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/N2 @@ -25,19 +25,23 @@ boundaryField { type calculated; } + sides { type calculated; } + base { type calculated; } + inlet { type calculated; } - frontBack + + frontAndBack { type empty; } diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/O2 b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/O2 index a5dab92ad1..ef57804773 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/O2 +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/O2 @@ -27,22 +27,26 @@ boundaryField inletValue $internalField; value $internalField; } + sides { type inletOutlet; inletValue $internalField; value $internalField; } + base { type zeroGradient; } + inlet { type fixedValue; value uniform 0; } - frontBack + + frontAndBack { type empty; } diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/T b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/T index 1997334fc8..8f4990b6ce 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/T +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/T @@ -27,22 +27,26 @@ boundaryField inletValue uniform 300; value uniform 300; } + sides { type inletOutlet; inletValue uniform 300; value uniform 300; } + base { type zeroGradient; } + inlet { type fixedValue; value uniform 300; } - frontBack + + frontAndBack { type empty; } diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/U b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/U index c0dc7624ca..9a1ad04efe 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/U +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/U @@ -29,23 +29,27 @@ boundaryField value uniform (0 0 0); } + sides { type pressureInletOutletVelocity; outletValue uniform (0 0 0); value uniform (0 0 0); } + base { type fixedValue; value uniform (0 0 0); } + inlet { type fixedValue; value uniform (0 0.05 0); } - frontBack + + frontAndBack { type empty; } diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/Ydefault b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/Ydefault index 003100fc05..74003f28f3 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/Ydefault +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/Ydefault @@ -45,7 +45,8 @@ boundaryField inletValue $internalField; value $internalField; } - frontBack + + frontAndBack { type empty; } diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/alphaSgs b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/alphaSgs index 330bea1037..5495f422b4 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/alphaSgs +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/alphaSgs @@ -25,19 +25,23 @@ boundaryField { type zeroGradient; } + sides { type zeroGradient; } + base { type zeroGradient; } + inlet { type zeroGradient; } - frontBack + + frontAndBack { type empty; } diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/k b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/k index 51ff80a123..377a3d0090 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/k +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/k @@ -27,23 +27,27 @@ boundaryField inletValue uniform 1e-4; value uniform 1e-4; } + sides { type inletOutlet; inletValue uniform 1e-4; value uniform 1e-4; } + base { type fixedValue; value uniform 1e-4; } + inlet { type fixedValue; value uniform 1e-4; } - frontBack + + frontAndBack { type empty; } diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/muSgs b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/muSgs index 0c723a78c8..203285a6c2 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/muSgs +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/muSgs @@ -25,19 +25,23 @@ boundaryField { type zeroGradient; } + sides { type zeroGradient; } + base { type zeroGradient; } + inlet { type zeroGradient; } - frontBack + + frontAndBack { type empty; } diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/omega b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/omega index fbe63fbec6..faca2f4dd8 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/omega +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/omega @@ -28,26 +28,29 @@ boundaryField inletValue uniform 0; value uniform 0; } + sides { type inletOutlet; inletValue uniform 0; value uniform 0; } + base { type zeroGradient; } + inlet { type fixedValue; value uniform 0.0; } - frontBack + + frontAndBack { type empty; } - } // ************************************************************************* // diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p index 3ae8577524..85a1f09dc4 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p @@ -32,17 +32,20 @@ boundaryField type calculated; value $internalField; } + base { type calculated; value $internalField; } + inlet { type calculated; value $internalField; } - frontBack + + frontAndBack { type empty; } diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p_rgh b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p_rgh index cdeb5e62da..54e8cbeb29 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p_rgh +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p_rgh @@ -26,6 +26,7 @@ boundaryField type buoyantPressure; value $internalField; } + sides { type totalPressure; @@ -37,17 +38,20 @@ boundaryField gamma 1; value $internalField; } + base { type buoyantPressure; value $internalField; } + inlet { type buoyantPressure; value $internalField; } - frontBack + + frontAndBack { type empty; } From 912a20b7a3e0fa01564085c4dd512cb50734bf6a Mon Sep 17 00:00:00 2001 From: andy Date: Fri, 2 Mar 2012 18:15:54 +0000 Subject: [PATCH 3/4] ENH: Propagated caching of HbyA across solvers --- applications/solvers/DNS/dnsFoam/dnsFoam.C | 17 +++++---- applications/solvers/combustion/XiFoam/pEqn.H | 20 ++++++----- .../solvers/combustion/engineFoam/pEqn.H | 19 ++++++---- .../solvers/combustion/fireFoam/pEqn.H | 20 ++++++----- .../solvers/combustion/reactingFoam/pEqn.H | 20 ++++++----- .../solvers/combustion/rhoReactingFoam/pEqn.H | 35 ++++++++++--------- .../solvers/compressible/rhoPimpleFoam/pEqn.H | 24 +++++++------ .../rhoPorousMRFPimpleFoam/pEqn.H | 27 ++++++++------ .../solvers/compressible/rhoSimpleFoam/pEqn.H | 21 +++++++---- .../rhoPorousMRFSimpleFoam/pEqn.H | 29 +++++++++------ .../compressible/rhoSimplecFoam/pEqn.H | 22 ++++++++---- .../solvers/compressible/sonicFoam/pEqn.H | 7 ++-- .../sonicFoam/sonicDyMFoam/pEqn.H | 7 ++-- .../electromagnetics/mhdFoam/mhdFoam.C | 18 ++++++---- .../buoyantBoussinesqPimpleFoam/pEqn.H | 21 ++++++----- .../buoyantBoussinesqSimpleFoam/pEqn.H | 23 +++++++----- .../heatTransfer/buoyantPimpleFoam/pEqn.H | 25 +++++++------ .../heatTransfer/buoyantSimpleFoam/pEqn.H | 23 +++++++----- .../chtMultiRegionFoam/fluid/pEqn.H | 24 +++++++------ .../lagrangian/LTSReactingParcelFoam/pEqn.H | 28 +++++++-------- .../lagrangian/coalChemistryFoam/pEqn.H | 20 ++++++----- .../pEqn.H | 28 +++++++-------- .../lagrangian/reactingParcelFilmFoam/pEqn.H | 20 ++++++----- .../lagrangian/reactingParcelFoam/pEqn.H | 20 ++++++----- .../sprayFoam/sprayEngineFoam/pEqn.H | 20 ++++++----- .../solvers/multiphase/cavitatingFoam/pEqn.H | 8 ++--- .../cavitatingFoam/resetPhiPatches.H | 15 -------- .../cavitatingFoam/resetPhivPatches.H | 19 ---------- 28 files changed, 321 insertions(+), 259 deletions(-) delete mode 100644 applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H delete mode 100644 applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C index 8eabdce99f..0733802b08 100644 --- a/applications/solvers/DNS/dnsFoam/dnsFoam.C +++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C @@ -86,23 +86,28 @@ int main(int argc, char *argv[]) for (int corr=1; corr<=1; corr++) { volScalarField rAU(1.0/UEqn.A()); + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); - U = rAU*UEqn.H(); - phi = (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, U, phi); + surfaceScalarField phiHbyA + ( + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, U, phi) + ); fvScalarMatrix pEqn ( - fvm::laplacian(rAU, p) == fvc::div(phi) + fvm::laplacian(rAU, p) == fvc::div(phiHbyA) ); pEqn.solve(); - phi -= pEqn.flux(); + phi = phiHbyA - pEqn.flux(); #include "continuityErrs.H" - U -= rAU*fvc::grad(p); + U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); } diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H index cb25e7029f..723df7af98 100644 --- a/applications/solvers/combustion/XiFoam/pEqn.H +++ b/applications/solvers/combustion/XiFoam/pEqn.H @@ -1,7 +1,8 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); if (pimple.transonic()) { @@ -10,7 +11,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -34,19 +35,22 @@ if (pimple.transonic()) } else { - phi = + surfaceScalarField phiHbyA + ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) ); @@ -54,7 +58,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -62,7 +66,7 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H index d1b7135649..00603dc7d9 100644 --- a/applications/solvers/combustion/engineFoam/pEqn.H +++ b/applications/solvers/combustion/engineFoam/pEqn.H @@ -1,7 +1,8 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); if (pimple.transonic()) { @@ -9,7 +10,7 @@ if (pimple.transonic()) ( "phid", fvc::interpolate(psi) - *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)) + *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U)) ); while (pimple.correctNonOrthogonal()) @@ -31,15 +32,19 @@ if (pimple.transonic()) } else { - phi = fvc::interpolate(rho) - *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)); + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho) + *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U)) + ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) ); @@ -47,7 +52,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -55,7 +60,7 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H index 872a0513a3..f5364ae314 100644 --- a/applications/solvers/combustion/fireFoam/pEqn.H +++ b/applications/solvers/combustion/fireFoam/pEqn.H @@ -2,25 +2,29 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU)); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); -surfaceScalarField phiU +surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); + +surfaceScalarField phiHbyA ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) + - phig ); -phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf(); while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( fvc::ddt(psi, rho)*gh - + fvc::div(phi) + + fvc::div(phiHbyA) + fvm::ddt(psi, p_rgh) - fvm::laplacian(rhorAUf, p_rgh) == @@ -32,7 +36,9 @@ while (pimple.correctNonOrthogonal()) if (pimple.finalNonOrthogonalIter()) { - phi += p_rghEqn.flux(); + phi = phiHbyA + p_rghEqn.flux(); + U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() - phig)/rhorAUf); + U.correctBoundaryConditions(); } } @@ -41,8 +47,6 @@ p = p_rgh + rho*gh; #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf); -U.correctBoundaryConditions(); K = 0.5*magSqr(U); dpdt = fvc::ddt(p); diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H index cb25e7029f..723df7af98 100644 --- a/applications/solvers/combustion/reactingFoam/pEqn.H +++ b/applications/solvers/combustion/reactingFoam/pEqn.H @@ -1,7 +1,8 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); if (pimple.transonic()) { @@ -10,7 +11,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -34,19 +35,22 @@ if (pimple.transonic()) } else { - phi = + surfaceScalarField phiHbyA + ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) ); @@ -54,7 +58,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -62,7 +66,7 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/combustion/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/rhoReactingFoam/pEqn.H index 4e94de243c..ef11c677ad 100644 --- a/applications/solvers/combustion/rhoReactingFoam/pEqn.H +++ b/applications/solvers/combustion/rhoReactingFoam/pEqn.H @@ -6,27 +6,25 @@ thermo.rho() -= psi*p; volScalarField rAU(1.0/UEqn.A()); - U = rAU*UEqn.H(); + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); if (pimple.transonic()) { - surfaceScalarField phiv + surfaceScalarField phiHbyA ( - (fvc::interpolate(U) & mesh.Sf()) + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - phi = fvc::interpolate(rho)*phiv; + surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA); - surfaceScalarField phid - ( - "phid", - fvc::interpolate(thermo.psi())*phiv - ); + phiHbyA *= fvc::interpolate(rho); fvScalarMatrix pDDtEqn ( - fvc::ddt(rho) + fvc::div(phi) + fvc::ddt(rho) + fvc::div(phiHbyA) + correction(fvm::ddt(psi, p) + fvm::div(phid, p)) ); @@ -42,23 +40,26 @@ if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } else { - phi = + surfaceScalarField phiHbyA + ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); fvScalarMatrix pDDtEqn ( fvc::ddt(rho) + psi*correction(fvm::ddt(p)) - + fvc::div(phi) + + fvc::div(phiHbyA) ); while (pimple.correctNonOrthogonal()) @@ -73,7 +74,7 @@ if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -84,7 +85,7 @@ #include "rhoEqn.H" #include "compressibleContinuityErrs.H" - U -= rAU*fvc::grad(p); + U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H index f7bb932fdc..cfe34f2873 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H @@ -3,7 +3,8 @@ rho = max(rho, rhoMin); rho = min(rho, rhoMax); rho.relax(); -U = rAU*UEqn().H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn().H(); if (pimple.nCorrPISO() <= 1) { @@ -17,7 +18,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -41,12 +42,15 @@ if (pimple.transonic()) } else { - phi = - fvc::interpolate(rho)* - ( - (fvc::interpolate(U) & mesh.Sf()) + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho) + *( + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); while (pimple.correctNonOrthogonal()) { @@ -54,7 +58,7 @@ else fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) ); @@ -62,7 +66,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -81,7 +85,7 @@ rho.relax(); Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H index 831f139e81..c5c7602a43 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H @@ -4,7 +4,8 @@ rho = min(rho, rhoMax); rho.relax(); volScalarField rAU(1.0/UEqn().A()); -U = rAU*UEqn().H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn().H(); if (pimple.nCorrPISO() <= 1) { @@ -18,7 +19,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -43,13 +44,17 @@ if (pimple.transonic()) } else { - phi = - fvc::interpolate(rho)* - ( - (fvc::interpolate(U) & mesh.Sf()) + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho) + *( + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); - mrfZones.relativeFlux(fvc::interpolate(rho), phi); + ) + ); + + mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA); while (pimple.correctNonOrthogonal()) { @@ -57,7 +62,7 @@ else fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) ); @@ -65,7 +70,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -83,7 +88,7 @@ rho.relax(); Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H index 4d3f1ac76e..56c444cdab 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H @@ -4,7 +4,9 @@ rho = min(rho, rhoMax); rho.relax(); volScalarField rAU(1.0/UEqn().A()); -U = rAU*UEqn().H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn().H(); + UEqn.clear(); bool closedVolume = false; @@ -14,7 +16,7 @@ if (simple.transonic()) surfaceScalarField phid ( "phid", - fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf()) + fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf()) ); while (simple.correctNonOrthogonal()) @@ -40,14 +42,19 @@ if (simple.transonic()) } else { - phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); - closedVolume = adjustPhi(phi, U, p); + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf()) + ); + + closedVolume = adjustPhi(phiHbyA, U, p); while (simple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( - fvm::laplacian(rho*rAU, p) == fvc::div(phi) + fvm::laplacian(rho*rAU, p) == fvc::div(phiHbyA) ); pEqn.setReference(pRefCell, pRefValue); @@ -56,7 +63,7 @@ else if (simple.finalNonOrthogonalIter()) { - phi -= pEqn.flux(); + phi = phiHbyA - pEqn.flux(); } } } @@ -67,7 +74,7 @@ else // Explicitly relax pressure for momentum corrector p.relax(); -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); // For closed-volume cases adjust the pressure and density levels diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H index b1f83304cd..c3e0a43a15 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H @@ -1,10 +1,12 @@ +volVectorField HbyA("HbyA", U); + if (pressureImplicitPorosity) { - U = trTU()&UEqn().H(); + HbyA = trTU() & UEqn().H(); } else { - U = trAU()*UEqn().H(); + HbyA = trAU()*UEqn().H(); } UEqn.clear(); @@ -16,7 +18,7 @@ if (simple.transonic()) surfaceScalarField phid ( "phid", - fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf()) + fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf()) ); mrfZones.relativeFlux(fvc::interpolate(psi), phid); @@ -45,10 +47,15 @@ if (simple.transonic()) } else { - phi = fvc::interpolate(rho*U) & mesh.Sf(); - mrfZones.relativeFlux(fvc::interpolate(rho), phi); + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho*HbyA) & mesh.Sf() + ); - closedVolume = adjustPhi(phi, U, p); + mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA); + + closedVolume = adjustPhi(phiHbyA, U, p); while (simple.correctNonOrthogonal()) { @@ -56,11 +63,11 @@ else if (pressureImplicitPorosity) { - tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phi)); + tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phiHbyA)); } else { - tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phi)); + tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phiHbyA)); } tpEqn().setReference(pRefCell, pRefValue); @@ -69,7 +76,7 @@ else if (simple.finalNonOrthogonalIter()) { - phi -= tpEqn().flux(); + phi = phiHbyA - tpEqn().flux(); } } } @@ -81,11 +88,11 @@ p.relax(); if (pressureImplicitPorosity) { - U -= trTU()&fvc::grad(p); + U = HbyA - (trTU() & fvc::grad(p)); } else { - U -= trAU()*fvc::grad(p); + U = HbyA - trAU()*fvc::grad(p); } U.correctBoundaryConditions(); diff --git a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H index ccc7a1b21e..a4d9325522 100644 --- a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H @@ -7,7 +7,10 @@ volScalarField p0(p); volScalarField AU(UEqn().A()); volScalarField AtU(AU - UEqn().H1()); -U = UEqn().H()/AU; + +volVectorField HbyA("HbyA", U); +HbyA = UEqn().H()/AU; + UEqn.clear(); bool closedVolume = false; @@ -19,7 +22,7 @@ if (simple.transonic()) surfaceScalarField phid ( "phid", - fvc::interpolate(psi*U) & mesh.Sf() + fvc::interpolate(psi*HbyA) & mesh.Sf() ); surfaceScalarField phic @@ -56,13 +59,18 @@ else { while (simple.correctNonOrthogonal()) { - phi = fvc::interpolate(rho*U) & mesh.Sf(); + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho*HbyA) & mesh.Sf() + ); + closedVolume = adjustPhi(phi, U, p); phi += fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf(); fvScalarMatrix pEqn ( - fvc::div(phi) + fvc::div(phiHbyA) //- fvm::laplacian(rho/AU, p) - fvm::laplacian(rho/AtU, p) ); @@ -73,7 +81,7 @@ else if (simple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -85,8 +93,8 @@ else // Explicitly relax pressure for momentum corrector p.relax(); -U -= (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU); -//U -= fvc::grad(p)/AU; +U = HbyA - (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU); +//U = HbyA - fvc::grad(p)/AU; U.correctBoundaryConditions(); diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H index a4311dfd2e..46cc36216c 100644 --- a/applications/solvers/compressible/sonicFoam/pEqn.H +++ b/applications/solvers/compressible/sonicFoam/pEqn.H @@ -1,14 +1,15 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); surfaceScalarField phid ( "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -33,5 +34,5 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H index 98cc22d15a..cad5bb7217 100644 --- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H +++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H @@ -1,14 +1,15 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = UEqn.H()/UEqn.A(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); surfaceScalarField phid ( "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U) ) ); @@ -29,5 +30,5 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); diff --git a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C index a91ee8885e..51f04820eb 100644 --- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C +++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C @@ -93,17 +93,21 @@ int main(int argc, char *argv[]) for (int corr=0; corr 0) + surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf()); + if (pZones.size() == 0) { - // ddtPhiCorr not well defined for cases with porosity - phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); - } - else - { - phi = - fvc::interpolate(rho) - *( - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + // ddtPhiCorr only used without porosity + phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi); } + phiHbyA *= fvc::interpolate(rho); + + fvScalarMatrix pDDtEqn ( fvc::ddt(rho) + psi*correction(fvm::ddt(p)) - + fvc::div(phi) + + fvc::div(phiHbyA) == parcels.Srho() + sources(psi, p, rho.name()) @@ -46,7 +42,7 @@ if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } @@ -58,7 +54,7 @@ #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent #include "compressibleContinuityErrs.H" - U -= rAU*fvc::grad(p); + U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); sources.correct(U); diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H index cb8553bbfb..3a61b28e60 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H +++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H @@ -1,7 +1,8 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*(UEqn == sources(rho, U))().H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*(UEqn == sources(rho, U))().H(); if (pimple.transonic()) { @@ -10,7 +11,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -39,19 +40,22 @@ if (pimple.transonic()) } else { - phi = + surfaceScalarField phiHbyA + ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) == coalParcels.Srho() @@ -64,7 +68,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -72,7 +76,7 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); sources.correct(U); diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H index 6ed7affe77..4c9b41bba4 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H @@ -6,27 +6,23 @@ thermo.rho() -= psi*p; volScalarField rAU(1.0/UEqn.A()); - U = rAU*(UEqn == sources(rho, U))().H(); + volVectorField HbyA("HbyA", U); + HbyA = rAU*(UEqn == sources(rho, U))().H(); - if (pZones.size() > 0) + surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf()); + if (pZones.size() == 0) { - // ddtPhiCorr not well defined for cases with porosity - phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); - } - else - { - phi = - fvc::interpolate(rho) - *( - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + // ddtPhiCorr only used without porosity + phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi); } + phiHbyA *= fvc::interpolate(rho); + + fvScalarMatrix pDDtEqn ( fvc::ddt(rho) + psi*correction(fvm::ddt(p)) - + fvc::div(phi) + + fvc::div(phiHbyA) == parcels.Srho() + sources(psi, p, rho.name()) @@ -47,7 +43,7 @@ if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } @@ -57,7 +53,7 @@ #include "rhoEqn.H" #include "compressibleContinuityErrs.H" - U -= rAU*fvc::grad(p); + U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); sources.correct(U); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H index 872a0513a3..f5364ae314 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H @@ -2,25 +2,29 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU)); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); -surfaceScalarField phiU +surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); + +surfaceScalarField phiHbyA ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) + - phig ); -phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf(); while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( fvc::ddt(psi, rho)*gh - + fvc::div(phi) + + fvc::div(phiHbyA) + fvm::ddt(psi, p_rgh) - fvm::laplacian(rhorAUf, p_rgh) == @@ -32,7 +36,9 @@ while (pimple.correctNonOrthogonal()) if (pimple.finalNonOrthogonalIter()) { - phi += p_rghEqn.flux(); + phi = phiHbyA + p_rghEqn.flux(); + U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() - phig)/rhorAUf); + U.correctBoundaryConditions(); } } @@ -41,8 +47,6 @@ p = p_rgh + rho*gh; #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf); -U.correctBoundaryConditions(); K = 0.5*magSqr(U); dpdt = fvc::ddt(p); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H index 5ee551a0ea..c504fdf2ea 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H @@ -1,7 +1,8 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); if (pimple.transonic()) { @@ -10,7 +11,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -36,19 +37,22 @@ if (pimple.transonic()) } else { - phi = + surfaceScalarField phiHbyA + ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) == parcels.Srho() @@ -58,7 +62,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -66,7 +70,7 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H index 5109788ea1..5c4fa0adc1 100644 --- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H +++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H @@ -1,7 +1,8 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); if (pimple.transonic()) { @@ -10,7 +11,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - ((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)) + ((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U)) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -36,19 +37,22 @@ if (pimple.transonic()) } else { - phi = + surfaceScalarField phiHbyA + ( + "phiHbyA", fvc::interpolate(rho) *( - ((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)) + ((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U)) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) == parcels.Srho() @@ -58,7 +62,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -66,7 +70,7 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H index 6fa1965a24..34898172ce 100644 --- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H +++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H @@ -9,11 +9,13 @@ )/psi; } - surfaceScalarField rhof(fvc::interpolate(rho, "rhof")); + surfaceScalarField rhof("rhof", fvc::interpolate(rho)); volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rAUf("rAUf", rhof*fvc::interpolate(rAU)); - volVectorField HbyA("HbyA", rAU*UEqn.H()); + + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); phiv = (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phiv); @@ -22,8 +24,6 @@ phiv -= phiGradp/rhof; - #include "resetPhivPatches.H" - while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn diff --git a/applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H b/applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H deleted file mode 100644 index e7d0c2f93e..0000000000 --- a/applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H +++ /dev/null @@ -1,15 +0,0 @@ -fvsPatchScalarFieldField& phiPatches = phi.boundaryField(); -const fvPatchScalarFieldField& rhoPatches = rho.boundaryField(); -const fvPatchVectorFieldField& Upatches = U.boundaryField(); -const fvsPatchVectorFieldField& SfPatches = mesh.Sf().boundaryField(); - -forAll(phiPatches, patchI) -{ - if (phi.boundaryField().types()[patchI] == "calculated") - { - calculatedFvsPatchScalarField& phiPatch = - refCast(phiPatches[patchI]); - - phiPatch == ((rhoPatches[patchI]*Upatches[patchI]) & SfPatches[patchI]); - } -} diff --git a/applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H b/applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H deleted file mode 100644 index 2f62dcfc37..0000000000 --- a/applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H +++ /dev/null @@ -1,19 +0,0 @@ -surfaceScalarField::GeometricBoundaryField& phivPatches = - phiv.boundaryField(); - -const volVectorField::GeometricBoundaryField& Upatches = - U.boundaryField(); - -const surfaceVectorField::GeometricBoundaryField& SfPatches = - mesh.Sf().boundaryField(); - -forAll(phivPatches, patchI) -{ - if (phiv.boundaryField().types()[patchI] == "calculated") - { - calculatedFvsPatchScalarField& phivPatch = - refCast(phivPatches[patchI]); - - phivPatch == (Upatches[patchI] & SfPatches[patchI]); - } -} From 7fb932f697912469d99049517614d302d9346237 Mon Sep 17 00:00:00 2001 From: laurence Date: Mon, 5 Mar 2012 09:23:38 +0000 Subject: [PATCH 4/4] BUG: cvMesh: revert the quick circumsphere test in buildParallelInterface --- .../conformalVoronoiMeshConformToSurface.C | 62 +++++++++++-------- 1 file changed, 36 insertions(+), 26 deletions(-) diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C index 79a79c3883..9bc2512f58 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C @@ -41,10 +41,10 @@ void Foam::conformalVoronoiMesh::conformToSurface() { reconformationMode reconfMode = reconformationControl(); - if (Pstream::parRun()) - { - seedProcessorBoundarySurfaces(true); - } +// if (Pstream::parRun()) +// { +// seedProcessorBoundarySurfaces(true); +// } if (reconfMode == rmNone) { @@ -74,15 +74,15 @@ void Foam::conformalVoronoiMesh::conformToSurface() storeSurfaceConformation(); } - if (Pstream::parRun()) - { - label nFarPoints = removeProcessorBoundarySeeds(true); - - reduce(nFarPoints, sumOp