From f500a857303a2bad026264f5b7077ea44403da58 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Tue, 22 Sep 2026 15:24:30 +0200 Subject: [PATCH 01/21] changes to .ini file --- avaframe/com1DFA/com1DFACfg.ini | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index 386282123..6b950cc0d 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -476,14 +476,23 @@ forestThMin = 0.6 forestVMax = 40. forestThMax = 10. +#++++++++++++ Entrainment Model +### Available Models: +#avaframe +#tjem +entrModel = avaframe #++++++++++++ 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 Avaframe Entrainment Model # Select entEroEnergy = 0 to enable frontal plowing; entEroEnergy > 0 means basal erosion (without frontal plowing) entEroEnergy = 5000 entShearResistance = 0 entDefResistance = 0 +#++++++++++++ Entrainment Shear Strength +# Used to determine the strength of the erodible snow cover when using TJEM Entrainment Model +entShearStrength = 500 + #++++++++++++ Deposition, Erosion and adaptive surface # adapt the topography every time step # activate with 1 From a5c69481d03dad372d808c8b12a38e6a707172d1 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Tue, 22 Sep 2026 15:34:20 +0200 Subject: [PATCH 02/21] add entrType parameter to com1DFA --- avaframe/com1DFA/com1DFA.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index e928523d3..c35b1bb76 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 = [ + "avaframe", + "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 = [ From 3351410c0da7b121ff988c8513ebe8ea6861618f Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Tue, 22 Sep 2026 15:35:48 +0200 Subject: [PATCH 03/21] change some functions to include parameter entrType --- avaframe/com1DFA/com1DFA.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index c35b1bb76..7fe58cad1 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -2198,6 +2198,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si entrModel = cfgGen["entrModel"].lower() entrType = entrModelsList.index(entrModel) + 1 log.debug("Entrainment Model used: %s, %s" % (entrModelsList[entrType - 1], entrType)) + print("Entrainment Model used: %s, %s" % (entrModelsList[entrType - 1], entrType)) # turn resistance model into integer ResModel = cfgGen["ResistanceModel"].lower() @@ -2311,6 +2312,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si dem, tCPU, frictType, + entrType, resistanceType, inputSimLines["reportAreaInfo"], ) @@ -2733,6 +2735,7 @@ def computeEulerTimeStep( dem, tCPU, frictType, + entrType, resistanceType, reportAreaInfo, ): From 209eb6e8470ed134eade9c11ab06f0ec9854e84e Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Tue, 22 Sep 2026 15:39:47 +0200 Subject: [PATCH 04/21] add TJEM entrainment model to Cython functions --- avaframe/com1DFA/DFAfunctionsCython.pxd | 2 + avaframe/com1DFA/DFAfunctionsCython.pyx | 55 ++++++++++++++++++++++++- 2 files changed, 55 insertions(+), 2 deletions(-) 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..96aee7266 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') @@ -393,7 +396,10 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType entrMassCell = entrMassRaster[indCellY, indCellX] if entrMassCell > 0.0: # compute entrained mass - dm, areaEntrPart = computeEntMassAndForce(dt, entrMassCell, areaPart, uMag, tau, entEroEnergy, rhoEnt) + 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, rhoEnt) # speed loss due to energy loss due to entrained mass dEnergyEntr = areaEntrPart * entShearResistance + dm * entDefResistance @@ -542,6 +548,51 @@ 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 rhoEnt): + """ 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 + + Returns + ------- + dm : float + entrained mass + areaEntrPart : float + Area for entrainement energy loss computation + entEroEnergy: float + erosion entrainment energy constant + rhoEnt: float + entrainement density + """ + cdef double width, ABotSwiped + # compute entrained mass + cdef double dm = 0.0 + cdef double areaEntrPart = 0.0 + + if entrMassCell > 0: + + # erosion: erode according to shear strength of snow cover + tau_c = entShearStrength + dm = max(0, tau - tau_c) * areaPart * dt / uMag + areaEntrPart = areaPart + tau = min(tau, tau_c) + + return dm, areaEntrPart, tau + cpdef double computeDetMass(double dt, double detCell, double areaPart, double uMag): From e2f52b5154218a18b3bd4fdbc6d9f348ee53a6b1 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 23 Sep 2026 09:35:48 +0200 Subject: [PATCH 05/21] add parameter to function call --- avaframe/com1DFA/com1DFA.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index 7fe58cad1..f4a52ed1e 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -2789,7 +2789,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) From 408c410333100813c1fdfdab8ee78c7aed4f4f0f Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 23 Sep 2026 09:40:00 +0200 Subject: [PATCH 06/21] update description of function computeEulerTimeStep --- avaframe/com1DFA/com1DFA.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index f4a52ed1e..992203dfe 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -2757,6 +2757,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 From 37fc73087600c86c692aa7376aba56e0ab75c175 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 23 Sep 2026 15:26:21 +0200 Subject: [PATCH 07/21] changes to Cython function computeEntMassAndForceTjem --- avaframe/com1DFA/DFAfunctionsCython.pyx | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 96aee7266..9b9c70030 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -572,9 +572,9 @@ cpdef (double, double, double) computeEntMassAndForceTjem(double dt, double entr dm : float entrained mass areaEntrPart : float - Area for entrainement energy loss computation - entEroEnergy: float - erosion entrainment energy constant + Area for entrainment energy loss computation + entShearStrength: float + shear strength of the entrainable snow mass rhoEnt: float entrainement density """ @@ -586,10 +586,13 @@ cpdef (double, double, double) computeEntMassAndForceTjem(double dt, double entr if entrMassCell > 0: # erosion: erode according to shear strength of snow cover + # in contrast to MoTVoellmy tau_c 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 tau_c = entShearStrength dm = max(0, tau - tau_c) * areaPart * dt / uMag areaEntrPart = areaPart - tau = min(tau, tau_c) + if dm > 0: + tau = min(tau, tau_c) return dm, areaEntrPart, tau From e4e6e2648eb52b9a0921582207bb6c827ad13e21 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Thu, 24 Sep 2026 11:01:26 +0200 Subject: [PATCH 08/21] move energy loss parameters to more suitable location --- avaframe/com1DFA/com1DFACfg.ini | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index 6b950cc0d..b5e9abeff 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -486,13 +486,16 @@ entrModel = avaframe # Used to determine speed loss via energy loss due to entrained mass when using Avaframe Entrainment Model # Select entEroEnergy = 0 to enable frontal plowing; entEroEnergy > 0 means basal erosion (without frontal plowing) entEroEnergy = 5000 -entShearResistance = 0 -entDefResistance = 0 #++++++++++++ Entrainment Shear Strength # Used to determine the strength of the erodible snow cover when using TJEM Entrainment Model entShearStrength = 500 +#++++++++++++Entrainment Energy loss +# Used to determine the energy loss because of entrainment +entShearResistance = 0 +entDefResistance = 0 + #++++++++++++ Deposition, Erosion and adaptive surface # adapt the topography every time step # activate with 1 From 6880e4efcf7583e460f8cf11b677fcb85770602a Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 28 Oct 2026 10:24:33 +0100 Subject: [PATCH 09/21] move code so the changed tau has an effect --- avaframe/com1DFA/DFAfunctionsCython.pyx | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 9b9c70030..d48b11a49 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -378,6 +378,15 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int entrType, int else: tau = 0.0 + # compute entrained mass and if wetSnow case compute total enthalpy of particles with entrainment + entrMassCell = entrMassRaster[indCellY, indCellX] + if entrMassCell > 0.0: + # compute entrained mass + 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, rhoEnt) + # adding bottom shear resistance contribution # norm of the friction force >=0 (if 0 -> detatchment) forceBotTang = - areaPart * tau @@ -392,14 +401,6 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int entrType, int 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] - if entrMassCell > 0.0: - # compute entrained mass - 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, rhoEnt) # speed loss due to energy loss due to entrained mass dEnergyEntr = areaEntrPart * entShearResistance + dm * entDefResistance From 7e5bef1c6de13f3b2b7cf4208414fbc94cff7bcf Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 28 Oct 2026 10:34:12 +0100 Subject: [PATCH 10/21] Revert "move code so the changed tau has an effect" This reverts commit 39c99bff7e2d122a88046a2ea502243a7e11537a. --- avaframe/com1DFA/DFAfunctionsCython.pyx | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index d48b11a49..9b9c70030 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -378,15 +378,6 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int entrType, int else: tau = 0.0 - # compute entrained mass and if wetSnow case compute total enthalpy of particles with entrainment - entrMassCell = entrMassRaster[indCellY, indCellX] - if entrMassCell > 0.0: - # compute entrained mass - 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, rhoEnt) - # adding bottom shear resistance contribution # norm of the friction force >=0 (if 0 -> detatchment) forceBotTang = - areaPart * tau @@ -401,6 +392,14 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int entrType, int 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] + if entrMassCell > 0.0: + # compute entrained mass + 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, rhoEnt) # speed loss due to energy loss due to entrained mass dEnergyEntr = areaEntrPart * entShearResistance + dm * entDefResistance From 0b0531a31a92a04654f184246e29b15480d18716 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 28 Oct 2026 15:34:13 +0100 Subject: [PATCH 11/21] move code so the changed tau has an effect --- avaframe/com1DFA/DFAfunctionsCython.pyx | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 9b9c70030..ab567b515 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -378,6 +378,14 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int entrType, int 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, rhoEnt) # adding bottom shear resistance contribution # norm of the friction force >=0 (if 0 -> detatchment) forceBotTang = - areaPart * tau @@ -392,17 +400,12 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int entrType, int 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 - 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, 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] From eb9211908afde588f5beaf719ba6c66d41b2a124 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 28 Oct 2026 15:47:38 +0100 Subject: [PATCH 12/21] minor corrections --- avaframe/com1DFA/DFAfunctionsCython.pyx | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index ab567b515..20bd326d0 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -569,6 +569,10 @@ cpdef (double, double, double) computeEntMassAndForceTjem(double dt, double entr particle speed (velocity magnitude) tau : float bottom shear stress + entShearStrength: float + shear strength of the entrainable snow mass + rhoEnt: float + entrainement density Returns ------- @@ -576,12 +580,7 @@ cpdef (double, double, double) computeEntMassAndForceTjem(double dt, double entr entrained mass areaEntrPart : float Area for entrainment energy loss computation - entShearStrength: float - shear strength of the entrainable snow mass - rhoEnt: float - entrainement density """ - cdef double width, ABotSwiped # compute entrained mass cdef double dm = 0.0 cdef double areaEntrPart = 0.0 @@ -589,13 +588,15 @@ cpdef (double, double, double) computeEntMassAndForceTjem(double dt, double entr if entrMassCell > 0: # erosion: erode according to shear strength of snow cover - # in contrast to MoTVoellmy tau_c here is not scaled with the flow density + + # 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 - tau_c = entShearStrength - dm = max(0, tau - tau_c) * areaPart * dt / uMag + tauC = entShearStrength + + dm = max(0.0, tau - tauC) * areaPart * dt / uMag areaEntrPart = areaPart if dm > 0: - tau = min(tau, tau_c) + tau = min(tau, tauC) return dm, areaEntrPart, tau From bd149e04db893694c10d1e60b7fd353f16801f26 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 28 Oct 2026 15:57:06 +0100 Subject: [PATCH 13/21] include velMagMin in order to avoid division by 0 --- avaframe/com1DFA/DFAfunctionsCython.pyx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 20bd326d0..365d82ec2 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -585,7 +585,8 @@ cpdef (double, double, double) computeEntMassAndForceTjem(double dt, double entr cdef double dm = 0.0 cdef double areaEntrPart = 0.0 - if entrMassCell > 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 From 19d2cb1eb5ad52a024d255902b887170b788f682 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 28 Oct 2026 16:02:23 +0100 Subject: [PATCH 14/21] remove print --- avaframe/com1DFA/com1DFA.py | 1 - 1 file changed, 1 deletion(-) diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index 992203dfe..343d2ea2e 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -2198,7 +2198,6 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si entrModel = cfgGen["entrModel"].lower() entrType = entrModelsList.index(entrModel) + 1 log.debug("Entrainment Model used: %s, %s" % (entrModelsList[entrType - 1], entrType)) - print("Entrainment Model used: %s, %s" % (entrModelsList[entrType - 1], entrType)) # turn resistance model into integer ResModel = cfgGen["ResistanceModel"].lower() From 967e26cd32b535d92a079970ea9c5c65ba308064 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 28 Oct 2026 16:12:49 +0100 Subject: [PATCH 15/21] change name of entrainment model and add units to description --- avaframe/com1DFA/com1DFACfg.ini | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index b5e9abeff..179ec1e08 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -478,17 +478,17 @@ forestThMax = 10. #++++++++++++ Entrainment Model ### Available Models: -#avaframe +#samosAT #tjem entrModel = avaframe #++++++++++++ Entrainment Erosion Energy -# Used to determine speed loss via energy loss due to entrained mass when using Avaframe Entrainment Model +# Used to determine speed loss via energy loss due to entrained mass when using Avaframe 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 +# Used to determine the strength of the erodible snow cover when using TJEM Entrainment Model [Pa] entShearStrength = 500 #++++++++++++Entrainment Energy loss From 7fd50ce054776ff4622d854b982a14f25f2e49e4 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Thu, 29 Oct 2026 15:43:17 +0100 Subject: [PATCH 16/21] rename to samosAT in com1DFA.py --- avaframe/com1DFA/com1DFA.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index 343d2ea2e..7a4dd1ffe 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -2192,7 +2192,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si # derive entrainment type # turn entrainment model into integer entrModelsList = [ - "avaframe", + "samosAT", "tjem", ] entrModel = cfgGen["entrModel"].lower() From 84c3057042ada42e8b33f4b83aa7f322fec604e9 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Thu, 29 Oct 2026 15:55:02 +0100 Subject: [PATCH 17/21] rename correctly --- avaframe/com1DFA/com1DFA.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index 7a4dd1ffe..1d40cfd48 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -2192,7 +2192,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si # derive entrainment type # turn entrainment model into integer entrModelsList = [ - "samosAT", + "samosat", "tjem", ] entrModel = cfgGen["entrModel"].lower() From c76607153cc5eb7d91810e009ec2fe32c7c68fab Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 4 Nov 2026 09:48:38 +0100 Subject: [PATCH 18/21] include velMagMin in function call --- avaframe/com1DFA/DFAfunctionsCython.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 365d82ec2..480f7df45 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -385,7 +385,7 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int entrType, int 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, rhoEnt) + 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 @@ -554,7 +554,7 @@ cpdef (double, double) computeEntMassAndForce(double dt, double entrMassCell, cpdef (double, double, double) computeEntMassAndForceTjem(double dt, double entrMassCell, double areaPart, double uMag, double tau, double entShearStrength, - double rhoEnt): + double velMagMin): """ compute force component due to entrained mass Parameters From d3738a59f12bf69d2962e48555f6910d97760e05 Mon Sep 17 00:00:00 2001 From: Felix Oesterle <6945681+fso42@users.noreply.github.com> Date: Thu, 5 Nov 2026 09:24:46 +0100 Subject: [PATCH 19/21] fix(com1DFA): limit entrained mass to available snow cover --- avaframe/com1DFA/DFAfunctionsCython.pyx | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 480f7df45..edb6d4259 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -571,8 +571,8 @@ cpdef (double, double, double) computeEntMassAndForceTjem(double dt, double entr bottom shear stress entShearStrength: float shear strength of the entrainable snow mass - rhoEnt: float - entrainement density + velMagMin : float + minimum velocity magnitude, avoids division by 0 for resting particles Returns ------- @@ -580,6 +580,8 @@ cpdef (double, double, double) computeEntMassAndForceTjem(double dt, double entr 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 @@ -595,6 +597,8 @@ cpdef (double, double, double) computeEntMassAndForceTjem(double dt, double entr 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) From 02df759ac69d02037485554aba28405b3d331351 Mon Sep 17 00:00:00 2001 From: Felix Oesterle <6945681+fso42@users.noreply.github.com> Date: Thu, 5 Nov 2026 12:48:15 +0100 Subject: [PATCH 20/21] fix(com1DFA): set default entrainment model to `samosAT` --- avaframe/com1DFA/com1DFACfg.ini | 4 +- avaframe/tests/test_DFAfunctionsCython.py | 58 +++++++++++++++++++++-- 2 files changed, 56 insertions(+), 6 deletions(-) diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index 179ec1e08..6144aa715 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -480,10 +480,10 @@ forestThMax = 10. ### Available Models: #samosAT #tjem -entrModel = avaframe +entrModel = samosAT #++++++++++++ Entrainment Erosion Energy -# Used to determine speed loss via energy loss due to entrained mass when using Avaframe Entrainment Model [m²/s²] +# 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 diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py index 123ede57c..ee2fd9b97 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) From cf104f153bcf21c3e4154fb913b26f0e14fc2811 Mon Sep 17 00:00:00 2001 From: Felix Oesterle <6945681+fso42@users.noreply.github.com> Date: Thu, 5 Nov 2026 13:55:48 +0100 Subject: [PATCH 21/21] test(DFAfunctionsCython): use tolerance for friction force assertion --- avaframe/tests/test_DFAfunctionsCython.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py index ee2fd9b97..c47d4e0fd 100644 --- a/avaframe/tests/test_DFAfunctionsCython.py +++ b/avaframe/tests/test_DFAfunctionsCython.py @@ -1312,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():