[com1]: Add input hydrograph (initial)

avaframe/com1DFA/DFAfunctionsCython.pyx

@@ -194,7 +194,6 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
       indCellY = indYDEM[k]
       # deduce area
       areaPart = m / (h * rho)
-
       # get cell and weights
       Lx0, Ly0, iCell, w[0], w[1], w[2], w[3] = DFAtlsC.getCellAndWeights(x, y, ncols, nrows, csz, interpOption)
 
@@ -2293,3 +2292,75 @@ def computeIniMovement(cfg, particles, dem, dT, fields):
   particles['m'] = np.asarray(mass)
 
   return particles, force
+
+
+def updateInitialVelocity(cfg, particles, dem, float velocityMag):
+  """
+  update particles' velocity in direction of the steepest descent
+
+  Parameters
+  ------------
+  cfg: configparser
+    configuration for DFA simulation
+  particles : dict
+    particles dictionary at t
+  dem : dict
+    dictionary with dem information
+  velocityMag: float
+    velocity of the particles
+
+  Returns
+  ------------
+  particles : dict
+    particles dictionary at t with updated velocity
+  """
+
+  cdef double[:] xArray = particles['x']
+  cdef double[:] yArray = particles['y']
+  cdef double[:] uxArray = particles['ux']
+  cdef double[:] uyArray = particles['uy']
+  cdef double[:] uzArray = particles['uz']
+  cdef double[:] velocityMagArray = particles['velocityMag']
+  cdef int nPart = particles['nPart']
+  cdef int nrows = dem['header']['nrows']
+  cdef int ncols = dem['header']['ncols']
+  cdef double[:, :] nxArray = dem['Nx']
+  cdef double[:, :] nyArray = dem['Ny']
+  cdef double[:, :] nzArray = dem['Nz']
+  cdef int interpOption = cfg.getint('interpOption')
+  cdef double csz = dem['header']['cellsize']
+
+  cdef double x, y, z
+  cdef double nx, ny, nz
+  cdef double tangentialTopoX, tangentialTopoY, tangentialTopoZ, tangentialTopoXNorm, tangentialTopoYNorm, tangentialTopoZNorm
+  cdef int Lx0, Ly0, iCell
+  cdef double w[4]
+  cdef int k
+
+  for k in range(nPart):
+    x = xArray[k]
+    y = yArray[k]
+
+    # get cell and weights
+    Lx0, Ly0, iCell, w[0], w[1], w[2], w[3] = DFAtlsC.getCellAndWeights(x, y, ncols, nrows, csz, interpOption)
+
+    # get normal at the particle location
+    nx, ny, nz = DFAtlsC.getVector(Lx0, Ly0, w[0], w[1], w[2], w[3], nxArray, nyArray, nzArray)
+    nx, ny, nz = DFAtlsC.normalize(nx, ny, nz)
+    # get vector in the tangent plane that points in the direction gravity would pull a particle along the surface
+    tangentialTopoX = nx * nz
+    tangentialTopoY = ny * nz
+    tangentialTopoZ = nz * nz -1
+    tangentialTopoXNorm, tangentialTopoYNorm, tangentialTopoZNorm = DFAtlsC.normalize(tangentialTopoX, tangentialTopoY, tangentialTopoZ)
+
+    uxArray[k] = tangentialTopoXNorm * velocityMag
+    uyArray[k] = tangentialTopoYNorm * velocityMag
+    uzArray[k] = tangentialTopoZNorm * velocityMag
+
+    velocityMagArray[k] = DFAtlsC.norm(uxArray[k], uyArray[k], uzArray[k])
+  particles['ux'] = np.asarray(uxArray)
+  particles['uy'] = np.asarray(uyArray)
+  particles['uz'] = np.asarray(uzArray)
+  particles['velocityMag'] = np.asarray(velocityMagArray)
+
+  return particles

avaframe/com1DFA/com1DFA.py

@@ -53,6 +53,7 @@
 from avaframe.com1DFA import checkCfg
 from avaframe.ana5Utils import distanceTimeAnalysis as dtAna
 import avaframe.out3Plot.outDistanceTimeAnalysis as dtAnaPlots
+import avaframe.com1DFA.debrisFunctions as debF
 import threading
 
 #######################################
