[com1DFA]: Add options for computing resistance (especially dealing with hEff)

avaframe/com1DFA/DFAfunctionsCython.pxd

@@ -10,7 +10,7 @@ cpdef (double, double) computeEntMassAndForce(double, double, double, double, do
 
 cpdef double computeDetMass(double, double, double, double)
 
-cpdef double computeResForce(double, double, double, double, double, double, int)
+cpdef double computeResForce(double, double, double, double, double, double, int, int)
 
 cdef (double, double, double) addArtificialViscosity(double, double, double, double, double, double, double, double,
                                                      int, int, double, double, double, double,

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):
+def computeForceC(cfg, particles, fields, dem, int frictType, 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):
       dictionary with dem information
   frictType: int
     identifier for friction law to be used
+  resistanceType: int
+    identifier for resistance model to be used
 
   Returns
   -------
@@ -369,7 +371,7 @@ def computeForceC(cfg, particles, fields, dem, int frictType):
 
       # adding resistance force due to obstacles
       cResCell = cResRaster[indCellY][indCellX]
-      cResPart = computeResForce(hRes, h, areaPart, rho, cResCell, uMag, explicitFriction)
+      cResPart = computeResForce(hRes, h, areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
       forceFrict[k] = forceFrict[k] - cResPart
 
       uxArray[k] = ux
@@ -500,8 +502,8 @@ cpdef double computeDetMass(double dt, double detCell,
   return dmDet
 
 
-cpdef double computeResForce(double hRes, double h, double areaPart, double rho,
-                             double cResCell, double uMag, int explicitFriction):
+cpdef double computeResForce(double hRes, double h, double areaPart, double rho, double cResCell,
+                             double uMag, int explicitFriction, int resistanceType):
   """ compute force component due to resistance
 
   Parameters
@@ -520,6 +522,8 @@ cpdef double computeResForce(double hRes, double h, double areaPart, double rho,
       particle speed (velocity magnitude)
   explicitFriction: int
     if 1 add resistance with an explicit method. Implicit otherwise
+  resistanceType: int
+    identifier for resistance model to be used
 
   Returns
   -------
@@ -530,11 +534,21 @@ cpdef double computeResForce(double hRes, double h, double areaPart, double rho,
   cdef double cRecResPart
   if(h < hRes):
       hResEff = h
+  # explicit formulation (explicitFriction == 1)
   if explicitFriction == 1:
-    # explicit formulation
-    cRecResPart = - rho * areaPart * hResEff * cResCell * uMag * uMag
+    if resistanceType == 1:
+      # cRes
+      cRecResPart = - rho * areaPart * hResEff * cResCell * uMag * uMag
+    elif resistanceType == 2:
+      # cResH
+      cRecResPart = - rho * areaPart * cResCell * uMag * uMag
   elif explicitFriction == 0:
-    cRecResPart = - rho * areaPart * hResEff * cResCell * uMag
+    if resistanceType == 1:
+      # cRes
+      cRecResPart = - rho * areaPart * hResEff * cResCell * uMag
+    elif resistanceType == 2:
+      # cResH
+      cRecResPart = - rho * areaPart * cResCell * uMag
   return cRecResPart
 
 

avaframe/com1DFA/com1DFA.py

@@ -1637,11 +1637,15 @@ def initializeResistance(cfg, dem, simTypeActual, resLine, reportAreaInfo, thres
     reportAreaInfo: dict
         simulation area information dictionary completed with entrainment area info
     """
-    cRes = cfg.getfloat("cRes")
     K = cfg.getfloat("detK")
     detrainment = cfg.getboolean("detrainment")
     detWithoutRes = cfg.getboolean("detWithoutRes")
 
+    ResModel = cfg["ResistanceModel"].lower()
+    if ResModel == "cres":
+        cRes = cfg.getfloat("cRes")
+    if ResModel == "cresh":
+        cRes = cfg.getfloat("cResH")
     # read dem header
     header = dem["originalHeader"]
     ncols = header["ncols"]
@@ -1754,6 +1758,15 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, simHash=""):
     frictType = frictModelsList.index(frictModel) + 1
     log.debug("Friction Model used: %s, %s" % (frictModelsList[frictType - 1], frictType))
 
+    # turn resistance model into integer
+    ResModel = cfgGen["ResistanceModel"].lower()
+    ResModelsList = [
+        "cres",
+        "cresh",
+    ]
+    resistanceType = ResModelsList.index(ResModel) + 1
+    log.debug("Resistance Model used: %s, %s" % (ResModelsList[resistanceType - 1], resistanceType))
+
     # Initialise Lists to save fields and add initial time step
     particlesList = []
     fieldsList = []
@@ -1809,7 +1822,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, simHash=""):
         log.debug("Computing time step t = %f s, dt = %f s" % (t, dt))
         # Perform computations
         particles, fields, zPartArray0, tCPU = computeEulerTimeStep(
-            cfgGen, particles, fields, zPartArray0, dem, tCPU, frictType
+            cfgGen, particles, fields, zPartArray0, dem, tCPU, frictType, resistanceType
         )
         # set max values of fields to dataframe
         if cfg["VISUALISATION"].getboolean("createRangeTimeDiagram"):
@@ -2126,7 +2139,7 @@ def writeMBFile(infoDict, avaDir, logName):
             )
 
 
