add dfBuoyancyFoam

Allwmake

@@ -24,5 +24,6 @@ wmake applications/solvers/df0DFoam
 wmake applications/solvers/dfLowMachFoam
 wmake applications/solvers/dfHighSpeedFoam
 wmake applications/solvers/dfSprayFoam
+wmake applications/solvers/dfBuoyancyFoam
 
 wmake applications/utilities/flameSpeed

applications/solvers/dfBuoyancyFoam/Make/files

@@ -0,0 +1,3 @@
+dfBuoyancyFoam.C
+
+EXE = $(DF_APPBIN)/dfBuoyancyFoam

applications/solvers/dfBuoyancyFoam/Make/options

@@ -0,0 +1,76 @@
+-include $(GENERAL_RULES)/mplibType
+
+EXE_INC = -std=c++14 \
+    -g \
+    -fopenmp \
+    -Wno-unused-variable \
+    -Wno-unused-but-set-variable \
+    -Wno-old-style-cast \
+    -I. \
+    $(PFLAGS) $(PINC) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionModels/pyrolysisModels/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(DF_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(DF_SRC)/lagrangian/spray/lnInclude \
+    -I$(LIB_SRC)/lagrangian/spray/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(DF_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+	-I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(DF_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/Pstream/mpi \
+    -I$(DF_SRC)/dfCanteraMixture/lnInclude \
+    -I$(DF_SRC)/dfChemistryModel/lnInclude \
+    -I$(DF_SRC)/dfCombustionModels/lnInclude \
+    -I$(CANTERA_ROOT)/include \
+    $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include,) \
+    $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include/torch/csrc/api/include,) \
+    $(PYTHON_INC_DIR) \
+    $(if $(AMGX_DIR), -I$(DF_ROOT)/src_gpu,) \
+    $(if $(AMGX_DIR), -I/usr/local/cuda/include,) \
+    $(if $(AMGX_DIR), -I$(AMGX_DIR)/include,) \
+    $(if $(ODE_GPU_SOLVER), -I$(OPENCC_PATH)/include,) \
+    $(if $(ODE_GPU_SOLVER), -DODE_GPU_SOLVER,)
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lturbulenceModels \
+    -lcompressibleTransportModels \
+    -llagrangian \
+    -lregionModels \
+    -L$(DF_LIBBIN) \
+    -ldfSurfaceFilmModels \
+    -ldfSLGThermo \
+    -ldfLagrangianIntermediate \
+    -ldfLagrangianTurbulence \
+    -ldfLagrangianSpray \
+    -ldfFluidThermophysicalModels \
+    -ldfCompressibleTurbulenceModels \
+    -ldfThermophysicalProperties \
+    -ldfCanteraMixture \
+    -ldfChemistryModel \
+    -ldfCombustionModels \
+    $(CANTERA_ROOT)/lib/libcantera.so \
+    $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libtorch.so,) \
+    $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libc10.so,) \
+    $(if $(LIBTORCH_ROOT),-rdynamic,) \
+    $(if $(LIBTORCH_ROOT),-lpthread,) \
+    $(if $(LIBTORCH_ROOT),$(DF_SRC)/dfChemistryModel/DNNInferencer/build/libDNNInferencer.so,) \
+    $(if $(PYTHON_LIB_DIR),$(PYTHON_LIB_DIR),) \
+    $(if $(AMGX_DIR), /usr/local/cuda/lib64/libcudart.so,) \
+	$(if $(AMGX_DIR), /usr/local/cuda/lib64/libnccl.so,) \
+    $(if $(AMGX_DIR), $(DF_ROOT)/src_gpu/build/libdfMatrix.so,) \
+    $(if $(AMGX_DIR), $(AMGX_DIR)/build/libamgxsh.so,) \
+    $(if $(ODE_GPU_SOLVER), $(ODE_GPU_SOLVER)/lib/libopencc.so,)

applications/solvers/dfBuoyancyFoam/UEqn.H

@@ -0,0 +1,34 @@
+    MRF.correctBoundaryVelocity(U);
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U) + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        parcels.SU(U)
+    //   + fvOptions(rho, U)
+    );
+
+    UEqn.relax();
+
+    // fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+          ==
+            fvc::reconstruct
+            (
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            )
+        );
+
+        // fvOptions.correct(U);
+        K = 0.5*magSqr(U);
+    }

applications/solvers/dfBuoyancyFoam/YEEqn.H

