Merge branch 'master' into molecularDynamics

.gitignore

@@ -7,6 +7,7 @@
 *.orig
 *.bak
 \#*\#
+.directory
 
 # CVS recovered versions - anywhere
 .#*
@@ -47,7 +48,9 @@ doc/[Dd]oxygen/man
 # source packages - anywhere
 *.tar.bz2
 *.tar.gz
+*.tar
 *.tgz
+*.gtgz
 
 # ignore the persistent .build tag in the main directory
 /.build

README

@@ -14,11 +14,10 @@
   General Public License terms under which you can copy the files.
 
 * System requirements
-  OpenFOAM is developed and tested on Linux, but should work with other Unix
-  style systems.  To check your system setup, execute the foamSystemCheck script
-  in the bin/ directory of the OpenFOAM installation. If no problems are
-  reported, proceed to "3. Installation"; otherwise contact your system
-  administrator.
+  OpenFOAM is developed and tested on Linux, but should work with other POSIX
+  systems.  To check your system setup, execute the foamSystemCheck script in
+  the bin/ directory of the OpenFOAM installation. If no problems are reported,
+  proceed to "3. Installation"; otherwise contact your system administrator.
 
   If the user wishes to run OpenFOAM in 32/64-bit mode they should consult the
   section "Running OpenFOAM in 32-bit mode".

applications/solvers/Lagrangian/kinematicParcelFoam/Make/files

@@ -1,3 +0,0 @@
-kinematicParcelFoam.C
-
-EXE = $(FOAM_APPBIN)/kinematicParcelFoam

applications/solvers/basic/potentialFoam/potentialFoam.C

@@ -97,7 +97,7 @@ int main(int argc, char *argv[])
     U.write();
     phi.write();
 
-    if (args.options().found("writep"))
+    if (args.optionFound("writep"))
     {
         p.write();
     }

applications/solvers/combustion/coalChemistryFoam/Make/files

@@ -0,0 +1,3 @@
+coalChemistryFoam.C
+
+EXE = $(FOAM_APPBIN)/coalChemistryFoam

applications/solvers/combustion/coalChemistryFoam/Make/options

@@ -0,0 +1,39 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/pdfs/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lcoalCombustion\
+    -lspecie \
+    -lbasicThermophysicalModels \
+    -lliquids \
+    -lliquidMixture \
+    -lsolids \
+    -lsolidMixture \
+    -lthermophysicalFunctions \
+    -lcombustionThermophysicalModels \
+    -lchemistryModel \
+    -lradiation \
+    -lODE

applications/solvers/combustion/coalChemistryFoam/UEqn.H

@@ -0,0 +1,17 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho.dimensionedInternalField()*g
+      + coalParcels.SU()
+      + limestoneParcels.SU()
+    );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }

applications/solvers/combustion/coalChemistryFoam/YEqn.H

@@ -0,0 +1,43 @@
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+
+{
+    label inertIndex = -1;
+    volScalarField Yt = 0.0*Y[0];
+
+    for (label i=0; i<Y.size(); i++)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+            solve
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+              ==
+                coalParcels.Srho(i)
+              + kappa*chemistry.RR(i)().dimensionedInternalField()
+            );
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}

applications/solvers/combustion/coalChemistryFoam/additionalOutput.H

@@ -0,0 +1,48 @@
+{
+    tmp<volScalarField> tdQ
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "dQ",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar
+            (
+                "zero",
+                dimensionSet(1, -3, -1, 0, 0, 0, 0),
+                0.0
+            )
+        )
+    );
+
+    scalarField& dQ = tdQ();
+
+    scalarField cp(dQ.size(), 0.0);
+
+    forAll(Y, i)
+    {
+        volScalarField RRi = chemistry.RR(i);
+
+        forAll(h, celli)
+        {
+            scalar Ti = T[celli];
+            cp[celli] += Y[i][celli]*chemistry.specieThermo()[i].Cp(Ti);
+            scalar hi = chemistry.specieThermo()[i].h(Ti);
+            scalar RR = RRi[celli];
+            dQ[celli] -= hi*RR;
+        }
+    }
+
+    forAll(dQ, celli)
+    {
+        dQ[celli] /= cp[celli];
+    }
+
+    tdQ().write();
+}

applications/solvers/combustion/coalChemistryFoam/chemistry.H

@@ -0,0 +1,25 @@
+{
+    Info << "Solving chemistry" << endl;
+
+    chemistry.solve
+    (
+        runTime.value() - runTime.deltaT().value(),
+        runTime.deltaT().value()
+    );
+
+    // turbulent time scale
+    if (turbulentReaction)
+    {
+        DimensionedField<scalar, volMesh> tk =
+            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+        DimensionedField<scalar, volMesh> tc =
+            chemistry.tc()().dimensionedInternalField();
+
+        // Chalmers PaSR model
+        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
+    }
+    else
+    {
+        kappa = 1.0;
+    }
+}

applications/solvers/combustion/coalChemistryFoam/coalChemistryFoam.C

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    coalChemistryFoam
+
+Description
+    Transient solver for compressible, turbulent flow with coal and
+    limestone parcel injections, and combustion.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "hCombustionThermo.H"
+#include "turbulenceModel.H"
+#include "basicThermoCloud.H"
+#include "CoalCloud.H"
+#include "chemistryModel.H"
+#include "chemistrySolver.H"
+#include "reactingThermoTypes.H"
+#include "timeActivatedExplicitCellSource.H"
+#include "radiationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readChemistryProperties.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "createClouds.H"
+    #include "createRadiationModel.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        rhoEffLagrangian = coalParcels.rhoEff() + limestoneParcels.rhoEff();
+        pDyn = 0.5*rho*magSqr(U);
+
+        coalParcels.evolve();
+
+        coalParcels.info();
+
+        Info<< endl;
+
+        limestoneParcels.evolve();
+
+        limestoneParcels.info();
+
+        Info<< endl;
+
+        #include "chemistry.H"
+        #include "rhoEqn.H"
+
+        // --- PIMPLE loop
+        for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
+        {
+            #include "UEqn.H"
+            #include "YEqn.H"
+            #include "hEqn.H"
+
+            // --- PISO loop
+            for (int corr=1; corr<=nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+
+            Info<< "T gas min/max   = " << min(T).value() << ", "
+                << max(T).value() << endl;
+        }
+
+        turbulence->correct();
+
+        rho = thermo->rho();
+
+        if (runTime.write())
+        {
+            #include "additionalOutput.H"
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/coalChemistryFoam/createClouds.H

@@ -0,0 +1,19 @@
+Info<< "\nConstructing coal cloud" << endl;
+CoalCloud<specieReactingProperties> coalParcels
+(
+    "coalCloud1",
+    rho,
+    U,
+    g,
+    thermo()
+);
+
+Info<< "\nConstructing limestone cloud" << endl;
+basicThermoCloud limestoneParcels
+(
+    "limestoneCloud1",
+    rho,
+    U,
+    g,
+    thermo()
+);

applications/solvers/combustion/coalChemistryFoam/createFields.H

@@ -0,0 +1,129 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<hCombustionThermo> thermo
+    (
+        hCombustionThermo::New(mesh)
+    );
+
+    combustionMixture& composition = thermo->composition();
+    PtrList<volScalarField>& Y = composition.Y();
+
+    word inertSpecie(thermo->lookup("inertSpecie"));
+
+    volScalarField& p = thermo->p();
+    volScalarField& h = thermo->h();
+    const volScalarField& T = thermo->T();
+    const volScalarField& psi = thermo->psi();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo->rho()
+    );
+
+    // lagrangian effective density field - used externally (optional)
+    volScalarField rhoEffLagrangian
+    (
+        IOobject
+        (
+            "rhoEffLagrangian",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimDensity, 0.0)
+    );
+
+    // dynamic pressure field - used externally (optional)
+    volScalarField pDyn
+    (
+        IOobject
+        (
+            "pDyn",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimPressure, 0.0)
+    );
+
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "compressibleCreatePhi.H"
+
+    DimensionedField<scalar, volMesh> kappa
+    (
+        IOobject
+        (
+            "kappa",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimless, 0.0)
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo()
+        )
+    );
+
+    Info<< "Creating field DpDt\n" << endl;
+    volScalarField DpDt =
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+    Info<< "\nConstructing explicit enthalpy cell source" << endl;
+    timeActivatedExplicitCellSource enthalpySource
+    (
+        "enthalpySource",
+        mesh,
+        h.dimensions()*phi.dimensions()/mesh.V().dimensions()
+    );
+
+    Info << "Constructing chemical mechanism" << endl;
+    chemistryModel chemistry
+    (
+        thermo(),
+        rho
+    );
+
+    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+    forAll (Y, i)
+    {
+        fields.add(Y[i]);
+    }
+    fields.add(h);

applications/solvers/combustion/coalChemistryFoam/hEqn.H

@@ -0,0 +1,22 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::ddt(rho, h)
+      + fvm::div(phi, h)
+      - fvm::laplacian(turbulence->alphaEff(), h)
+     ==
+        DpDt
+     +  coalParcels.Sh()
+     +  limestoneParcels.Sh()
+     +  enthalpySource.Su()
+     +  radiation->Sh(thermo())
+    );
+
+    hEqn.relax();
+
+    hEqn.solve();
+
+    thermo->correct();
+
+    radiation->correct();
+}

applications/solvers/combustion/coalChemistryFoam/pEqn.H

@@ -0,0 +1,72 @@
+rho = thermo->rho();
+
+volScalarField rUA = 1.0/UEqn.A();
+U = rUA*UEqn.H();
+
+if (transonic)
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(thermo->psi())
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho*rUA, p)
+         ==
+            coalParcels.Srho()
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    phi =
+        fvc::interpolate(rho)*
+        (
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rUA, p)
+         ==
+            coalParcels.Srho()
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

applications/solvers/combustion/coalChemistryFoam/readChemistryProperties.H

@@ -0,0 +1,22 @@
+Info<< "Reading chemistry properties\n" << endl;
+
+IOdictionary chemistryProperties
+(
+    IOobject
+    (
+        "chemistryProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction"));
+
+dimensionedScalar Cmix("Cmix", dimless, 1.0);
+
+if (turbulentReaction)
+{
+    chemistryProperties.lookup("Cmix") >> Cmix;
+}

applications/solvers/combustion/coalChemistryFoam/rhoEqn.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2007 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,17 +22,22 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
 \*---------------------------------------------------------------------------*/
 