-def computeEulerTimeStep(cfg, particles, fields, zPartArray0, dem, tCPU, frictType):
+def computeEulerTimeStep(cfg, particles, fields, zPartArray0, dem, tCPU, frictType, resistanceType):
     """compute next time step using an euler forward scheme
 
     Parameters
@@ -2145,6 +2158,8 @@ def computeEulerTimeStep(cfg, particles, fields, zPartArray0, dem, tCPU, frictTy
         computation time dictionary
     frictType: int
         indicator for chosen type of friction model
+    resistanceType: int
+        identifier for chosen type of resistance model
 
     Returns
     -------
@@ -2159,7 +2174,7 @@ def computeEulerTimeStep(cfg, particles, fields, zPartArray0, dem, tCPU, frictTy
     startTime = time.time()
     # loop version of the compute force
     log.debug("Compute Force C")
-    particles, force, fields = DFAfunC.computeForceC(cfg, particles, fields, dem, frictType)
+    particles, force, fields = DFAfunC.computeForceC(cfg, particles, fields, dem, frictType, resistanceType)
     tCPUForce = time.time() - startTime
     tCPU["timeForce"] = tCPU["timeForce"] + tCPUForce
 

avaframe/com1DFA/com1DFACfg.ini

@@ -350,10 +350,16 @@ entTempRef = -10.0
 cpIce = 2050.0
 
 #++++++++++++ Resistance force parameters
-# height of obstacles [m]
+# Choose one resistance model:
+# cRes : F = - cRes *hEff * rho * A* u²; hEff = min(hRes,h) [cRes] = 1/m
+# cResH : F = - cResH * rho * A* u², [cResH] =[]
+ResistanceModel = cRes
+# height of the obstacles [m]
 hRes = 10
-# characteristic coefficient
+# resistance parameters for different models
 cRes = 0.003
+cResH = 0.003
+
 #++++++++++++ Detrainment parameter
 # Flag if snow detrainment at obstacles (in resistance area) is computed
 # detrainment can be added to resistance force (DetAndRes = True) or computed solely (DetWithoutRes = True)
@@ -362,6 +368,7 @@ detWithoutRes = False
 # detrainment parameter defined by Feistl et al. (2014)
 # value need to be refined
 detK = 20  
+
 #++++++++++++ Entrainment Erosion Energy
 # Used to determine speed loss via energy loss due to entrained mass
 entEroEnergy = 5000

avaframe/tests/test_DFAfunctionsCython.py

@@ -133,23 +133,39 @@ def test_computeResForce(capfd):
     cResCell = 1
     uMag = 10
     explicitFriction = 0
+    #cRes
+    resistanceType = 1
     cResPart = DFAfunC.computeResForce(
-        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction)
-#    print(cResPart)
+        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
+    print(cResPart)
     assert cResPart == -2000
 
     h = 3
     cResPart = DFAfunC.computeResForce(
-        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction)
-#    print(cResPart)
+        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
+    print(cResPart)
     assert cResPart == -4000
 
     explicitFriction = 1
     cResPart = DFAfunC.computeResForce(
-        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction)
-#    print(cResPart)
+        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
+    print(cResPart)
     assert cResPart == -40000
 
+    #cResH
+    explicitFriction = 0
+    resistanceType = 2
+    cResPart = DFAfunC.computeResForce(
+        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
+    print(cResPart)
+    assert cResPart == -2000
+
+    h = 1
+    cResPart = DFAfunC.computeResForce(
+        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
+    print(cResPart)
+    assert cResPart == -2000
+
 
 def test_account4FrictionForce(capfd):
     """ Test the account4FrictionForce function"""

avaframe/tests/test_com1DFA.py

@@ -690,7 +690,11 @@ def test_initializeResistance():
 
     # setup required input
     cfg = configparser.ConfigParser()
-    cfg["GENERAL"] = {"cRes": 0.003, "detK": 10, "detrainment": False, "detWithoutRes": False}
+    cfg["GENERAL"] = {"cRes": 0.003,
+                      "ResistanceModel": "cRes",
+                       "detK": 10, 
+                      "detrainment": False, 
+                      "detWithoutRes": False}
 
     nrows = 11
     ncols = 15
@@ -737,7 +741,11 @@ def test_initializeResistance():
     assert reportAreaInfo["resistance"] == "Yes"
     assert reportAreaInfo["detrainment"] == "No"
 
-    cfg["GENERAL"] = {"cRes": 0.003, "detK": 10.0, "detrainment": True, "detWithoutRes": False}
+    cfg["GENERAL"] = {"cRes": 0.003,
+                      "ResistanceModel": "cRes",
+                       "detK": 10.0, 
+                      "detrainment": True, 
+                      "detWithoutRes": False}
     cResRaster, detRaster, reportAreaInfo = com1DFA.initializeResistance(
         cfg["GENERAL"], dem, simTypeActual, resLine, reportAreaInfo, thresholdPointInPoly
     )
@@ -751,7 +759,11 @@ def test_initializeResistance():
     assert reportAreaInfo["resistance"] == "Yes"
     assert reportAreaInfo["detrainment"] == "Yes"
 
-    cfg["GENERAL"] = {"cRes": 0.003, "detK": 10.0, "detrainment": True, "detWithoutRes": True}
+    cfg["GENERAL"] = {"cRes": 0.003,
+                      "ResistanceModel": "cRes",
+                       "detK": 10.0, 
+                      "detrainment": True, 
+                      "detWithoutRes": True}
     cResRaster, detRaster, reportAreaInfo = com1DFA.initializeResistance(
         cfg["GENERAL"], dem, simTypeActual, resLine, reportAreaInfo, thresholdPointInPoly
     )
@@ -765,6 +777,18 @@ def test_initializeResistance():
     assert reportAreaInfo["resistance"] == "No"
     assert reportAreaInfo["detrainment"] == "Yes"
 
+    cfg["GENERAL"] = {"cResH": 0.003,
+                      "ResistanceModel": "cResH",
+                       "detK": 10.0, 
+                      "detrainment": False, 
+                      "detWithoutRes": False}
+    cResRaster, detRaster, reportAreaInfo = com1DFA.initializeResistance(
+        cfg["GENERAL"], dem, simTypeActual, resLine, reportAreaInfo, thresholdPointInPoly
+    )
+    assert np.array_equal(cResRaster, testArray)
+    assert np.sum(cResRaster) == 0.363
+    assert reportAreaInfo["resistance"] == "Yes"
+
 
 def test_setDEMOriginToZero():
     """test if origin is set to zero"""
@@ -1841,6 +1865,7 @@ def test_initializeSimulation(tmp_path):
         "entTempRef": "-10.",
         "cpIce": "2050.",
         "TIni": "-10.",
+        "ResistanceModel": "cRes"
     }
     # setup dem input
     demHeader = {}

docs/theoryCom1DFA.rst

@@ -217,6 +217,14 @@ coefficient :math:`c_{\text{res}}` that depends on the structure of the obstacle
 
 Note: in previous versions, this formula included information about tree diameter, tree spacing, etc. Please check out previous documentation versions for details.
 
+There is also the option to compute the resistance force :math:`F_i^{\text{res}}`  independent of the effective height :math:`h^{\text{eff}}`:
+
+- .. math::
+   F_i^{\text{res}} = -c_{\text{resH}}\,\rho_0\,A\,
+    \overline{u}^2\,
+    \frac{\overline{u}_i}{\|\overline{u}\|}
+
+with characteristic coefficient :math:`c_{\text{resH}}`.
 
 
 Surface integral forces: