v11legacy

.vscode/settings.json

@@ -2,6 +2,13 @@
     "files.associations": {
         "MonodGrowth.C": "cpp",
         "phreeqcRM.C": "cpp",
-        "simpleFirstOrderKineticMole.C": "cpp"
+        "simpleFirstOrderKineticMole.C": "cpp",
+        "darcyFoam.C": "cpp",
+        "constantVelocityFoam.C": "cpp",
+        "impesFoam.C": "cpp",
+        "dbsFoam.C": "cpp",
+        "PDRFoam.C": "cpp",
+        "ostream": "cpp",
+        "KozenyCarman.C": "cpp"
     }
 }

Allwmake

@@ -18,7 +18,7 @@ wmake libraries/porousBoundaryConditions
 # Compiling solvers
 wmake solvers/darcyFoam
 wmake solvers/impesFoam
-wmake solvers/dbsSimpleFoam
+#wmake solvers/dbsSimpleFoam
 wmake solvers/dbsFoam
 wmake solvers/constantVelocityFoam
 #wmake solvers/anisoImpesFoam

libraries/porousModels/Make/files

@@ -22,14 +22,6 @@ porousModels/absolutePermeabilityModels/heterogeneousScalarConstant/heterogeneou
 porousModels/absolutePermeabilityModels/ColebrookWhite/ColebrookWhite.C
 porousModels/absolutePermeabilityModels/none/none.C
 
-porousModels/dispersionModels/dispersionModel/dispersionModel.C
-porousModels/dispersionModels/dispersionModel/dispersionModelNew.C
-porousModels/dispersionModels/linearDispersion/linearDispersion.C
-porousModels/dispersionModels/archiesLaw/archiesLaw.C
-porousModels/dispersionModels/taylorAris/taylorAris.C
-porousModels/dispersionModels/diffusionOnly/diffusionOnly.C
-porousModels/dispersionModels/dispersionNone/dispersionNone.C
-
 
 porousModels/dispersionTensorModels/dispersionTensorModel/dispersionTensorModel.C
 porousModels/dispersionTensorModels/dispersionTensorModel/dispersionTensorModelNew.C

libraries/porousModels/geochemicalModels/phreeqcRM/phreeqcRM.C