-#include "upwindCECStencilObject.H"
-
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
 {
-    defineTypeNameAndDebug(upwindCECStencilObject, 0);
+    solve
+    (
+        fvm::ddt(rho)
+      + fvc::div(phi)
+      ==
+        coalParcels.Srho()
+    );
 }
 
-
 // ************************************************************************* //

applications/solvers/combustion/reactingFoam/createFields.H

@@ -37,7 +37,6 @@ volVectorField U
 
 volScalarField& p = thermo->p();
 const volScalarField& psi = thermo->psi();
-const volScalarField& T = thermo->T();
 volScalarField& h = thermo->h();
 
 
@@ -99,5 +98,5 @@ volScalarField dQ
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("zero", dimensionSet(1,-3,-1,0,0,0,0), 0.0)
+    dimensionedScalar("zero", dimensionSet(1, -3, -1, 0, 0, 0, 0), 0.0)
 );

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/buoyantBoussinesqPisoFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/heatTransfer/buoyantPisoFoam/UEqn.H

@@ -15,12 +15,10 @@
         (
             UEqn
          ==
-           -fvc::reconstruct
+            fvc::reconstruct
             (
-                (
-                    fvc::snGrad(pd)
-                  + ghf*fvc::snGrad(rho)
-                ) * mesh.magSf()
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              - fvc::snGrad(p)*mesh.magSf()
             )
         );
     }

applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H

@@ -59,27 +59,6 @@
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
     );
 
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-    dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
-
-    Info<< "Creating field pd\n" << endl;
-    volScalarField pd
-    (
-        IOobject
-        (
-            "pd",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    p = pd + rho*gh + pRef;
     thermo->correct();
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);

applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H

@@ -1,5 +1,5 @@
 {
-    bool closedVolume = pd.needReference();
+    bool closedVolume = p.needReference();
 
     rho = thermo->rho();
 
@@ -17,38 +17,35 @@
         )
     );
 
-    phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
+    phi = phiU + rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, pd)
-          + fvc::ddt(psi)*pRef
-          + fvc::ddt(psi, rho)*gh
+            fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rhorUAf, pd)
+          - fvm::laplacian(rhorUAf, p)
         );
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi += pdEqn.flux();
+            phi += pEqn.flux();
         }
     }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
     U.correctBoundaryConditions();
 
-    p == pd + rho*gh + pRef;
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
     #include "rhoEqn.H"
@@ -62,6 +59,4 @@
             /fvc::domainIntegrate(thermo->psi());
         rho = thermo->rho();
     }
-
-    pd == p - (rho*gh + pRef);
 }

applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H

@@ -13,14 +13,11 @@
     (
         UEqn()
      ==
-       -fvc::reconstruct
+        fvc::reconstruct
         (
-            (
-                fvc::snGrad(pd)
-              + ghf*fvc::snGrad(rho)
-            ) * mesh.magSf()
+            fvc::interpolate(rho)*(g & mesh.Sf())
+          - fvc::snGrad(p)*mesh.magSf()
         )
     ).initialResidual();
 
     maxResidual = max(eqnResidual, maxResidual);
-

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C

@@ -39,15 +39,14 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
 
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -55,17 +54,17 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
-#       include "initConvergenceCheck.H"
+        #include "readSIMPLEControls.H"
+        #include "initConvergenceCheck.H"
 
-        pd.storePrevIter();
+        p.storePrevIter();
         rho.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
         {
-#           include "UEqn.H"
-#           include "hEqn.H"
-#           include "pEqn.H"
+            #include "UEqn.H"
+            #include "hEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();
@@ -76,7 +75,7 @@ int main(int argc, char *argv[])
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
 
-#       include "convergenceCheck.H"
+        #include "convergenceCheck.H"
     }
 
     Info<< "End\n" << endl;

applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H

@@ -51,38 +51,16 @@
         )
     );
 
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-    dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
-
-    Info<< "Creating field pd\n" << endl;
-    volScalarField pd
-    (
-        IOobject
-        (
-            "pd",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    p = pd + rho*gh + pRef;
     thermo->correct();
 
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
     setRefCell
     (
-        pd,
+        p,
         mesh.solutionDict().subDict("SIMPLE"),
-        pdRefCell,
-        pdRefValue
+        pRefCell,
+        pRefValue
     );
 
 

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -1,4 +1,6 @@
 {
+    rho = thermo->rho();
+
     volScalarField rUA = 1.0/UEqn().A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
 
@@ -7,60 +9,57 @@
 
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     bool closedVolume = adjustPhi(phi, U, p);
-    surfaceScalarField buoyancyPhi = ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
-    phi -= buoyancyPhi;
+    surfaceScalarField buoyancyPhi =
+        rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
+    phi += buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::laplacian(rhorUAf, pd) == fvc::div(phi)
+            fvm::laplacian(rhorUAf, p) == fvc::div(phi)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, p[pRefCell]);
 
         // retain the residual from the first iteration
         if (nonOrth == 0)
         {
-            eqnResidual = pdEqn.solve().initialResidual();
+            eqnResidual = pEqn.solve().initialResidual();
             maxResidual = max(eqnResidual, maxResidual);
         }
         else
         {
-            pdEqn.solve();
+            pEqn.solve();
         }
 
         if (nonOrth == nNonOrthCorr)
         {
+            // For closed-volume cases adjust the pressure and density levels
+            // to obey overall mass continuity
+            if (closedVolume)
+            {
+                p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
+                    /fvc::domainIntegrate(thermo->psi());
+            }
+
             // Calculate the conservative fluxes
-            phi -= pdEqn.flux();
+            phi -= pEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
-            pd.relax();
+            p.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rhorUAf);
+            U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rhorUAf);
             U.correctBoundaryConditions();
         }
     }
 
     #include "continuityErrs.H"
 
-    p == pd + rho*gh + pRef;
-
-    // For closed-volume cases adjust the pressure and density levels
-    // to obey overall mass continuity
-    if (closedVolume)
-    {
-        p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
-            /fvc::domainIntegrate(thermo->psi());
-    }
-
     rho = thermo->rho();
     rho.relax();
     Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
         << endl;
-
-    pd == p - (rho*gh + pRef);
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files

@@ -1,10 +1,6 @@
 regionProperties/regionProperties.C
 
-coupleManager/coupleManager.C
-
 derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
-derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
-derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C
 derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
 
 fluid/compressibleCourantNo.C

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options

@@ -2,9 +2,7 @@ EXE_INC = \
     -Ifluid \
     -Isolid \
     -IregionProperties \
-    -IcoupleManager \
-    -IderivedFvPatchFields/solidWallTemperatureCoupled \
-    -IderivedFvPatchFields/solidWallHeatFluxTemperatureCoupled \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -67,6 +67,7 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
 #       include "readTimeControls.H"
+#       include "readPIMPLEControls.H"
 
         if (fluidRegions.size())
         {
@@ -78,22 +79,36 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        forAll(fluidRegions, i)
+        if (nOuterCorr != 1)
         {
-            Info<< "\nSolving for fluid region "
-                << fluidRegions[i].name() << endl;
-#           include "setRegionFluidFields.H"
-#           include "readFluidMultiRegionPISOControls.H"
-#           include "solveFluid.H"
+            forAll(fluidRegions, i)
+            {
+#               include "setRegionFluidFields.H"
+#               include "storeOldFluidFields.H"
+            }
         }
 
-        forAll(solidRegions, i)
+
+        // --- PIMPLE loop
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
-            Info<< "\nSolving for solid region "
-                << solidRegions[i].name() << endl;
-#           include "setRegionSolidFields.H"
-#           include "readSolidMultiRegionPISOControls.H"
-#           include "solveSolid.H"
+            forAll(fluidRegions, i)
+            {
+                Info<< "\nSolving for fluid region "
+                    << fluidRegions[i].name() << endl;
+#               include "setRegionFluidFields.H"
+#               include "readFluidMultiRegionPIMPLEControls.H"
+#               include "solveFluid.H"
+            }
+
+            forAll(solidRegions, i)
+            {
+                Info<< "\nSolving for solid region "
+                    << solidRegions[i].name() << endl;
+#               include "setRegionSolidFields.H"
+#               include "readSolidMultiRegionPIMPLEControls.H"
+#               include "solveSolid.H"
+            }
         }
 
         runTime.write();

applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.C

