add ini files, update runScripts

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/ana1Tests/damBreakTestCfg.ini

@@ -0,0 +1,140 @@
+### 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
+
+# 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|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 = 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)
+# NOTE: initial and last time step are always saved!
+tSteps = 0:1
+
+#++++++++++++++++ particle Initialisation +++++++++
+# particle initialisation step - set iniStep to True to activate redistribution of particles to reduce SPH force
+# this is in a development stage - hence parameters are set for development and will be adjusted after extensive testing
+iniStep = True
+# max number of iterations - high number might cause significant increase in computational time
+maxIterations = 30
+# buffer zone factor multiplied with sphKernelRadius
+bufferZoneFactor = 4
+
+#+++++++++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
+# End time [s]
+tEnd = 20
+# to use a variable time step (time step depends on kernel radius)
+sphKernelRadiusTimeStepping = True
+# 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) 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
+# 2) Ata artificial viscosity
+viscOption = 0
+
+#++++++++++++++++ Particles
+# 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 = 15
+# variation of nppK exponent
+aPPK = 0|-0.5|-1|-2|-3
+
+# if splitOption=0
+# threshold for splitting particles, split if: mPart > massPerPart x thresholdMassSplit
+thresholdMassSplit = 5
+
+
+#+++++++++++++Mesh and interpolation
+# remesh the input DEM
+# expected mesh size [m]
+meshCellSize = 3
+
+#++++++++++++Friction model
+# add the friction using an explicit formulation (1)
+# 0 use an implicit method
+explicitFriction = 1
+# friction type (samosAT, Coulomb)
+frictModel = Coulomb
+#+++++++++++++SamosAt friction model
+mu = 0.3838640350354158

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]')
@@ -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,3 +93,5 @@ frictModel = Coulomb
 #+++++++++++++General Friction parameters
 # tan of bed friction angle used for: samosAT, Coulomb, Voellmy
 mu = 0.4
+#+++++++++++++Voellmy friction model
+xsi = 10000.