@@ -0,0 +1,90 @@
+tmp<fv::convectionScheme<scalar>> mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+{
+    // radiation->correct();
+    combustion->correct();
+
+    HRR = combustion->Qdot();
+
+    volScalarField Yt(0.0*Y[0]);
+
+    forAll(Y, i)
+    {
+        if (i != inertIndex)
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YiEqn
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->alphaEff(), Yi)
+              ==
+                parcels.SYi(i, Yi)
+            //   + surfaceFilm.Srho(i)
+              + combustion->R(Yi)
+            //   + fvOptions(rho, Yi)
+            );
+            YiEqn.relax();
+
+            // fvOptions.constrain(YiEqn);
+
+            YiEqn.solve("Yi");
+
+            // fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), U),
+                p,
+                "div(phiv,p)"
+            )
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        // combustion->Qdot() // ha Eqn has no reaction source
+    //   + radiation->Sh(thermo, he)
+         parcels.Sh(he)
+    //   + surfaceFilm.Sh()
+    //   + fvOptions(rho, he)
+    );
+
+    EEqn.relax();
+
+    // fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    // fvOptions.correct(he);
+
+    // thermo.correct();
+    chemistry->correctThermo();
+
+    Info<< "min/max(T) = "
+        << min(T).value() << ", " << max(T).value() << endl;
+}

applications/solvers/dfBuoyancyFoam/createClouds.H

@@ -0,0 +1,9 @@
+Info<< "\nConstructing reacting cloud" << endl;
+basicReactingCloud parcels
+(
+    "reactingCloud1",
+    rho,
+    U,
+    g,
+    slgThermo
+);

applications/solvers/dfBuoyancyFoam/createFieldRefs.H

@@ -0,0 +1,5 @@
+const volScalarField& psi = thermo.psi();
+const volScalarField& T = thermo.T();
+// regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
+const word inertSpecie(chemistry->lookup("inertSpecie"));
+const label inertIndex(chemistry->species()[inertSpecie]);

applications/solvers/dfBuoyancyFoam/createFields.H

@@ -0,0 +1,168 @@
+Info<< "Reading thermophysical properties\n" << endl;
+//autoPtr<psiReactionThermo> pThermo(psiReactionThermo::New(mesh));
+//psiReactionThermo& thermo = pThermo();
+//thermo.validate(args.executable(), "h", "e");
+
+CanteraMixture::setEnergyName("ha");
+// fluidThermo* pThermo = new heRhoThermo<rhoThermo, CanteraMixture>(mesh, word::null);
+fluidThermo*   pThermo = new heRhoThermo<rhoThermo, dfSingleStepReactingMixture>(mesh, word::null);
+fluidThermo& thermo = *pThermo;
+SLGThermo slgThermo(mesh, thermo);
+
+// basicSpecieMixture& composition = thermo.composition();
+// PtrList<volScalarField>& Y = composition.Y();
+
+// const word inertSpecie(thermo.lookup("inertSpecie"));
+// if (!composition.species().found(inertSpecie))
+// {
+//     FatalIOErrorIn(args.executable().c_str(), thermo)
+//         << "Inert specie " << inertSpecie << " not found in available species "
+//         << composition.species()
+//         << exit(FatalIOError);
+// }
+
+Info<< "Creating field rho\n" << endl;
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
+
+volScalarField& p = thermo.p();
+
+Info<< "\nReading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "compressibleCreatePhi.H"
+
+#include "createMRF.H"
+
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+Info<< "Creating combustion model\n" << endl;
+autoPtr<CombustionModel<basicThermo>> combustion
+(
+    CombustionModel<basicThermo>::New(thermo, turbulence())
+);
+
+dfChemistryModel<basicThermo>* chemistry = combustion->chemistry();
+PtrList<volScalarField>& Y = chemistry->Y();
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+#include "readpRef.H"
+
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+mesh.setFluxRequired(p_rgh.name());
+
+#include "phrghEqn.H"
+
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll(Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(thermo.he());
+
+volScalarField HRR
+(
+    IOobject
+    (
+        "HRR",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(dimEnergy/dimVolume/dimTime, 0)
+);
+
+IOdictionary additionalControlsDict
+(
+    IOobject
+    (
+        "additionalControls",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+Switch solvePrimaryRegion
+(
+    additionalControlsDict.lookup("solvePrimaryRegion")
+);
+
+Switch solvePyrolysisRegion
+(
+    additionalControlsDict.lookupOrDefault<bool>("solvePyrolysisRegion", true)
+);
+
+
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar(p.dimensions()/dimTime, 0)
+);
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+#include "createClouds.H"
+#include "createSurfaceFilmModel.H"
+// #include "createPyrolysisModel.H"
+// #include "createRadiationModel.H"
+// #include "createFvOptions.H"