@@ -1,186 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "coupleManager.H"
-#include "OFstream.H"
-#include "regionProperties.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::coupleManager::coupleManager(const fvPatch& patch)
-:
-    patch_(patch),
-    neighbourRegionName_("undefined-neighbourRegionName"),
-    neighbourPatchName_("undefined-neighbourPatchName"),
-    neighbourFieldName_("undefined-neighbourFieldName"),
-    localRegion_(patch_.boundaryMesh().mesh())
-{}
-
-
-Foam::coupleManager::coupleManager
-(
-    const fvPatch& patch,
-    const dictionary& dict
-)
-:
-    patch_(patch),
-    neighbourRegionName_(dict.lookup("neighbourRegionName")),
-    neighbourPatchName_(dict.lookup("neighbourPatchName")),
-    neighbourFieldName_(dict.lookup("neighbourFieldName")),
-    localRegion_(patch_.boundaryMesh().mesh())
-{}
-
-
-Foam::coupleManager::coupleManager
-(
-    const coupleManager& cm
-)
-:
-    patch_(cm.patch()),
-    neighbourRegionName_(cm.neighbourRegionName()),
-    neighbourPatchName_(cm.neighbourPatchName()),
-    neighbourFieldName_(cm.neighbourFieldName()),
-    localRegion_(patch_.boundaryMesh().mesh())
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::coupleManager::~coupleManager()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::coupleManager::regionOwner() const
-{
-    const fvMesh& nbrRegion = neighbourRegion();
-
-    const regionProperties& props =
-        localRegion_.objectRegistry::parent().lookupObject<regionProperties>
-        (
-            "regionProperties"
-        );
-
-    label myIndex = findIndex(props.fluidRegionNames(), localRegion_.name());
-    if (myIndex == -1)
-    {
-        label i = findIndex(props.solidRegionNames(), localRegion_.name());
-
-        if (i == -1)
-        {
-            FatalErrorIn("coupleManager::regionOwner() const")
-                << "Cannot find region " << localRegion_.name()
-                << " neither in fluids " << props.fluidRegionNames()
-                << " nor in solids " << props.solidRegionNames()
-                << exit(FatalError);
-        }
-        myIndex = props.fluidRegionNames().size() + i;
-    }
-    label nbrIndex = findIndex(props.fluidRegionNames(), nbrRegion.name());
-    if (nbrIndex == -1)
-    {
-        label i = findIndex(props.solidRegionNames(), nbrRegion.name());
-
-        if (i == -1)
-        {
-            FatalErrorIn("coupleManager::regionOwner() const")
-                << "Cannot find region " << nbrRegion.name()
-                << " neither in fluids " << props.fluidRegionNames()
-                << " nor in solids " << props.solidRegionNames()
-                << exit(FatalError);
-        }
-        nbrIndex = props.fluidRegionNames().size() + i;
-    }
-
-    return myIndex < nbrIndex;
-}
-
-
-void Foam::coupleManager::checkCouple() const
-{
-    Info<< "neighbourRegionName_ = " << neighbourRegionName_ << endl;
-    Info<< "neighbourPatchName_ = " << neighbourPatchName_ << endl;
-    Info<< "neighbourFieldName_ = " << neighbourFieldName_ << endl;
-
-    const fvPatch& nPatch = neighbourPatch();
-
-    if (patch_.size() != nPatch.size())
-    {
-        FatalErrorIn("Foam::coupleManager::checkCouple()")
-            << "Unequal patch sizes:" << nl
-            << "    patch name (size) = " << patch_.name()
-            << "(" << patch_.size() << ")" << nl
-            << "    neighbour patch name (size) = "
-            << nPatch.name() << "(" << patch_.size() << ")" << nl
-            << abort(FatalError);
-    }
-}
-
-
-void Foam::coupleManager::coupleToObj() const
-{
-    const fvPatch& nPatch = neighbourPatch();
-
-    OFstream obj
-        (
-             patch_.name() + "_to_" + nPatch.name() + "_couple.obj"
-        );
-    const vectorField& c1 = patch_.Cf();
-    const vectorField& c2 = neighbourPatch().Cf();
-
-    if (c1.size() != c2.size())
-    {
-        FatalErrorIn("coupleManager::coupleToObj() const")
-            << "Coupled patches are of unequal size:" << nl
-            << "    patch0 = " << patch_.name()
-            << "(" << patch_.size() <<  ")" << nl
-            << "    patch1 = " << nPatch.name()
-            << "(" << nPatch.size() <<  ")" << nl
-            << abort(FatalError);
-    }
-
-    forAll(c1, i)
-    {
-        obj << "v " << c1[i].x() << " " << c1[i].y() << " " << c1[i].z() << nl
-            << "v " << c2[i].x() << " " << c2[i].y() << " " << c2[i].z() << nl
-            << "l " << (2*i + 1) << " " << (2*i + 2) << endl;
-    }
-}
-
-
-void Foam::coupleManager::writeEntries(Ostream& os) const
-{
-    os.writeKeyword("neighbourRegionName");
-    os << neighbourRegionName_ << token::END_STATEMENT << nl;
-    os.writeKeyword("neighbourPatchName");
-    os << neighbourPatchName_ << token::END_STATEMENT << nl;
-    os.writeKeyword("neighbourFieldName");
-    os << neighbourFieldName_ << token::END_STATEMENT << nl;
-}
-
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.H

