From 6d7300367943403086d3354388b247ef2f71ec0e Mon Sep 17 00:00:00 2001 From: paulasp Date: Fri, 18 Jul 2025 11:24:40 +0200 Subject: [PATCH 1/2] Allow hydrograph as input, additionally to release area first try to get source in different time steps first try for particles with initial velocity hydrographs works with csv table Array containing z - location of particles is updated with hydrograph particles comment FSO: move specific hydrograph functions in a separate debris file adapt tests add tests for new functions add check for requirements of hydrograph csv file update velocity magnitude check if existing particles are withon hydrograph polygon --- avaframe/com1DFA/DFAfunctionsCython.pyx | 73 ++++- avaframe/com1DFA/com1DFA.py | 58 ++-- avaframe/com1DFA/com1DFACfg.ini | 15 + avaframe/com1DFA/debrisFunctions.py | 261 +++++++++++++++ avaframe/in1Data/getInput.py | 50 ++- avaframe/in3Utils/geoTrans.py | 46 ++- avaframe/tests/test_DFAfunctionsCython.py | 67 ++++ avaframe/tests/test_com1DFA.py | 4 +- avaframe/tests/test_debrisFunctions.py | 298 ++++++++++++++++++ avaframe/tests/test_geoTrans.py | 37 +++ avaframe/tests/test_particleInitialisation.py | 1 + 11 files changed, 875 insertions(+), 35 deletions(-) create mode 100644 avaframe/com1DFA/debrisFunctions.py create mode 100644 avaframe/tests/test_debrisFunctions.py diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 7e8b8698d..f6fc61ff5 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/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 \ No newline at end of file diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index 55b9740a6..1d01e5406 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/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,26 +1666,26 @@ 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: @@ -1681,15 +1693,15 @@ def initializeFields(cfg, dem, particles, releaseLine): 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 diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index e00963570..121e0ef4d 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/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 diff --git a/avaframe/com1DFA/debrisFunctions.py b/avaframe/com1DFA/debrisFunctions.py new file mode 100644 index 000000000..b0cb6eb96 --- /dev/null +++ b/avaframe/com1DFA/debrisFunctions.py @@ -0,0 +1,261 @@ +""" +Functions that are used specifically for modeling debris flows (within DebrisFrame). +""" + +# Load modules +import logging +import numpy as np +import copy + +import avaframe.com1DFA.DFAfunctionsCython as DFAfunC +import avaframe.com1DFA.particleTools as particleTools +import avaframe.in3Utils.geoTrans as geoTrans +import avaframe.com1DFA.com1DFA as com1DFA +import avaframe.in2Trans.shpConversion as shpConv +from avaframe.in1Data import getInput as gI + +# create local logger +# change log level in calling module to DEBUG to see log messages +log = logging.getLogger(__name__) + + +def releaseHydrograph(cfg, inputSimLines, particles, fields, dem, zPartArray0, t, atol=1e-05): + """ + Update particles with "new" particles initialised by a hydrograph. + + Parameters + --------- + cfg: configparser + configuration settings + inputSimLines : dict + dictionary with input data dictionaries (releaseLine, hydrographAreaLine,...) + particles : dict + particles dictionary at t that are in the flow already + fields: dict + fields dictionary at t + dem: dict + dictionary with info on DEM data + zPartArray0: dict + dictionary containing z - value of particles at timestep 0 + t: float + timestep of iteration + atol: float + look for matching time steps with atol tolerance - default is atol=1.e-5 + + Returns + --------- + particles: dict + particles dictionary at t including the hydrograph particles + fields: dict + updated fields dictionary at t including the hydrograph particles + zPartArray0: dict + dictionary containing z - value of particles at timestep 0 + """ + hydrValues = inputSimLines["hydrographAreaLine"]["values"] + if np.isclose(t, hydrValues["timeStep"], atol=atol, rtol=0).any(): + iTup = np.where(np.isclose(t, hydrValues["timeStep"], atol=atol, rtol=0)) + # iTup is a tuple containing an array with one value in the first position, so we can extract the index: + i = iTup[0].item() + log.info( + "add hydrograph at timestep: %f s with thickness %s m and velocity %s m/s" + % (t, hydrValues["thickness"][i], hydrValues["velocity"][i]) + ) + # similar workflow to secondary release! + particles, zPartArray0 = addHydrographParticles( + cfg, + particles, + inputSimLines, + hydrValues["thickness"][i], + hydrValues["velocity"][i], + dem, + zPartArray0, + ) + particles = DFAfunC.getNeighborsC(particles, dem) + # update fields (compute grid values) + if fields["computeTA"]: + particles = DFAfunC.computeTrajectoryAngleC(particles, zPartArray0) + particles, fields = DFAfunC.updateFieldsC(cfg["GENERAL"], particles, dem, fields) + + return particles, fields, zPartArray0 + + +def addHydrographParticles(cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0): + """ + add new particles initialized by a hydrograph to particles that are in the flow already + + Parameters + --------- + cfg: configparser object + configuration settings + particles : dict + particles dictionary at t that are in the flow already + inputSimLines : dict + dictionary with input data dictionaries (releaseLine, hydrographAreaLine,...) + thickness: float + thickness of incoming hydrograph + velocityMag: float + velocity of incoming hydrograph + dem: dict + dictionary with info on DEM data + zPartArray0: numpy array + z - value of particles at timestep 0 + + Returns + --------- + particles: dict + particles dictionary at t including the hydrograph particles + zPartArray0: dict + dictionary containing z - value of particles at timestep 0 + """ + hydrLine = inputSimLines["hydrographAreaLine"] + hydrLine["header"] = dem["originalHeader"].copy() + hydrLine = geoTrans.prepareArea( + hydrLine, + dem, + np.sqrt(2), + thList=[thickness], + combine=True, + checkOverlap=False, + ) + + # check if already existing particles are within the hydrograph polygon + # it's possible that there are still a few particles in the polygon with low velocities + # TODO: could think of a threshold of number of particles that are still allowed in the polygons? + mask = geoTrans.getParticlesInPolygon( + particles, hydrLine, cfg["GENERAL"].getfloat("thresholdPointInHydr") + ) + if np.sum(mask) > 0: + # if there is at least one particle within the polygon (including the buffer): + message = ( + "Already existing particles are within the hydrograph polygon, which can cause numerical instabilities (at timestep: %02f s)" + % (particles["t"] + particles["dt"]) + ) + # timestep in particles is not updated yet + log.error(message) + raise ValueError(message) + + particlesHydrograph = com1DFA.initializeParticles( + cfg["GENERAL"], + hydrLine, + dem, + ) + particlesHydrograph = DFAfunC.updateInitialVelocity( + cfg["GENERAL"], particlesHydrograph, dem, velocityMag + ) + + particles = particleTools.mergeParticleDict(particles, particlesHydrograph) + # save initial z position for travel angle computation + zPartArray0 = np.append(zPartArray0, copy.deepcopy(particlesHydrograph["z"])) + return particles, zPartArray0 + + +def checkHydrograph(hydrographValues, hydrCsv): + """ + check if hydrograph satisfied the following requirements: + - hydrograph-timesteps are unique + - provided release-thickness is larger than zero + - the hydrograph-timesteps are not too close (that the particle density becomes too high) + + Parameters + ----------- + hydrCsv: str + directory to csv table containing hydrograph values + hydrographValues: dict + contains hydrograph values: timestep, thickness, velocity + """ + # check if timesteps are unique + timeStepUnique = np.unique(hydrographValues["timeStep"]) + if timeStepUnique.ndim == 0: + if timeStepUnique != hydrographValues["timeStep"]: + message = "The provided hydrograph time steps in %s are not unique" % (hydrCsv) + log.error(message) + raise ValueError(message) + elif len(timeStepUnique) != len(hydrographValues["timeStep"]): + message = "The provided hydrograph timesteps in %s are not unique" % (hydrCsv) + log.error(message) + raise ValueError(message) + + # check that hydrograph thickness > 0 + for th in hydrographValues["thickness"]: + if th <= 0: + message = "For every release time step a thickness > 0 needs to be provided in %s" % (hydrCsv) + log.error(message) + raise ValueError(message) + + +def preparehydrographAreaLine(inputSimFiles, demOri, cfg): + """ + read hydrograph polygon and values + + Parameters + ---------- + inputSimFiles : dict + dictionary containing + - hydrographFile: str, path to hydrograph polygon file + - hydrographCsv: str, path to hydrograph values (csv-)file + cfg: configparser object + configuration for simType + demOri : dict + dictionary with dem info (header original origin), raster data correct mesh cell size + + Returns + ------- + hydrLine: dict + contains hydrograph outline and values, among other things: + - x, y, z + - values: timeStep, thickness, velocity + """ + try: + hydrFile = inputSimFiles["hydrographFile"] + hydrLine = shpConv.readLine(hydrFile, "", demOri) + hydrLine["fileName"] = hydrFile + hydrLine["type"] = "Hydrograph" + gI.checkForMultiplePartsShpArea( + cfg["GENERAL"]["avalancheDir"], hydrLine, "com1DFA", type="hydrograph" + ) + except: + message = "No hydrograph shp file found" + log.error(message) + raise FileNotFoundError(message) + + try: + hydrLine["values"] = gI.getHydrographCsv(inputSimFiles["hydrographCsv"]) + hydrLine["thicknessSource"] = ["csv file"] + except: + message = "No hydrograph csv file found" + log.error(message) + raise FileNotFoundError(message) + + checkHydrograph(hydrLine["values"], inputSimFiles["hydrographCsv"]) + + return hydrLine + + +def checkTravelledDistance(cfgGen, hydrographValues, hydrCsv): + """ + not used at the moment (related to timeStepDistance in the configuration file)! + check if time steps of hydrograph are not to close that the particle density becomes too high + check that particles moved out of hydrograph area before new particles are initialized + time between hydrograph time steps + first timestep is skipped since this is always ok. + + Parameters + ----------- + hydrCsv: str + directory to csv table containing hydrograph values + cfgGen: configparser + configuration settings, section "GENERAL" + hydrographValues: dict + contains hydrograph values: timestep, thickness, velocity + """ + timeStepUnique = np.unique(hydrographValues["timeStep"]) + if timeStepUnique.ndim > 0: + hydrDT = np.append(hydrographValues["timeStep"], 0) - np.append(0, hydrographValues["timeStep"]) + vel = np.where(np.array(hydrographValues["velocity"]) > 0, np.array(hydrographValues["velocity"]), 1) + distance = vel[:-1] * hydrDT[1:-1] + + if np.any(distance < cfgGen.getfloat("timeStepDistance")): + message = "Please select timesteps with greater spacing in %s." % (hydrCsv) + # TODO: error or warning? + log.error(message) + raise ValueError(message) diff --git a/avaframe/in1Data/getInput.py b/avaframe/in1Data/getInput.py index 1d165c9d7..40f7dc8db 100644 --- a/avaframe/in1Data/getInput.py +++ b/avaframe/in1Data/getInput.py @@ -296,6 +296,13 @@ def getInputDataCom1DFA(avaDir): entResInfo["resRemeshed"] = "No" entResInfo["bhdRemeshed"] = "No" + hydrographFile, entResInfo["hydrograph"], _ = getAndCheckInputFiles( + inputDir, "HYDR", "Hydrograph", fileExt="shp" + ) + hydrographCsv, entResInfo["hydrographCsv"], _ = getAndCheckInputFiles( + inputDir, "HYDR", "Hydrograph parameters (csv)", fileExt="csv" + ) + # return DEM, first item of release, entrainment and resistance areas inputSimFiles = { "demFile": demFile, @@ -310,6 +317,8 @@ def getInputDataCom1DFA(avaDir): "kFile": kFile, "tauCFile": tauCFile, "bhdFile": bhdFile, + "hydrographFile": hydrographFile, + "hydrographCsv": hydrographCsv, } for thFile in ["rel", "secondaryRel", "ent"]: @@ -356,7 +365,7 @@ def getAndCheckInputFiles(inputDir, folder, inputType, fileExt="shp", fileSuffix """ available = "No" - supportedFileFormats = [".shp", ".asc", ".tif"] + supportedFileFormats = [".shp", ".asc", ".tif", ".csv"] # Define the directory to search and the extensions if fileExt == "": @@ -395,7 +404,8 @@ def getAndCheckInputFiles(inputDir, folder, inputType, fileExt="shp", fileSuffix if OutputFile.suffix not in supportedFileFormats: message = ( - "Unsupported file format found for OutputFile %s; shp, asc, tif are allowed" % OutputFile + "Unsupported file format found for OutputFile %s; shp, asc, tif, csv are allowed" + % OutputFile ) log.error(message) raise AssertionError(message) @@ -1131,3 +1141,39 @@ def deriveLineRaster( log.info("%s read from %s" % (rasterNameStr[rasterType], str(rasterPath))) return rasterPath, lineDict + + +def getHydrographCsv(hydrCsv): + """ + get hydrograph values from the csv table + TODO: now the first column is defined as timestep, the second as thickness, third as velocity + see DebrisFrame Issue #18 + + Parameters + ----------- + hydrCsv: str + path to csv file containing hydrograph values + + Returns + ----------- + hydrographValues: dict + contains hydrograph values: timestep, thickness, velocity + """ + hydrParameters = np.genfromtxt(hydrCsv, delimiter=",", filling_values=0, skip_header=1) + + if hydrParameters.ndim == 2: + # sort the columns according to the first column (timestep) + hydrParameters = hydrParameters[np.argsort(hydrParameters[:, 0])] + hydrographValues = { + "timeStep": hydrParameters[:, 0], + "thickness": hydrParameters[:, 1], + "velocity": hydrParameters[:, 2], + } + else: + hydrographValues = { + "timeStep": [hydrParameters[0]], + "thickness": [hydrParameters[1]], + "velocity": [hydrParameters[2]], + } + + return hydrographValues diff --git a/avaframe/in3Utils/geoTrans.py b/avaframe/in3Utils/geoTrans.py index a2464d084..81f221dfe 100644 --- a/avaframe/in3Utils/geoTrans.py +++ b/avaframe/in3Utils/geoTrans.py @@ -1323,6 +1323,37 @@ def checkParticlesInRelease(particles, line, radius): particles : dict particles dictionary where particles outside of the polygon have been removed """ + Mask = getParticlesInPolygon(particles, line, radius) + # also remove particles with negative mass + mask = np.where(particles["m"] <= 0, False, True) + Mask = np.logical_and(Mask, mask) + nRemove = len(Mask) - np.sum(Mask) + if nRemove > 0: + particles = particleTools.removePart(particles, Mask, nRemove, "") + log.debug("removed %s particles because they are not within the release polygon" % (nRemove)) + + return particles + + +def getParticlesInPolygon(particles, line, radius): + """ + check which particles are within a polygon (including a buffer) + + Parameters + ---------- + particles : dict + particles dictionary + line: dict + line dictionary of polygon + radius: float + this value is added to the polygon coordinates to decide whether a particles + is located within or outside of the polygon + + Returns + ------- + Mask: numpy array + is True at particle indices that are inside the polygon, and False for outside the polygon + """ NameRel = line["Name"] StartRel = line["Start"] LengthRel = line["Length"] @@ -1336,18 +1367,11 @@ def checkParticlesInRelease(particles, line, radius): "y": line["y"][int(start) : int(end)], "Name": name, } + # get an array that have at the same indices of a particle a True if the particle is within the polygon + # and a False if it is outside the polygon (including a buffer (radius)) mask = pointInPolygon(line["header"], particles, avapath, radius) Mask = np.logical_or(Mask, mask) - - # also remove particles with negative mass - mask = np.where(particles["m"] <= 0, False, True) - Mask = np.logical_and(Mask, mask) - nRemove = len(Mask) - np.sum(Mask) - if nRemove > 0: - particles = particleTools.removePart(particles, Mask, nRemove, "") - log.debug("removed %s particles because they are not within the release polygon" % (nRemove)) - - return particles + return Mask def pointInPolygon(demHeader, points, Line, radius): @@ -1368,7 +1392,7 @@ def pointInPolygon(demHeader, points, Line, radius): Returns ------- Mask : 1D numpy array - Mask of particles to keep + Mask of particles that are within the polygon """ xllc = demHeader["xllcenter"] yllc = demHeader["yllcenter"] diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py index 488ccf91b..387f909b0 100644 --- a/avaframe/tests/test_DFAfunctionsCython.py +++ b/avaframe/tests/test_DFAfunctionsCython.py @@ -10,6 +10,7 @@ import avaframe.com1DFA.DFAfunctionsCython as DFAfunC import avaframe.com1DFA.DFAtools as DFAtls import avaframe.com1DFA.DFAToolsCython as DFAtlsC +import avaframe.in3Utils.geoTrans as geoTrans def test_getNeighborsC(capfd): @@ -1005,6 +1006,72 @@ def test_computeForceC(): assert (forceFrictArray == test_forceFrictArray).all() +def test_updateInitialVelocity(): + particles = { + "x": np.array([10.0]), + "y": np.array([10.0]), + "ux": np.array([0.0]), + "uy": np.array([0.0]), + "uz": np.array([0.0]), + "nPart": 1, + "velocityMag": np.array([0.0]), + } + + # incline plane + dem = { + "header": {"nrows": 5, "ncols": 5, "cellsize": 5}, + "rasterData": np.array( + [ + [100, 100, 100, 100, 100], + [95, 95, 95, 95, 95], + [90, 90, 90, 90, 90], + [85, 85, 85, 85, 85], + [80, 80, 80, 80, 80], + ] + ), + } + + cfg = configparser.ConfigParser() + cfg["GENERAL"] = {"interpOption": "2"} + + velocityMag = 10.0 + dem = geoTrans.getNormalMesh(dem, num=1) + + particlesTest = { + "x": np.array([10.0]), + "y": np.array([10.0]), + "ux": np.array([0.0]), + "uy": np.array([np.sqrt(50)]), + "uz": np.array([-np.sqrt(50)]), + "nPart": 1, + "velocityMag": np.array([velocityMag]), + } + + particlesVelocity = DFAfunC.updateInitialVelocity(cfg["GENERAL"], particles, dem, velocityMag) + for key in particlesTest: + assert np.isclose(particlesTest[key], particlesVelocity[key], atol=1e-4) + assert ( + DFAtls.norm(particlesVelocity["ux"], particlesVelocity["uy"], particlesVelocity["uz"]) == velocityMag + ) + + # flat plane + dem["rasterData"] = np.ones((dem["header"]["ncols"], dem["header"]["nrows"])) + dem = geoTrans.getNormalMesh(dem, num=1) + particlesTest = { + "x": np.array([10.0]), + "y": np.array([10.0]), + "ux": np.array([0.0]), + "uy": np.array([0.0]), + "uz": np.array([0.0]), + "nPart": 1, + "velocityMag": np.array([0.0]), + } + + particlesVelocity = DFAfunC.updateInitialVelocity(cfg["GENERAL"], particles, dem, velocityMag) + for key in particlesTest: + assert np.isclose(particlesTest[key], particlesVelocity[key], atol=1e-3) + + """ TODO: When calling pytest, the following function raises an error ("Fatal Python error: Aborted") (see issue #1002?) diff --git a/avaframe/tests/test_com1DFA.py b/avaframe/tests/test_com1DFA.py index 6ac6f78cc..feec8bffd 100644 --- a/avaframe/tests/test_com1DFA.py +++ b/avaframe/tests/test_com1DFA.py @@ -1286,6 +1286,7 @@ def test_releaseSecRelArea(): "thickness": [0.5, 1.0, 0.5], "rasterData": [secRelRaster1, secRelRaster2, secRelRaster3], "initializedFrom": "shapefile", + "type": "secondary release", } secondaryReleaseInfo["header"] = demHeader secondaryReleaseInfo["header"]["xllcenter"] = dem["originalHeader"]["xllcenter"] @@ -1364,7 +1365,7 @@ def test_getRelThFromPart(): # setup required input cfg = configparser.ConfigParser() cfg["GENERAL"] = {"relThFromFile": "True", "relTh": ""} - inputSimLines = {"releaseLine": {"thickness": ["1.2", "1.5"], "id": ["0", "1"]}} + inputSimLines = {"releaseLine": {"thickness": ["1.2", "1.5"], "id": ["0", "1"], "type": "Release"}} relThField = "" # call function to be tested @@ -1441,6 +1442,7 @@ def test_initializeParticles(): "Name": [""], "thickness": [1.0], "rasterData": relRaster, + "type": "Release", } releaseLine["header"] = demHeader diff --git a/avaframe/tests/test_debrisFunctions.py b/avaframe/tests/test_debrisFunctions.py new file mode 100644 index 000000000..0d6369069 --- /dev/null +++ b/avaframe/tests/test_debrisFunctions.py @@ -0,0 +1,298 @@ +""" Tests for module debrisFunctions """ + +import numpy as np +import pytest +import configparser + +import avaframe.com1DFA.debrisFunctions as debF + + +def test_addHydrographParticles(): + inputSimLines = { + "hydrographAreaLine": { + "Name": ["testHydr"], + "Start": np.asarray([0.0]), + "Length": np.asarray([5]), + "type": "Hydrograph", + "x": np.asarray( + [ + 0, + 10.0, + 10.0, + 0.0, + 0.0, + ] + ) + - 2.5, + "y": np.asarray([0.0, 0.0, 10.0, 10.0, 0.0]) - 2.5, + "thicknessSource": ["csv file"], + "thickness": 1, + } + } + thickness = inputSimLines["hydrographAreaLine"]["thickness"] + velocityMag = 0 + + demHeader = {} + demHeader["xllcenter"] = 0 + demHeader["yllcenter"] = 0 + demHeader["cellsize"] = 5.0 + demHeader["nodata_value"] = -9999 + demHeader["nrows"] = 7 + demHeader["ncols"] = 7 + dem = {"header": demHeader} + dem["rasterData"] = np.ones((demHeader["nrows"], demHeader["ncols"])) + dem["originalHeader"] = dem["header"] + dem["areaRaster"] = np.ones((demHeader["nrows"], demHeader["ncols"])) + dem["Nx"] = np.zeros_like(dem["rasterData"]) + dem["Ny"] = np.zeros_like(dem["rasterData"]) + dem["Nz"] = np.zeros_like(dem["rasterData"]) + + cfg = configparser.ConfigParser() + cfg["GENERAL"] = { + "resType": "ppr|pft|pfv", + "rho": "1000.", + "gravAcc": "9.81", + "cpIce": "2050", + "TIni": "-10", + "avalancheDir": "data/avaKotHYDR", + "massPerParticleDeterminationMethod": "MPPDH", + "interpOption": "2", + "initialiseParticlesFromFile": "False", + "iniStep": "False", + "seed": "12345", + "sphKernelRadius": "1", + "deltaTh": "1", + "initPartDistType": "uniform", + "thresholdPointInPoly": "0.001", + "massPerPart": "1000", + "thresholdPointInHydr": "0", + } + + particles = { + "nPart": 3, + "x": np.array([12, 20, 30]), + "y": np.array([5, 10, 30]), + "z": np.array([1, 1, 1]), + "m": np.array([1000, 1000, 1000]), + "idFixed": np.array([0, 0, 0]), + "t": 1.0, + "dt": 0.1, + } + nPart = particles["nPart"] + particles["totalEnthalpy"] = ( + cfg["GENERAL"].getfloat("TIni") * cfg["GENERAL"].getfloat("cpIce") + + cfg["GENERAL"].getfloat("gravAcc") * particles["z"] + ) + particles["massPerPart"] = 1000 + particles["mTot"] = np.sum(particles["m"]) + particles["tPlot"] = 0 + particles["h"] = np.ones(nPart) + particles["ux"] = np.zeros(nPart) + particles["uy"] = np.zeros(nPart) + particles["uz"] = np.zeros(nPart) + particles["uAcc"] = np.zeros(nPart) + particles["velocityMag"] = np.zeros(nPart) + particles["trajectoryLengthXY"] = np.zeros(nPart) + particles["trajectoryLengthXYCor"] = np.zeros(nPart) + particles["trajectoryLengthXYZ"] = np.zeros(nPart) + particles["trajectoryAngle"] = np.zeros(nPart) + particles["stoppCriteria"] = False + particles["peakForceSPH"] = 0.0 + particles["forceSPHIni"] = 0.0 + particles["peakMassFlowing"] = 0 + particles["xllcenter"] = dem["originalHeader"]["xllcenter"] + particles["yllcenter"] = dem["originalHeader"]["yllcenter"] + particles["nExitedParticles"] = 0.0 + particles["dmDet"] = np.zeros(nPart) + particles["dmEnt"] = np.zeros(nPart) + + zPartArray0 = np.array([1, 1, 1]) + newParticleNumber = 4 + particlesTest = { + "nPart": newParticleNumber + 3, + "mTot": 7000, + "x": np.append(particles["x"], np.array([0, 5, 0, 5])), + "y": np.append(particles["y"], np.array([0, 0, 5, 5])), + "z": np.append(particles["z"], np.ones([newParticleNumber])), + "m": np.append(particles["m"], np.ones([newParticleNumber]) * 1000), + } + zPartArray0Test = np.ones(particlesTest["nPart"]) + + particlesHydr, zPartArray0Hydr = debF.addHydrographParticles( + cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0 + ) + + assert np.all(np.equal(zPartArray0Hydr, zPartArray0Test)) + for key in particlesTest: + if key in ["nPart", "mTot"]: + assert particlesTest[key] == particlesHydr[key] + else: + assert np.all(np.equal(particlesTest[key], particlesHydr[key])) + for key in ["ux", "uy", "uz", "velocityMag"]: + assert np.all(np.equal(np.zeros(particlesTest["nPart"]), particlesHydr[key])) + + cfg["GENERAL"]["deltaTh"] = "0.25" + cfg["GENERAL"]["initPartDistType"] = "random" + cfg["GENERAL"]["thresholdMassSplit"] = "1.5" + + particlesHydr, zPartArray0Hydr = debF.addHydrographParticles( + cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0 + ) + assert particlesHydr["nPart"] == 16 + 3 + for key in ["ux", "uy", "uz", "velocityMag", "x", "y", "z"]: + assert len(particlesHydr[key]) == particlesHydr["nPart"] + assert particlesHydr["mTot"] == 7000 + + particles["x"] = np.array([4, 10, 30]) + particles["y"] = np.array([5, 3, 30]) + + with pytest.raises(ValueError): + debF.addHydrographParticles(cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0) + +""" +TODO: When calling pytest, the following function raises an error ("Fatal Python error: Aborted") +(see issue #1002?) + +def test_releaseHydrograph(): + inputSimLines = { + "hydrographAreaLine": { + "Name": ["testHydr"], + "Start": np.asarray([0.0]), + "Length": np.asarray([5]), + "type": "Hydrograph", + "x": np.asarray( + [ + 0, + 10.0, + 10.0, + 0.0, + 0.0, + ] + ) + - 2.5, + "y": np.asarray([0.0, 0.0, 10.0, 10.0, 0.0]) - 2.5, + "thicknessSource": ["csv file"], + "thickness": 1, + "values": { + "timeStep": np.array([0, 1]), + "thickness": np.array([1, 1]), + "velocity": np.array([0, 0]), + }, + } + } + thickness = inputSimLines["hydrographAreaLine"]["thickness"] + velocityMag = 0 + + demHeader = {} + demHeader["xllcenter"] = 0 + demHeader["yllcenter"] = 0 + demHeader["cellsize"] = 5.0 + demHeader["nodata_value"] = -9999 + demHeader["nrows"] = 7 + demHeader["ncols"] = 7 + dem = {"header": demHeader} + dem["rasterData"] = np.ones((demHeader["nrows"], demHeader["ncols"])) + dem["originalHeader"] = dem["header"] + dem["areaRaster"] = np.ones((demHeader["nrows"], demHeader["ncols"])) + dem["Nx"] = np.zeros_like(dem["rasterData"]) + dem["Ny"] = np.zeros_like(dem["rasterData"]) + dem["Nz"] = np.zeros_like(dem["rasterData"]) + dem["headerNeighbourGrid"] = demHeader + + cfg = configparser.ConfigParser() + cfg["GENERAL"] = { + "resType": "ppr|pft|pfv", + "rho": "1000.", + "rhoEnt": "1000.", + "gravAcc": "9.81", + "cpIce": "2050", + "TIni": "-10", + "avalancheDir": "data/avaKotHYDR", + "massPerParticleDeterminationMethod": "MPPDH", + "interpOption": "2", + "initialiseParticlesFromFile": "False", + "iniStep": "False", + "seed": "12345", + "sphKernelRadius": "1", + "deltaTh": "1", + "initPartDistType": "uniform", + "thresholdPointInPoly": "0.001", + "massPerPart": "1000", + } + + particles = { + "nPart": 3, + "x": np.array([10, 20, 30]), + "y": np.array([5, 10, 30]), + "z": np.array([1, 1, 1]), + "m": np.array([1000, 1000, 1000]), + "idFixed": np.array([0, 0, 0]), + "dmDet": np.array([0, 0, 0]), + "dmEnt": np.array([0, 0, 0]), + "ux": np.array([0, 0, 0]), + "uy": np.array([0, 0, 0]), + "uz": np.array([0, 0, 0]), + "trajectoryAngle": np.array([0, 0, 0]), + "stoppedParticles": { + "m": np.empty(0), + "x": np.empty(0), + "y": np.empty(0), + }, + } + nPart = particles["nPart"] + particles["totalEnthalpy"] = ( + cfg["GENERAL"].getfloat("TIni") * cfg["GENERAL"].getfloat("cpIce") + + cfg["GENERAL"].getfloat("gravAcc") * particles["z"] + ) + particles["massPerPart"] = 1000 + particles["mTot"] = np.sum(particles["m"]) + particles["tPlot"] = 0 + particles["h"] = np.ones(nPart) + particles["uAcc"] = np.zeros(nPart) + particles["velocityMag"] = np.zeros(nPart) + particles["trajectoryLengthXY"] = np.zeros(nPart) + particles["trajectoryLengthXYCor"] = np.zeros(nPart) + particles["trajectoryLengthXYZ"] = np.zeros(nPart) + particles["stoppCriteria"] = False + particles["peakForceSPH"] = 0.0 + particles["forceSPHIni"] = 0.0 + particles["peakMassFlowing"] = 0 + particles["xllcenter"] = dem["originalHeader"]["xllcenter"] + particles["yllcenter"] = dem["originalHeader"]["yllcenter"] + particles["nExitedParticles"] = 0.0 + + fields = { + "computeTA": 0, + "computeKE": 0, + "computeP": 0, + "pfv": np.zeros_like(dem["rasterData"]), + "ppr": np.zeros_like(dem["rasterData"]), + "pft": np.zeros_like(dem["rasterData"]), + "pta": np.zeros_like(dem["rasterData"]), + "pke": np.zeros_like(dem["rasterData"]), + "dmDet": np.zeros_like(dem["rasterData"]), + } + fields["pft"][[2, 4, 6], [1, 2, 6]] = 1 + + zPartArray0 = np.array([1, 1, 1]) + + t = 0.00 + + newParticleNumber = 4 + particlesTest = { + "nPart": newParticleNumber + 3, + "mTot": 7000, + "x": np.append(particles["x"], np.array([0, 5, 0, 5])), + "y": np.append(particles["y"], np.array([0, 0, 5, 5])), + "z": np.append(particles["z"], np.ones([newParticleNumber])), + "m": np.append(particles["m"], np.ones([newParticleNumber]) * 1000), + } + zPartArray0Test = np.ones(particlesTest["nPart"]) + + debF.updateParticlesHydrograph(cfg, inputSimLines, particles, fields, dem, zPartArray0, t) + +""" + + +if __name__ == "__main__": + test_addHydrographParticles() diff --git a/avaframe/tests/test_geoTrans.py b/avaframe/tests/test_geoTrans.py index 30bdad6ba..6f20facad 100644 --- a/avaframe/tests/test_geoTrans.py +++ b/avaframe/tests/test_geoTrans.py @@ -510,6 +510,43 @@ def test_checkParticlesInRelease(): assert particles2["nPart"] == 4 +def test_getParticlesInPolygon(): + """test if particles are within release polygon""" + # setup required input + hydrographLine = { + "Name": ["testRel"], + "Start": np.asarray([0.0]), + "Length": np.asarray([5]), + "type": "Release", + "x": np.asarray([0, 10.0, 10.0, 0.0, 0.0]), + "y": np.asarray([0.0, 0.0, 10.0, 10.0, 0.0]), + } + demHeader = {} + demHeader["xllcenter"] = 0.0 + demHeader["yllcenter"] = 0.0 + demHeader["cellsize"] = 1.0 + demHeader["noDataValue"] = -9999 + demHeader["nrows"] = 5 + demHeader["ncols"] = 5 + hydrographLine["header"] = demHeader + particles = { + "x": np.asarray([2.4, 9.7, 10.02, 11.5]), + "y": np.asarray([2.4, 9.7, 10.2, 11.5]), + "nPart": 4, + "m": np.asarray([1.4, 1.7, 1.4, 1.8]), + } + radius = 0 + + MaskTest = np.array([True, True, False, False]) + Mask = geoTrans.getParticlesInPolygon(particles, hydrographLine, radius) + assert np.array_equal(MaskTest, Mask) + + radius = 0.5 + MaskTest = np.array([True, True, True, False]) + Mask = geoTrans.getParticlesInPolygon(particles, hydrographLine, radius) + assert np.array_equal(MaskTest, Mask) + + def test_remeshDEM(tmp_path): """test size of interpolated data onto new mesh""" diff --git a/avaframe/tests/test_particleInitialisation.py b/avaframe/tests/test_particleInitialisation.py index dbb631afa..bd5c51b64 100644 --- a/avaframe/tests/test_particleInitialisation.py +++ b/avaframe/tests/test_particleInitialisation.py @@ -163,6 +163,7 @@ def test_getIniPosition(tmp_path): "header": demOri["header"], "y": np.asarray([5.0, 5.0, 10.0, 10.0, 5.0]), "rasterData": relRaster, + "type": "Release", } } From 69ddad383ecd3af358362afe75b7d0ee0d70dbde Mon Sep 17 00:00:00 2001 From: Felix Oesterle Date: Wed, 19 Nov 2025 21:33:44 +0100 Subject: [PATCH 2/2] refactor(runStandardTestsCom1DFA): switch benchmark filtering back to tags --- avaframe/runStandardTestsCom1DFA.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/avaframe/runStandardTestsCom1DFA.py b/avaframe/runStandardTestsCom1DFA.py index 117020395..5c29c7c1e 100644 --- a/avaframe/runStandardTestsCom1DFA.py +++ b/avaframe/runStandardTestsCom1DFA.py @@ -254,10 +254,10 @@ def runSingleTest( # filter benchmarks for tag standardTest # filterType = "TAGS" # valuesList = ["resistance"] -# filterType = "TAGS" -# valuesList = ["standardTest", "standardTestSnowGlide"] -filterType = "NAME" -valuesList = ["avaInclinedPlaneEntresTest"] +filterType = "TAGS" +valuesList = ["standardTest", "standardTestSnowGlide"] +# filterType = "NAME" +# valuesList = ["avaInclinedPlaneEntresTest"] testList = tU.filterBenchmarks(testDictList, filterType, valuesList, condition="or")