From 7feba96f60d1e76eea5a5be0189b1eb19d251044 Mon Sep 17 00:00:00 2001 From: Anna Wirbel Date: Fri, 8 May 2026 14:59:26 +0200 Subject: [PATCH 1/4] add option to only save certain particle props --- avaframe/com1DFA/DFAfunctionsCython.pyx | 9 +++++++ avaframe/com1DFA/com1DFA.py | 35 +++++++++++++++++++------ avaframe/com1DFA/com1DFACfg.ini | 3 +++ 3 files changed, 39 insertions(+), 8 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 919707a3c..42a595ef9 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -1299,6 +1299,9 @@ def updateFieldsC(cfg, particles, dem, fields): cdef double[:, :] FTEntBilinear = np.zeros((nrows, ncols)) # declare intermediate step variables cdef double[:] hBB = np.zeros((nPart)) + cdef double[:] indxPart = np.zeros((nPart)) + cdef double[:] indyPart = np.zeros((nPart)) + cdef double[:] iCellPart = np.zeros((nPart)) cdef double m, dmDet, demEnt, h, x, y, z, s, ux, uy, uz, nx, ny, nz, hbb, hLim, areaPart, trajectoryAngle cdef double mStop, xStop, yStop cdef int k, i, l @@ -1328,6 +1331,9 @@ def updateFieldsC(cfg, particles, dem, fields): # for the travel angle we simply do a nearest interpolation indx = math.round(x / csz) indy = math.round(y / csz) + indxPart[k] = indx + indyPart[k] = indy + iCellPart[k] = iCell if computeTA: trajectoryAngle = trajectoryAngleArray[k] travelAngleField[indy, indx] = max(travelAngleField[indy, indx], trajectoryAngle) @@ -1429,6 +1435,9 @@ def updateFieldsC(cfg, particles, dem, fields): hBB[k] = hbb particles['h'] = np.asarray(hBB) + particles['indXPart'] = np.asarray(indxPart) + particles['indYPart'] = np.asarray(indyPart) + particles['iCell'] = np.asarray(iCellPart) return particles, fields diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index 5b60a9911..ab4fcd7ad 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -1586,6 +1586,9 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel particles["massEntrained"] = 0.0 particles["massDetrained"] = 0.0 particles["massStopped"] = 0.0 + particles["indXPart"] = np.zeros(np.shape(hPartArray)) + particles["indYPart"] = np.zeros(np.shape(hPartArray)) + particles["iCell"] = np.zeros(np.shape(hPartArray)) # remove particles that might lay outside of the release polygon if ( not cfg.getboolean("iniStep") @@ -2159,7 +2162,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si exportFields(cfg, t, fields, dem, outDir, cuSimName, TSave="initial") if "particles" in resTypes: - savePartToPickle(particles, outDirData, cuSimName) + savePartToPickle(particles, outDirData, cuSimName, cfg["EXPORTS"]["particleProperties"]) # Update dtSave to remove the initial timestep we just saved dtSave = updateSavingTimeStep(dtSaveOriginal, cfgGen, t) @@ -2284,7 +2287,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si # export particles dictionaries of saving time steps if "particles" in resTypes: - savePartToPickle(particles, outDirData, cuSimName) + savePartToPickle(particles, outDirData, cuSimName, cfg["EXPORTS"]["particleProperties"]) # export particles properties for visulation if cfg["VISUALISATION"].getboolean("writePartToCSV"): @@ -2416,7 +2419,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si # export particles dictionaries of saving time steps if "particles" in resTypes: - savePartToPickle(particles, outDirData, cuSimName) + savePartToPickle(particles, outDirData, cuSimName, cfg["EXPORTS"]["particleProperties"]) # save contour line for each sim only if the field is properly computed (not a dummy array) contourResType = cfg["VISUALISATION"]["contourResType"] @@ -2834,7 +2837,7 @@ def releaseSecRelArea(cfg, particles, fields, dem, zPartArray0, reportAreaInfo): return particles, zPartArray0, reportAreaInfo -def savePartToPickle(dictList, outDir, logName): +def savePartToPickle(dictList, outDir, logName, particleProperties=""): """Save each dictionary from a list to a pickle in outDir; works also for one dictionary instead of list Note: particle coordinates are still in com1DFA reference system with origin 0,0 @@ -2846,16 +2849,32 @@ def savePartToPickle(dictList, outDir, logName): path to output directory logName : str simulation Id + particleProperties: str + particle properties to be saved, if empty all particle properties are saved, t (time info) always appended """ + # create list of particle properties and append t (time info) + if particleProperties != "": + particleProperties = particleProperties.split("|") + if "t" not in particleProperties: + particleProperties.append("t") + if isinstance(dictList, list): for dict in dictList: - fi = open(outDir / ("particles_%s_%09.4f.pickle" % (logName, dict["t"])), "wb") - pickle.dump(dict, fi) - fi.close() + if particleProperties != "": + particlesToSave = {key: dict[key] for key in particleProperties} + else: + particlesToSave = dict + fi = open(outDir / ("particles_%s_%09.4f.pickle" % (logName, dict["t"])), "wb") + pickle.dump(particlesToSave, fi) + fi.close() else: + if particleProperties != "": + particlesToSave = {key: dictList[key] for key in particleProperties} + else: + particlesToSave = dictList fi = open(outDir / ("particles_%s_%09.4f.pickle" % (logName, dictList["t"])), "wb") - pickle.dump(dictList, fi) + pickle.dump(particlesToSave, fi) fi.close() diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index 6c305d86b..666eebc2b 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -570,4 +570,7 @@ exportData = True exportRasters = False # use LZW compression when writing TIFF raster files useCompression = True +# particle properties - list all properties that shall be saved, t is always added +particleProperties = + From a43c799ced23b58aff79ece6d7a9de28e4523623 Mon Sep 17 00:00:00 2001 From: Anna Wirbel Date: Mon, 11 May 2026 15:16:27 +0200 Subject: [PATCH 2/4] fix bug --- avaframe/com1DFA/com1DFA.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index ab4fcd7ad..963b89269 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -2865,9 +2865,9 @@ def savePartToPickle(dictList, outDir, logName, particleProperties=""): particlesToSave = {key: dict[key] for key in particleProperties} else: particlesToSave = dict - fi = open(outDir / ("particles_%s_%09.4f.pickle" % (logName, dict["t"])), "wb") - pickle.dump(particlesToSave, fi) - fi.close() + fi = open(outDir / ("particles_%s_%09.4f.pickle" % (logName, dict["t"])), "wb") + pickle.dump(particlesToSave, fi) + fi.close() else: if particleProperties != "": particlesToSave = {key: dictList[key] for key in particleProperties} From 3e2135d00f13856d8eac67a4f3b2acc7629eaae1 Mon Sep 17 00:00:00 2001 From: Anna Wirbel Date: Tue, 12 May 2026 16:27:44 +0200 Subject: [PATCH 3/4] check for particle properties, rename key --- avaframe/com1DFA/com1DFA.py | 86 +++++++++++++++++++++++++++++---- avaframe/com1DFA/com1DFACfg.ini | 2 +- 2 files changed, 77 insertions(+), 11 deletions(-) diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index 963b89269..bd148bbc7 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -2162,7 +2162,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si exportFields(cfg, t, fields, dem, outDir, cuSimName, TSave="initial") if "particles" in resTypes: - savePartToPickle(particles, outDirData, cuSimName, cfg["EXPORTS"]["particleProperties"]) + savePartToPickle(particles, outDirData, cuSimName, cfg=cfg) # Update dtSave to remove the initial timestep we just saved dtSave = updateSavingTimeStep(dtSaveOriginal, cfgGen, t) @@ -2287,7 +2287,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si # export particles dictionaries of saving time steps if "particles" in resTypes: - savePartToPickle(particles, outDirData, cuSimName, cfg["EXPORTS"]["particleProperties"]) + savePartToPickle(particles, outDirData, cuSimName, cfg=cfg) # export particles properties for visulation if cfg["VISUALISATION"].getboolean("writePartToCSV"): @@ -2419,7 +2419,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si # export particles dictionaries of saving time steps if "particles" in resTypes: - savePartToPickle(particles, outDirData, cuSimName, cfg["EXPORTS"]["particleProperties"]) + savePartToPickle(particles, outDirData, cuSimName, cfg=cfg) # save contour line for each sim only if the field is properly computed (not a dummy array) contourResType = cfg["VISUALISATION"]["contourResType"] @@ -2837,7 +2837,7 @@ def releaseSecRelArea(cfg, particles, fields, dem, zPartArray0, reportAreaInfo): return particles, zPartArray0, reportAreaInfo -def savePartToPickle(dictList, outDir, logName, particleProperties=""): +def savePartToPickle(dictList, outDir, logName, cfg=""): """Save each dictionary from a list to a pickle in outDir; works also for one dictionary instead of list Note: particle coordinates are still in com1DFA reference system with origin 0,0 @@ -2849,15 +2849,81 @@ def savePartToPickle(dictList, outDir, logName, particleProperties=""): path to output directory logName : str simulation Id - particleProperties: str - particle properties to be saved, if empty all particle properties are saved, t (time info) always appended + cfg: str or configparser object + ['EXPORTS'] and ['GENERAL'] settings to provide particle properties to be saved, + if empty str all particle properties are saved, t (time info) always appended """ + dictKeys = [ + "nPart", + "x", + "y", + "trajectoryLengthXY", + "trajectoryLengthXYCor", + "trajectoryLengthXYZ", + "z", + "m", + "dmDet", + "massPerPart", + "nPPK", + "mTot", + "h", + "ux", + "uy", + "uz", + "uAcc", + "stoppCriteria", + "kineticEne", + "trajectoryAngle", + "potentialEne", + "peakKinEne", + "peakMassFlowing", + "simName", + "xllcenter", + "yllcenter", + "ID", + "nID", + "parentID", + "t", + "inCellDEM", + "indXDEM", + "indYDEM", + "indPartInCell", + "partInCell", + "secondaryReleaseInfo", + "iterate", + "idFixed", + "peakForceSPH", + "forceSPHIni", + "totalEnthalpy", + "velocityMag", + "nExitedParticles", + "tPlot", + "dmEnt", + "stoppedParticles", + "massInitialized", + "massEntrained", + "massDetrained", + "massStopped", + ] + # create list of particle properties and append t (time info) - if particleProperties != "": - particleProperties = particleProperties.split("|") - if "t" not in particleProperties: - particleProperties.append("t") + if isinstance(cfg, configparser.ConfigParser): + # first check if particle properties are valid + nonExisting = [ + item for item in cfg["EXPORTS"]["exportParticleProperties"].split("|") if item not in dictKeys + ] + if len(nonExisting) > 0: + message = "These particle properties are not available %s" % nonExisting + log.error(message) + raise AttributeError(message) + + particleProperties = list(set(["t"] + cfg["EXPORTS"]["exportParticleProperties"].split("|"))) + if cfg["TRACKPARTICLES"].getboolean("trackParticles"): + trackParticleProperties = cfg["TRACKPARTICLES"]["particleProperties"].split("|") + particleProperties = set( + ["x", "y", "z", "ux", "uy", "uz", "m", "h"] + particleProperties + trackParticleProperties + ) if isinstance(dictList, list): for dict in dictList: diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index 666eebc2b..d33a0bea9 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -571,6 +571,6 @@ exportRasters = False # use LZW compression when writing TIFF raster files useCompression = True # particle properties - list all properties that shall be saved, t is always added -particleProperties = +exportParticleProperties = From 085842a2c98f41013a0e39823493d8173c7b84d1 Mon Sep 17 00:00:00 2001 From: Anna Wirbel Date: Wed, 13 May 2026 11:05:36 +0200 Subject: [PATCH 4/4] add pytest --- avaframe/com1DFA/DFAfunctionsCython.pyx | 9 ---- avaframe/com1DFA/com1DFA.py | 43 ++++++++------- avaframe/com1DFA/com1DFACfg.ini | 2 +- avaframe/tests/test_com1DFA.py | 69 +++++++++++++++++++++++++ 4 files changed, 94 insertions(+), 29 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 42a595ef9..919707a3c 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -1299,9 +1299,6 @@ def updateFieldsC(cfg, particles, dem, fields): cdef double[:, :] FTEntBilinear = np.zeros((nrows, ncols)) # declare intermediate step variables cdef double[:] hBB = np.zeros((nPart)) - cdef double[:] indxPart = np.zeros((nPart)) - cdef double[:] indyPart = np.zeros((nPart)) - cdef double[:] iCellPart = np.zeros((nPart)) cdef double m, dmDet, demEnt, h, x, y, z, s, ux, uy, uz, nx, ny, nz, hbb, hLim, areaPart, trajectoryAngle cdef double mStop, xStop, yStop cdef int k, i, l @@ -1331,9 +1328,6 @@ def updateFieldsC(cfg, particles, dem, fields): # for the travel angle we simply do a nearest interpolation indx = math.round(x / csz) indy = math.round(y / csz) - indxPart[k] = indx - indyPart[k] = indy - iCellPart[k] = iCell if computeTA: trajectoryAngle = trajectoryAngleArray[k] travelAngleField[indy, indx] = max(travelAngleField[indy, indx], trajectoryAngle) @@ -1435,9 +1429,6 @@ def updateFieldsC(cfg, particles, dem, fields): hBB[k] = hbb particles['h'] = np.asarray(hBB) - particles['indXPart'] = np.asarray(indxPart) - particles['indYPart'] = np.asarray(indyPart) - particles['iCell'] = np.asarray(iCellPart) return particles, fields diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index bd148bbc7..694548625 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -10,7 +10,6 @@ import os import pathlib import pickle -import platform import re import time from datetime import datetime @@ -1586,9 +1585,6 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel particles["massEntrained"] = 0.0 particles["massDetrained"] = 0.0 particles["massStopped"] = 0.0 - particles["indXPart"] = np.zeros(np.shape(hPartArray)) - particles["indYPart"] = np.zeros(np.shape(hPartArray)) - particles["iCell"] = np.zeros(np.shape(hPartArray)) # remove particles that might lay outside of the release polygon if ( not cfg.getboolean("iniStep") @@ -2909,21 +2905,30 @@ def savePartToPickle(dictList, outDir, logName, cfg=""): # create list of particle properties and append t (time info) if isinstance(cfg, configparser.ConfigParser): - # first check if particle properties are valid - nonExisting = [ - item for item in cfg["EXPORTS"]["exportParticleProperties"].split("|") if item not in dictKeys - ] - if len(nonExisting) > 0: - message = "These particle properties are not available %s" % nonExisting - log.error(message) - raise AttributeError(message) - - particleProperties = list(set(["t"] + cfg["EXPORTS"]["exportParticleProperties"].split("|"))) - if cfg["TRACKPARTICLES"].getboolean("trackParticles"): - trackParticleProperties = cfg["TRACKPARTICLES"]["particleProperties"].split("|") - particleProperties = set( - ["x", "y", "z", "ux", "uy", "uz", "m", "h"] + particleProperties + trackParticleProperties - ) + if cfg["EXPORTS"]["exportParticleProperties"] == "": + particleProperties = "" + else: + # first check if particle properties are valid + nonExisting = [ + item + for item in cfg["EXPORTS"]["exportParticleProperties"].split("|") + if item not in dictKeys + ] + if len(nonExisting) > 0: + message = "These particle properties are not available %s" % nonExisting + log.error(message) + raise AttributeError(message) + + particleProperties = list(set(["t"] + cfg["EXPORTS"]["exportParticleProperties"].split("|"))) + if cfg["TRACKPARTICLES"].getboolean("trackParticles"): + trackParticleProperties = cfg["TRACKPARTICLES"]["particleProperties"].split("|") + particleProperties = set( + ["x", "y", "z", "ux", "uy", "uz", "m", "h"] + + particleProperties + + trackParticleProperties + ) + else: + particleProperties = "" if isinstance(dictList, list): for dict in dictList: diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index d33a0bea9..0927cca80 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -524,7 +524,7 @@ thresholdPointInPoly = 0.001 [TRACKPARTICLES] # if particles should be tracked - don't forget to specify the "tSteps" you want to -# save further up (for example tStep = 0:1 will lead to tracking patiles every 1 second) +# save further up (for example tStep = 0:1 will lead to tracking particles every 1 second) trackParticles = False # centerTrackPartPoint of the location of the particles to track (x|y coordinates) centerTrackPartPoint = 2933|-4010 diff --git a/avaframe/tests/test_com1DFA.py b/avaframe/tests/test_com1DFA.py index 63bb20ee6..abf772987 100644 --- a/avaframe/tests/test_com1DFA.py +++ b/avaframe/tests/test_com1DFA.py @@ -1821,6 +1821,75 @@ def test_savePartToPickle(tmp_path): assert np.array_equal(particlesRead3["m"], particles1["m"]) assert particlesRead3["t"] == 0.0 + # call function to be tested + logName = "simNameTest4" + cfg = configparser.ConfigParser() + cfg["EXPORTS"] = {"exportParticleProperties": "x|m"} + cfg["TRACKPARTICLES"] = {"trackParticles": False} + com1DFA.savePartToPickle(particles1, outDir, logName, cfg=cfg) + + # read pickle + picklePath4 = outDir / "particles_simNameTest4_0000.0000.pickle" + particlesRead4 = pickle.load(open(picklePath4, "rb")) + + assert np.array_equal(particlesRead4["x"], particles1["x"]) + assert "y" not in particlesRead4.keys() + assert np.array_equal(particlesRead4["m"], particles1["m"]) + assert particlesRead4["t"] == 0.0 + + # call function to be tested + logName = "simNameTest5" + cfg = configparser.ConfigParser() + cfg["EXPORTS"] = {"exportParticleProperties": "x|m"} + cfg["TRACKPARTICLES"] = {"trackParticles": True, "particleProperties": "iCell"} + particles1["ux"] = np.asarray([1.0, 2.0, 3.0]) + particles1["uy"] = np.asarray([1.0, 4.0, 5.0]) + particles1["uz"] = np.asarray([10.0, 11.0, 11.0]) + particles1["iCell"] = np.asarray([10.0, 11.0, 11.0]) + particles2["ux"] = np.asarray([1.0, 2.0, 3.0]) + particles2["uy"] = np.asarray([1.0, 4.0, 5.0]) + particles2["uz"] = np.asarray([10.0, 11.0, 11.0]) + particles2["iCell"] = np.asarray([10.0, 11.0, 11.0]) + particles1["z"] = np.asarray([1.0, 2.0, 3.0]) + particles2["z"] = np.asarray([1.0, 4.0, 5.0]) + particles1["h"] = np.asarray([1.0, 2.0, 3.0]) + particles2["h"] = np.asarray([1.0, 4.0, 5.0]) + com1DFA.savePartToPickle(particles1, outDir, logName, cfg=cfg) + + # read pickle + picklePath5 = outDir / "particles_simNameTest5_0000.0000.pickle" + particlesRead5 = pickle.load(open(picklePath5, "rb")) + + assert np.array_equal(particlesRead5["x"], particles1["x"]) + assert "y" in particlesRead5.keys() + assert "ux" in particlesRead5.keys() + assert np.array_equal(particlesRead5["iCell"], particles1["iCell"]) + assert np.array_equal(particlesRead5["m"], particles1["m"]) + assert particlesRead5["t"] == 0.0 + + # call function to be tested + logName = "simNameTest6" + cfg = configparser.ConfigParser() + cfg["EXPORTS"] = {"exportParticleProperties": "x|m|hallo"} + cfg["TRACKPARTICLES"] = {"trackParticles": False} + + with pytest.raises(AttributeError) as e: + com1DFA.savePartToPickle(particles1, outDir, logName, cfg=cfg) + assert ("These particle properties are not available") in str(e.value) + + # call function to be tested + logName = "simNameTest7" + cfg = configparser.ConfigParser() + cfg["EXPORTS"] = {"exportParticleProperties": ""} + cfg["TRACKPARTICLES"] = {"trackParticles": False} + com1DFA.savePartToPickle(particles1, outDir, logName, cfg=cfg) + # read pickle + picklePath7 = outDir / "particles_simNameTest7_0000.0000.pickle" + particlesRead7 = pickle.load(open(picklePath7, "rb")) + + for pProp in particlesRead7: + assert pProp in ['ux', 'uy', 'uz', 'iCell', 'z', 'x', 'y', 'm', 'h', 't'] + def test_exportFields(tmp_path): """test exporting fields to ascii files"""