@@ -562,6 +563,8 @@ def prepareInputData(inputSimFiles, cfg):
         - secondaryRelFile : str, path to secondaryRelease file
         - entFiles : str, path to entrainment file
         - resFile : str, path to resistance file
+        - hydrographFile: str, path to hydrograph polygon file
+        - hydrographCsv: str, path to hydrograph values (csv-)file
         - entResInfo : flag dict
         flag if Yes entrainment and/or resistance areas found and used for simulation
         flag True if a Secondary Release file found and activated
@@ -585,6 +588,7 @@ def prepareInputData(inputSimFiles, cfg):
         - resLine : dict, resistance line dictionary
         - entrainmentArea : str, entrainment file name
         - resistanceArea : str, resistance file name
+        - hydrographAreaLine: dict, hydrograph line dictionary
         - entResInfo : flag dict
         flag if Yes entrainment and/or resistance areas found and used for simulation
         flag True if a Secondary Release file found and activated
@@ -741,6 +745,11 @@ def prepareInputData(inputSimFiles, cfg):
     else:
         damLine = None
 
+    if cfg["GENERAL"].getboolean("hydrograph"):
+        hydrLine = debF.preparehydrographAreaLine(inputSimFiles, demOri, cfg)
+    else:
+        hydrLine = None
+
     inputSimLines = {
         "releaseLine": releaseLine,
         "secondaryReleaseLine": secondaryReleaseLine,
@@ -756,6 +765,7 @@ def prepareInputData(inputSimFiles, cfg):
         "xiFile": inputSimFiles["xiFile"],
         "kFile": inputSimFiles["kFile"],
         "tauCFile": inputSimFiles["tauCFile"],
+        "hydrographAreaLine": hydrLine,
     }
 
     return demOri, inputSimLines
@@ -1613,6 +1623,8 @@ def getRelThFromPart(cfg, releaseLine, relThField, thName):
 
     if len(relThField) != 0:
         relThForPart = np.amax(relThField)
+    elif releaseLine["type"] == "Hydrograph":
+        relThForPart = releaseLine["thickness"]
     elif cfg.getboolean("%sThFromFile" % thName):
         relThForPart = np.amax(np.asarray(releaseLine["thickness"], dtype=float))
     else:
@@ -1654,42 +1666,42 @@ def initializeFields(cfg, dem, particles, releaseLine):
     nrows = header["nrows"]
     # initialize fields
     fields = {}
-    fields["pfv"] = np.zeros((nrows, ncols))
-    fields["pft"] = np.zeros((nrows, ncols))
-    fields["FV"] = np.zeros((nrows, ncols))
-    fields["FT"] = np.zeros((nrows, ncols))
-    fields["FM"] = np.zeros((nrows, ncols))
-    fields["Vx"] = np.zeros((nrows, ncols))
-    fields["Vy"] = np.zeros((nrows, ncols))
-    fields["Vz"] = np.zeros((nrows, ncols))
-    fields["dmDet"] = np.zeros((nrows, ncols))
-    fields["FTStop"] = np.zeros((nrows, ncols))
-    fields["FTDet"] = np.zeros((nrows, ncols))
-    fields["FTEnt"] = np.zeros((nrows, ncols))
-    fields["sfcChange"] = np.zeros((nrows, ncols))
-    fields["sfcChangeTotal"] = np.zeros((nrows, ncols))
-    fields["demAdapted"] = np.zeros((nrows, ncols))
+    fields["pfv"] = np.zeros((nrows, ncols))  # peak flow velocity [m/s]
+    fields["pft"] = np.zeros((nrows, ncols))  # peal flow thickness [m]
+    fields["FV"] = np.zeros((nrows, ncols))  # flow velocity [m/s]
+    fields["FT"] = np.zeros((nrows, ncols))  # flow thickness [m]
+    fields["FM"] = np.zeros((nrows, ncols))  # flow mass [kg]
+    fields["Vx"] = np.zeros((nrows, ncols))  # velocity in x direction [m/s]
+    fields["Vy"] = np.zeros((nrows, ncols))  # velocity in y direction [m/s]
+    fields["Vz"] = np.zeros((nrows, ncols))  # velocity in z direction [m/s]
+    fields["dmDet"] = np.zeros((nrows, ncols))  # flowing mass change due to detrainment [kg]
+    fields["FTStop"] = np.zeros((nrows, ncols))  # flow thickness that is stopped [m]
+    fields["FTDet"] = np.zeros((nrows, ncols))  # flow thickness that is detrained [m]
+    fields["FTEnt"] = np.zeros((nrows, ncols))  # flow thickness that is entrained [m]
+    fields["sfcChange"] = np.zeros((nrows, ncols))  # depth that changes the surface [m]
+    fields["sfcChangeTotal"] = np.zeros((nrows, ncols))  # total depth that changed the surface [m]
+    fields["demAdapted"] = np.zeros((nrows, ncols))  # adapted topography [m]
     # for optional fields, initialize with dummys (minimum size array). The cython functions then need something
     # even if it is empty to run properly
     if ("TA" in resTypesLast) or ("pta" in resTypesLast):
