From e32bf4b7a3b7c635a6cf7bdfe30f680c5ac62201 Mon Sep 17 00:00:00 2001 From: Anna Wirbel Date: Wed, 3 Dec 2025 11:38:17 +0100 Subject: [PATCH] (com1) add timeInfo resType add info in docs add discrete colormap option and units to colorbar add pytests fix small issues add label to plot colorbar update docu update pytest remove unused import add timeInfo to resTpyes --- avaframe/com1DFA/com1DFA.py | 4 + avaframe/com1DFA/com1DFACfg.ini | 4 +- avaframe/com1DFA/com1DFATools.py | 24 ++++ avaframe/com1DFA/deriveParameterSet.py | 13 +++ .../Inputs/simiSol_com1DFACfg.ini | 2 +- avaframe/out1Peak/outPlotAllPeak.py | 16 ++- avaframe/out3Plot/plotUtils.py | 28 +++++ avaframe/out3Plot/plotUtilsCfg.ini | 6 +- .../data/testSimiSol/simiSol_com1DFACfg.ini | 108 ------------------ avaframe/tests/test_com1DFA.py | 7 +- avaframe/tests/test_com1DFATools.py | 24 ++++ avaframe/tests/test_deriveParameterSet.py | 9 +- avaframe/tests/test_simiSol.py | 20 +++- docs/moduleCom1DFA.rst | 4 +- 14 files changed, 142 insertions(+), 127 deletions(-) delete mode 100644 avaframe/tests/data/testSimiSol/simiSol_com1DFACfg.ini diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index 1d01e5406..9285a99b7 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -1681,6 +1681,7 @@ def initializeFields(cfg, dem, particles, releaseLine): 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] + fields["timeInfo"] = np.zeros((nrows, ncols)) # first time # 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): @@ -2674,6 +2675,9 @@ def computeEulerTimeStep( particles = DFAfunC.computeTrajectoryAngleC(particles, zPartArray0) particles, fields = DFAfunC.updateFieldsC(cfg, particles, dem, fields) + # update field that indicates when cell was first affected by mass + fields = com1DFATools.updateTimeField(fields, particles["t"]) + # adapt DEM considering erosion and deposition # only adapt DEM when in one grid cell the changing height > threshold thresholdAdaptSfc = cfg.getfloat("thresholdAdaptSfc") diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index e2b3cd77d..c9e9cd5e9 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -11,8 +11,8 @@ simTypeList = available modelType = dfa #+++++++++++++ Output++++++++++++ -# desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA, dmDet, sfcChange, demAdapted, particles) - separated by | -resType = ppr|pft|pfv +# desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA, dmDet, sfcChange, demAdapted, timeInfo, particles) - separated by | +resType = ppr|pft|pfv|timeInfo # 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) diff --git a/avaframe/com1DFA/com1DFATools.py b/avaframe/com1DFA/com1DFATools.py index 5fde3a4a7..a1b52b2b5 100644 --- a/avaframe/com1DFA/com1DFATools.py +++ b/avaframe/com1DFA/com1DFATools.py @@ -481,3 +481,27 @@ def chooseDemPlot(dem, adaptedDemBackground=False): demPlot["rasterData"] = dem["originalRasterData"] log.info("The original DEM is used in the plots.") return demPlot + + +def updateTimeField(fields, timeStep): + """update filed indicating first time step mass entered a cell + + Parameters + ----------- + fields: dict + dictionary with fields + timeStep: float + actual time step + + Returns + --------- + fields: dict + updated timeInfo field + """ + + FT = fields["FT"] + + # set time step to previously not affected cells + fields["timeInfo"] = np.where(((fields["timeInfo"] == 0) & (FT != 0)), timeStep, fields["timeInfo"]) + + return fields diff --git a/avaframe/com1DFA/deriveParameterSet.py b/avaframe/com1DFA/deriveParameterSet.py index c463fa5f3..c6848a722 100644 --- a/avaframe/com1DFA/deriveParameterSet.py +++ b/avaframe/com1DFA/deriveParameterSet.py @@ -153,6 +153,7 @@ def checkResType(fullCfg, section, key, value): "FTDet", "sfcChange", "demAdapted", + "timeInfo", ] message = "The parameter % s is not a valid resType. It will not be saved" newResType = [] @@ -411,6 +412,18 @@ def checkThicknessSettings(cfg, thName, inputSimFiles): ) log.error(message) raise AssertionError(message) + # Check: If raster input file but thFromFile = False, error + elif ( + not cfg["GENERAL"].getboolean(thFlag) + and inputSimFiles["entResInfo"][thName + "FileType"] in [".asc", ".tif"] + and inputSimFiles["entResInfo"]["flag" + nameTypes[thName]] == "Yes" + ): + message = "If %s input file is of type .asc or .tif - %s needs to be set to True" % ( + nameStrings[thName], + thFlag, + ) + log.error(message) + raise AssertionError(message) else: message = "Check %s - needs to be True or False" % thFlag log.error(message) diff --git a/avaframe/data/avaSimilaritySol/Inputs/simiSol_com1DFACfg.ini b/avaframe/data/avaSimilaritySol/Inputs/simiSol_com1DFACfg.ini index e3e16cea2..9075d8989 100644 --- a/avaframe/data/avaSimilaritySol/Inputs/simiSol_com1DFACfg.ini +++ b/avaframe/data/avaSimilaritySol/Inputs/simiSol_com1DFACfg.ini @@ -19,7 +19,7 @@ initPartDistType = random #+++++++++SNOW properties # True if release thickness should be read from shapefile file; if False - relTh read from ini file -relThFromFile = False +relThFromFile = True #++++++++++++Time stepping parameters # End time [s] diff --git a/avaframe/out1Peak/outPlotAllPeak.py b/avaframe/out1Peak/outPlotAllPeak.py index a54151eb8..f3e392726 100644 --- a/avaframe/out1Peak/outPlotAllPeak.py +++ b/avaframe/out1Peak/outPlotAllPeak.py @@ -10,7 +10,6 @@ matplotlib.use("agg") from matplotlib import pyplot as plt import pathlib -from mpl_toolkits.axes_grid1 import make_axes_locatable import avaframe.out3Plot.plotUtils as pU import avaframe.in1Data.getInput as gI @@ -306,9 +305,18 @@ def addConstrainedDataField(fileName, resType, demField, ax, cellSize, alpha=1.0 if len(np.nonzero(data)[0]) > 0.0: # add Colorbar fig = ax.get_figure() - divider = make_axes_locatable(ax) - cax = divider.append_axes("right", size="5%", pad=0.05) - fig.colorbar(im1, cax=cax) + cax = ax.inset_axes([1.05, 0.075, 0.05, 0.925], transform=ax.transAxes) + fig.colorbar(im1, ax=ax, cax=cax) + cax2 = ax.inset_axes([1.0, 0.0, 0.05, 0.075], transform=ax.transAxes) + plt.text( + 0.5, + 0.5, + ("[%s]" % unit), + horizontalalignment="left", + verticalalignment="center", + transform=cax2.transAxes, + ) + cax2.set_visible(False) return ax, rowsMinPlot, colsMinPlot, extentCellCorners diff --git a/avaframe/out3Plot/plotUtils.py b/avaframe/out3Plot/plotUtils.py index d0987c692..5396d361a 100644 --- a/avaframe/out3Plot/plotUtils.py +++ b/avaframe/out3Plot/plotUtils.py @@ -214,6 +214,31 @@ ] cmapSfcChange = copy.copy(cmapCrameri.nuuk.reversed()) +# colormap for timeInfo +levtimeInfo = list(fU.splitIniValueToArraySteps(cfgPlotUtils["timeInfoColorLevels"])) +# lipari reversed color map +colorstimInfo = [ + "#fdf5da", + "#f2debb", + "#e9c99f", + "#e5b58a", + "#e7a37a", + "#ea906d", + "#e57b62", + "#d46b5e", + "#bc6461", + "#a36267", + "#8e616c", + "#7c6071", + "#6b5f76", + "#5b5d79", + "#45587a", + "#294b70", + "#13385a", + "#082540", + "#031326", +] +cmaptimeInfo = copy.copy(cmapCrameri.lipari.reversed()) ############################################### # Set colormaps to use @@ -237,6 +262,8 @@ cmapSfcChange = {"cmap": cmapSfcChange, "colors": colorsSfcC, "levels": levSfcC} +cmapTime = {"cmap": cmaptimeInfo, "colors": colorstimInfo, "levels": levtimeInfo} + # for zdelta # Remark FSO: the try except comes from cmcrameri v1.5 not having lipari, but it is still # widely used (Okt 2024). TODO: remove in future versions @@ -269,6 +296,7 @@ "dmDet": cmapDmDet, "demAdapted": cmapGreys, "sfcChange": cmapSfcChange, + "timeInfo": cmapTime, } cmapDEM = cmapGreys diff --git a/avaframe/out3Plot/plotUtilsCfg.ini b/avaframe/out3Plot/plotUtilsCfg.ini index 8371e55a0..5d085f30d 100644 --- a/avaframe/out3Plot/plotUtilsCfg.ini +++ b/avaframe/out3Plot/plotUtilsCfg.ini @@ -86,6 +86,7 @@ unitzdelta = m unitdmDet = kg unitsfcChange = m unitdemAdapted = m +unittimeInfo = s # threshold levels elevMaxppr = 100 elevMaxpft = 1 @@ -112,6 +113,8 @@ travelAngleColorLevels = 28|29|30|31|32|33|34 probaColorLevels = 0|0.25|0.50|0.75|1. # for thicknessChange surfaceChangeLevels = -1|-0.3|-0.1|0|0.1|0.3|1 +# for time Info +timeInfoColorLevels = 10|20|30|40|50|60|70|80|90|100|110|120|130|140|150|160|170|180|200 # contour levels (when adding contour lines on a plot) contourLevelsppr = 1|3|5|10|25|50|100|250 contourLevelspft = 0.1|0.25|0.5|0.75|1 @@ -144,4 +147,5 @@ nameVy = y velocity nameVz = z velocity nameTA = travel angle namezdelta = energy line height -namedmDet = detrained mass \ No newline at end of file +namedmDet = detrained mass +nametimeInfo = firstTime \ No newline at end of file diff --git a/avaframe/tests/data/testSimiSol/simiSol_com1DFACfg.ini b/avaframe/tests/data/testSimiSol/simiSol_com1DFACfg.ini deleted file mode 100644 index d83c27c21..000000000 --- a/avaframe/tests/data/testSimiSol/simiSol_com1DFACfg.ini +++ /dev/null @@ -1,108 +0,0 @@ -### Config File - This file contains the main settings for the similarity -### solution run with com1DFAPy - - -[GENERAL] -#+++++++++++++ Output++++++++++++ -# desired result Parameters (ppr, pft, pfv, FT, FV, P, particles) - separated by | -resType = FT|FV|Vx|Vy|Vz -tSteps = 0:1 - -#+++++++++SNOW properties -# True if release thickness should be read from shapefile file; if False - relTh read from ini file -relThFromFile = False - -#++++++++++++Time stepping parameters -tEnd = 5 -# to use a variable time step (time step depends on kernel radius) -sphKernelRadiusTimeStepping = True -# courant number -cMax = 0.04 - -#+++++++++++++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 = 8 -# Choice of artificial viscosity -# 0) No artificial viscosity -# 1) SAMOS artificial viscosity -# 2) Ata artificial viscosity -viscOption = 1 - -#++++++++++++++++ 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 -# reference kernel radius -sphKR0 = 5 -# reference number of particles per kernel radius -nPPK0 = 15|20 -# variation of appK exponent -aPPK = -0.5|-1 - -# threshold for splitting particles, split if: mPart > massPerPart x thresholdMassSplit -thresholdMassSplit = 5 - - -#+++++++++++++Mesh and interpolation -# remesh the input DEM -# expected mesh size [m] -meshCellSize = 8 - -#+++++++++++++Flow model parameters+++++++++ -# subgridMixingFactor -subgridMixingFactor = 10 - -#++++++++++++Friction model -# friction type (samosAT, Coulomb) -frictModel = Coulomb -# add the friction using an explicit formulation (1) -# 0 use an implicit method -explicitFriction = 1 -#+++++++++++++SamosAt friction model -mu = 0.466307658 - -[SIMISOL] -# which time step should be saved -tSave = 5 -# dimensioning parameters -L_x = 80. -L_y = 80. - -# release thickness -relTh = 4. - -# bed friction angle -bedFrictionAngle = 25. -# internal friction angle -internalFrictionAngle = 25. -# plane inclination angle -planeinclinationAngle = 35. - -# Flag earthpressure coefficients -# if false takes 1 -flagEarth = False - -# save analysis plots at time step dtSol -dtSol = 50. - -#++++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 = True -# 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.02 diff --git a/avaframe/tests/test_com1DFA.py b/avaframe/tests/test_com1DFA.py index feec8bffd..ee2fa9713 100644 --- a/avaframe/tests/test_com1DFA.py +++ b/avaframe/tests/test_com1DFA.py @@ -1915,13 +1915,14 @@ def test_initializeFields(): # print("compute TA", fields["computeTA"]) # print("compute P", fields["computeP"]) - assert len(fields) == 23 + assert len(fields) == 24 assert fields["computeTA"] is False assert fields["computeKE"] is False assert fields["computeP"] assert np.sum(fields["pfv"]) == 0.0 assert np.sum(fields["pta"]) == 0.0 assert np.sum(fields["ppr"]) == 0.0 + assert np.sum(fields["pke"]) == 0.0 assert np.sum(fields["FV"]) == 0.0 assert np.sum(fields["P"]) == 0.0 assert np.sum(fields["TA"]) == 0.0 @@ -1937,6 +1938,8 @@ def test_initializeFields(): assert np.sum(fields["FTDet"]) == 0.0 assert np.sum(fields["FTStop"]) == 0.0 assert np.sum(fields["FTEnt"]) == 0.0 + assert np.sum(fields["timeInfo"]) == 0.0 + assert np.sum(fields["sfcChangeTotal"]) == 0.0 cfg["REPORT"] = {"plotFields": "pft|pfv"} cfg["GENERAL"] = { @@ -1947,7 +1950,7 @@ def test_initializeFields(): } # call function to be tested particles, fields = com1DFA.initializeFields(cfg, dem, particles, "") - assert len(fields) == 23 + assert len(fields) == 24 assert fields["computeTA"] assert fields["computeKE"] assert fields["computeP"] is False diff --git a/avaframe/tests/test_com1DFATools.py b/avaframe/tests/test_com1DFATools.py index d34df7867..98cc5074e 100644 --- a/avaframe/tests/test_com1DFATools.py +++ b/avaframe/tests/test_com1DFATools.py @@ -150,3 +150,27 @@ def test_updateResCoeffFields(tmp_path): assert fields["detRaster"][0, 0] == 3 assert fields["cResRaster"][2, 3] == 0 assert fields["detRaster"][2, 3] == 3 + + +def test_updateTimeField(): + """test updating the field of timeInfo""" + + nrows = 5 + ncols = 6 + + timeInfo = np.zeros((nrows, ncols)) + timeInfo[0, 1:3] = 1.0 + + FT = np.zeros((nrows, ncols)) + FT[0, 1:3] = 0.5 + FT[1, 2:4] = 2.5 + timeStep = 5.0 + fields = {"timeInfo": timeInfo, "FT": FT} + + fields = com1DFATools.updateTimeField(fields, timeStep) + + assert fields["timeInfo"][0, 1] == 1.0 + assert fields["timeInfo"][0, 2] == 1.0 + assert fields["timeInfo"][1, 2] == 5.0 + assert fields["timeInfo"][1, 3] == 5.0 + assert not fields["timeInfo"][2:6, :].any() diff --git a/avaframe/tests/test_deriveParameterSet.py b/avaframe/tests/test_deriveParameterSet.py index 289819d38..a6d6ab749 100644 --- a/avaframe/tests/test_deriveParameterSet.py +++ b/avaframe/tests/test_deriveParameterSet.py @@ -391,11 +391,14 @@ def test_checkThicknessSettings(): thName = "relTh" - thicknessSettingsCorrect = dP.checkThicknessSettings(cfg, thName, inputSimFiles) - - assert thicknessSettingsCorrect + with pytest.raises(AssertionError) as e: + assert dP.checkThicknessSettings(cfg, "relTh", inputSimFiles) + assert str(e.value) == ( + "If Release area input file is of type .asc or .tif - relThFromFile needs to be set to True" + ) cfg["GENERAL"]["relThRangeVariation"] = "50$4" + cfg["GENERAL"]["relThFromFile"] = "True" with pytest.raises(AssertionError) as e: assert dP.checkThicknessSettings(cfg, "relTh", inputSimFiles) diff --git a/avaframe/tests/test_simiSol.py b/avaframe/tests/test_simiSol.py index 0daa4cc5c..aeed221ec 100644 --- a/avaframe/tests/test_simiSol.py +++ b/avaframe/tests/test_simiSol.py @@ -18,9 +18,7 @@ def test_mainCompareSimSolCom1DFA(tmp_path): dirname = pathlib.Path(__file__).parents[0] - simiSolCfg = ( - dirname / ".." / "tests" / "data" / "testSimiSol" / "simiSol_com1DFACfg.ini" - ) + simiSolCfg = dirname / ".." / "data" / "avaSimilaritySol" / "Inputs" / "simiSol_com1DFACfg.ini" sourceDir = dirname / ".." / "data" / "avaSimilaritySol" / "Inputs" destDir = tmp_path / "avaSimilaritySol" / "Inputs" avalancheDir = tmp_path / "avaSimilaritySol" @@ -33,15 +31,27 @@ def test_mainCompareSimSolCom1DFA(tmp_path): cfgMain = cfgUtils.getGeneralConfig() cfgMain['MAIN']['avalancheDir'] = str(avalancheDir) cfg = cfgUtils.getModuleConfig(com1DFA, simiSolCfg) + # adjust settings for faster computation times + cfg["GENERAL"]["cMax"] = "0.04" + cfg["GENERAL"]["sphKernelRadius"] = "8" + cfg["GENERAL"]["nPPK0"] = "15|20" + cfg["GENERAL"]["aPPK"] = "-0.5|-1" + cfg["GENERAL"]["meshCellSize"] = "8" + + # write updated cfg to file + cfgFile = pathlib.Path(destDir, "%s.ini" % ("simiSolUpdated_com1DFACfg")) + with open(cfgFile, "w") as conf: + cfg.write(conf) + # Define release thickness distribution demFile = gI.getDEMPath(avalancheDir) relDict = simiSolTest.getReleaseThickness(avalancheDir, cfg, demFile) # call com1DFA to perform simulations - provide configuration file and release thickness function # (may be multiple sims) - _, _, _, simDF = com1DFA.com1DFAMain(cfgMain, cfgInfo=simiSolCfg) + _, _, _, simDF = com1DFA.com1DFAMain(cfgMain, cfgInfo=cfgFile) simDF, _ = cfgUtils.readAllConfigurationInfo(avalancheDir) - solSimi = simiSolTest.mainSimilaritySol(simiSolCfg) + solSimi = simiSolTest.mainSimilaritySol(cfgFile) simDF = simiSolTest.postProcessSimiSol( avalancheDir, cfgMain, cfg, simDF, solSimi, outDirTest diff --git a/docs/moduleCom1DFA.rst b/docs/moduleCom1DFA.rst index c5f0222e6..77cb853c8 100644 --- a/docs/moduleCom1DFA.rst +++ b/docs/moduleCom1DFA.rst @@ -122,7 +122,8 @@ Release-, entrainment thickness settings Thickness is unambiguous: it is measured normal to the slope. Release, entrainment and secondary release thickness can be specified in two different ways, 1) directly from the provided -input file (shape file or raster file) or 2) through the :py:mod:`com1DFA` configuration file: +input file (shape file or raster file) or 2) through the :py:mod:`com1DFA` configuration file +(only available if input file is of type shape file): 1. **Read from input file**: @@ -306,6 +307,7 @@ The result types that can be chosen to be exported are (all correspond to fields * FTDet - thickness of detrained mass computed based on dmDet / (rho * area of cell) * sfcChange - flow depth that changed the surface topography due to detrainment, stopping and entrainment * demAdapted - adapted DEM considering stopping/ detrainment/ entrainment +* timeInfo - time step at which a cell was first affected by flow * particles (:ref:`com1DFAAlgorithm:Particle properties`) Have a look at the designated subsection Output in ``com1DFA/com1DFACfg.ini``.