From 8d5de304b9a7e938041478f147f8b71bc6eda8c5 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Tue, 22 Sep 2026 15:24:30 +0200 Subject: [PATCH 1/8] 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 d57746d1d96ac6cbc01d5502e8de81254f5e4c3f Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Tue, 22 Sep 2026 15:34:20 +0200 Subject: [PATCH 2/8] 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 1bcf0af8d6e2052cc44c9e63335b49342c589138 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Tue, 22 Sep 2026 15:35:48 +0200 Subject: [PATCH 3/8] 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 54664a2455d4bed4585ebcd229d57956647f6e09 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Tue, 22 Sep 2026 15:39:47 +0200 Subject: [PATCH 4/8] 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 2873a91d8a9f29a517c2bcb01effbfe4d1a63f1c Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 23 Sep 2026 09:35:48 +0200 Subject: [PATCH 5/8] 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 1e9ab762efb7e268db227e2c50c40ea4ca8b80ba Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 23 Sep 2026 09:40:00 +0200 Subject: [PATCH 6/8] 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 7077a5d9488e9fdce08aee415321159b3d70ab41 Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Wed, 23 Sep 2026 15:26:21 +0200 Subject: [PATCH 7/8] 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 167255c75408288fd9b2892e5325de9a8288069f Mon Sep 17 00:00:00 2001 From: Michael Pallua Date: Thu, 24 Sep 2026 11:01:26 +0200 Subject: [PATCH 8/8] 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