-        fields["pta"] = np.zeros((nrows, ncols))
-        fields["TA"] = np.zeros((nrows, ncols))
+        fields["pta"] = np.zeros((nrows, ncols))  # peak travel angle [°]
+        fields["TA"] = np.zeros((nrows, ncols))  # travel angle [°]
         fields["computeTA"] = True
         log.debug("Computing Travel Angle")
     else:
         fields["pta"] = np.zeros((1, 1))
         fields["TA"] = np.zeros((1, 1))
         fields["computeTA"] = False
     if "pke" in resTypesLast:
-        fields["pke"] = np.zeros((nrows, ncols))
+        fields["pke"] = np.zeros((nrows, ncols))  # peak kinetic energy [kJ/m²]
         fields["computeKE"] = True
         log.debug("Computing Kinetic energy")
     else:
         fields["pke"] = np.zeros((1, 1))
         fields["computeKE"] = False
     if ("P" in resTypesLast) or ("ppr" in resTypesLast):
-        fields["P"] = np.zeros((nrows, ncols))
-        fields["ppr"] = np.zeros((nrows, ncols))
+        fields["P"] = np.zeros((nrows, ncols))  # pressure [kPa]
+        fields["ppr"] = np.zeros((nrows, ncols))  # peak pressure [kPa]
         fields["computeP"] = True
         log.debug("Computing Pressure")
     else:
@@ -2131,6 +2143,11 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
     while t <= tEnd * (1.0 + 1.0e-13) and particles["iterate"]:
         startTime = time.time()
         log.debug("Computing time step t = %f s, dt = %f s" % (t, dt))
+
+        if cfgGen.getboolean("hydrograph"):
+            particles, fields, zPartArray0 = debF.releaseHydrograph(
+                cfg, inputSimLines, particles, fields, dem, zPartArray0, t
+            )
         # Perform computations
         particles, fields, zPartArray0, tCPU, dem = computeEulerTimeStep(
             cfgGen,
@@ -3594,6 +3611,7 @@ def adaptDEM(dem, fields, cfg):
     """adapt topography in respect to erosion and deposition
 
     Parameters
+    ---------
     dem: dict
         dictionary with info on DEM data
     fields : dict

avaframe/com1DFA/com1DFACfg.ini

@@ -132,6 +132,21 @@ entThDistVariation =
 # entrainment thickness (only considered if ENT file is shapefile and entThFromFile=False) [m]
 entTh =
 
+#+++++++++++++Hydrograph
+# if hydrograph is True, add hydrograph, provide the hydrograph values in a csv-table in the HYDR folder
+hydrograph = False
+# when checking if an already existing particle is within a hydrograph polygon
+# (that would cause numerical instabilities and rise an error), one can decide to add a buffer
+# around the polygon (0 means take strictly the points inside, a very small value
+# will include the points located on the polygon line)
+thresholdPointInHydr = 0.01
+# disabled at the moment:
+# distance that particles need to travel before new particles are initialized through the hydrograph
+# the distance is computed out of the hydrograph timesteps and the hydrograph velocity,
+#  if the velocity is 0, it's set to 1 m/s:
+# (distance = (timestep[i] - timestep[i-1]) * velocity)
+timeStepDistance = 5
+
 #++++++++++++Time stepping parameters
 # fixed time step (also used as first time step when using CFL) [s]
 dt = 0.1