@@ -1,170 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    coupleManager
-
-Description
-    Object to handle the coupling of region patches. It can be queried to
-    return the neighbour information.
-
-SourceFiles
-    coupleManager.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef coupleManager_H
-#define coupleManager_H
-
-#include "Ostream.H"
-#include "dictionary.H"
-#include "fvPatch.H"
-#include "fvMesh.H"
-#include "Time.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                         Class coupleManager Declaration
-\*---------------------------------------------------------------------------*/
-
-class coupleManager
-{
-    // Private data
-
-        //- Reference to the local fvPatch
-        const fvPatch& patch_;
-
-        //- Name of neighbour region
-        word neighbourRegionName_;
-
-        //- Name of patch on the neighbour region
-        word neighbourPatchName_;
-
-        //- Name of field on the neighbour region
-        word neighbourFieldName_;
-
-        //- Reference to the local region
-        const fvMesh& localRegion_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise assignment
-        void operator=(const coupleManager&);
-
-
-public:
-
-    // Constructors
-
-        //- Construct from fvPatch
-        coupleManager(const fvPatch& patch);
-
-        //- Construct from fvPatch and dictionary
-        coupleManager(const fvPatch& patch, const dictionary& dict);
-
-        //- Copy constructor
-        coupleManager(const coupleManager& cm);
-
-
-    // Destructor
-
-        ~coupleManager();
-
-
-    // Member Functions
-
-        // Access
-
-            //- Return a reference to the local patch
-            inline const fvPatch& patch() const;
-
-            //- Return the name of the neighbour region
-            inline const word& neighbourRegionName() const;
-
-            //- Return the name of the patch on the neighbour region
-            inline const word& neighbourPatchName() const;
-
-            //- Return the name of the field on the neighbour region
-            inline const word& neighbourFieldName() const;
-
-            //- Return a reference to the neighbour mesh
-            inline const fvMesh& neighbourRegion() const;
-
-            //- Return the neighbour patch ID
-            inline label neighbourPatchID() const;
-
-            //- Return a reference to the neighbour patch
-            inline const fvPatch& neighbourPatch() const;
-
-            //- Return a reference to the neighbour patch field
-            template<class Type>
-            inline const fvPatchField<Type>& neighbourPatchField() const;
-
-            //- Am I owner (= first to evaluate) of this region interface? 
-            bool regionOwner() const;
-
-            //- Check that the couple is valid
-            void checkCouple() const;
-
-
-        // Edit
-
-            //- Return the name of the neighbour region
-            word& neighbourRegionName();
-
-            //- Return the name of the patch on the neighbour region
-            word& neighbourPatchName();
-
-            //- Return the name of the field on the neighbour region
-            word& neighbourFieldName();
-
-
-        // Write
-
-           //- Write couple to obj file
-           void coupleToObj() const;
-
-           //- Write entries for patches
-           void writeEntries(Ostream& os) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "coupleManagerI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManagerI.H

@@ -1,98 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-inline const Foam::fvPatch& Foam::coupleManager::patch() const
-{
-    return patch_;
-};
-
-
-inline const Foam::word& Foam::coupleManager::neighbourRegionName() const
-{
-    return neighbourRegionName_;
-};
-
-
-inline const Foam::word& Foam::coupleManager::neighbourPatchName() const
-{
-    return neighbourPatchName_;
-};
-
-
-inline const Foam::word& Foam::coupleManager::neighbourFieldName() const
-{
-     return neighbourFieldName_;
-};
-
-
-inline const Foam::fvMesh& Foam::coupleManager::neighbourRegion() const
-{
-    return localRegion_.objectRegistry::parent()
-        .lookupObject<fvMesh>(neighbourRegionName_);
-}
-
-
-inline Foam::label Foam::coupleManager::neighbourPatchID() const
-{
-    return neighbourRegion().boundaryMesh().findPatchID(neighbourPatchName_);
-}
-
-
-inline const Foam::fvPatch& Foam::coupleManager::neighbourPatch() const
-{
-    return neighbourRegion().boundary()[neighbourPatchID()];
-}
-
-
-template<class Type>
-inline const Foam::fvPatchField<Type>&
-Foam::coupleManager::neighbourPatchField() const
-{
-    return neighbourPatch().lookupPatchField
-        <GeometricField<Type, fvPatchField, volMesh>, Type>
-            (neighbourFieldName_);
-}
-
-
-inline Foam::word& Foam::coupleManager::neighbourRegionName()
-{
-    return neighbourRegionName_;
-};
-
-
-inline Foam::word& Foam::coupleManager::neighbourPatchName()
-{
-    return neighbourPatchName_;
-};
-
-
-inline Foam::word& Foam::coupleManager::neighbourFieldName()
-{
-    return neighbourFieldName_;
-};
-
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C

@@ -1,150 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "solidWallHeatFluxTemperatureCoupledFvPatchScalarField.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-
-#include "solidWallTemperatureCoupledFvPatchScalarField.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
-solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(p, iF),
-    coupleManager_(p),
-    KName_("undefined-K")
-{}
-
-
-Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
-solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-(
-    const solidWallHeatFluxTemperatureCoupledFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
-    coupleManager_(ptf.coupleManager_),
-    KName_(ptf.KName_)
-{}
-
-
-Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
-solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedGradientFvPatchScalarField(p, iF),
-    coupleManager_(p, dict),
-    KName_(dict.lookup("K"))
-{
-    if (dict.found("value"))
-    {
-        fvPatchField<scalar>::operator=
-        (
-            scalarField("value", dict, p.size())
-        );
-    }
-    else
-    {
-        evaluate();
-    }
-}
-
-
-Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
-solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-(
-    const solidWallHeatFluxTemperatureCoupledFvPatchScalarField& whftcsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(whftcsf, iF),
-    coupleManager_(whftcsf.coupleManager_),
-    KName_(whftcsf.KName_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const fvPatchField<scalar>& neighbourField =
-        coupleManager_.neighbourPatchField<scalar>();
-
-    const fvPatchField<scalar>& K =
-        patch().lookupPatchField<volScalarField, scalar>(KName_);
-
-    gradient() = -refCast<const solidWallTemperatureCoupledFvPatchScalarField>
-        (neighbourField).flux()/K;
-
-    fixedGradientFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    fvPatchScalarField::write(os);
-    coupleManager_.writeEntries(os);
-    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
-    writeEntry("value", os);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-makePatchTypeField
-(
-    fvPatchScalarField,
-    solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-);
-
-} // End namespace Foam
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.H

@@ -1,165 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-
-Description
-    Fixed heat-flux boundary condition for temperature, to be used by the
-    conjugate heat transfer solver.
-
-    Example usage:
-        myInterfacePatchName
-        {
-            type                solidWallHeatFluxTemperatureCoupled;
-            neighbourRegionName fluid;
-            neighbourPatchName  fluidSolidInterface;
-            neighbourFieldName  T;
-            K                   K;
-            value               uniform 300;
-        }
-
-SourceFiles
-    solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef solidWallHeatFluxTemperatureCoupledFvPatchScalarField_H
-#define solidWallHeatFluxTemperatureCoupledFvPatchScalarField_H
-
-#include "fvPatchFields.H"
-#include "fixedGradientFvPatchFields.H"
-#include "coupleManager.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-   Class solidWallHeatFluxTemperatureCoupledFvPatchScalarField Declaration
-\*---------------------------------------------------------------------------*/
-
-class solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-:
-    public fixedGradientFvPatchScalarField
-{
-    // Private data
-
-        //- Couple manager object
-        coupleManager coupleManager_;
-
-        //- Name of thermal conductivity field
-        word KName_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("solidWallHeatFluxTemperatureCoupled");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given
-        // solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        //  onto a new patch
-        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallHeatFluxTemperatureCoupledFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-                    (
-                        *this
-                    )
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallHeatFluxTemperatureCoupledFvPatchScalarField&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-                    (
-                        *this,
-                        iF
-                    )
-            );
-        }
-
-
-    // Member functions
-
-        //- Update the coefficients associated with the patch field
-        virtual void updateCoeffs();
-
-        //- Write
-        virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C

@@ -28,8 +28,62 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "fvPatchFieldMapper.H"
 #include "volFields.H"
+#include "directMappedPatchBase.H"
 #include "regionProperties.H"
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::interfaceOwner
+(
+    const polyMesh& nbrRegion
+) const
+{
+    const fvMesh& myRegion = patch().boundaryMesh().mesh();
+
+    const regionProperties& props =
+        myRegion.objectRegistry::parent().lookupObject<regionProperties>
+        (
+            "regionProperties"
+        );
+
+    label myIndex = findIndex(props.fluidRegionNames(), myRegion.name());
+    if (myIndex == -1)
+    {
+        label i = findIndex(props.solidRegionNames(), myRegion.name());
+
+        if (i == -1)
+        {
+            FatalErrorIn
+            (
+                "solidWallMixedTemperatureCoupledFvPatchScalarField"
+                "::interfaceOwner(const polyMesh&) const"
+            )   << "Cannot find region " << myRegion.name()
+                << " neither in fluids " << props.fluidRegionNames()
+                << " nor in solids " << props.solidRegionNames()
+                << exit(FatalError);
+        }
+        myIndex = props.fluidRegionNames().size() + i;
+    }
+    label nbrIndex = findIndex(props.fluidRegionNames(), nbrRegion.name());
+    if (nbrIndex == -1)
+    {
+        label i = findIndex(props.solidRegionNames(), nbrRegion.name());
+
+        if (i == -1)
+        {
+            FatalErrorIn("coupleManager::interfaceOwner(const polyMesh&) const")
+                << "Cannot find region " << nbrRegion.name()
+                << " neither in fluids " << props.fluidRegionNames()
+                << " nor in solids " << props.solidRegionNames()
+                << exit(FatalError);
+        }
+        nbrIndex = props.fluidRegionNames().size() + i;
+    }
+
+    return myIndex < nbrIndex;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
@@ -40,7 +94,7 @@ solidWallMixedTemperatureCoupledFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(p, iF),
-    coupleManager_(p),
+    neighbourFieldName_("undefined-neighbourFieldName"),
     KName_("undefined-K")
 {
     this->refValue() = 0.0;
@@ -60,7 +114,7 @@ solidWallMixedTemperatureCoupledFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(ptf, p, iF, mapper),
-    coupleManager_(ptf.coupleManager_),
+    neighbourFieldName_(ptf.neighbourFieldName_),
     KName_(ptf.KName_),
     fixesValue_(ptf.fixesValue_)
 {}
@@ -75,14 +129,46 @@ solidWallMixedTemperatureCoupledFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(p, iF),
-    coupleManager_(p, dict),
+    neighbourFieldName_(dict.lookup("neighbourFieldName")),
     KName_(dict.lookup("K"))
 {
+    if (!isA<directMappedPatchBase>(this->patch().patch()))
+    {
+        FatalErrorIn
+        (
+            "solidWallMixedTemperatureCoupledFvPatchScalarField::"
+            "solidWallMixedTemperatureCoupledFvPatchScalarField\n"
+            "(\n"
+            "    const fvPatch& p,\n"
+            "    const DimensionedField<scalar, volMesh>& iF,\n"
+            "    const dictionary& dict\n"
+            ")\n"
+        )   << "\n    patch type '" << p.type()
+            << "' not type '" << directMappedPatchBase::typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << dimensionedInternalField().name()
+            << " in file " << dimensionedInternalField().objectPath()
+            << exit(FatalError);
+    }
+
     fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
-    refValue() = static_cast<scalarField>(*this);
-    refGrad() = 0.0;
-    valueFraction() = 1.0;
-    fixesValue_ = true;
+
+    if (dict.found("refValue"))
+    {
+        // Full restart
+        refValue() = scalarField("refValue", dict, p.size());
+        refGrad() = scalarField("refGradient", dict, p.size());
+        valueFraction() = scalarField("valueFraction", dict, p.size());
+        fixesValue_ = readBool(dict.lookup("fixesValue"));
+    }
+    else
+    {
+        // Start from user entered data. Assume fixedValue.
+        refValue() = *this;
+        refGrad() = 0.0;
+        valueFraction() = 1.0;
+        fixesValue_ = true;
+    }
 }
 
 
@@ -94,7 +180,7 @@ solidWallMixedTemperatureCoupledFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(wtcsf, iF),
-    coupleManager_(wtcsf.coupleManager_),
+    neighbourFieldName_(wtcsf.neighbourFieldName_),
     KName_(wtcsf.KName_),
     fixesValue_(wtcsf.fixesValue_)
 {}
@@ -116,33 +202,111 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
         return;
     }
 
+    // Get the coupling information from the directMappedPatchBase
+    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
+    (
+        patch().patch()
+    );
+    const polyMesh& nbrMesh = mpp.sampleMesh();
+
     tmp<scalarField> intFld = patchInternalField();
 
+    if (interfaceOwner(nbrMesh))
+    {
+        // Note: other side information could be cached - it only needs
+        // to be updated the first time round the iteration (i.e. when
+        // switching regions) but unfortunately we don't have this information.
+
+        const mapDistribute& distMap = mpp.map();
+        const fvPatch& nbrPatch = refCast<const fvMesh>
+        (
+            nbrMesh
+        ).boundary()[mpp.samplePolyPatch().index()];
+
+
+        // Calculate the temperature by harmonic averaging
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        const solidWallMixedTemperatureCoupledFvPatchScalarField& nbrField =
+        refCast<const solidWallMixedTemperatureCoupledFvPatchScalarField>
+        (
+            nbrPatch.lookupPatchField<volScalarField, scalar>
+            (
+                neighbourFieldName_
+            )
+        );
+
+        // Swap to obtain full local values of neighbour internal field
+        scalarField nbrIntFld = nbrField.patchInternalField();
+        mapDistribute::distribute
+        (
+            Pstream::defaultCommsType,
+            distMap.schedule(),
+            distMap.constructSize(),
+            distMap.subMap(),           // what to send
+            distMap.constructMap(),     // what to receive
+            nbrIntFld
+        );
+
+        // Swap to obtain full local values of neighbour K*delta
+        scalarField nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
+        mapDistribute::distribute
+        (
+            Pstream::defaultCommsType,
+            distMap.schedule(),
+            distMap.constructSize(),
+            distMap.subMap(),           // what to send
+            distMap.constructMap(),     // what to receive
+            nbrKDelta
+        );
+
+
+        tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
+
+        // Calculate common wall temperature. Reuse *this to store common value.
+        scalarField Twall
+        (
+            (myKDelta()*intFld() + nbrKDelta*nbrIntFld)
+          / (myKDelta() + nbrKDelta)
+        );
+        // Assign to me
+        fvPatchScalarField::operator=(Twall);
+        // Distribute back and assign to neighbour
+        mapDistribute::distribute
+        (
+            Pstream::defaultCommsType,
+            distMap.schedule(),
+            nbrField.size(),
+            distMap.constructMap(),     // reverse : what to send
+            distMap.subMap(),
+            Twall
+        );
+        const_cast<solidWallMixedTemperatureCoupledFvPatchScalarField&>
+        (
+            nbrField
+        ).fvPatchScalarField::operator=(Twall);
+    }
+
+
+    // Switch between fixed value (of harmonic avg) or gradient
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     label nFixed = 0;
 
     // Like snGrad but bypass switching on refValue/refGrad.
-    tmp<scalarField> normalGradient =
-        (*this-intFld())
-      * patch().deltaCoeffs();
+    tmp<scalarField> normalGradient = (*this-intFld())*patch().deltaCoeffs();
 
     if (debug)
     {
+        scalar Q = gSum(K()*patch().magSf()*normalGradient());
+
         Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
             << "updateCoeffs() :"
             << " patch:" << patch().name()
+            << " heatFlux:" << Q
             << " walltemperature "
-            << " min:"
-            <<  returnReduce
-                (
-                    (this->size() > 0 ? min(*this) : VGREAT),
-                    minOp<scalar>()
-                )
-            << " max:"
-            <<  returnReduce
-                (
-                    (this->size() > 0 ? max(*this) : -VGREAT),
-                    maxOp<scalar>()
-                )
+            << " min:" << gMin(*this)
+            << " max:" << gMax(*this)
             << " avg:" << gAverage(*this)
             << endl;
     }
@@ -150,7 +314,7 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
     forAll(*this, i)
     {
         // if outgoing flux use fixed value.
-        if (intFld()[i] > operator[](i))
+        if (normalGradient()[i] < 0.0)
         {
             this->refValue()[i] = operator[](i);
             this->refGrad()[i] = 0.0;   // not used
@@ -185,80 +349,16 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
 }
 
 
-void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::evaluate
-(
-    const Pstream::commsTypes
-)
-{
-    if (!this->updated())
-    {
-        this->updateCoeffs();
-    }
-
-    if (!coupleManager_.regionOwner())
-    {
-        // I am the last one to evaluate.
-
-        tmp<scalarField> intFld = patchInternalField();
-
-        const fvPatch& nbrPatch = coupleManager_.neighbourPatch();
-
-        solidWallMixedTemperatureCoupledFvPatchScalarField& nbrField =
-        refCast<solidWallMixedTemperatureCoupledFvPatchScalarField>
-        (
-            const_cast<fvPatchField<scalar>&>
-            (
-                coupleManager_.neighbourPatchField<scalar>()
-            )
-        );
-        tmp<scalarField> nbrIntFld = nbrField.patchInternalField();
-        tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
-        tmp<scalarField> nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
-
-        // Calculate common wall temperature and assign to both sides
-        scalarField::operator=
-        (
-            (myKDelta()*intFld + nbrKDelta()*nbrIntFld)
-          / (myKDelta() + nbrKDelta())
-        );
-
-        nbrField.scalarField::operator=(*this);
-
-        if (debug)
-        {
-            Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
-                << "updateCoeffs() :"
-                << " patch:" << patch().name()
-                << " setting master and slave to wall temperature "
-                << " min:"
-                <<  returnReduce
-                    (
-                        (this->size() > 0 ? min(*this) : VGREAT),
-                        minOp<scalar>()
-                    )
-                << " max:"
-                <<  returnReduce
-                    (
-                        (this->size() > 0 ? max(*this) : -VGREAT),
-                        maxOp<scalar>()
-                    )
-                << " avg:" << gAverage(*this)
-                << endl;
-        }
-    }
-
-    fvPatchScalarField::evaluate();
-}
-
-
 void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::write
 (
     Ostream& os
 ) const
 {
     mixedFvPatchScalarField::write(os);
-    coupleManager_.writeEntries(os);
+    os.writeKeyword("neighbourFieldName")<< neighbourFieldName_
+        << token::END_STATEMENT << nl;
     os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("fixesValue") << fixesValue_ << token::END_STATEMENT << nl;
 }
 
 

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H

@@ -37,13 +37,22 @@ Description
         myInterfacePatchName
         {
             type                solidWallMixedTemperatureCoupled;
-            neighbourRegionName fluid;
-            neighbourPatchName  fluidSolidInterface;
             neighbourFieldName  T;
             K                   K;
             value               uniform 300;
         }
 
+    Needs to be on underlying directMapped(Wall)FvPatch.
+
+    Note: runs in parallel with arbitrary decomposition. Uses directMapped
+    functionality to calculate exchange.
+
+    Note: lags interface data so both sides use same data.
+    - problem: schedule to calculate average would interfere
+    with standard processor swaps.
+    - so: updateCoeffs sets both to same Twall. Only need to do
+    this for last outer iteration but don't have access to this.
+
 SourceFiles
     solidWallMixedTemperatureCoupledFvPatchScalarField.C
 
@@ -54,7 +63,6 @@ SourceFiles
 
 #include "fvPatchFields.H"
 #include "mixedFvPatchFields.H"
-#include "coupleManager.H"
 #include "fvPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -72,14 +80,21 @@ class solidWallMixedTemperatureCoupledFvPatchScalarField
 {
     // Private data
 
-        //- Couple manager object
-        coupleManager coupleManager_;
-
+        //- Name of field on the neighbour region
+        const word neighbourFieldName_;
+        
         //- Name of thermal conductivity field
-        word KName_;
+        const word KName_;
 
         bool fixesValue_;
 
+
+    // Private Member Functions
+
+        //- Am I or neighbour owner of interface
+        bool interfaceOwner(const polyMesh& nbrRegion) const;
+
+
 public:
 
     //- Runtime type information
@@ -162,12 +177,6 @@ public:
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
 
-        //- Evaluate the patch field
-        virtual void evaluate
-        (
-            const Pstream::commsTypes commsType=Pstream::blocking
-        );
-
         //- Write
         virtual void write(Ostream&) const;
 };

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C

@@ -1,156 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "solidWallTemperatureCoupledFvPatchScalarField.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::solidWallTemperatureCoupledFvPatchScalarField::
-solidWallTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(p, iF),
-    coupleManager_(p),
-    KName_("undefined-K")
-{}
-
-
-Foam::solidWallTemperatureCoupledFvPatchScalarField::
-solidWallTemperatureCoupledFvPatchScalarField
-(
-    const solidWallTemperatureCoupledFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
-    coupleManager_(ptf.coupleManager_),
-    KName_(ptf.KName_)
-{}
-
-
-Foam::solidWallTemperatureCoupledFvPatchScalarField::
-solidWallTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedValueFvPatchScalarField(p, iF),
-    coupleManager_(p, dict),
-    KName_(dict.lookup("K"))
-{
-    if (dict.found("value"))
-    {
-        fvPatchField<scalar>::operator=
-        (
-            scalarField("value", dict, p.size())
-        );
-    }
-    else
-    {
-        evaluate();
-    }
-}
-
-
-Foam::solidWallTemperatureCoupledFvPatchScalarField::
-solidWallTemperatureCoupledFvPatchScalarField
-(
-    const solidWallTemperatureCoupledFvPatchScalarField& wtcsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(wtcsf, iF),
-    coupleManager_(wtcsf.coupleManager_),
-    KName_(wtcsf.KName_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::solidWallTemperatureCoupledFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const fvPatchField<scalar>& neighbourField =
-        coupleManager_.neighbourPatchField<scalar>();
-
-    operator==(neighbourField);
-
-    fixedValueFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::solidWallTemperatureCoupledFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    fvPatchScalarField::write(os);
-    coupleManager_.writeEntries(os);
-    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
-    writeEntry("value", os);
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::solidWallTemperatureCoupledFvPatchScalarField::flux() const
-{
-    const fvPatchScalarField& Kw =
-        patch().lookupPatchField<volScalarField, scalar>(KName_);
-
-    const fvPatchScalarField& Tw = *this;
-
-    return Tw.snGrad()*Kw;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-makePatchTypeField
-(
-    fvPatchScalarField,
-    solidWallTemperatureCoupledFvPatchScalarField
-);
-
-} // End namespace Foam
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.H

@@ -1,160 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    solidWallHeatFluxCoupledFvPatchScalarField
-
-Description
-    Fixed value boundary condition for temperature, to be used by the
-    conjugate heat transfer solver.
-
-    Example usage:
-        myInterfacePatchName
-        {
-            type                solidWallTemperatureCoupled;
-            neighbourRegionName fluid;
-            neighbourPatchName  fluidSolidInterface;
-            neighbourFieldName  T;
-            K                   K;
-            value               uniform 300;
-        }
-
-SourceFiles
-    solidWallTemperatureCoupledFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef solidWallTemperatureCoupledFvPatchScalarField_H
-#define solidWallTemperatureCoupledFvPatchScalarField_H
-
-#include "fvPatchFields.H"
-#include "fixedValueFvPatchFields.H"
-#include "coupleManager.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-        Class solidWallTemperatureCoupledFvPatchScalarField Declaration
-\*---------------------------------------------------------------------------*/
-
-class solidWallTemperatureCoupledFvPatchScalarField
-:
-    public fixedValueFvPatchScalarField
-{
-    // Private data
-
-        //- Couple manager object
-        coupleManager coupleManager_;
-
-        //- Name of thermal conductivity field
-        word KName_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("solidWallTemperatureCoupled");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        solidWallTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        solidWallTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given solidWallTemperatureCoupledFvPatchScalarField
-        //  onto a new patch
-        solidWallTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallTemperatureCoupledFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallTemperatureCoupledFvPatchScalarField(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        solidWallTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallTemperatureCoupledFvPatchScalarField&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallTemperatureCoupledFvPatchScalarField(*this, iF)
-            );
-        }
-
-
-    // Member functions
-
-        //- (intensive) flux
-        tmp<scalarField> flux() const;
-
-        //- Update the coefficients associated with the patch field
-        virtual void updateCoeffs();
-
-        //- Write
-        virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H

@@ -14,12 +14,10 @@
          (
             UEqn()
          ==
-           -fvc::reconstruct
+            fvc::reconstruct
             (
-                (
-                    fvc::snGrad(pd)
-                  + ghf*fvc::snGrad(rho)
-                ) * mesh.magSf()
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              - fvc::snGrad(p)*mesh.magSf()
             )
         );
     }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H

@@ -6,43 +6,15 @@
     PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
     PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
     PtrList<volScalarField> DpDtFluid(fluidRegions.size());
-    PtrList<volScalarField> ghFluid(fluidRegions.size());
-    PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
-    PtrList<volScalarField> pdFluid(fluidRegions.size());
 
     List<scalar> initialMassFluid(fluidRegions.size());
 
-    dimensionedScalar pRef
-    (
-        "pRef",
-        dimensionSet(1, -1, -2, 0, 0),
-        rp.lookup("pRef")
-    );
-
     // Populate fluid field pointer lists
     forAll(fluidRegions, i)
     {
         Info<< "*** Reading fluid mesh thermophysical properties for region "
             << fluidRegions[i].name() << nl << endl;
 
-        Info<< "    Adding to pdFluid\n" << endl;
-        pdFluid.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "pd",
-                    runTime.timeName(),
-                    fluidRegions[i],
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                fluidRegions[i]
-            )
-        );
-
         Info<< "    Adding to thermoFluid\n" << endl;
         thermoFluid.set
         (
@@ -145,6 +117,7 @@
             i,
             new volScalarField
             (
+                "DpDt",
                 fvc::DDt
                 (
                     surfaceScalarField
@@ -157,36 +130,5 @@
             )
         );
 
-        const dictionary& environmentalProperties =
-            fluidRegions[i].lookupObject<IOdictionary>
-                ("environmentalProperties");
-        dimensionedVector g(environmentalProperties.lookup("g"));
-
-        Info<< "    Adding to ghFluid\n" << endl;
-        ghFluid.set
-        (
-            i,
-            new volScalarField
-            (
-                "gh",
-                 g & fluidRegions[i].C()
-            )
-        );
-        ghfFluid.set
-        (
-            i,
-            new surfaceScalarField
-            (
-                "ghf",
-                 g & fluidRegions[i].Cf()
-            )
-        );
-
-        Info<< "    Updating p from pd\n" << endl;
-        thermoFluid[i].p() == pdFluid[i] + rhoFluid[i]*ghFluid[i] + pRef;
-        thermoFluid[i].correct();
-
         initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
     }
-
-

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H

@@ -1,5 +1,5 @@
 {
-    tmp<fvScalarMatrix> hEqn
+    fvScalarMatrix hEqn
     (
         fvm::ddt(rho, h)
       + fvm::div(phi, h)
@@ -7,8 +7,16 @@
      ==
         DpDt
     );
-    hEqn().relax();
-    hEqn().solve();
+    if (oCorr == nOuterCorr-1)
+    {
+        hEqn.relax();
+        hEqn.solve(mesh.solver("hFinal"));
+    }
+    else
+    {
+        hEqn.relax();
+        hEqn.solve();
+    }
 
     thermo.correct();
 

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

@@ -1,5 +1,5 @@
 {
-    bool closedVolume = pd.needReference();
+    bool closedVolume = p.needReference();
 
     rho = thermo.rho();
 
@@ -17,31 +17,34 @@
         )
     );
 
-    phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
+    phi = phiU + fvc::interpolate(rho)*(g & mesh.Sf())*rhorUAf;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, pd)
-          + fvc::ddt(psi)*pRef
-          + fvc::ddt(psi, rho)*gh
+            fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, pd)
+          - fvm::laplacian(rhorUAf, p)
         );
 
-        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+        if
+        (
+            oCorr == nOuterCorr-1
+         && corr == nCorr-1
+         && nonOrth == nNonOrthCorr
+        )
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi += pdEqn.flux();
+            phi += pEqn.flux();
         }
     }
 
