Prepare com1DFA code according to theory paper

avaframe/ana1Tests/energyLineTest.py

@@ -1,6 +1,5 @@
 """
 Energy line test
-
 This module runs a DFA simulation extracts the center of mass path
 and compares it to the analytic geometric/alpha line solution
 """
@@ -29,7 +28,6 @@
 
 def mainEnergyLineTest(avalancheDir, energyLineTestCfg, com1DFACfg, simName, dem):
     """This is the core function of the energyLineTest module
-
     This module extracts the center of mass path from a DFA simulation
     and compares it to he analytic geometric/alpha line solution
     """
@@ -49,7 +47,6 @@ def mainEnergyLineTest(avalancheDir, energyLineTestCfg, com1DFACfg, simName, dem
 
 def generateCom1DFAEnergyPlot(avalancheDir, energyLineTestCfg, com1DFACfg, avaProfileMass, dem, fieldsList, simName):
     """ Make energy test analysis and plot results
-
     Parameters
     -----------
     avalancheDir: pathlib
@@ -130,6 +127,9 @@ def generateCom1DFAEnergyPlot(avalancheDir, energyLineTestCfg, com1DFACfg, avaPr
 
     # add alpha line
     ax2.plot(sGeomL, zGeomL, 'b-', label=r'$\alpha$ line (%.2f°)' % alphaDeg)
+    if energyLineTestCfg['energyLineTest'].getboolean('shiftedAlphaLine'):
+        ax2.plot(avaProfileMass['s'], avaProfileMass['z'][-1] - (avaProfileMass['s']-avaProfileMass['s'][-1]) * mu,
+                 'b-.', label=r'shifted $\alpha$ line (%.2f°)' % alphaDeg)
     zLim = ax2.get_ylim()
     sLim = ax2.get_xlim()
     zMin = zLim[0]
@@ -171,6 +171,9 @@ def generateCom1DFAEnergyPlot(avalancheDir, energyLineTestCfg, com1DFACfg, avaPr
 
     # add alpha line
     ax3.plot(sGeomL, zGeomL-z0, 'b-', label=r'$\alpha$ line (%.2f°)' % alphaDeg)
+    if energyLineTestCfg['energyLineTest'].getboolean('shiftedAlphaLine'):
+        ax3.plot(avaProfileMass['s'], avaProfileMass['z'][-1]-z0 - (avaProfileMass['s']-avaProfileMass['s'][-1]) * mu,
+                 'b-.', label=r'shifted $\alpha$ line (%.2f°)' % alphaDeg)
 
     # add horizontal line at the final mass averaged position
     # compute plot limits
@@ -186,6 +189,8 @@ def generateCom1DFAEnergyPlot(avalancheDir, energyLineTestCfg, com1DFACfg, avaPr
     zMax = avaProfileMass['z'][0] - sMin*np.tan(min(runOutAngleRad, alphaRad))-z0
     if avaProfileMass['z'][-1] == zIntersection:
         zMin = zMin - (zMax-zMin)*0.1
+    if errorZ < 1e-3:
+        zMin = zMin - (zMax-zMin)*0.1
     ax3.vlines(x=avaProfileMass['s'][-1], ymin=zMin, ymax=avaProfileMass['z'][-1]-z0,
                color='r', linestyle='--')
     ax3.vlines(x=sIntersection, color='b', ymin=zMin, ymax=zIntersection-z0, linestyle='--')
@@ -222,7 +227,6 @@ def generateCom1DFAEnergyPlot(avalancheDir, energyLineTestCfg, com1DFACfg, avaPr
 
 def getRunOutAngle(avaProfileMass, indStart=0, indEnd=-1):
     """Compute the center of mass runout angle
-
     Parameters
     -----------
     avaProfileMass: dict
@@ -232,7 +236,6 @@ def getRunOutAngle(avaProfileMass, indStart=0, indEnd=-1):
         0 by default - so full mass averaged path from top
     indEnd: int
         index of the start of the mass averaged pass (to discard the bottom extension). -1 by default
-
     Returns
     --------
     runOutAngleRad: float
@@ -252,10 +255,8 @@ def getRunOutAngle(avaProfileMass, indStart=0, indEnd=-1):
 def getAlphaProfileIntersection(energyLineTestCfg, avaProfileMass, mu, csz):
     """Extend the  profile path and compute the intersection
     between the theoretical energy line and the path profile
-
     The profile is extended by a line. The line slope is computed
     from the slope of the regression on the las points of the profile
-
     Parameters
     -----------
     energyLineTestCfg: configParser
@@ -266,7 +267,6 @@ def getAlphaProfileIntersection(energyLineTestCfg, avaProfileMass, mu, csz):
         friction coefficient
     csz: float
         dem cell size
-
     Returns
     --------
     slopeExt: float
@@ -321,7 +321,6 @@ def getAlphaProfileIntersection(energyLineTestCfg, avaProfileMass, mu, csz):
 
 def getEnergyInfo(avaProfileMass, g, mu, sIntersection, zIntersection, runOutAngleDeg, alphaDeg):
     """Compute energy dots and errors
-
     Parameters
     -----------
     avaProfileMass: dict
@@ -340,7 +339,6 @@ def getEnergyInfo(avaProfileMass, g, mu, sIntersection, zIntersection, runOutAng
         center of mass runout angle in radians
     runOutAngleDeg: float
         center of mass runout angle in degrees
-
     Returns
     --------
     zEne: numpy 1D array

avaframe/ana1Tests/energyLineTestCfg.ini

@@ -17,6 +17,9 @@ plotScor = False
 # if False, FT|FV|FM need to be saved in resType
 pathFromPart = True
 
+# for plotting the results, add the shifted alpha line
+shiftedAlphaLine = False
+
 [com1DFA_override]
 # use default com1DFA config as base configuration (True) and override following parameters
 # if False and local is available use local
@@ -44,6 +47,8 @@ relTh = 1
 #++++++++++++Time stepping parameters
 # to use a variable time step (time step depends on kernel radius)
 sphKernelRadiusTimeStepping = True
+# Upper time step limit coefficient if option sphKernelRadiusTimeStepping is chosen.
+cMax = 0.02
 # stopCriterion (stops when massFlowing<0.01*peakMassFlowing or ke<0.01*pke)
 stopCrit = 0.0
 
@@ -65,15 +70,19 @@ viscOption = 0
 # number of particles defined by: MPPDIR= mass per particle direct, MPPDH= mass per particles through release thickness,
 # MPPKR= mass per particles through number of particles per kernel radius
 massPerParticleDeterminationMethod = MPPKR
-nPPK0 = 1
 
 #+++++++++++++Flow model parameters+++++++++
 # subgridMixingFactor
 subgridMixingFactor = 100.
-# curvature acceleration coefficient
-# take curvature term into account in the gravity acceleration term
+# take curvature term into account in the gravity acceleration term for computing the friction force
+# (and gradient if curvAccInGradient=1)
+# 0 if deactivated, 1 if activated
+curvAccInFriction = 0
+# take curvature term into account in the tangential momentum equation
 # 0 if deactivated, 1 if activated
-curvAcceleration = 0
+curvAccInTangent = 0
+# if curvAccInFriction=1, take curvature term into account in the pressure gradient (if curvAccInGradient=1)
+curvAccInGradient = 0
 
 #++++++++++++Friction model
 # add the friction using an explicit formulation (1)