Theory paper code

avaframe/ana1Tests/damBreak.py → avaframe/ana1Tests/damBreakTest.py

@@ -21,7 +21,7 @@
 log = logging.getLogger(__name__)
 
 
-def damBreakSol(avaDir, cfg, cfgC, outDirTest):
+def damBreakSol(avaDir, cfg, damBreakCfg, com1DFACfg, outDirTest):
     """ Compute analytical Flow thickness and velocity for dam break test case
 
         for a granular flow over a dry rough sloping bed with the Savage Hutter model
@@ -32,8 +32,10 @@ def damBreakSol(avaDir, cfg, cfgC, outDirTest):
             path to avalanche directory
         cfg: configParser object
             main avaframe configuration - here showPlot flag used
-        cfgC: configParser object
-            configuration setting for avalanche simulation including DAMBREAK section
+        damBreakCfg: configParser object
+            dam break configuration
+        com1DFACfg: configParser object
+            com1DFA configuration
         outDirTest: pathlib path
             path to output directory
 
@@ -58,21 +60,22 @@ def damBreakSol(avaDir, cfg, cfgC, outDirTest):
 
     # Set Parameters
     # Coordinate system chosen in the direction of the inclined plane
-    g = cfgC['GENERAL'].getfloat('gravAcc')       # acceleration due to gravity [ms-2]
-    phiDeg = cfgC['DAMBREAK'].getfloat('phi')       # slope angle [°]
-    deltaDeg = cfgC['DAMBREAK'].getfloat('delta')       # friction angle [°]
+    cfgGen = com1DFACfg['GENERAL']
+    g = cfgGen.getfloat('gravAcc')       # acceleration due to gravity [ms-2]
+    phiDeg = damBreakCfg.getfloat('phi')       # slope angle [°]
+    deltaDeg = damBreakCfg.getfloat('delta')       # friction angle [°]
     phi = np.radians(phiDeg)                          # slope angle [rad]
     delta = np.radians(deltaDeg)                        # bed friction angle [rad]
     gz = g * np.cos(phi)                          # projection of g perpendicular to the inclined plane
     m0 = gz * (np.tan(phi) - np.tan(delta))       # constant x-acceleration resulting from gravity and friction force
-    hL = cfgC['GENERAL'].getfloat('relTh')        # initial height [m] in Riemann problem in state 1 (x<0), hR (x>0)=0
+    hL = cfgGen.getfloat('relTh')        # initial height [m] in Riemann problem in state 1 (x<0), hR (x>0)=0
     cL = np.sqrt(gz * hL)
-    x0R = cfgC['DAMBREAK'].getfloat('xBack') / np.cos(phi)
-    tSave = cfgC['DAMBREAK'].getfloat('tSave')
+    x0R = damBreakCfg.getfloat('xBack') / np.cos(phi)
+    tSave = damBreakCfg.getfloat('tSave')
     # Define time [0-1] seconds and space [-2,2] meters domains multiplied times 100
-    tAna = np.arange(0, cfgC['DAMBREAK'].getfloat('tEnd'), cfgC['DAMBREAK'].getfloat('dt'))
-    x = np.arange(cfgC['DAMBREAK'].getfloat('xStart') / np.cos(phi) , cfgC['DAMBREAK'].getfloat('xEnd') / np.cos(phi),
-                  cfgC['DAMBREAK'].getfloat('dx'))
+    tAna = np.arange(0, damBreakCfg.getfloat('tEnd'), damBreakCfg.getfloat('dt'))
+    x = np.arange(damBreakCfg.getfloat('xStart') / np.cos(phi), damBreakCfg.getfloat('xEnd') / np.cos(phi),
+                  damBreakCfg.getfloat('dx'))
     nt = len(tAna)
     nx = len(x)
     # Initialise Flow thickness solution and velocity
@@ -227,13 +230,19 @@ def analyzeResults(avalancheDir, fieldsList, timeList, solDam, fieldHeader, cfg,
     xAna = solDam['xAna']
     hAna = solDam['hAna']
     uAna = solDam['uAna']
-    # get plots limits
-    limits = outAna1Plots.getPlotLimits(cfgDam, fieldsList, fieldHeader)
     hErrorL2Array = np.zeros((len(fieldsList)))
     vhErrorL2Array = np.zeros((len(fieldsList)))
     hErrorLMaxArray = np.zeros((len(fieldsList)))
     vhErrorLMaxArray = np.zeros((len(fieldsList)))
     count = 0
+    tSave = cfgDam.getfloat('tSave')
+    indT = min(np.searchsorted(timeList, tSave), min(len(timeList)-1, len(fieldsList)-1))
+    tSave = timeList[indT]
+    # get plots limits
+    if cfgDam.getboolean('plotSequence'):
+        limits = outAna1Plots.getPlotLimits(cfgDam, fieldsList, fieldHeader)
+    else:
+        limits = outAna1Plots.getPlotLimits(cfgDam, [fieldsList[indT]], fieldHeader)
     # run the comparison routine for each saved time step
     for t, field in zip(timeList, fieldsList):
         # get similartiy solution h, u at required time step
@@ -251,9 +260,9 @@ def analyzeResults(avalancheDir, fieldsList, timeList, solDam, fieldHeader, cfg,
 
         # make analytical results 2D
         nrows, ncols = np.shape(hNumPlus)
-        O = np.ones((nrows, ncols))
-        hDamPlus = O*hDamPlus
-        uDamPlus = O*uDamPlus
+        OnesRaster = np.ones((nrows, ncols))
+        hDamPlus = OnesRaster*hDamPlus
+        uDamPlus = OnesRaster*uDamPlus
 
         hEL2Plus, hL2RelPlus, hLmaxPlus, hLmaxRelPlus = anaTools.normL2Scal(hDamPlus, hNumPlus, cellSize,
                                                                             np.cos(phiRad))
@@ -273,15 +282,12 @@ def analyzeResults(avalancheDir, fieldsList, timeList, solDam, fieldHeader, cfg,
         log.debug("L2 %s error on the Flow thickness at t=%.2f s is : %.4f" % (title, t, hEL2Plus))
         log.debug("L2 %s error on the momentum at t=%.2f s is : %.4f" % (title, t, vhL2Plus))
         # Make all individual time step comparison plot
-        if cfgDam.getboolean('plotSequence'):
+        if cfgDam.getboolean('plotSequence') or t == tSave:
             outAna1Plots.plotComparisonDam(cfgDam, simHash, fieldsList[0], field, fieldHeader, solDam, t, limits,
                                            outDirTest)
         count = count + 1
 
-    tSave = cfgDam.getfloat('tSave')
-    indT = min(np.searchsorted(timeList, tSave), min(len(timeList)-1, len(fieldsList)-1))
-    tSave = timeList[indT]
-    if cfgDam.getboolean('plotErrorTime') and len(timeList)>1:
+    if cfgDam.getboolean('plotErrorTime') and len(timeList) > 1:
         # Create result plots
 
         outAna1Plots.plotErrorTime(timeList, hErrorL2Array, hErrorLMaxArray, vhErrorL2Array, vhErrorLMaxArray,

avaframe/data/avaDamBreak/Inputs/damBreak_com1DFACfg.ini → avaframe/ana1Tests/damBreakTestCfg.ini

@@ -1,15 +1,67 @@
-### Config File - This file contains the main settings for the simulation run for the dam break problem
-## This file is part of Avaframe.
+### Config File - This file contains the main settings for the dam break test module
+## Set your parameters
+# This file will be overridden by local_damBreakTestCfg.ini if it exists
+# So copy this file to local_damBreakTestCfg.ini, adjust your variables there
 
+[DAMBREAK]
+#sphKernel radius for the com1DFA simulations
+sphKernelRadius = 6|5|4|3
+
+# slope angle [°]
+phi = 22
+# bed friction angle [°]
+delta = 21
+u0 = 0.
+# time end and resolution for the analytic solution
+tEnd = 30
+dt = 0.1
+# space resolution
+dx = 0.5
+# domain extend for error computation
+# in x direction
+xStart = -200
+xEnd = 220
+# in y direction
+yStart = -50
+yEnd = 50
+# start x position of the dam
+xBack = -120
+# xFront = 0
 
-[GENERAL]
+# set time step of analysis
+tSave = 15
+# use only the component of the velocity/momentum in the flow direction (vector going down the inclined plane in x dir)
+projectVelocity = False
+
+#++++Plotting parameters++++++
+# list of parameters to display in the summary plot (list of parameters separated by |)
+paramInfo = sphKernelRadius|aPPK|nPPK0|cMax|nPart
+# list of parameters to display in the convergence plot (list of parameters separated by |)
+paramConvInfo = sphKR0|nPPK0|cMax
+# plotting flags
+# only analyze and plot the tSave time step
+onlyLast = False
+# plot error function of time for each simulation
+plotErrorTime = False
+# plot individual figure for the h, hu and u error for each saved time step
+plotSequence = False
+# use relative difference
+relativError = True
+# when plotting, the domain extent is scaleCoef times bigger than the data extent
+scaleCoef = 1.05
+
+
+
+[com1DFA_override]
+# use default com1DFA config as base configuration (True) and override following parameters
+# if False and local is available use local
+defaultConfig = True
 #++++++++++++++++ Simulation type
 # list of simulations that shall be performed (null, ent, res, entres, available (use all available input data))
 simTypeList = null
-
 #+++++++++++++ Output++++++++++++
 # desired result Parameters (ppr, pfd, pfv, FD, FV, P, particles) - separated by |
-resType = FD|FV|Vx|Vy|Vz
+resType = FT|FV|Vx|Vy|Vz
 # saving time step, i.e.  time in seconds (first and last time step are always saved)
 # option 1: give an interval with start:interval in seconds (tStep = 0:5 - this will save desired results every 5 seconds for the full simulation)
 # option 2: explicitly list all desired time steps (closest to actual computational time step) separated by | (example tSteps = 1|50.2|100)
@@ -24,29 +76,30 @@ iniStep = True
 maxIterations = 30
 # buffer zone factor multiplied with sphKernelRadius
 bufferZoneFactor = 4
-fdOptionIni = False
 
 #+++++++++SNOW properties
 # True if release thickness should be read from shapefile file; if False - relTh read from ini file
 relThFromShp = False
 relTh = 1
 
 #++++++++++++Time stepping parameters
-tEnd = 5
+# End time [s]
+tEnd = 20
 # to use a variable time step (time step depends on kernel radius)
 sphKernelRadiusTimeStepping = True
-# courant number if option cflTimeStepping is chosen.
-# Upper time step limit coefficient if option sphKernelRadiusTimeStepping is chosen.
+# courant number
 cMax = 0.01
 
 #+++++++++++++SPH parameters
 # SPH gradient option
 # 0) No pressure gradients computed
 # 1) SamosAT style (no reprojecion on the surface, dz = 0 and gz is used)
-# 2) SamosAT done in the cartesian coord system (reprojecion on the surface, dz != 0 and g3 is used)
-# 3) SamosAT but done in the local coord system (will hopefully allow us to add the earth pressure coef)
+# 2) SamostAt but done in the cartesian coord system (will hopefully allow us to add the earth pressure coef)
+# 3) SamostAt but done in the local coord system (will hopefully allow us to add the earth pressure coef)
 # and this time reprojecion on the surface, dz not 0 and g3 is used
 sphOption = 2
+# sph kernel smoothing length [m]
+sphKernelRadius = 3
 # Choice of artificial viscosity
 # 0) No artificial viscosity
 # 1) SAMOS artificial viscosity
@@ -57,25 +110,25 @@ 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
+# number of particles per kernel radius (if MPPKR is used)
 # is computed with: nPPK = nPPK0 * (sphKR/sphKR0)^aPPK
 # where sphKR is the sphKernelRadius specified further up
 # reference kernel radius
 sphKR0 = 5
 # reference number of particles per kernel radius
-nPPK0 = 10|20|30
+nPPK0 = 15
 # variation of nppK exponent
-aPPK = -0.5|-1|-1.5|-2
+aPPK = 0|-0.5|-1|-2|-3
 
 # if splitOption=0
 # threshold for splitting particles, split if: mPart > massPerPart x thresholdMassSplit
-thresholdMassSplit = 10
+thresholdMassSplit = 5
 
 
+#+++++++++++++Mesh and interpolation
 # remesh the input DEM
 # expected mesh size [m]
 meshCellSize = 3
-# sph kernel smoothing length [m]
-sphKernelRadius = 3
 
 #++++++++++++Friction model
 # add the friction using an explicit formulation (1)
@@ -85,46 +138,3 @@ explicitFriction = 1
 frictModel = Coulomb
 #+++++++++++++SamosAt friction model
 mu = 0.3838640350354158
-
-
-[DAMBREAK]
-# slope angle [°]
-phi = 22
-# bed friction angle [°]
-delta = 21
-u0 = 0.
-# time end and resolution for the analytic solution
-tEnd = 30
-dt = 0.1
-# space resolution
-dx = 0.5
-# domain extend for error computation
-# in x direction
-xStart = -200
-xEnd = 220
-# in y direction
-yStart = -50
-yEnd = 50
-# start x position of the dam
-xBack = -120
-# xFront = 0
-
-# set time step of analysis
-tSave = 5
-# use only the component of the velocity/momentum in the flow direction (vector going down the inclined plane in x dir)
-projectVelocity = False
-
-#++++Plotting parameters++++++
-# list of parameters to display in the summary plot (list of parameters separated by |)
-paramInfo = sphKernelRadius|aPPK|nPPK0|nPart
-# plotting flags
-# only analyze and plot the tSave time step
-onlyLast = False
-# plot error function of time for each simulation
-plotErrorTime = False
-# plot individual figure for the h, hu and u error for each saved time step
-plotSequence = False
-# use relative difference
-relativError = True
-# when plotting, the domain extent is scaleCoef times bigger than the data extent
-scaleCoef = 1.05

avaframe/ana1Tests/energyLineTest.py

@@ -135,11 +135,12 @@ def generateCom1DFAEnergyPlot(avalancheDir, energyLineTestCfg, com1DFACfg, avaPr
     zMin = zLim[0]
     ax2.vlines(x=avaProfileMass['s'][-1], ymin=zMin, ymax=avaProfileMass['z'][-1],
                color='r', linestyle='--')
-    ax2.vlines(x=sIntersection, color='b', ymin=zMin, ymax=zIntersection, linestyle='--')
+    if ~np.isnan(sIntersection):
+        ax2.vlines(x=sIntersection, color='b', ymin=zMin, ymax=zIntersection, linestyle='--')
+        ax2.hlines(y=zIntersection, color='b', xmin=0, xmax=sIntersection, linestyle='--')
 
     ax2.hlines(y=avaProfileMass['z'][-1], xmin=0, xmax=avaProfileMass['s'][-1],
                color='r', linestyle='--')
-    ax2.hlines(y=zIntersection, color='b', xmin=0, xmax=sIntersection, linestyle='--')
 
     ax2.set_xlabel('s [m]')
     ax2.set_ylabel('z [m]')
@@ -183,8 +184,8 @@ def generateCom1DFAEnergyPlot(avalancheDir, energyLineTestCfg, com1DFACfg, avaPr
     rmseVelocityElevation = resultEnergyTest['rmseVelocityElevation']
     errorS = abs(runOutSError)
     errorZ = abs(runOutZError)
-    sMin = min(avaProfileMass['s'][-1], sIntersection) - max(errorS, 0)
-    sMax = max(avaProfileMass['s'][-1], sIntersection) + max(errorS, 0)
+    sMin = np.nanmin(np.array([avaProfileMass['s'][-1], sIntersection])) - np.nanmax(np.array([errorS, 0]))
+    sMax = np.nanmax(np.array([avaProfileMass['s'][-1], sIntersection])) + np.nanmax(np.array([errorS, 0]))
     zMin = avaProfileMass['z'][-1] + min(slopeExt*(sMax - avaProfileMass['s'][-1]), 0 - 2*errorZ)-z0
     zMax = avaProfileMass['z'][0] - sMin*np.tan(min(runOutAngleRad, alphaRad))-z0
     if avaProfileMass['z'][-1] == zIntersection:
@@ -193,11 +194,12 @@ def generateCom1DFAEnergyPlot(avalancheDir, energyLineTestCfg, com1DFACfg, avaPr
         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='--')
+    if ~np.isnan(sIntersection):
+        ax3.vlines(x=sIntersection, color='b', ymin=zMin, ymax=zIntersection-z0, linestyle='--')
+        ax3.hlines(y=zIntersection-z0, color='b', xmin=sMin, xmax=sIntersection, linestyle='--')
 
     ax3.hlines(y=avaProfileMass['z'][-1]-z0, xmin=sMin, xmax=avaProfileMass['s'][-1],
                color='r', linestyle='--')
-    ax3.hlines(y=zIntersection-z0, color='b', xmin=sMin, xmax=sIntersection, linestyle='--')
 
     ax3.set_xlabel('s [m]')
     ax3.set_ylabel('$\Delta z$ [m]')
@@ -292,7 +294,7 @@ def getAlphaProfileIntersection(energyLineTestCfg, avaProfileMass, mu, csz):
     # find the intersection segment
     idx = np.argwhere(np.diff(np.sign(z - alphaLine))).flatten()
     idx = idx[-1]
-    coefExt = 0
+    coefExt = 1
     # did not find the intersection, look further
     while s[idx] == 0 and coefExt < 4:
         s = np.append(avaProfileMass['s'], (1+coefExt)*avaProfileMass['s'][-1])
@@ -310,8 +312,14 @@ def getAlphaProfileIntersection(energyLineTestCfg, avaProfileMass, mu, csz):
     zP1 = z[idx+1]
     zA0 = alphaLine[idx]
     zA1 = alphaLine[idx+1]
-    sIntersection = s0 + (s1-s0)*(zA0-zP0)/((zP1-zP0)-(zA1-zA0))
-    zIntersection = zP0 + (sIntersection-s0) * (zP1-zP0) / (s1-s0)
+    if s0 == 0:
+        sIntersection = np.nan
+        zIntersection = np.nan
+        s[-1] = avaProfileMass['s'][-1]
+        z[-1] = avaProfileMass['z'][-1]
+    else:
+        sIntersection = s0 + (s1-s0)*(zA0-zP0)/((zP1-zP0)-(zA1-zA0))
+        zIntersection = zP0 + (sIntersection-s0) * (zP1-zP0) / (s1-s0)
     if coefExt > 0:
         s[-1] = sIntersection
         z[-1] = zIntersection
@@ -363,7 +371,10 @@ def getEnergyInfo(avaProfileMass, g, mu, sIntersection, zIntersection, runOutAng
     GK = sIntersection * mu
     z_end = avaProfileMass['z'][0] - GK
     zEneTarget = avaProfileMass['z'][0] - avaProfileMass['s'] * mu
-    sGeomL = [avaProfileMass['s'][0], sIntersection]
+    if ~np.isnan(sIntersection):
+        sGeomL = [avaProfileMass['s'][0], sIntersection]
+    else:
+        sGeomL = [avaProfileMass['s'][0], avaProfileMass['s'][-1]]
     zGeomL = [avaProfileMass['z'][0], z_end]
     # compute errors
     # rmse of the energy height

avaframe/ana1Tests/energyLineTestCfg.ini

@@ -93,5 +93,20 @@ frictModel = Coulomb
 #+++++++++++++General Friction parameters
 # tan of bed friction angle used for: samosAT, Coulomb, Voellmy
 mu = 0.4
+#+++++++++++++Voellmy friction model
+xsi = 10000.
 
-dam = True
+#+++++++++++++++++ Dam Parameters
+# the dam foot print is given in Intputs/DAM as a shape file line
+# the dam is located on the left side of the dam line
+# (when one travels from the first point to the last point of the shapefile dam line)
+# use dam?
+dam = False
+# behavior of the particles when they hit the dam (are they fully reflected: restitutionCoefficient = 1 or only reflected
+# along the tangent component of the dam: restitutionCoefficient = 0 or something in between: restitutionCoefficient = 0-1)
+restitutionCoefficient = 1
+# number of iterations allowed for the interactions between a particle and the dam during one time step:
+# 1: particle can only bounce once (may cross the dam at another location then)
+nIterDam = 1
+# 1 to activate energy loss due to snow flowing over dam
+dissDam = 1