diff --git a/avaframe/com1DFA/DFAfunctionsCython.pxd b/avaframe/com1DFA/DFAfunctionsCython.pxd index 7987a1417..22611f0f4 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pxd +++ b/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) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 8738b850e..edb6d4259 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/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): diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index e928523d3..1d40cfd48 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/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) diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index 386282123..6144aa715 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/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 diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py index 123ede57c..c47d4e0fd 100644 --- a/avaframe/tests/test_DFAfunctionsCython.py +++ b/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,6 +1155,9 @@ 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) @@ -1115,7 +1165,7 @@ def test_computeForceC(): # 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():