@@ -49,27 +52,24 @@
     U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
     U.correctBoundaryConditions();
 
-    // Update pressure field (including bc)
-    p == pd + rho*gh + pRef;
+    // Update pressure substantive derivative
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
     // Solve continuity
-#   include "rhoEqn.H"
+    #include "rhoEqn.H"
 
     // Update continuity errors
-#   include "compressibleContinuityErrors.H"
+    #include "compressibleContinuityErrors.H"
 
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
     {
-        p += (massIni - fvc::domainIntegrate(psi*p))/fvc::domainIntegrate(psi);
+        p += (massIni - fvc::domainIntegrate(psi*p))
+            /fvc::domainIntegrate(psi);
         rho = thermo.rho();
     }
 
     // Update thermal conductivity
     K = thermoFluid[i].Cp()*turb.alphaEff();
-
-    // Update pd (including bc)
-    pd == p - (rho*gh + pRef);
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H

@@ -0,0 +1,9 @@
+    const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
+
+    int nCorr(readInt(pimple.lookup("nCorrectors")));
+
+    int nNonOrthCorr =
+        pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
+
+    bool momentumPredictor =
+        pimple.lookupOrDefault<Switch>("momentumPredictor", true);

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H

@@ -1 +0,0 @@
-    solve(fvm::ddt(rho) + fvc::div(phi));

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -7,12 +7,15 @@
     surfaceScalarField& phi = phiFluid[i];
     compressible::turbulenceModel& turb = turbulence[i];
     volScalarField& DpDt = DpDtFluid[i];
-    const volScalarField& gh = ghFluid[i];
-    const surfaceScalarField& ghf = ghfFluid[i];
-    volScalarField& pd = pdFluid[i];
 
     volScalarField& p = thermo.p();
     const volScalarField& psi = thermo.psi();
     volScalarField& h = thermo.h();
 
     const dimensionedScalar massIni("massIni", dimMass, initialMassFluid[i]);
+
+    const dictionary& environmentalProperties =
+        fluidRegions[i].lookupObject<IOdictionary>
+        ("environmentalProperties");
+
+    const dimensionedVector g(environmentalProperties.lookup("g"));

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H

@@ -1,15 +1,18 @@
-#   include "rhoEqn.H"
-    for (int ocorr=0; ocorr<nOuterCorr; ocorr++)
-    {
-    #   include "UEqn.H"
+if (oCorr == 0)
+{
+    #include "rhoEqn.H"
+}
 
-    #   include "hEqn.H"
+#include "UEqn.H"
 
-     // --- PISO loop
+#include "hEqn.H"
 
-        for (int corr=0; corr<nCorr; corr++)
-        {
-    #       include "pEqn.H"
-        }
-    }
-    turb.correct();
+// --- PISO loop
+for (int corr=0; corr<nCorr; corr++)
+{
+    #include "pEqn.H"
+}
+
+turb.correct();
+
+rho = thermo.rho();

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/storeOldFluidFields.H

@@ -0,0 +1,2 @@
+    p.storePrevIter();
+    rho.storePrevIter();

applications/solvers/heatTransfer/chtMultiRegionFoam/readPIMPLEControls.H

@@ -0,0 +1,7 @@
+    // We do not have a top-level mesh. Construct the fvSolution for
+    // the runTime instead.
+    fvSolution solutionDict(runTime);
+
+    const dictionary& pimple = solutionDict.subDict("PIMPLE");
+
+    int nOuterCorr(readInt(pimple.lookup("nOuterCorrectors")));

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPIMPLEControls.H

@@ -0,0 +1,4 @@
+    const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
+
+    int nNonOrthCorr =
+        pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H

@@ -1,4 +1,4 @@
-    dictionary piso = solidRegions[i].solutionDict().subDict("PISO");
+    const dictionary& piso = solidRegions[i].solutionDict().subDict("PISO");
 
     int nNonOrthCorr = 0;
     if (piso.found("nNonOrthogonalCorrectors"))

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H

@@ -1,4 +1,4 @@
-//    fvMesh& mesh = solidRegions[i];
+    fvMesh& mesh = solidRegions[i];
 
     volScalarField& rho = rhos[i];
     volScalarField& cp = cps[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H

@@ -1,10 +1,13 @@
 {
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        solve
+        tmp<fvScalarMatrix> TEqn
         (
-            fvm::ddt(rho*cp, T) - fvm::laplacian(K, T)
+            fvm::ddt(rho*cp, T)
+          - fvm::laplacian(K, T)
         );
+        TEqn().relax();
+        TEqn().solve();
     }
 
     Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl;

applications/solvers/lagrangian/reactingParcelFoam/Make/files

@@ -0,0 +1,3 @@
+reactingParcelFoam.C
+
+EXE = $(FOAM_APPBIN)/reactingParcelFoam

applications/solvers/lagrangian/reactingParcelFoam/Make/options

@@ -0,0 +1,37 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/pdfs/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lspecie \
+    -lbasicThermophysicalModels \
+    -lliquids \
+    -lliquidMixture \
+    -lsolids \
+    -lsolidMixture \
+    -lthermophysicalFunctions \
+    -lcombustionThermophysicalModels \
+    -lchemistryModel \
+    -lradiation \
+    -lODE

applications/solvers/lagrangian/reactingParcelFoam/UEqn.H

@@ -0,0 +1,16 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho.dimensionedInternalField()*g
+      + parcels.SU()
+    );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }

applications/solvers/lagrangian/reactingParcelFoam/YEqn.H

@@ -0,0 +1,43 @@
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+
+{
+    label inertIndex = -1;
+    volScalarField Yt = 0.0*Y[0];
+
+    for (label i=0; i<Y.size(); i++)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+            solve
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+              ==
+                parcels.Srho(i)
+              + kappa*chemistry.RR(i)().dimensionedInternalField()
+            );
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}

applications/solvers/lagrangian/reactingParcelFoam/additionalOutput.H

@@ -0,0 +1,48 @@
+{
+    tmp<volScalarField> tdQ
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "dQ",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar
+            (
+                "zero",
+                dimensionSet(1, -3, -1, 0, 0, 0, 0),
+                0.0
+            )
+        )
+    );
+
+    scalarField& dQ = tdQ();
+
+    scalarField cp(dQ.size(), 0.0);
+
+    forAll(Y, i)
+    {
+        volScalarField RRi = chemistry.RR(i);
+
+        forAll(h, celli)
+        {
+            scalar Ti = T[celli];
+            cp[celli] += Y[i][celli]*chemistry.specieThermo()[i].Cp(Ti);
+            scalar hi = chemistry.specieThermo()[i].h(Ti);
+            scalar RR = RRi[celli];
+            dQ[celli] -= hi*RR;
+        }
+    }
+
+    forAll(dQ, celli)
+    {
+        dQ[celli] /= cp[celli];
+    }
+
+    tdQ().write();
+}

applications/solvers/lagrangian/reactingParcelFoam/chemistry.H

@@ -0,0 +1,25 @@
+{
+    Info << "Solving chemistry" << endl;
+
+    chemistry.solve
+    (
+        runTime.value() - runTime.deltaT().value(),
+        runTime.deltaT().value()
+    );
+
+    // turbulent time scale
+    if (turbulentReaction)
+    {
+        DimensionedField<scalar, volMesh> tk =
+            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+        DimensionedField<scalar, volMesh> tc =
+            chemistry.tc()().dimensionedInternalField();
+
+        // Chalmers PaSR model
+        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
+    }
+    else
+    {
+        kappa = 1.0;
+    }
+}

applications/solvers/lagrangian/reactingParcelFoam/createClouds.H

@@ -0,0 +1,9 @@
+Info<< "\nConstructing reacting cloud" << endl;
+BasicReactingCloud<specieReactingProperties> parcels
+(
+    "reactingCloud1",
+    rho,
+    U,
+    g,
+    thermo()
+);

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -0,0 +1,90 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<hCombustionThermo> thermo
+    (
+        hCombustionThermo::New(mesh)
+    );
+
+    combustionMixture& composition = thermo->composition();
+    PtrList<volScalarField>& Y = composition.Y();
+
+    word inertSpecie(thermo->lookup("inertSpecie"));
+
+    volScalarField& p = thermo->p();
+    volScalarField& h = thermo->h();
+    const volScalarField& T = thermo->T();
+    const volScalarField& psi = thermo->psi();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo->rho()
+    );
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "compressibleCreatePhi.H"
+
+    DimensionedField<scalar, volMesh> kappa
+    (
+        IOobject
+        (
+            "kappa",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimless, 0.0)
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo()
+        )
+    );
+
+    Info<< "Creating field DpDt\n" << endl;
+    volScalarField DpDt =
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+    Info << "Constructing chemical mechanism" << endl;
+    chemistryModel chemistry
+    (
+        thermo(),
+        rho
+    );
+
+    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+    forAll (Y, i)
+    {
+        fields.add(Y[i]);
+    }
+    fields.add(h);

applications/solvers/lagrangian/reactingParcelFoam/hEqn.H

@@ -0,0 +1,20 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::ddt(rho, h)
+      + fvm::div(phi, h)
+      - fvm::laplacian(turbulence->alphaEff(), h)
+     ==
+        DpDt
+     +  parcels.Sh()
+     +  radiation->Sh(thermo())
+    );
+
+    hEqn.relax();
+
+    hEqn.solve();
+
+    thermo->correct();
+
+    radiation->correct();
+}

applications/solvers/lagrangian/reactingParcelFoam/pEqn.H

