diff --git a/applications/solvers/faTwoLayerAvalancheFoam/Make/files b/applications/solvers/faTwoLayerAvalancheFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..e1f1c5cf1be61c81168869fd1194502c5ea145e1
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/Make/files
@@ -0,0 +1,3 @@
+faTwoLayerAvalancheFoam.C
+
+EXE = $(FOAM_MODULE_APPBIN)/faTwoLayerAvalancheFoam
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/Make/options b/applications/solvers/faTwoLayerAvalancheFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..8e8eb40484cd32347e7dd50a96ef6126c9ec9051
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/Make/options
@@ -0,0 +1,21 @@
+sinclude $(GENERAL_RULES)/module-path-user
+
+/* Failsafe - user location */
+ifeq (,$(strip $(FOAM_MODULE_APPBIN)))
+ FOAM_MODULE_APPBIN = $(FOAM_USER_APPBIN)
+endif
+ifeq (,$(strip $(FOAM_MODULE_LIBBIN)))
+ FOAM_MODULE_LIBBIN = $(FOAM_USER_LIBBIN)
+endif
+
+EXE_INC = \
+ -I$(LIB_SRC)/finiteArea/lnInclude \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/cfdTools/general/lnInclude \
+ -I../../../src/avalanche/lnInclude
+
+EXE_LIBS = \
+ -lfiniteArea \
+ -lfiniteVolume \
+ -L$(FOAM_LIBBIN) -lfaAvalanche
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/PFcreateFaFields.H b/applications/solvers/faTwoLayerAvalancheFoam/PFcreateFaFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..51a1b7e4cd5c3ca32970d5b3fee9d58e6f83ea7f
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/PFcreateFaFields.H
@@ -0,0 +1,166 @@
+ Info << "Reading field h2" << endl;
+ areaScalarField h2
+ (
+ IOobject
+ (
+ "h2",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ aMesh
+ );
+
+
+ Info << "Reading field Us2" << endl;
+ areaVectorField Us2
+ (
+ IOobject
+ (
+ "Us2",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ aMesh
+ );
+
+ areaVectorField UsDiff
+ (
+ IOobject
+ (
+ "UsDiff",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ Us2-Us1
+ );
+
+ Info << "Reading field c2" << endl;
+ areaScalarField c2
+ (
+ IOobject
+ (
+ "c2",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ aMesh
+ );
+
+ areaScalarField powderLayerAlone
+ (
+ IOobject
+ (
+ "powderLayerAlone",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ aMesh,
+ dimensionedScalar(dimless)
+ );
+
+ edgeScalarField phis2
+ (
+ IOobject
+ (
+ "phis2",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ fac::interpolate(Us2) & aMesh.Le()
+ );
+
+ edgeScalarField phi2s2
+ (
+ IOobject
+ (
+ "phi2s2",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ fac::interpolate(h2*Us2) & aMesh.Le()
+ );
+
+ areaScalarField geff
+ (
+ IOobject
+ (
+ "geff",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ gn
+ );
+
+ areaScalarField Sdp
+ (
+ IOobject
+ (
+ "Sdp",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ aMesh,
+ dimensionedScalar(dimVelocity)
+ );
+
+ areaScalarField Spd
+ (
+ IOobject
+ (
+ "Spd",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ aMesh,
+ dimensionedScalar(dimVelocity)
+ );
+
+ areaScalarField Sap
+ (
+ IOobject
+ (
+ "Sap",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ aMesh,
+ dimensionedScalar(dimVelocity)
+ );
+
+ areaVectorField tau2
+ (
+ IOobject
+ (
+ "tau2",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ aMesh,
+ dimensionedVector(dimPressure/dimDensity)
+ );
+
+ Info << endl;
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/PFcreateSubmodels.H b/applications/solvers/faTwoLayerAvalancheFoam/PFcreateSubmodels.H
new file mode 100644
index 0000000000000000000000000000000000000000..9440ca5546956dfb5575988a6f4408389f13fa9d
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/PFcreateSubmodels.H
@@ -0,0 +1,24 @@
+autoPtr<suspensionFrictionModel> friction2
+(
+ suspensionFrictionModel::New(transportProperties2, UsDiff, h2, c2)
+);
+
+autoPtr<suspensionEntrainmentModel> entrainment2
+(
+ suspensionEntrainmentModel::New(transportProperties2, Us2, h2, hentrain, c2, tau2)
+);
+
+autoPtr<ambientEntrainmentModel> ambientEntrainment2
+(
+ ambientEntrainmentModel::New(transportProperties2, Us2, h2, c2)
+);
+
+autoPtr<suspensionDepositionModel> sedimentation2
+(
+ suspensionDepositionModel::New(transportProperties2, Us2, h2, c2, tau2)
+);
+
+autoPtr<couplingModel> coupling
+(
+ couplingModel::New(transportProperties, Us1, h1, pb1, Us2, h2, c2)
+);
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/PFreadTransportProperties.H b/applications/solvers/faTwoLayerAvalancheFoam/PFreadTransportProperties.H
new file mode 100644
index 0000000000000000000000000000000000000000..8e50736c86e4c71ba37c228fd14a802fc7fccabc
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/PFreadTransportProperties.H
@@ -0,0 +1,40 @@
+dictionary transportProperties2 = transportProperties.subDict("powderCloud");
+
+dimensionedScalar hmin2("hmin", dimLength, transportProperties2);
+
+dimensionedScalar cmin2("cmin", dimless, transportProperties2);
+
+dimensionedScalar h02("h0", dimLength, transportProperties2);
+
+dimensionedScalar u02("u0", dimVelocity, transportProperties2);
+
+dimensionedScalar xi2("xi", dimless, transportProperties2);
+
+dimensionedScalar R2("R", dimless, transportProperties2);
+
+dimensionedScalar Ds2("Ds", dimLength, transportProperties2);
+
+dimensionedScalar nu2("nu", dimViscosity, transportProperties2);
+
+dimensionedScalar packingd2("phi0", dimless, transportProperties2);
+
+dimensionedScalar hwem2("hwem", dimLength, transportProperties2);
+
+dimensionedScalar hentmin2("hentmin", dimLength, transportProperties2);
+
+Switch waterSink2(transportProperties2.getOrDefault<Switch>("waterSink", false));
+
+Switch bindHeight2(transportProperties.getOrDefault<Switch>("bindHeight", true));
+
+Switch correctMomentum2(transportProperties.getOrDefault<Switch>("correctMomentum", true));
+
+Info << "Running powder cloud with" << endl
+ << " xi2 " << xi1 << endl
+ << " hwem " << hwem2 << endl
+ << " hmin " << hmin2 << endl
+ << " cmin " << cmin2 << endl
+ << " h0 " << h02 << endl
+ << " u0 " << u02 << endl
+ << " waterSink is " << waterSink2 << endl
+ << " bindHeight is " << bindHeight2 << endl << endl;
+
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/SHcreateFaFields.H b/applications/solvers/faTwoLayerAvalancheFoam/SHcreateFaFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..8daefd0309c8a47542096bc50ec02185f653dea1
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/SHcreateFaFields.H
@@ -0,0 +1,128 @@
+ Info << "Reading field h1" << endl;
+ areaScalarField h1
+ (
+ IOobject
+ (
+ "h1",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ aMesh
+ );
+
+
+ Info << "Reading field Us1" << endl;
+ areaVectorField Us1
+ (
+ IOobject
+ (
+ "Us1",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ aMesh
+ );
+
+
+ edgeScalarField phis1
+ (
+ IOobject
+ (
+ "phis1",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ fac::interpolate(Us1) & aMesh.Le()
+ );
+
+
+ edgeScalarField phi2s1
+ (
+ IOobject
+ (
+ "phi2s1",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ fac::interpolate(h1*Us1) & aMesh.Le()
+ );
+
+
+ areaVectorField n = aMesh.faceAreaNormals();
+
+ areaScalarField gn
+ (
+ IOobject
+ (
+ "gn",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ g & n
+ );
+
+ areaVectorField gs
+ (
+ IOobject
+ (
+ "gs",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ g - gn*n
+ );
+
+ areaScalarField pb1
+ (
+ IOobject
+ (
+ "pb",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ aMesh,
+ dimensionedScalar(dimPressure)
+ );
+
+ areaVectorField tau1
+ (
+ IOobject
+ (
+ "tau",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ aMesh,
+ dimensionedVector(dimPressure)
+ );
+
+ //entrainment height
+ areaScalarField hentrain
+ (
+ IOobject
+ (
+ "hentrain",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ aMesh,
+ dimensionedScalar(dimLength)
+ );
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/SHcreateSubmodels.H b/applications/solvers/faTwoLayerAvalancheFoam/SHcreateSubmodels.H
new file mode 100644
index 0000000000000000000000000000000000000000..e615796ffc083e7876f3cf60423306fd0b74aa6c
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/SHcreateSubmodels.H
@@ -0,0 +1,14 @@
+autoPtr<frictionModel> friction1
+(
+ frictionModel::New(transportProperties1, Us1, h1, pb1)
+);
+
+autoPtr<entrainmentModel> entrainment1
+(
+ entrainmentModel::New(transportProperties1, Us1, h1, hentrain, pb1, tau1)
+);
+
+autoPtr<depositionModel> deposition1
+(
+ depositionModel::New(transportProperties1, Us1, h1, hentrain, pb1, tau1)
+);
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/SHreadTransportProperties.H b/applications/solvers/faTwoLayerAvalancheFoam/SHreadTransportProperties.H
new file mode 100644
index 0000000000000000000000000000000000000000..3c2dd0c681a3ec3f89da3f84a3c68531aea6c0f1
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/SHreadTransportProperties.H
@@ -0,0 +1,33 @@
+
+IOdictionary transportProperties
+(
+ IOobject
+ (
+ "transportProperties",
+ runTime.constant(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ )
+);
+
+dictionary transportProperties1 = transportProperties.subDict("denseFlow");
+
+dimensionedScalar hmin1("hmin", dimLength, transportProperties1);
+
+dimensionedScalar xi1("xi", dimless, transportProperties1);
+
+dimensionedScalar xit1("xit", dimless, transportProperties1);
+
+dimensionedScalar packingd("phi", dimless, transportProperties1);
+
+Switch pressureFeedback1(transportProperties1.get<Switch>("pressureFeedback"));
+
+Switch bindHeight1(transportProperties1.getOrDefault<Switch>("bindHeight", true));
+
+Info << "Running dense flow with" << endl
+ << " xi " << xi1 << endl
+ << " xi_t" << xit1 << endl
+ << " hmin " << hmin1 << endl
+ << " pressureFeedback is " << pressureFeedback1 << endl
+ << " bindHeight is " << bindHeight1 << endl << endl;
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/advanceDenseFlowModel.H b/applications/solvers/faTwoLayerAvalancheFoam/advanceDenseFlowModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..629bef84005f3363e0e0d55fc96483e7f3f13106
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/advanceDenseFlowModel.H
@@ -0,0 +1,75 @@
+#include "calcBasalstress.H"
+
+const areaVectorField & tauSc1 = friction1->tauSc();
+const areaScalarField & tauSp1 = friction1->tauSp();
+
+faVectorMatrix Us1Eqn
+(
+ fam::ddt(h1, Us1)
+ + xi1*fam::div(phi2s1, Us1)
+ + tauSc1
+ + fam::Sp(tauSp1, Us1)
+ ==
+ gs*h1
+ - fac::grad(pb1*h1/(2.*friction1->rho()))
+ - fam::Sp
+ (
+ xit1*Sdp/packingd,
+ Us1
+ )
+ + Spd*Us2
+);
+
+
+if (!final)
+ Us1Eqn.relax();
+
+SolverPerformance<vector> Us1Residual = solve(Us1Eqn);
+
+tau1 = (tauSc1 + tauSp1*Us1)*friction1->rho();
+phis1 = (fac::interpolate(Us1) & aMesh.Le());
+
+const areaScalarField & Sm1 = entrainment1->Sm();
+const areaScalarField & Sd1 = deposition1->Sd();
+
+
+h1.storePrevIter();
+
+faScalarMatrix h1Eqn
+(
+ fam::ddt(h1)
+ + fam::div(phis1, h1)
+ ==
+ Sm1/(hentrain+dimensionedScalar("hentrainmin", dimLength, 1e-5))*hentrain
+ - fam::Sp
+ (
+ Sd1/(h1 + dimensionedScalar("small", dimLength, 1e-5)),
+ h1
+ )
+ + Spd/(c2+dimensionedScalar("cmin", dimless, 1e-7))*c2/packingd
+ - fam::Sp
+ (
+ Sdp/(h1 + dimensionedScalar("small", dimLength, 1e-5))/packingd,
+ h1
+ )
+);
+
+solverPerformance h1Residual = h1Eqn.solve();
+
+phi2s1 = h1Eqn.flux();
+
+if (!final)
+ h1.relax();
+
+
+// Bind h
+if (bindHeight1)
+{
+ h1 = max
+ (
+ h1,
+ hmin1
+ );
+}
+h1.correctBoundaryConditions();
+Us1.correctBoundaryConditions();
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/advancePowderCloudModel.H b/applications/solvers/faTwoLayerAvalancheFoam/advancePowderCloudModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..1ac657ee164586021e0716a3fec4e8a01c5b18e4
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/advancePowderCloudModel.H
@@ -0,0 +1,114 @@
+
+UsDiff = Us2-Us1;
+
+const areaVectorField & tauSc2 = friction2->tauSc();
+const areaScalarField & tauSp2 = friction2->tauSp();
+
+geff = max(gn - (fac::ndiv(phis2, Us2)&n), dimensionedScalar(dimVelocity/dimTime));
+
+//Parker et al. (1986), Eq. (5)
+faVectorMatrix Us2Eqn
+(
+ fam::ddt(h2, Us2)
+ + xi2*fam::div(phi2s2, Us2)
+ + tauSc2 // ustar**2 expl. part
+ - tauSp2*Us1 // ustar**2 expl. part
+ + fam::Sp(tauSp2, Us2) // ustar**2 impl. part
+ ==
+ R2*gs*c2*h2
+ - 1./2.*R2*fac::grad(geff*c2*h2*h2)
+ + xit1*Sdp*Us1*(1.+R2)
+ - fam::Sp
+ (
+ Spd,
+ Us2
+ )
+);
+
+
+if (correctMomentum2)
+{
+ Us2Eqn += R2*c2*fam::ddt(h2, Us2)
+ + R2*h2*Us2*fac::ddt(c2)
+ + xi2*R2*fam::div(fac::interpolate(c2)*phi2s2, Us2);
+}
+
+if (!final)
+ Us2Eqn.relax();
+
+SolverPerformance<vector> Us2Residual = solve(Us2Eqn);
+
+Us2.correctBoundaryConditions();
+
+phis2 = (fac::interpolate(Us2) & aMesh.Le());
+tau2 = tauSc2 + tauSp2*Us2;
+
+//Calculate sediment entrainment
+areaScalarField Sm2 = entrainment2->Sm();
+
+//no entrainment if there is not a minimal flow thickness and not if the dense flow layer is below
+Sm2 = Sm2*pos(h2-hentmin2)*powderLayerAlone;
+Sm2 = min(Sm2, hentrain/runTime.deltaT());
+
+//Calculate ambient fluid entrainment
+Sap = ambientEntrainment2->Sm();
+Sap = Sap*pos(hwem2-h2); //clip entrainment if maximum water height is reached
+
+//Additional sink term for water
+const areaScalarField vs(R2*geff*Ds2*Ds2/18./nu2); //Terminal velocity
+const areaScalarField dw(pos(h2-hmin2*2.)*vs);
+
+
+//Parker et al. (1986), Eq. (3)
+faScalarMatrix h2Eqn
+(
+ fam::ddt(h2)
+ + fam::div(phis2, h2)
+ ==
+ Sap
+ + Sdp/(h1 + dimensionedScalar("small", dimLength, 1e-5)) * h1 / packingd2
+ - fam::Sp
+ (
+ Spd/(h2 + dimensionedScalar("small", dimLength, 1e-5)) / packingd2,
+ h2
+ )
+);
+
+if (waterSink2)
+{
+ h2Eqn -= fam::Sp(dw/(h2+dimensionedScalar("hmin", dimLength, 1e-5)),h2);
+}
+
+solverPerformance h2Residual = h2Eqn.solve();
+
+phi2s2 = h2Eqn.flux();
+
+if (!final)
+{
+ h2.relax();
+}
+
+if (bindHeight2)
+{
+ h2 = max(h2, hmin2);
+}
+
+
+faScalarMatrix c2Eqn
+(
+ fam::ddt(h2, c2)
+ + fam::div(phi2s2, c2)
+ ==
+ Sm2/(hentrain+dimensionedScalar("hentrainmin", dimLength, 1e-5))*hentrain
+ + Sdp/(h1 + dimensionedScalar("small", dimLength, 1e-5)) * h1
+ - fam::Sp(Spd/(c2+dimensionedScalar("cmin", dimless, 1e-7)), c2)
+);
+
+if(!final)
+ c2Eqn.relax();
+
+solverPerformance c2Residual = c2Eqn.solve();
+
+if (bindHeight2)
+ c2 = max(c2, cmin2);
+
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/calcBasalstress.H b/applications/solvers/faTwoLayerAvalancheFoam/calcBasalstress.H
new file mode 100644
index 0000000000000000000000000000000000000000..233375bc4433accdd567c15cfe1233e9ae9ca689
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/calcBasalstress.H
@@ -0,0 +1,13 @@
+
+pb1 = -friction1->rho()*xi1*fac::ndiv(phi2s1, Us1)&n;
+
+if (pressureFeedback1)
+{
+ pb1 -= fac::ngrad(pb1*h1/2.)&n;
+}
+
+pb1 += friction1->rho()*(g*h1)&n;
+
+// Clipping
+pb1 = Foam::max(pb1, dimensionedScalar(dimPressure));
+
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/faTwoLayerAvalancheFoam.C b/applications/solvers/faTwoLayerAvalancheFoam/faTwoLayerAvalancheFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..97e693bcc9c111bc63274a286b6f47c74aa1284d
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/faTwoLayerAvalancheFoam.C
@@ -0,0 +1,245 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | avalanche module
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Application
+ faSavageHutterFoam
+
+Description
+ A depth-integrated solver for mixed snow avalanches.
+ The dense flow is described with the Savage-Hutter model.
+ The powder cloud is described with the Parker-Fukushima model.
+ The coupling model is new.
+
+ Publication in preparation.
+
+Author
+ Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+#include "frictionModel.H"
+#include "entrainmentModel.H"
+#include "depositionModel.H"
+#include "suspensionFrictionModel.H"
+#include "suspensionEntrainmentModel.H"
+#include "ambientEntrainmentModel.H"
+#include "suspensionDepositionModel.H"
+#include "couplingModel.H"
+#include "SolverPerformance.H"
+
+#define CREATE_FIELDS SHcreateFaFields.H
+#define CREATE_FIELDS_2 PFcreateFaFields.H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+ argList::addNote
+ (
+ "(avalanche)\n"
+ "A depth-integrated solver for mixed snow avalanches using"
+ " Finite Area Methods."
+ );
+
+ #include "postProcessFA.H"
+
+ #include "setRootCase.H"
+ #include "createTime.H"
+ #include "createMesh.H"
+ #include "createFaMesh.H"
+ #include "readGravitationalAcceleration.H"
+
+ #include "SHreadTransportProperties.H"
+ #include "PFreadTransportProperties.H"
+
+ #include "SHcreateFaFields.H"
+ #include "PFcreateFaFields.H"
+
+ #include "SHcreateSubmodels.H"
+ #include "PFcreateSubmodels.H"
+
+ #include "createTimeControls.H"
+
+ Info<< "\nStarting time loop\n" << endl;
+
+ #include "readSolutionControls.H"
+
+ Info<< nl
+ << "Numerical settings" << nl
+ << " max number of iterations " << nCorr << nl
+ << " min number of iterations " << minCorr << nl
+ << " TOL h " << hResidualMax << nl
+ << " TOL Us " << UsResidualMax << nl
+ << " TOL c " << cResidualMax << nl << endl;
+
+ const bool initDeltaT = runTime.controlDict().get<bool>("initDeltaT");
+
+ if (initDeltaT)
+ {
+ Info<< "Initializing Delta T" << endl;
+ #include "readTimeControls.H"
+ #include "surfaceCourantNo1.H"
+ #include "surfaceCourantNo2.H"
+
+ runTime.setDeltaT
+ (
+ min(maxCo/(max(CoNum1, CoNum2) + SMALL)*runTime.deltaT().value(), maxDeltaT)
+ );
+ }
+
+ bool final = false;
+
+ while (runTime.run())
+ {
+ #include "readSolutionControls.H"
+ #include "readTimeControls.H"
+ #include "surfaceCourantNo1.H"
+ #include "surfaceCourantNo2.H"
+
+ scalar CoNum = max(CoNum1, CoNum2);
+ #include "setDeltaT.H"
+
+ runTime++;
+
+ Info<< "Time = " << runTime.timeName() << nl << endl;
+ final = false;
+
+ for (int iCorr = 0; ; iCorr++)
+ {
+
+ //Calculate coupling fluxes
+ Sdp = coupling->Sdp();
+ Sdp = min(Sdp, h1/runTime.deltaT());
+
+ //Calculate sedimation rate
+ Spd = sedimentation2->Sd();
+ Spd = min(Spd, h2*c2/runTime.deltaT());
+
+ //Info << "Upward mass flux = [" << min(Sdp).value() << ", " << max(Sdp).value() << "]" << endl;
+ //Info << "Upward momentum flux = [" << min(Sdp*mag(Us1)).value() << ", " << max(Sdp*mag(Us1)).value() << "]" << endl;
+
+ powderLayerAlone = pos(hmin1*3-h1);
+
+ #include "advanceDenseFlowModel.H"
+
+ #include "advancePowderCloudModel.H"
+
+ //Keep track of intact snow cover
+ faScalarMatrix hentrainEqn
+ (
+ fam::ddt(hentrain)
+ ==
+ Sd1
+ - fam::Sp
+ (
+ Sm1/(hentrain + dimensionedScalar("small", dimLength, SMALL)),
+ hentrain
+ )
+ - fam::Sp(
+ Sm2/(hentrain+dimensionedScalar("small", dimLength, SMALL))/packingd,
+ hentrain)
+
+ );
+
+ hentrainEqn.solve();
+
+
+
+ if (final)
+ {
+ if
+ (
+ h1Residual.initialResidual() < hResidualMax
+ && mag(Us1Residual.initialResidual()) < UsResidualMax
+ && h2Residual.initialResidual() < hResidualMax
+ && mag(Us2Residual.initialResidual()) < UsResidualMax
+ && c2Residual.initialResidual() < cResidualMax
+ )
+ {
+ Info<< "Layer 1: Reached residual in h1 = "
+ << h1Residual.initialResidual()
+ << " < " << hResidualMax
+ << " and in Us1 = "
+ << Us1Residual.initialResidual()
+ << " < " << UsResidualMax << endl
+ << "Layer 2: Reached residual in h2 = "
+ << h2Residual.initialResidual()
+ << " < " << hResidualMax
+ << " and in Us2 = "
+ << Us2Residual.initialResidual()
+ << " < " << UsResidualMax
+ << ", in c2 = "
+ << c2Residual.initialResidual()
+ << " < " << cResidualMax << ", " << endl
+ << "stopping loop!" << endl;
+ }
+ else
+ {
+ Info<< "Reached maximum numbers of iterations, "
+ << "stopping loop!" << endl;
+ }
+
+ const scalar Vsl = gSum((hentrain*aMesh.S())())/packingd.value();
+ const scalar Vdl = gSum((h1*aMesh.S())())/packingd.value();
+ const scalar Vpl = gSum((h2*c2*aMesh.S())());
+ Info << "Grain vol. in staionary layer = " << Vsl << endl;
+ Info << "Grain vol. in dense flow layer = " << Vdl << endl;
+ Info << "Grain vol. in powder cloud layer = " << Vpl << endl;
+ Info << "Grain vol. in all layers = " << Vsl + Vdl + Vpl << endl;
+ break;
+ }
+
+ if
+ (
+ (
+ h1Residual.initialResidual() < hResidualMax
+ && mag(Us1Residual.initialResidual()) < UsResidualMax
+ && iCorr >= minCorr
+ )
+ || iCorr >= nCorr
+ )
+ {
+ final = true;
+ }
+ }
+
+ if (runTime.outputTime())
+ {
+ runTime.write();
+ }
+
+ runTime.printExecutionTime(Info);
+ }
+
+
+ Info<< nl << "End" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/readSolutionControls.H b/applications/solvers/faTwoLayerAvalancheFoam/readSolutionControls.H
new file mode 100644
index 0000000000000000000000000000000000000000..5e537e3e073219af07575ab47ffc2b0d1094f5d2
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/readSolutionControls.H
@@ -0,0 +1,16 @@
+
+ scalar nCorr =
+ aMesh.solutionDict().getOrDefault<int>("nOuterCorrectors", 50);
+
+ scalar minCorr =
+ aMesh.solutionDict().getOrDefault<int>("minCorrectors", 3);
+
+ scalar hResidualMax =
+ aMesh.solutionDict().getOrDefault<scalar>("hResidualMax", 0.);
+
+ scalar UsResidualMax =
+ aMesh.solutionDict().getOrDefault<scalar>("UsResidualMax", 1);
+
+ scalar cResidualMax =
+ aMesh.solutionDict().getOrDefault<scalar>("cResidualMax", 1);
+
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/surfaceCourantNo1.H b/applications/solvers/faTwoLayerAvalancheFoam/surfaceCourantNo1.H
new file mode 100644
index 0000000000000000000000000000000000000000..abc0b6985ccf40d55f78f3e808a4c4495db12c48
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/surfaceCourantNo1.H
@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | foam-extend: Open Source CFD
+ \\ / O peration | Version: 3.2
+ \\ / A nd | Web: http://www.foam-extend.org
+ \\/ M anipulation | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+ This file is part of foam-extend.
+
+ foam-extend is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by the
+ Free Software Foundation, either version 3 of the License, or (at your
+ option) any later version.
+
+ foam-extend is distributed in the hope that it will be useful, but
+ WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
+
+Global
+ surfaceCourantNo
+
+Author
+ Hrvoje Jasak, Wikki Ltd. All rights reserved.
+
+Description
+ Calculates and outputs the mean and maximum Courant Numbers for the
+ Finite Area method.
+
+\*---------------------------------------------------------------------------*/
+
+scalar CoNum1 = 0.0;
+scalar meanCoNum1 = 0.0;
+scalar velMag1 = 0.0;
+if (aMesh.nInternalEdges())
+{
+ edgeScalarField edgeSpeed
+ (
+ (xi1*fac::interpolate(Us1) & aMesh.Le()/aMesh.magLe())
+ + Foam::sqrt(max(dimensionedScalar("gh_min", dimensionSet(0, 2, -2, 0, 0, 0, 0), 0), fac::interpolate(h1)*fac::interpolate(g&n)))
+ );
+
+ velMag1 = max(edgeSpeed.internalField()).value();
+
+ CoNum1 = max(aMesh.edgeInterpolation::deltaCoeffs()*edgeSpeed.internalField()).value()*runTime.deltaT().value();
+
+ meanCoNum1 = sum(aMesh.edgeInterpolation::deltaCoeffs()*edgeSpeed.internalField()*aMesh.magLe()).value()/gSum(aMesh.magLe())*runTime.deltaT().value();
+}
+
+
+Info<< "Layer 1: Courant Number mean: " << meanCoNum1
+ << " max: " << CoNum1
+ << " velocity max magnitude: " << velMag1 << endl;
+
+// ************************************************************************* //
diff --git a/applications/solvers/faTwoLayerAvalancheFoam/surfaceCourantNo2.H b/applications/solvers/faTwoLayerAvalancheFoam/surfaceCourantNo2.H
new file mode 100644
index 0000000000000000000000000000000000000000..77e64e3807d40eaf3d0cfcc2120fe8041cd6e599
--- /dev/null
+++ b/applications/solvers/faTwoLayerAvalancheFoam/surfaceCourantNo2.H
@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | foam-extend: Open Source CFD
+ \\ / O peration | Version: 3.2
+ \\ / A nd | Web: http://www.foam-extend.org
+ \\/ M anipulation | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+ This file is part of foam-extend.
+
+ foam-extend is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by the
+ Free Software Foundation, either version 3 of the License, or (at your
+ option) any later version.
+
+ foam-extend is distributed in the hope that it will be useful, but
+ WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
+
+Global
+ surfaceCourantNo
+
+Author
+ Hrvoje Jasak, Wikki Ltd. All rights reserved.
+
+Description
+ Calculates and outputs the mean and maximum Courant Numbers for the
+ Finite Area method.
+
+\*---------------------------------------------------------------------------*/
+
+scalar CoNum2 = 0.0;
+scalar meanCoNum2 = 0.0;
+scalar velMag2 = 0.0;
+if (aMesh.nInternalEdges())
+{
+ edgeScalarField edgeSpeed
+ (
+ (xi2*fac::interpolate(Us2) & aMesh.Le()/aMesh.magLe())
+ + Foam::sqrt(max(dimensionedScalar("gh_min", dimensionSet(0, 2, -2, 0, 0, 0, 0), 0), fac::interpolate(R2*c2*h2)*fac::interpolate(gn)))
+ );
+
+ velMag2 = max(edgeSpeed.internalField()).value();
+
+ CoNum2 = max(aMesh.edgeInterpolation::deltaCoeffs()*edgeSpeed.internalField()).value()*runTime.deltaT().value();
+
+ meanCoNum2 = sum(aMesh.edgeInterpolation::deltaCoeffs()*edgeSpeed.internalField()*aMesh.magLe()).value()/gSum(aMesh.magLe())*runTime.deltaT().value();
+}
+
+
+Info<< "Layer 2: Courant Number mean: " << meanCoNum2
+ << " max: " << CoNum2
+ << " velocity max magnitude: " << velMag2 << endl;
+
+// ************************************************************************* //
diff --git a/src/avalanche/Make/files b/src/avalanche/Make/files
index 52f8cdddcc3f1edd0451f77a5ba8462a9f3fc961..a19a1cba4750c24260cf03ff0cc779b5f6dcd587 100644
--- a/src/avalanche/Make/files
+++ b/src/avalanche/Make/files
@@ -1,3 +1,8 @@
+coupling = coupling
+$(coupling)/couplingModel/couplingModel.C
+$(coupling)/couplingModel/couplingModelNew.C
+$(coupling)/couplingInertial/couplingInertial.C
+
deposition = deposition
$(deposition)/depositionModel/depositionModel.C
$(deposition)/depositionModel/depositionModelNew.C
diff --git a/src/avalanche/coupling/couplingInertial/couplingInertial.C b/src/avalanche/coupling/couplingInertial/couplingInertial.C
new file mode 100644
index 0000000000000000000000000000000000000000..ece4e981aee82a83fc2ceedb90d281d60e6b9763
--- /dev/null
+++ b/src/avalanche/coupling/couplingInertial/couplingInertial.C
@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | avalanche module
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Author
+ Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+#include "couplingInertial.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace couplingModels
+{
+ defineTypeNameAndDebug(couplingInertial, 0);
+ addToRunTimeSelectionTable(couplingModel, couplingInertial, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::couplingModels::couplingInertial::couplingInertial
+(
+ const dictionary& couplingProperties,
+ const areaVectorField& Us1,
+ const areaScalarField& h1,
+ const areaScalarField& pb1,
+ const areaVectorField& Us2,
+ const areaScalarField& h2,
+ const areaScalarField& c2
+)
+:
+ couplingModel(type(), couplingProperties, Us1, h1, pb1, Us2, h2, c2),
+ I0_("I0", dimless, coeffDict_),
+ u0_("u0", dimless, coeffDict_),
+ d_("d", dimLength, coeffDict_),
+ rhos_("rhos", dimDensity, coeffDict_),
+ I_
+ (
+ IOobject
+ (
+ "I",
+ Us1_.time().timeName(),
+ Us1_.db(),
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ Us1_.mesh(),
+ dimensionedScalar(dimless)
+ )
+{
+ Info << " " << u0_ << nl
+ << " " << d_ << nl
+ << " " << rhos_ << nl << endl;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+const Foam::areaScalarField&
+Foam::couplingModels::couplingInertial::Sdp() const
+{
+ //Inertial number at the base, assuming it is constant along flow depth
+ I_ = 5./2.*mag(Us1_)/(h1_+dimensionedScalar(dimLength, 1e-2))*d_*sqrt(rhos_/(pb1_+dimensionedScalar(dimPressure, 1e-2)));
+ Sdp_ = mag(Us1_)*u0_*mag(I_-I0_)*h1_/(h1_+dimensionedScalar(dimLength, 1e-2));
+ return Sdp_;
+}
+
+
+bool Foam::couplingModels::couplingInertial::read
+(
+ const dictionary& couplingProperties
+)
+{
+ readDict(type(), couplingProperties);
+
+ coeffDict_.readEntry("u0", u0_);
+ coeffDict_.readEntry("I0", I0_);
+ coeffDict_.readEntry("d", d_);
+ coeffDict_.readEntry("rhos", rhos_);
+
+ return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/avalanche/coupling/couplingInertial/couplingInertial.H b/src/avalanche/coupling/couplingInertial/couplingInertial.H
new file mode 100644
index 0000000000000000000000000000000000000000..cc8fdb73a2534c3645a75df54ef09d6be02a1b5e
--- /dev/null
+++ b/src/avalanche/coupling/couplingInertial/couplingInertial.H
@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | avalanche module
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+ Foam::couplingModels::couplingInertial
+Description
+ A coupling model for mixed snow avalanches following the scaling with the
+ Inertial number I.
+
+ Publication in preparation.
+
+SourceFiles
+ couplingInertial.C
+
+Author
+ Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef couplingInertial_H
+#define couplingInertial_H
+
+#include "couplingModel.H"
+#include "dimensionedScalar.H"
+#include "volFields.H"
+#include "IOdictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "dimensionedScalar.H"
+#include "tmp.H"
+#include "autoPtr.H"
+#include "faMatrices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace couplingModels
+{
+
+/*---------------------------------------------------------------------------*\
+ Class Erosionenergy Declaration
+\*---------------------------------------------------------------------------*/
+
+class couplingInertial
+:
+ public couplingModel
+{
+ // Private data
+
+
+ //- reference inertial number
+ dimensionedScalar I0_;
+
+ //- coupling volume flux at I-I0 = 1
+ dimensionedScalar u0_;
+
+ //- particle diameter for calculation I
+ dimensionedScalar d_;
+
+ //- rhos for calculating I
+ dimensionedScalar rhos_;
+
+ //- Inertial number
+ mutable areaScalarField I_;
+public:
+
+ //- Runtime type information
+ TypeName("couplingInertial");
+
+
+ // Constructors
+
+ //- Construct from components
+ couplingInertial
+ (
+ const dictionary& couplingProperties,
+ const areaVectorField& Us1,
+ const areaScalarField& h1,
+ const areaScalarField& pb1,
+ const areaVectorField& Us2,
+ const areaScalarField& h2,
+ const areaScalarField& c2
+ );
+
+
+ //- Destructor
+ virtual ~couplingInertial() = default;
+
+
+ // Member Functions
+
+ //- Return the volume flux
+ virtual const areaScalarField& Sdp() const;
+
+ //- Read couplingProperties dictionary
+ virtual bool read(const dictionary& couplingProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace frictionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/avalanche/coupling/couplingModel/couplingModel.C b/src/avalanche/coupling/couplingModel/couplingModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..5f980c9cb103b0c8b21a2924d1841ca32914d36b
--- /dev/null
+++ b/src/avalanche/coupling/couplingModel/couplingModel.C
@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | avalanche module
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Author
+ Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "couplingModel.H"
+#include "fvCFD.H"
+#include "faCFD.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(couplingModel, 0);
+ defineRunTimeSelectionTable(couplingModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+void Foam::couplingModel::readDict
+(
+ const word& type,
+ const dictionary& dict
+)
+{
+ couplingProperties_ = dict;
+ coeffDict_ = couplingProperties_.optionalSubDict(type + "Coeffs");
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::couplingModel::couplingModel
+(
+ const word& type,
+ const dictionary& couplingProperties,
+ const areaVectorField& Us1,
+ const areaScalarField& h1,
+ const areaScalarField& pb1,
+ const areaVectorField& Us2,
+ const areaScalarField& h2,
+ const areaScalarField& c2
+)
+:
+ couplingProperties_(couplingProperties),
+ coeffDict_
+ (
+ couplingProperties_.optionalSubDict(type + "Coeffs")
+ ),
+ Us1_(Us1),
+ h1_(h1),
+ pb1_(pb1),
+ Us2_(Us2),
+ h2_(h2),
+ c2_(c2),
+ Sdp_
+ (
+ IOobject
+ (
+ "Sdp",
+ Us1_.time().timeName(),
+ Us1_.db(),
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ Us1_.mesh(),
+ dimensionedScalar(dimVelocity)
+ )
+{
+ Info<< " with " << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/src/avalanche/coupling/couplingModel/couplingModel.H b/src/avalanche/coupling/couplingModel/couplingModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..74eb4a8e164c82fed95be49a487cf26fda6ab293
--- /dev/null
+++ b/src/avalanche/coupling/couplingModel/couplingModel.H
@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | avalanche module
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Namespace
+ Foam::couplingModel
+
+Description
+ A namespace for various coupling model implementations for mixed snow
+ avalanches. This model couples two shallow flow layers.
+
+Class
+ Foam::couplingModel
+
+Description
+ An abstract base class for coupling models.
+
+SourceFiles
+ couplingModel.C
+ couplingModelNew.C
+
+Author
+ Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef couplingModel_H
+#define couplingModel_H
+
+#include "IOdictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "dimensionedScalar.H"
+#include "tmp.H"
+#include "autoPtr.H"
+#include "faMatrices.H"
+#include "areaFieldsFwd.H"
+#include "FieldFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class couplingModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class couplingModel
+{
+protected:
+
+ // Protected data
+
+ dictionary couplingProperties_;
+
+ //- Model coefficients dictionary
+ dictionary coeffDict_;
+
+ //- Reference to the surface velocity field
+ const areaVectorField& Us1_;
+
+ //- Reference to the flow height field
+ const areaScalarField& h1_;
+
+ //- Reference to the basal pressure field
+ const areaScalarField& pb1_;
+
+ //- Reference to the surface velocity field
+ const areaVectorField& Us2_;
+
+ //- Reference to the flow height field
+ const areaScalarField& h2_;
+
+ //- Reference to the phase fraction field
+ const areaScalarField& c2_;
+
+ //- Flux term
+ mutable areaScalarField Sdp_;
+
+
+ // Protected Member Functions
+
+ //- Read/update the couplingProperties and coeffDict dictionaries
+ void readDict(const word& type, const dictionary& dict);
+
+
+ //- Disallow copy construct
+ couplingModel(const couplingModel&) = delete;
+
+ //- Disallow default bitwise assignment
+ void operator=(const couplingModel&) = delete;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("couplingModel");
+
+
+ // Declare run-time constructor selection table
+
+#ifndef SWIG
+ declareRunTimeSelectionTable
+ (
+ autoPtr,
+ couplingModel,
+ dictionary,
+ (
+ const dictionary& couplingProperties,
+ const areaVectorField& Us1,
+ const areaScalarField& h1,
+ const areaScalarField& pb1,
+ const areaVectorField& Us2,
+ const areaScalarField& h2,
+ const areaScalarField& c2
+ ),
+ (couplingProperties, Us1, h1, pb1, Us2, h2, c2)
+ );
+#endif
+
+
+ // Selectors
+
+ //- Return a reference to the selected coupling model
+ static autoPtr<couplingModel> New
+ (
+ const dictionary& couplingProperties,
+ const areaVectorField& Us1,
+ const areaScalarField& h1,
+ const areaScalarField& pb1,
+ const areaVectorField& Us2,
+ const areaScalarField& h2,
+ const areaScalarField& c2
+ );
+
+
+ // Constructors
+
+ //- Construct from components
+ couplingModel
+ (
+ const word& type,
+ const dictionary& couplingProperties,
+ const areaVectorField& Us1,
+ const areaScalarField& h1,
+ const areaScalarField& pb1,
+ const areaVectorField& Us2,
+ const areaScalarField& h2,
+ const areaScalarField& c2
+ );
+
+
+ //- Destructor
+ virtual ~couplingModel() = default;
+
+
+ // Member Functions
+
+ //- Read couplingProperties dictionary
+ virtual bool read(const dictionary& couplingProperties) = 0;
+
+ //- Return the coupling properties dictionary
+ const dictionary& couplingProperties() const
+ {
+ return couplingProperties_;
+ }
+
+ //- Const access to the model coefficients dictionary
+ virtual const dictionary& coeffDict() const
+ {
+ return coeffDict_;
+ }
+
+ //- Return the Source by coupling
+ virtual const areaScalarField& Sdp() const = 0;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/avalanche/coupling/couplingModel/couplingModelNew.C b/src/avalanche/coupling/couplingModel/couplingModelNew.C
new file mode 100644
index 0000000000000000000000000000000000000000..84694c2e6619c30ac03d84c33d7bf43213792a9c
--- /dev/null
+++ b/src/avalanche/coupling/couplingModel/couplingModelNew.C
@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | avalanche module
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Author
+ Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "couplingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::couplingModel> Foam::couplingModel::New
+(
+ const dictionary& dict,
+ const areaVectorField& Us1,
+ const areaScalarField& h1,
+ const areaScalarField& pb1,
+ const areaVectorField& Us2,
+ const areaScalarField& h2,
+ const areaScalarField& c2
+)
+{
+ const word modelName(dict.get<word>("couplingModel"));
+
+ Info<< "Selecting coupling model " << modelName << endl;
+
+ auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelName);
+
+ if (!cstrIter.found())
+ {
+ FatalIOErrorInLookup
+ (
+ dict,
+ "couplingModel",
+ modelName,
+ *dictionaryConstructorTablePtr_
+ ) << exit(FatalIOError);
+ }
+
+ return autoPtr<couplingModel>
+ (
+ cstrIter()(dict, Us1, h1, pb1, Us2, h2, c2)
+ );
+}
+
+
+// ************************************************************************* //