@@ -949,7 +949,7 @@ void Foam::geochemicalModels::phreeqcRM::updateMineralDistribution()
         {
             forAll(mineralList_,s)
             {
-                Ys_[s].storeOldTime();
+                Ys_[s].storeOldTimes();
                 if(activatePhaseEquilibrium_[s] == true)
                 {
                     Ys_[s][i] = so[2*s*nxyz_ + i]*Vm_[s].value()*1e3;

solvers/constantVelocityFoam/constantVelocityFoam.C

@@ -31,18 +31,24 @@ Developers
     - Cyprien Soulaine
 \*---------------------------------------------------------------------------*/
 
-#include "fvCFD.H"
+#include "argList.H"
 #include "incompressiblePhase.H"
 #include "geochemicalModel.H"
 
+#include "fvcFlux.H"
+#include "fvcSurfaceIntegrate.H"
+
+
+using namespace Foam;
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
 
 //    argList::addOption("phase","a","specify the phase name");
 //    Foam::argList args(argc,argv);
-    #include "setRootCaseLists.H"
+    #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"

solvers/constantVelocityFoam/createFields.H

@@ -1,6 +1,22 @@
 //word phaseName("");
 //if (args.optionFound("phase")) phaseName = args.option("phase");
 
+Info << nl << "Reading field p\n" << endl;
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("p",dimPressure,1e5)
+);
+
+
 Info << nl << "Reading field U\n" << endl;
 volVectorField U
 (

solvers/darcyFoam/darcyFoam.C

@@ -30,11 +30,26 @@ Description
 Developers
     - Pierre Horgue
     - 10/02/2020 - CS : modified with the porousModel class
+    - 17/11/2023 - CS : upgrade to v11
+
 \*---------------------------------------------------------------------------*/
 
-#include "fvCFD.H"
+
+#include "argList.H"
 #include "incompressiblePhase.H"
 #include "geochemicalModel.H"
+
+#include "uniformDimensionedFields.H"
+//#include "fvcDdt.H"
+//#include "fvcGrad.H"
+#include "fvcFlux.H"
+#include "fvcReconstruct.H"
+
+#include "fvmDdt.H"
+//#include "fvmDiv.H"
+#include "fvmLaplacian.H"
+
+using namespace Foam;
 //#include "porousModel.H"
 //#include "sourceEventFile.H"
 
@@ -45,7 +60,7 @@ int main(int argc, char *argv[])
 
 //    argList::addOption("phase","a","specify the phase name");
 //    Foam::argList args(argc,argv);
-    #include "setRootCaseLists.H"
+    #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
@@ -88,7 +103,7 @@ Info<< "init OK " << nl << endl; //
         U.correctBoundaryConditions();
 
         runTime.write();
-        if(runTime.outputTime())
+        if(runTime.writeTime())
         {
             sourceTerm.write();
         }

solvers/dbsFoam/Make/options

@@ -1,20 +1,29 @@
 EXE_INC = \
+    -I$(LIB_SRC)/physicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
+    -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/ThermophysicalTransportModels/thermophysicalTransportModel/lnInclude \
+    -I$(LIB_SRC)/ThermophysicalTransportModels/fluid/lnInclude \
+    -I$(LIB_SRC)/ThermophysicalTransportModels/fluidThermo/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
     -I../../libraries/porousModels/lnInclude \
     -I../../libraries/otherModels/lnInclude \
     -I../../ThirdParty/phreeqcrm/include \
 
 EXE_LIBS = \
+    -lmeshTools \
+    -lspecie \
+    -lmomentumTransportModels \
+    -lcompressibleMomentumTransportModels \
     -lfiniteVolume \
     -lfvModels \
     -lfvConstraints \
-    -lsampling \
-    -lmeshTools \
     -lspecie \
     -L$(FOAM_USER_LIBBIN) \
     -lporousModels \
     -lotherModels \
-    -L../../ThirdParty/phreeqcrm/lib/ -lphreeqcrm-3.6.2
+    -L../../ThirdParty/phreeqcrm/lib/ -lphreeqcrm-3.6.2

solvers/dbsFoam/StCourantNo.H

@@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Global
+    CourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    const scalarField sumPhi
+    (
+        fvc::surfaceSum(mag(phiSt))().primitiveField()
+       /rho.primitiveField()
+    );
+
+    StCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    const scalar meanStCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+
+    Info<< "St courant Number mean: " << meanStCoNum
+        << " max: " << StCoNum << endl;
+}
+
+// ************************************************************************* //

solvers/dbsFoam/UEqn.H

@@ -1,26 +1,21 @@
-// Solve the Momentum equation
+    fvVectorMatrix UEqn
+    (
+        1./(eps+SMALL)*(fvm::ddt(U) + fvm::div(phi, U))
+    //  - fvm::laplacian(turbulence->nu()/(eps+SMALL),U, "laplacian(nu,U)")
+    //  + fvm::Sp(turbulence->nu()*invK,U)
+    - fvm::laplacian(nu/(eps+SMALL),U, "laplacian(nu,U)")
+    + fvm::Sp(nu*invK,U)
+    ==
+        fvModels.source(U)
+    );
 
-MRF.correctBoundaryVelocity(U);
+    UEqn.relax();
 
-tmp<fvVectorMatrix> tUEqn
-(
-    1./(eps+SMALL)*(fvm::ddt(U) + fvm::div(phi, U))
-//  - fvm::laplacian(turbulence->nu()/(eps+SMALL),U, "laplacian(nu,U)")
-//  + fvm::Sp(turbulence->nu()*invK,U)
-  - fvm::laplacian(nu/(eps+SMALL),U, "laplacian(nu,U)")
-  + fvm::Sp(nu*invK,U)
- ==
-    fvModels.source(U)
-);
-fvVectorMatrix& UEqn = tUEqn.ref();
+    fvConstraints.constrain(UEqn);
 
-UEqn.relax();
+    if (pimple.momentumPredictor())
+    {
+        solve(UEqn == -fvc::grad(p));
 
-fvConstraints.constrain(UEqn);
-
-if (pimple.momentumPredictor())
-{
-    solve(UEqn == -fvc::grad(p));
-
-    fvConstraints.constrain(U);
-}
+        fvConstraints.constrain(U);
+    }

solvers/dbsFoam/createFieldRefs.H

@@ -0,0 +1 @@
+const volScalarField& psi = thermo.psi();

solvers/dbsFoam/createFields.H

@@ -1,4 +1,3 @@
-#include "createRDeltaT.H"
 
 Info<< "Reading field p\n" << endl;
 volScalarField p
@@ -14,13 +13,13 @@ volScalarField p
     mesh
 );
 
-Info<< "Reading field U\n" << endl;
+Info<< "\nReading field U\n" << endl;
 volVectorField U
 (
     IOobject
     (
         "U",
-        runTime.timeName(),
+        runTime.name(),
         mesh,
         IOobject::MUST_READ,
         IOobject::AUTO_WRITE
@@ -30,8 +29,6 @@ volVectorField U
 
 #include "createPhi.H"
 
-pressureReference pressureReference(p, pimple.dict());
-
 mesh.schemes().setFluxRequired(p.name());
 
 
@@ -51,7 +48,27 @@ volScalarField sourceTerm
 );
 
 
-IOdictionary transportProperties
+
+/*
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.name(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar(p.dimensions()/dimTime, 0)
+);
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+*/
+
+
+IOdictionary physicalProperties
 (
     IOobject
     (
@@ -63,10 +80,7 @@ IOdictionary transportProperties
     )
 );
 
-//autoPtr<geochemicalModel> geochemistry
-//  = geochemicalModel::New(mesh,transportProperties);
-
-geochemicalModel geochemistry(mesh,transportProperties);
+geochemicalModel geochemistry(mesh,physicalProperties);
 
 
 Info << nl << "Reading/generating porosity field eps" << endl;
@@ -78,35 +92,22 @@ const volScalarField & invK = geochemistry.inversePermeability();
 
 Switch activateHeleShaw
 (
-    transportProperties.lookupOrDefault<Switch>("activateHeleShaw",false)
+    physicalProperties.lookupOrDefault<Switch>("activateHeleShaw",false)
 );
 
 
 volScalarField HeleShawCorrection ("HeleShawCorrection", 0.0*invK);
 
 if (activateHeleShaw)
 {
-    Foam::absolutePermeabilityModels::HeleShaw HeleShaw(mesh,transportProperties);
+    Foam::absolutePermeabilityModels::HeleShaw HeleShaw(mesh,physicalProperties);
     HeleShawCorrection = HeleShaw.inversePermeability();
     Info << "HeleShaw correction is activated" << nl << endl;
 }
 
 const volScalarField & nu = geochemistry.nu();
 
 
-/*
-Info << nl << "Reading/generating momentumTransportModel";
-//autoPtr<viscosityModel> viscosity(viscosityModel::New(mesh));
-
-
-autoPtr<incompressible::momentumTransportModel> turbulence
-(
-    incompressible::momentumTransportModel::New(U, phi, geochemistry.viscosity())
-);
-Info << " ... OK" << endl;
-*/
-
-
 #include "createMRF.H"
 #include "createFvModels.H"
 #include "createFvConstraints.H"