@@ -0,0 +1,72 @@
+rho = thermo->rho();
+
+volScalarField rUA = 1.0/UEqn.A();
+U = rUA*UEqn.H();
+
+if (transonic)
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(thermo->psi())
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho*rUA, p)
+         ==
+            parcels.Srho()
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    phi =
+        fvc::interpolate(rho)*
+        (
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rUA, p)
+         ==
+            parcels.Srho()
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    reactingParcelFoam
+
+Description
+    Transient PISO solver for compressible, laminar or turbulent flow with
+    reacting Lagrangian parcels.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "hCombustionThermo.H"
+#include "turbulenceModel.H"
+#include "BasicReactingCloud.H"
+#include "chemistryModel.H"
+#include "chemistrySolver.H"
+#include "reactingThermoTypes.H"
+#include "radiationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readChemistryProperties.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "createClouds.H"
+    #include "createRadiationModel.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        parcels.evolve();
+
+        parcels.info();
+
+        #include "chemistry.H"
+        #include "rhoEqn.H"
+
+        // --- PIMPLE loop
+        for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
+        {
+            #include "UEqn.H"
+            #include "YEqn.H"
+
+            // --- PISO loop
+            for (int corr=1; corr<=nCorr; corr++)
+            {
+                #include "hEqn.H"
+                #include "pEqn.H"
+            }
+
+            Info<< "T gas min/max   = " << min(T).value() << ", "
+                << max(T).value() << endl;
+        }
+
+        turbulence->correct();
+
+        rho = thermo->rho();
+
+        if (runTime.write())
+        {
+            #include "additionalOutput.H"
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/reactingParcelFoam/readChemistryProperties.H

@@ -0,0 +1,22 @@
+Info<< "Reading chemistry properties\n" << endl;
+
+IOdictionary chemistryProperties
+(
+    IOobject
+    (
+        "chemistryProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction"));
+
+dimensionedScalar Cmix("Cmix", dimless, 1.0);
+
+if (turbulentReaction)
+{
+    chemistryProperties.lookup("Cmix") >> Cmix;
+}

applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2007 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,17 +22,22 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
 \*---------------------------------------------------------------------------*/
 
-#include "upwindCFCStencilObject.H"
-
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
 {
-    defineTypeNameAndDebug(upwindCFCStencilObject, 0);
+    solve
+    (
+        fvm::ddt(rho)
+      + fvc::div(phi)
+      ==
+        parcels.Srho()
+    );
 }
 
-
 // ************************************************************************* //

applications/solvers/lagrangian/trackedReactingParcelFoam/Make/files

@@ -0,0 +1,3 @@
+trackedReactingParcelFoam.C
+
+EXE = $(FOAM_APPBIN)/trackedReactingParcelFoam

applications/solvers/lagrangian/trackedReactingParcelFoam/Make/options

@@ -0,0 +1,38 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/pdfs/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lspecie \
+    -lbasicThermophysicalModels \
+    -lliquids \
+    -lliquidMixture \
+    -lsolids \
+    -lsolidMixture \
+    -lthermophysicalFunctions \
+    -lcombustionThermophysicalModels \
+    -lchemistryModel \
+    -lradiation \
+    -lODE

applications/solvers/lagrangian/trackedReactingParcelFoam/UEqn.H

@@ -0,0 +1,47 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho.dimensionedInternalField()*g
+      + reactingParcels.SU()
+    );
+
+    UEqn.relax();
+
+    tmp<volScalarField> trAU;
+    tmp<volTensorField> trTU;
+
+    if (pressureImplicitPorosity)
+    {
+        tmp<volTensorField> tTU = tensor(I)*UEqn.A();
+        pZones.addResistance(UEqn, tTU());
+        trTU = inv(tTU());
+        trTU().rename("rAU");
+
+        volVectorField gradp = fvc::grad(p);
+
+        for (int UCorr=0; UCorr<nUCorr; UCorr++)
+        {
+            U = trTU() & (UEqn.H() - gradp);
+        }
+        U.correctBoundaryConditions();
+    }
+    else
+    {
+        pZones.addResistance(UEqn);
+
+        solve
+        (
+            UEqn == -fvc::grad(p)
+        );
+
+        trAU = 1.0/UEqn.A();
+        trAU().rename("rAU");
+    }
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }

applications/solvers/lagrangian/trackedReactingParcelFoam/YEqn.H

@@ -0,0 +1,46 @@
+
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+
+{
+    label inertIndex = -1;
+    volScalarField Yt = 0.0*Y[0];
+
+    for (label i=0; i<Y.size(); i++)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+            solve
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+              ==
+                reactingParcels.Srho(i)
+              + kappa*chemistry.RR(i)().dimensionedInternalField()
+              + pointMassSources.Su(i)
+            );
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+
+}

applications/solvers/lagrangian/trackedReactingParcelFoam/additionalOutput.H

@@ -0,0 +1,48 @@
+{
+    tmp<volScalarField> tdQ
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "dQ",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar
+            (
+                "zero",
+                dimensionSet(1, -3, -1, 0, 0, 0, 0),
+                0.0
+            )
+        )
+    );
+
+    scalarField& dQ = tdQ();
+
+    scalarField cp(dQ.size(), 0.0);
+
+    forAll(Y, i)
+    {
+        volScalarField RRi = chemistry.RR(i);
+
+        forAll(h, celli)
+        {
+            scalar Ti = T[celli];
+            cp[celli] += Y[i][celli]*chemistry.specieThermo()[i].Cp(Ti);
+            scalar hi = chemistry.specieThermo()[i].h(Ti);
+            scalar RR = RRi[celli];
+            dQ[celli] -= hi*RR;
+        }
+    }
+
+    forAll(dQ, celli)
+    {
+        dQ[celli] /= cp[celli];
+    }
+
+    tdQ().write();
+}

applications/solvers/lagrangian/trackedReactingParcelFoam/chemistry.H

@@ -0,0 +1,25 @@
+{
+    Info << "Solving chemistry" << endl;
+
+    chemistry.solve
+    (
+        runTime.value() - runTime.deltaT().value(),
+        runTime.deltaT().value()
+    );
+
+    // turbulent time scale
+    if (turbulentReaction)
+    {
+        DimensionedField<scalar, volMesh> tk =
+            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+        DimensionedField<scalar, volMesh> tc =
+            chemistry.tc()().dimensionedInternalField();
+
+        // Chalmers PaSR model
+        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
+    }
+    else
+    {
+        kappa = 1.0;
+    }
+}

applications/solvers/lagrangian/trackedReactingParcelFoam/createClouds.H

@@ -0,0 +1,9 @@
+Info<< "\nConstructing reacting cloud" << endl;
+BasicTrackedReactingCloud<specieReactingProperties> reactingParcels
+(
+    "reactingCloud1",
+    rho,
+    U,
+    g,
+    thermo()
+);

applications/solvers/lagrangian/trackedReactingParcelFoam/createFields.H

@@ -0,0 +1,111 @@
+    Info<< "Reading thermophysical properties" << nl << endl;
+
+    autoPtr<hCombustionThermo> thermo
+    (
+        hCombustionThermo::New(mesh)
+    );
+
+    combustionMixture& composition = thermo->composition();
+    PtrList<volScalarField>& Y = composition.Y();
+
+    word inertSpecie(thermo->lookup("inertSpecie"));
+
+    volScalarField& p = thermo->p();
+    volScalarField& h = thermo->h();
+    const volScalarField& T = thermo->T();
+    const volScalarField& psi = thermo->psi();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo->rho()
+    );
+
+    Info<< "Reading field U" << nl << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "compressibleCreatePhi.H"
+
+    DimensionedField<scalar, volMesh> kappa
+    (
+        IOobject
+        (
+            "kappa",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimless, 0.0)
+    );
+
+    Info<< "Creating turbulence model" << nl << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo()
+        )
+    );
+
+    Info<< "Creating field DpDt" << nl << endl;
+    volScalarField DpDt =
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+    Info << "Constructing chemical mechanism" << nl << endl;
+    chemistryModel chemistry
+    (
+        thermo(),
+        rho
+    );
+
+    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+    forAll (Y, i)
+    {
+        fields.add(Y[i]);
+    }
+    fields.add(h);
+
+    Info<< "Creating porous zones" << nl << endl;
+    porousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);
+
+    label nUCorr = 0;
+    if (pZones.size())
+    {
+        // nUCorrectors for pressureImplicitPorosity
+        if (mesh.solutionDict().subDict("PISO").found("nUCorrectors"))
+        {
+            mesh.solutionDict().subDict("PISO").lookup("nUCorrectors")
+                >> nUCorr;
+        }
+
+        if (nUCorr > 0)
+        {
+            pressureImplicitPorosity = true;
+        }
+    }
+

applications/solvers/lagrangian/trackedReactingParcelFoam/createMulticomponentPointSources.H

@@ -0,0 +1,8 @@
+Info<< "\nConstructing multi-component mass flow rate point sources" << endl;
+timeActivatedExplicitMulticomponentPointSource pointMassSources
+(
+    "pointMassSources",
+    mesh,
+    Y,
+    dimMass/dimVolume/dimTime
+);

applications/solvers/lagrangian/trackedReactingParcelFoam/hEqn.H

@@ -0,0 +1,20 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::ddt(rho, h)
+      + fvm::div(phi, h)
+      - fvm::laplacian(turbulence->alphaEff(), h)
+     ==
+        DpDt
+     +  reactingParcels.Sh()
+     +  radiation->Sh(thermo())
+    );
+
+    hEqn.relax();
+
+    hEqn.solve();
+
+    thermo->correct();
+
+    radiation->correct();
+}

applications/solvers/lagrangian/trackedReactingParcelFoam/pEqn.H

@@ -0,0 +1,110 @@
+rho = thermo->rho();
+
+if (pressureImplicitPorosity)
+{
+    U = trTU()&UEqn.H();
+}
+else
+{
+    U = trAU()*UEqn.H();
+}
+
+if (transonic)
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(thermo->psi())
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+//          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        tmp<fvScalarMatrix> lapTerm;
+
+        if (pressureImplicitPorosity)
+        {
+            lapTerm = fvm::laplacian(rho*trTU(), p);
+        }
+        else
+        {
+            lapTerm = fvm::laplacian(rho*trAU(), p);
+        }
+
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - lapTerm()
+         ==
+            reactingParcels.Srho()
+          + pointMassSources.Su()
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    phi =
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(trAU(), rho, U, phi)
+        );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        tmp<fvScalarMatrix> lapTerm;
+
+        if (pressureImplicitPorosity)
+        {
+            lapTerm = fvm::laplacian(rho*trTU(), p);
+        }
+        else
+        {
+            lapTerm = fvm::laplacian(rho*trAU(), p);
+        }
+
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phi)
+          - lapTerm()
+         ==
+            reactingParcels.Srho()
+          + pointMassSources.Su()
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+if (pressureImplicitPorosity)
+{
+    U -= trTU()&fvc::grad(p);
+}
+else
+{
+    U -= trAU()*fvc::grad(p);
+}
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+

