[com1DFA] Add TJEM entrainment model

avaframe/com1DFA/DFAfunctionsCython.pxd

@@ -8,6 +8,8 @@
 
 cpdef (double, double) computeEntMassAndForce(double, double, double, double, double, double, double)
 
+cpdef (double, double, double) computeEntMassAndForceTjem(double, double, double, double, double, double, double)
+
 cpdef double computeDetMass(double, double, double, double)
 
 cpdef double computeResForce(double, double, double, double, int, int)

avaframe/com1DFA/DFAfunctionsCython.pyx

@@ -27,7 +27,7 @@ import avaframe.in3Utils.geoTrans as geoTrans
 log = logging.getLogger(__name__)
 
 
-def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType):
+def computeForceC(cfg, particles, fields, dem, int frictType, int entrType, int resistanceType):
   """ compute forces acting on the particles (without the SPH component)
 
   Cython implementation implementation
@@ -44,6 +44,8 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
       dictionary with dem information
   frictType: int
     identifier for friction law to be used
+  entrType: int
+    identifier for entrainment law to be used
   resistanceType: int
     identifier for resistance model to be used
 
@@ -73,6 +75,7 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
   cdef double BSamosAtMedium = cfg.getfloat('Bsamosatmedium')
   cdef double RSamosAtMedium = cfg.getfloat('Rsamosatmedium')
   cdef double entEroEnergy = cfg.getfloat('entEroEnergy')
+  cdef double entShearStrength = cfg.getfloat('entShearStrength')
   cdef double entShearResistance = cfg.getfloat('entShearResistance')
   cdef double entDefResistance = cfg.getfloat('entDefResistance')
   cdef double rho = cfg.getfloat('rho')
@@ -375,6 +378,14 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
           else:
             tau = 0.0
 
+      # select entrainment model
+      # compute entrained mass and update tau (in case of TJEM)
+      entrMassCell = entrMassRaster[indCellY, indCellX]
+      if entrMassCell > 0.0:
+        if entrType == 1:
+            dm, areaEntrPart = computeEntMassAndForce(dt, entrMassCell, areaPart, uMag, tau, entEroEnergy, rhoEnt)
+        if entrType == 2:
+            dm, areaEntrPart, tau = computeEntMassAndForceTjem(dt, entrMassCell, areaPart, uMag, tau, entShearStrength, velMagMin)
       # adding bottom shear resistance contribution
       # norm of the friction force >=0 (if 0 -> detatchment)
       forceBotTang = - areaPart * tau
@@ -389,14 +400,12 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
           uMagRes = uMag
         forceFrict[k] = forceFrict[k] - forceBotTang/uMagRes
 
-      # compute entrained mass and if wetSnow case compute total enthalpy of particles with entrainment
-      entrMassCell = entrMassRaster[indCellY, indCellX]
+      # continue calculating entrainment
       if entrMassCell > 0.0:
-        # compute entrained mass
-        dm, areaEntrPart = computeEntMassAndForce(dt, entrMassCell, areaPart, uMag, tau, entEroEnergy, rhoEnt)
         # speed loss due to energy loss due to entrained mass
         dEnergyEntr = areaEntrPart * entShearResistance + dm * entDefResistance
 
+        # if wetSnow case compute total enthalpy of particles with entrainment
         if frictType == 4 and dm > 0.0:
           # enthalpy change due to entrained mass
           entrEnthCell = entrEnthRaster[indCellY, indCellX]
@@ -542,6 +551,60 @@ cpdef (double, double) computeEntMassAndForce(double dt, double entrMassCell,
 
   return dm, areaEntrPart
 
+cpdef (double, double, double) computeEntMassAndForceTjem(double dt, double entrMassCell,
+                                              double areaPart, double uMag,
+                                              double tau, double entShearStrength,
+                                              double velMagMin):
+  """ compute force component due to entrained mass
+
+  Parameters
+  ----------
+  dt: float
+    time step
+  entrMassCell : float
+      available mass for entrainement
+  areaPart : float
+      particle area
+  uMag : float
+      particle speed (velocity magnitude)
+  tau : float
+      bottom shear stress
+  entShearStrength: float
+    shear strength of the entrainable snow mass
+  velMagMin : float
+    minimum velocity magnitude, avoids division by 0 for resting particles
+
+  Returns
+  -------
+  dm : float
+      entrained mass
+  areaEntrPart : float
+      Area for entrainment energy loss computation
+  tau : float
+      bottom shear stress, limited to the shear strength of the snow cover if entrainment occurs
+  """
+  # compute entrained mass
+  cdef double dm = 0.0
+  cdef double areaEntrPart = 0.0
+
+  # avoid very small uMag so as not to divide by 0
+  if entrMassCell > 0 and uMag > velMagMin:
+
+      # erosion: erode according to shear strength of snow cover
+
+      # in contrast to MoTVoellmy tauC here is not scaled with the flow density
+      # because com1DFA does not scale shear stress at the bottom of the avalanche with the flow density while MoTVoellmy does
+      tauC = entShearStrength
+
+      dm = max(0.0, tau - tauC) * areaPart * dt / uMag
+      # cannot entrain more than the mass available under the particle
+      dm = min(dm, entrMassCell * areaPart)
+      areaEntrPart = areaPart
+      if dm > 0:
+        tau = min(tau, tauC)
+
+  return dm, areaEntrPart, tau
+
 
 cpdef double computeDetMass(double dt, double detCell,
                                               double areaPart, double uMag):

avaframe/com1DFA/com1DFA.py

@@ -2189,6 +2189,16 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
     frictType = frictModelsList.index(frictModel) + 1
     log.debug("Friction Model used: %s, %s" % (frictModelsList[frictType - 1], frictType))
 
+    # derive entrainment type
+    # turn entrainment model into integer
+    entrModelsList = [
+        "samosat",
+        "tjem",
+    ]
+    entrModel = cfgGen["entrModel"].lower()
+    entrType = entrModelsList.index(entrModel) + 1
+    log.debug("Entrainment Model used: %s, %s" % (entrModelsList[entrType - 1], entrType))
+
     # turn resistance model into integer
     ResModel = cfgGen["ResistanceModel"].lower()
     ResModelsList = [
@@ -2301,6 +2311,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
             dem,
             tCPU,
             frictType,
+            entrType,
             resistanceType,
             inputSimLines["reportAreaInfo"],
         )
@@ -2723,6 +2734,7 @@ def computeEulerTimeStep(
     dem,
     tCPU,
     frictType,
+    entrType,
     resistanceType,
     reportAreaInfo,
 ):
@@ -2744,6 +2756,8 @@ def computeEulerTimeStep(
         computation time dictionary
     frictType: int
         indicator for chosen type of friction model
+    entrType: int
+        indicator for chosen type of entrainment model
     resistanceType: int
         identifier for chosen type of resistance model
 
@@ -2776,7 +2790,7 @@ def computeEulerTimeStep(
     startTime = time.time()
     # loop version of the compute force
     log.debug("Compute Force C")
-    particles, force, fields = DFAfunC.computeForceC(cfg, particles, fields, dem, frictType, resistanceType)
+    particles, force, fields = DFAfunC.computeForceC(cfg, particles, fields, dem, frictType, entrType, resistanceType)
     tCPUForce = time.time() - startTime
     tCPU["timeForce"] = tCPU["timeForce"] + tCPUForce
     # compute lateral force (SPH component of the calculation)

avaframe/com1DFA/com1DFACfg.ini

@@ -476,11 +476,23 @@ forestThMin = 0.6
 forestVMax = 40.
 forestThMax = 10.
 
+#++++++++++++ Entrainment Model
+### Available Models:
+#samosAT
+#tjem
+entrModel = samosAT
 
 #++++++++++++ Entrainment Erosion Energy
-# Used to determine speed loss via energy loss due to entrained mass
+# Used to determine speed loss via energy loss due to entrained mass when using samosAT Entrainment Model [m²/s²]
 # Select entEroEnergy = 0 to enable frontal plowing; entEroEnergy > 0 means basal erosion (without frontal plowing)
 entEroEnergy = 5000
+
+#++++++++++++ Entrainment Shear Strength
+# Used to determine the strength of the erodible snow cover when using TJEM Entrainment Model [Pa]
+entShearStrength = 500
+
+#++++++++++++Entrainment Energy loss
+# Used to determine the energy loss because of entrainment
 entShearResistance = 0
 entDefResistance = 0
 

avaframe/tests/test_DFAfunctionsCython.py

@@ -129,6 +129,52 @@ def test_computeEntMassAndForce(capfd):
     assert areaEntrPart == 1
 
 
+def test_computeEntMassAndForceTjem():
+    """Test the computeEntMassAndForceTjem function"""
+    dt = 0.1
+    areaPart = 1
+    uMag = 10
+    entShearStrength = 500
+    velMagMin = 1e-3
+
+    # no entrainable mass
+    dm, areaEntrPart, tau = DFAfunC.computeEntMassAndForceTjem(
+        dt, 0, areaPart, uMag, 1000, entShearStrength, velMagMin
+    )
+    assert dm == 0
+    assert areaEntrPart == 0
+    assert tau == 1000
+
+    # tau above shear strength: entrainment and tau limited to shear strength
+    dm, areaEntrPart, tau = DFAfunC.computeEntMassAndForceTjem(
+        dt, 200, areaPart, uMag, 1000, entShearStrength, velMagMin
+    )
+    assert dm == pytest.approx(5.0)
+    assert areaEntrPart == 1
+    assert tau == entShearStrength
+
+    # tau below shear strength: no entrainment, tau unchanged
+    dm, areaEntrPart, tau = DFAfunC.computeEntMassAndForceTjem(
+        dt, 200, areaPart, uMag, 400, entShearStrength, velMagMin
+    )
+    assert dm == 0
+    assert tau == 400
+
+    # resting particle: no division by 0, no entrainment
+    dm, areaEntrPart, tau = DFAfunC.computeEntMassAndForceTjem(
+        dt, 200, areaPart, 0, 1000, entShearStrength, velMagMin
+    )
+    assert dm == 0
+    assert areaEntrPart == 0
+    assert tau == 1000
+
+    # entrained mass is capped by the mass available under the particle
+    dm, areaEntrPart, tau = DFAfunC.computeEntMassAndForceTjem(
+        dt, 1, areaPart, uMag, 1000, entShearStrength, velMagMin
+    )
+    assert dm == pytest.approx(1.0)
+
+
 def test_computeDetMass(capfd):
     """Test the computeDetMass function"""
     dt = 0.1
@@ -1018,6 +1064,7 @@ def test_computeForceC():
         "Bsamosatmedium": "4.13",
         "Rsamosatmedium": "0.05",
         "entEroEnergy": "5000",
+        "entShearStrength": "500",
         "entShearResistance": "0",
         "entDefResistance": "0",
         "xsivoellmy": "4000.",
@@ -1108,14 +1155,17 @@ def test_computeForceC():
     # resistanceType
     resistanceType = 1  # default; since cResRaster = 0 there is no resistance
 
+    # entrainment model (1 = samosAT)
+    entrType = 1
+
     # apply function
     # initialize array function fricition force
     forceFrictArray = np.zeros_like(test_tauArray, dtype=np.float64)
     # looping over ObrienAndJulien, HerschelAndBulkley and Bingham with initial config settings (see above)
     # i = index of result in array, j = frictType
     for i, j in enumerate(range(10, 13)):
         particles, force, *_ = DFAfunC.computeForceC(
-            cfg["GENERAL"], particles, fields, dem, j, resistanceType
+            cfg["GENERAL"], particles, fields, dem, j, entrType, resistanceType
         )
         forceFrictArray[i] = force["forceFrict"][0]
 
@@ -1177,7 +1227,7 @@ def test_computeForceC():
     cfg["GENERAL"]["etaObrienAndJulien"] = "0"
     cfg["GENERAL"]["tauyObrienAndJulien"] = "0"
 
-    particles, force, *_ = DFAfunC.computeForceC(cfg["GENERAL"], particles, fields, dem, 10, resistanceType)
+    particles, force, *_ = DFAfunC.computeForceC(cfg["GENERAL"], particles, fields, dem, 10, entrType, resistanceType)
     forceFrictArray[3] = force["forceFrict"][0]
     # compute tauy and eta via alpha, beta and cv
     test_etaObrienAndJulien = test_alpha1EtaObrienAndJulien * math.exp(
@@ -1212,7 +1262,7 @@ def test_computeForceC():
     # set tauyCfg = 0
     cfg["GENERAL"]["tauyHerschelAndBulkley"] = "0"
 
-    particles, force, *_ = DFAfunC.computeForceC(cfg["GENERAL"], particles, fields, dem, 11, resistanceType)
+    particles, force, *_ = DFAfunC.computeForceC(cfg["GENERAL"], particles, fields, dem, 11, entrType, resistanceType)
     forceFrictArray[4] = force["forceFrict"][0]
     # compute tauy via alpha, beta and cv
     test_tauyHerschelandBulkley = test_alpha2TauyHerschelAndBulkley * math.exp(
@@ -1240,7 +1290,7 @@ def test_computeForceC():
     cfg["GENERAL"]["etaBingham"] = "0"
     cfg["GENERAL"]["tauyBingham"] = "0"
 
-    particles, force, *_ = DFAfunC.computeForceC(cfg["GENERAL"], particles, fields, dem, 12, resistanceType)
+    particles, force, *_ = DFAfunC.computeForceC(cfg["GENERAL"], particles, fields, dem, 12, entrType, resistanceType)
     forceFrictArray[5] = force["forceFrict"][0]
     # compute tauy and eta via alpha, beta and cv
     test_etaBingham = test_alpha1EtaBingham * math.exp(test_beta1EtaBingham * test_cvSediment)
@@ -1262,7 +1312,7 @@ def test_computeForceC():
         test_uMagRes = test_uMag
     test_forceFrictArray = test_forceFrict - test_forceBotTang / test_uMagRes
 
-    assert (forceFrictArray == test_forceFrictArray).all()
+    assert np.allclose(forceFrictArray, test_forceFrictArray, rtol=1e-12, atol=0)
 
 
 def test_updateInitialVelocity():