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``.