applications/solvers/lagrangian/trackedReactingParcelFoam/readChemistryProperties.H

@@ -0,0 +1,22 @@
+Info<< "Reading chemistry properties\n" << endl;
+
+IOdictionary chemistryProperties
+(
+    IOobject
+    (
+        "chemistryProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction"));
+
+dimensionedScalar Cmix("Cmix", dimless, 1.0);
+
+if (turbulentReaction)
+{
+    chemistryProperties.lookup("Cmix") >> Cmix;
+}

applications/solvers/lagrangian/trackedReactingParcelFoam/rhoEqn.H

@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    solve
+    (
+        fvm::ddt(rho)
+      + fvc::div(phi)
+      ==
+        reactingParcels.Srho()
+      + pointMassSources.Su()
+    );
+}
+
+// ************************************************************************* //

applications/solvers/lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "hCombustionThermo.H"
+#include "turbulenceModel.H"
+#include "BasicTrackedReactingCloud.H"
+#include "chemistryModel.H"
+#include "chemistrySolver.H"
+#include "reactingThermoTypes.H"
+#include "radiationModel.H"
+#include "porousZones.H"
+#include "timeActivatedExplicitMulticomponentPointSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readChemistryProperties.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "createRadiationModel.H"
+    #include "createClouds.H"
+    #include "createMulticomponentPointSources.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        reactingParcels.evolve();
+
+        reactingParcels.info();
+
+        #include "chemistry.H"
+        #include "rhoEqn.H"
+
+        // --- PIMPLE loop
+        for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
+        {
+            #include "UEqn.H"
+            #include "YEqn.H"
+
+            // --- PISO loop
+            for (int corr=1; corr<=nCorr; corr++)
+            {
+                #include "hEqn.H"
+                #include "pEqn.H"
+            }
+
+            Info<< "T gas min/max   = " << min(T).value() << ", "
+                << max(T).value() << endl;
+        }
+
+        turbulence->correct();
+
+        rho = thermo->rho();
+
+        if (runTime.write())
+        {
+            #include "additionalOutput.H"
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/Make/files

@@ -0,0 +1,3 @@
+uncoupledKinematicParcelFoam.C
+
+EXE = $(FOAM_APPBIN)/uncoupledKinematicParcelFoam

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/Make/options

@@ -6,7 +6,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
 
 EXE_LIBS = \
     -llagrangian \

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/createFields.H

@@ -47,11 +47,7 @@
     );
 
     word kinematicCloudName("kinematicCloud");
-
-    if (args.options().found("cloudName"))
-    {
-        kinematicCloudName = args.options()["cloudName"];
-    }
+    args.optionReadIfPresent("cloudName", kinematicCloudName);
 
     Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
     basicKinematicCloud kinematicCloud

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H

@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
+                  - fvc::snGrad(p)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H

@@ -12,7 +12,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 

applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -88,24 +88,6 @@
 
     dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
 
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        )
-    );
-
     volScalarField rho1 = rho10 + psi1*p;
     volScalarField rho2 = rho20 + psi2*p;
 
@@ -152,11 +134,11 @@
     );
 
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }

applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H

@@ -2,17 +2,17 @@
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rUAf = fvc::interpolate(rUA);
 
-    tmp<fvScalarMatrix> pdEqnComp;
+    tmp<fvScalarMatrix> pEqnComp;
 
     if (transonic)
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
     }
     else
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
     }
 
 
@@ -26,16 +26,16 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqnIncomp
+        fvScalarMatrix pEqnIncomp
         (
             fvc::div(phi)
-          - fvm::laplacian(rUAf, pd)
+          - fvm::laplacian(rUAf, p)
         );
 
         if
@@ -51,9 +51,9 @@
                     max(alpha1, scalar(0))*(psi1/rho1)
                   + max(alpha2, scalar(0))*(psi2/rho2)
                 )
-               *pdEqnComp()
-              + pdEqnIncomp,
-                mesh.solver(pd.name() + "Final")
+               *pEqnComp()
+              + pEqnIncomp,
+                mesh.solver(p.name() + "Final")
             );
         }
         else
@@ -64,8 +64,8 @@
                     max(alpha1, scalar(0))*(psi1/rho1)
                   + max(alpha2, scalar(0))*(psi2/rho2)
                 )
-               *pdEqnComp()
-              + pdEqnIncomp
+               *pEqnComp()
+              + pEqnIncomp
             );
         }
 
@@ -73,26 +73,21 @@
         {
             dgdt =
                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(pdEqnComp & pd);
-            phi += pdEqnIncomp.flux();
+               *(pEqnComp & p);
+            phi += pEqnIncomp.flux();
         }
     }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 
-    p = max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        );
+    p.max(pMin);
 
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(pd) " << min(pd).value() << endl;
+    Info<< "min(p) " << min(p).value() << endl;
 
     // Make the fluxes relative to the mesh motion
     fvc::makeRelative(phi, U);

applications/solvers/multiphase/compressibleInterFoam/UEqn.H

@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
+                  - fvc::snGrad(p)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -46,11 +46,6 @@
     #include "createPhi.H"
 
 
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi);
 
@@ -88,24 +83,6 @@
 
     dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
 
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        )
-    );
-
     volScalarField rho1 = rho10 + psi1*p;
     volScalarField rho2 = rho20 + psi2*p;
 

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -2,17 +2,17 @@
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rUAf = fvc::interpolate(rUA);
 
-    tmp<fvScalarMatrix> pdEqnComp;
+    tmp<fvScalarMatrix> pEqnComp;
 
     if (transonic)
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
     }
     else
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
     }
 
 
@@ -26,16 +26,16 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf;
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqnIncomp
+        fvScalarMatrix pEqnIncomp
         (
             fvc::div(phi)
-          - fvm::laplacian(rUAf, pd)
+          - fvm::laplacian(rUAf, p)
         );
 
         solve
@@ -44,31 +44,27 @@
                 max(alpha1, scalar(0))*(psi1/rho1)
               + max(alpha2, scalar(0))*(psi2/rho2)
             )
-           *pdEqnComp()
-          + pdEqnIncomp
+           *pEqnComp()
+          + pEqnIncomp
         );
 
         if (nonOrth == nNonOrthCorr)
         {
             dgdt =
                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(pdEqnComp & pd);
-            phi += pdEqnIncomp.flux();
+               *(pEqnComp & p);
+            phi += pEqnIncomp.flux();
         }
     }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 
-    p = max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        );
+    p.max(pMin);
 
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(pd) " << min(pd).value() << endl;
+    Info<< "min(p) " << min(p).value() << endl;
 }

applications/solvers/multiphase/interDyMFoam/correctPhi.H

@@ -12,7 +12,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 
@@ -27,7 +27,7 @@
             fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
         );
 
-        pcorrEqn.setReference(pdRefCell, pdRefValue);
+        pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
         if (nonOrth == nNonOrthCorr)

applications/solvers/multiphase/interDyMFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -92,47 +92,17 @@
         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }
     }
 
 
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rho*(g & mesh.C())
-    );
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
-
+    label pRefCell = 0;
     scalar pRefValue = 0.0;
-
-    if (pd.needReference())
-    {
-        pRefValue = readScalar
-        (
-            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
-        );
-
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pdRefCell)
-        );
-    }
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);

applications/solvers/multiphase/interDyMFoam/interDyMFoam.C

@@ -114,18 +114,6 @@ int main(int argc, char *argv[])
             #include "pEqn.H"
         }
 
-        p = pd + rho*gh;
-
-        if (pd.needReference())
-        {
-            p += dimensionedScalar
-            (
-                "p",
-                p.dimensions(),
-                pRefValue - getRefCellValue(p, pdRefCell)
-            );
-        }
-
         turbulence->correct();
 
         runTime.write();

applications/solvers/multiphase/interDyMFoam/pEqn.H

@@ -7,38 +7,38 @@
 
     phi = phiU +
     (
-        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-      - ghf*fvc::snGrad(rho)
-    )*rAUf*mesh.magSf();
+        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+      + fvc::interpolate(rho)*(g & mesh.Sf())
+    )*rAUf;
 
-    if (pd.needReference())
+    if (p.needReference())
     {
         fvc::makeRelative(phi, U);
-        adjustPhi(phi, U, pd);
+        adjustPhi(phi, U, p);
         fvc::makeAbsolute(phi, U);
     }
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::laplacian(rAUf, pd) == fvc::div(phi)
+            fvm::laplacian(rAUf, p) == fvc::div(phi)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, pRefValue);
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi -= pdEqn.flux();
+            phi -= pEqn.flux();
         }
     }
 

applications/solvers/multiphase/interFoam/UEqn.H

@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
+                  - fvc::snGrad(p)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/interFoam/correctPhi.H

@@ -1,11 +1,11 @@
 {
 #   include "continuityErrs.H"
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }
@@ -22,7 +22,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 
@@ -37,7 +37,7 @@
             fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
         );
 
-        pcorrEqn.setReference(pdRefCell, pdRefValue);
+        pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
         if (nonOrth == nNonOrthCorr)

applications/solvers/multiphase/interFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -83,45 +83,9 @@
     );
 
 
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rho*gh
-    );
-
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
-
+    label pRefCell = 0;
     scalar pRefValue = 0.0;
-
-    if (pd.needReference())
-    {
-        pRefValue = readScalar
-        (
-            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
-        );
-
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pdRefCell)
-        );
-    }
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
 
 
     // Construct interface from alpha1 distribution

applications/solvers/multiphase/interFoam/interFoam.C

@@ -89,18 +89,6 @@ int main(int argc, char *argv[])
 
         #include "continuityErrs.H"
 
-        p = pd + rho*gh;
-
-        if (pd.needReference())
-        {
-            p += dimensionedScalar
-            (
-                "p",
-                p.dimensions(),
-                pRefValue - getRefCellValue(p, pdRefCell)
-            );
-        }
-
         turbulence->correct();
 
         runTime.write();

applications/solvers/multiphase/interFoam/pEqn.H

@@ -12,33 +12,33 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf;
 
-    adjustPhi(phi, U, pd);
+    adjustPhi(phi, U, p);
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::laplacian(rUAf, pd) == fvc::div(phi)
+            fvm::laplacian(rUAf, p) == fvc::div(phi)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, pRefValue);
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi -= pdEqn.flux();
+            phi -= pEqn.flux();
         }
     }