[com1]: make deposition erodable (initially)

avaframe/com1DFA/DFAfunctionsCython.pyx

@@ -1393,6 +1393,8 @@ def updateFieldsC(cfg, particles, dem, fields):
             PKE[j, i] = kineticEnergy[j, i]
 
       # thickness change due to detrainment, stopping and entrainment
+      # MassDetBilinear, MassStopBilinear <= 0, MassEntBilinear >= 0,
+      # FTDetBilinear, FTStopBilinear >= 0, FTEntBilinear <= 0
       FTDetBilinear[j, i] = - MassDetBilinear[j, i] / (areaRaster[j, i] * rho)  # / m * FTBilinear[j, i]
       FTStopBilinear[j, i] = - MassStopBilinear[j, i] / (areaRaster[j, i] * rho)  # / m * FTBilinear[j, i] 
       FTEntBilinear[j, i] = - MassEntBilinear[j, i] / (areaRaster[j, i] * rhoEnt)
@@ -1409,6 +1411,7 @@ def updateFieldsC(cfg, particles, dem, fields):
   fields['FTStop'] = np.asarray(FTStopBilinear).copy()
   fields['FTDet'] = np.asarray(FTDetBilinear).copy()
   fields['FTEnt'] = np.asarray(FTEntBilinear).copy()
+  fields['mStop'] = np.asarray(MassStopBilinear).copy()
 
   if computeP:
     fields['ppr'] = np.asarray(PP).copy()

avaframe/com1DFA/com1DFA.py

@@ -38,6 +38,7 @@
 import avaframe.com1DFA.DFAToolsCython as DFAtllsC
 import avaframe.com1DFA.damCom1DFA as damCom1DFA
 import avaframe.in2Trans.rasterUtils as IOf
+import avaframe.in2Trans.transformFields as transField
 import avaframe.in3Utils.fileHandlerUtils as fU
 from avaframe.in3Utils import cfgUtils
 import avaframe.out3Plot.outDebugPlots as debPlot
@@ -552,7 +553,7 @@ def setThickness(cfg, lineTh, typeTh):
     if cfg["GENERAL"].getboolean(thFlag):
         if cfg["INPUT"]["thFromIni"] != "" and typeTh in cfg["INPUT"]["thFromIni"]:
             lineTh["thicknessSource"] = ["ini file"] * len(lineTh["thickness"])
-        elif cfg["GENERAL"].getboolean("timeDependentRelease"):
+        elif cfg["GENERAL"].getboolean("timeDependentRelease") and typeTh == "relTh":
             lineTh["thicknessSource"] = ["csv file"] * len(lineTh["thickness"])
         else:
             lineTh["thicknessSource"] = ["shp file"] * len(lineTh["thickness"])
@@ -561,7 +562,7 @@ def setThickness(cfg, lineTh, typeTh):
 
     # set thickness value info read from ini file that has been updated from shp or ini previously
     for count, id in enumerate(lineTh["id"]):
-        if cfg["GENERAL"].getboolean("timeDependentRelease"):
+        if cfg["GENERAL"].getboolean("timeDependentRelease") and typeTh == "relTh":
             lineTh["thickness"][count] = float(lineTh["thickness"][count].item())
 
         elif cfg["GENERAL"].getboolean(thFlag):
@@ -1314,7 +1315,7 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     # initialize entrainment and resistance
     # get info of simType and whether or not to initialize resistance and entrainment
     simTypeActual = cfgGen["simTypeActual"]
-    entrMassRaster, entrEnthRaster, reportAreaInfo = initializeMassEnt(
+    entrMassRaster, entrEnthRaster, entrDepthRaster, reportAreaInfo = initializeMassEnt(
         dem,
         simTypeActual,
         inputSimLines["entLine"],
@@ -1326,6 +1327,12 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     # check if entrainment and release overlap
     entrMassRaster = geoTrans.checkOverlap(entrMassRaster, relRaster, "Entrainment", "Release", crop=True)
     entrEnthRaster = geoTrans.checkOverlap(entrEnthRaster, relRaster, "Entrainment", "Release", crop=True)
+
+    if cfg["GENERAL"].getboolean("entrainableDeposition"):
+        demNotErodableRaster = dem["rasterData"] - entrDepthRaster
+    else:
+        demNotErodableRaster = dem["rasterData"]
+
     # check for overlap with the secondary release area
     if secondaryReleaseInfo["flagSecondaryRelease"] == "Yes":
         for secIndex, secRelRaster in enumerate(secondaryReleaseInfo["rasterData"]):
@@ -1385,6 +1392,8 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     # surfacic entrainment mass available (unit kg/m²)
     fields["entrMassRaster"] = entrMassRaster
     fields["entrEnthRaster"] = entrEnthRaster
+    fields["entrDepth"] = entrDepthRaster
+    fields["demNotErodableRaster"] = demNotErodableRaster
     entreainableMass = np.nansum(fields["entrMassRaster"] * dem["areaRaster"])
     log.info("Mass available for entrainment: %.2f kg" % (entreainableMass))
 
@@ -1963,9 +1972,14 @@ def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPoin
     -------
     entrMassRaster : 2D numpy array
         raster of available mass for entrainment
+    entrEnthRaster: 2D numpy array
+        raster of available entrainment enthalpie
+    entrDepthDict["rasterData"]: 2D numpy array
+        raster of available depth for entrainment
     reportAreaInfo: dict
         simulation area information dictionary completed with entrainment area info
     """
+
     # read dem header
     header = dem["originalHeader"]
     ncols = header["ncols"]
@@ -1976,26 +1990,29 @@ def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPoin
         log.info("Entrainment area features: %s" % (entLine["Name"]))
         if entLine["initializedFrom"] == "shapefile":
             entLine = geoTrans.prepareArea(entLine, dem, thresholdPointInPoly, thList=entLine["thickness"])
-        entrMassRaster = entLine["rasterData"]
+        entrThicknessRaster = entLine["rasterData"]
         # ToDo: not used in samos but implemented
         # tempRaster = cfg['GENERAL'].getfloat('entTempRef') + (dem['rasterData'] - cfg['GENERAL'].getfloat('entMinZ'))
         # * cfg['GENERAL'].getfloat('entTempGrad')
         # entrEnthRaster = np.where(tempRaster < 0, tempRaster*cfg['GENERAL'].getfloat('cpIce'),
         #                           tempRaster*cfg['GENERAL'].getfloat('cpWtr')/cfg['GENERAL'].getfloat('hFusion'))
         entrEnthRaster = np.where(
-            entrMassRaster > 0,
+            entrThicknessRaster > 0,
             cfg["GENERAL"].getfloat("entTempRef") * cfg["GENERAL"].getfloat("cpIce"),
             0,
         )
         reportAreaInfo["entrainment"] = "Yes"
     else:
-        entrMassRaster = np.zeros((nrows, ncols))
+        entrThicknessRaster = np.zeros((nrows, ncols))
         entrEnthRaster = np.zeros((nrows, ncols))
         reportAreaInfo["entrainment"] = "No"
 
-    entrMassRaster = entrMassRaster * cfg["GENERAL"].getfloat("rhoEnt")
+    entrMassRaster = entrThicknessRaster * cfg["GENERAL"].getfloat("rhoEnt")
+    # convert entrainable thickness into entrainable depth
+    entrThDict = {"rasterData": entrThicknessRaster, "header": dem["header"]}
+    entrDepthDict, _, _ = transField.convertDepthThickness(entrThDict, dem, typeOfInput="thickness")
 
-    return entrMassRaster, entrEnthRaster, reportAreaInfo
+    return entrMassRaster, entrEnthRaster, entrDepthDict["rasterData"], reportAreaInfo
 
 
 def initializeResistance(cfg, dem, simTypeActual, resLine, reportAreaInfo, thresholdPointInPoly):
@@ -3499,6 +3516,20 @@ def prepareVarSimDict(standardCfg, inputSimFiles, variationDict, simNameExisting
                 cfgSim["INPUT"]["entThFile"] = pathToEnt
                 inputSimFiles["entResInfo"]["entRemeshed"] = remeshedEnt
             cfgSim["INPUT"]["entrainmentScenario"] = str(pathlib.Path("ENT", inputSimFiles["entFile"].name))
+            if cfgSim["GENERAL"]["entrainableDeposition"] == "True" and (
+                cfgSim["GENERAL"]["adaptSfcEntrainment"] == "0"
+                or (
+                    cfgSim["GENERAL"]["adaptSfcStopped"] == "0"
+                    and cfgSim["GENERAL"]["adaptSfcDetrainment"] == "0"
+                )
+            ):
+                # set entrainable Deposition to False if entrainment or adapt Sfc are switched off
+                message = f"When entrainable deposition is True, adaptSfcEntrainment and adaptSfcDeposition (or adaptSfcDetrainment) need to be 1."
+                log.error(message)
+                raise ValueError(message)
+        else:
+            cfgSim["GENERAL"]["entrainableDeposition"] = "False"
+            log.info("Deposition is not entrainable when it's not an ent sim type.")
 
         # add info about resistance file path to the cfg
         if "res" in row._asdict()["simTypeList"] and inputSimFiles["resFile"] is not None:
@@ -3907,7 +3938,15 @@ def adaptDEM(dem, fields, cfg):
     dem: dict
         dictionary with info on DEM data
     fields : dict
-        fields dictionary
+        fields dictionary containing:
+            FTDet (flow thickness change due to detrainment >= 0)
+            FTStop (flow thickness change due to stopping >= 0)
+            FTEnt (flow thickness change due to entrainment <= 0)
+            mStop (stopped mass <=0)
+            sfcChangeTotal (so far changed surface)
+            entrDepth (entrainable depth that is on top of the not erodable topography)
+            entrMassRaster (entrainable mass)
+            demNotErodableRaster (topography that is not erodable)
     cfg: dict
         configuration settings
 
@@ -3916,17 +3955,23 @@ def adaptDEM(dem, fields, cfg):
     dem: dict
         dictionary with info on DEM data containing adapted topography
     fields : dict
-        fields dictionary containing adapted DEM
+        fields dictionary containing adapted DEM, adapted entrainable mass and depth raster
     """
 
     ZDEM = dem["rasterData"].copy()
-    FTDet = fields["FTDet"]
-    FTStop = fields["FTStop"]
-    FTEnt = fields["FTEnt"]
+    FTDet = fields["FTDet"]  # >= 0
+    FTStop = fields["FTStop"]  # >= 0
+    FTEnt = fields["FTEnt"]  # <= 0
+    mStop = fields["mStop"]  # <= 0
+
     sfcChangeTotal = fields["sfcChangeTotal"]
     sfcChange = np.zeros_like(FTDet)
     ZDEMadapt = ZDEM
 
+    depthStop = np.zeros_like(ZDEM)
+    depthDet = np.zeros_like(ZDEM)
+    depthEnt = np.zeros_like(ZDEM)
+
     _, _, NzNormed = DFAtls.normalize(dem["Nx"].copy(), dem["Ny"].copy(), dem["Nz"].copy())
 
     if cfg.getboolean("adaptSfcStopped"):
@@ -3944,6 +3989,14 @@ def adaptDEM(dem, fields, cfg):
         depthEnt = FTEnt / NzNormed
         ZDEMadapt += depthEnt
         sfcChange += depthEnt
+    if cfg.getboolean("entrainableDeposition"):
+        hErodable = fields["entrDepth"]
+        hErodable += depthStop
+        hErodable += depthDet
+        hErodable += depthEnt
+        ZDEMadapt = fields["demNotErodableRaster"] + hErodable
+        fields["entrDepth"] = hErodable
+        fields["entrMassRaster"] -= mStop
 
     dem["rasterData"] = ZDEMadapt
     fields["demAdapted"] = ZDEMadapt

avaframe/com1DFA/com1DFACfg.ini

@@ -11,7 +11,7 @@ simTypeList = available
 modelType = dfa
 
 #+++++++++++++ Output++++++++++++
-# desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA, dmDet, sfcChange, demAdapted, timeInfo, particles) - separated by |
+# desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA, dmDet, sfcChange, demAdapted, timeInfo, entrDepth, particles) - separated by |
 resType = ppr|pft|pfv|timeInfo
 # saving time step, i.e.  time in seconds (last time step is always saved; initial time step only if explicitly specified)
 # 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)
@@ -489,6 +489,9 @@ entDefResistance = 0
 adaptSfcStopped = 0
 adaptSfcDetrainment = 0
 adaptSfcEntrainment = 0
+# deposition is only entrainable if it is an entrainment simulation and if adaptSfcStopped = 1 and adaptSfcEntrainment = 1
+# then mass that is stopped can be entrained (also if the entrainment area is not in this stopped area)
+entrainableDeposition = False
 # only adapt topography if changing height in at least one cell is > thresholdAdaptSfc [m]
 thresholdAdaptSfc = 0.1
 # use the adapted topography as background in the report plots

avaframe/com1DFA/deriveParameterSet.py

@@ -154,6 +154,7 @@ def checkResType(fullCfg, section, key, value):
             "sfcChange",
             "demAdapted",
             "timeInfo",
+            "entrDepth",
         ]
         message = "The parameter % s is not a valid resType. It will not be saved"
         newResType = []
@@ -285,7 +286,7 @@ def getThicknessValue(cfg, inputSimFiles, fName, thType):
 
     # if thickness should be read from shape file
     if cfg["GENERAL"].getboolean(thFlag):
-        if cfg["GENERAL"].getboolean("timeDependentRelease"):
+        if cfg["GENERAL"].getboolean("timeDependentRelease") and thType == "relTh":
             cfg["INPUT"]["relThInfo"] = "timeDependent"
         # if at least one but not all features in a shapefile have a thickness value - error
         elif ("None" in thicknessList) and thType != "entTh":
@@ -902,8 +903,11 @@ def appendThicknessToCfg(cfg):
             for count, id in enumerate(idList):
                 thNameId = thType + id
                 if thNameId in cfg["GENERAL"].keys():
-                    log.info("Thickness value for %s already set in initial config file, \
-                              read from there not from shp file" % thNameId)
+                    log.info(
+                        "Thickness value for %s already set in initial config file, \
+                              read from there not from shp file"
+                        % thNameId
+                    )
                 else:
                     cfgGen[thNameId] = str(float(thicknessList[count]))
 

avaframe/in2Trans/transformFields.py

@@ -1,5 +1,5 @@
 """
-    ASCII file reader and handler
+ASCII file reader and handler
 
 """
 
@@ -16,14 +16,14 @@
 log = logging.getLogger(__name__)
 
 
-def convertDepthThickness(inputDict, demDict, typeOfInput='depth'):
+def convertDepthThickness(inputDict, demDictOrig, typeOfInput="depth"):
     """convert depthField to thicknessField or the other way around depending on type, using a DEM to compute the slope
     angle required for transformation
     also writes field to a new directory called transformed where depthFile/thicknessFile is located
 
     Parameters
     -----------
-    demDict: dict
+    demDictOrig: dict
         dictionary with dem header and rasterData (numpy nd array of z values)
     inputDict: dict
         dictionary with depthField/thicknessField header and rasterData (numpy nd array of depth values)
@@ -37,6 +37,7 @@ def convertDepthThickness(inputDict, demDict, typeOfInput='depth'):
 
     """
     # get normal vector of the grid mesh
+    demDict = {"header": demDictOrig["header"].copy(), "rasterData": demDictOrig["rasterData"]}
     demDict = gT.getNormalMesh(demDict)
     _, _, NzNormed = DFAtls.normalize(demDict["Nx"], demDict["Ny"], demDict["Nz"])
 
@@ -45,41 +46,41 @@ def convertDepthThickness(inputDict, demDict, typeOfInput='depth'):
     demDict["header"]["nodata_value"] = inputDict["header"]["nodata_value"]
 
     # multiply depth with cos(slopeAngle)
-    if typeOfInput == 'depth':
+    if typeOfInput == "depth":
         outField = inputRasterNew * NzNormed
-    elif typeOfInput == 'thickness':
+    elif typeOfInput == "thickness":
         outField = inputRasterNew / NzNormed
     else:
-        message = 'Type for thickness/depth conversion is: %s - not valid' % typeOfInput
+        message = "Type for thickness/depth conversion is: %s - not valid" % typeOfInput
         log.error(message)
         raise AssertionError(message)
 
     slopeAngleField = np.rad2deg(np.arccos(NzNormed))
 
     # create thickness dict
-    outputDict = {'header': demDict['header'], 'rasterData': outField}
+    outputDict = {"header": demDict["header"], "rasterData": outField}
 
     return outputDict, inputRasterNew, slopeAngleField
 
 
-def fetchPointValuesFromField(dataDF, xyPoints, resType, interpMethod='bilinear'):
-    """ derive field values at xyPoints using a interpMethod (options: nearest and bilinear)
+def fetchPointValuesFromField(dataDF, xyPoints, resType, interpMethod="bilinear"):
+    """derive field values at xyPoints using a interpMethod (options: nearest and bilinear)
 
-        Parameters
-        -----------
-        dataDF: pandas dataFrame
-            dataframe with info on simulations (including result file paths,configuration)
-        xyPoints: dict
-            dictionary with keys x, y and point of interest coordinates
-        resType: str
-            name of result type (pfv, pft, ...)
-        interpMethod: str
-            interpolation method to derive values of field at xyPoints
+    Parameters
+    -----------
+    dataDF: pandas dataFrame
+        dataframe with info on simulations (including result file paths,configuration)
+    xyPoints: dict
+        dictionary with keys x, y and point of interest coordinates
+    resType: str
+        name of result type (pfv, pft, ...)
+    interpMethod: str
+        interpolation method to derive values of field at xyPoints
 
-        Returns
-        ---------
-        dataDF: pandas dataFrame
-            updated pandas dataFrame with new column pointValues_resType providing value of resType at xyPoint
+    Returns
+    ---------
+    dataDF: pandas dataFrame
+        updated pandas dataFrame with new column pointValues_resType providing value of resType at xyPoint
 
     """
 
@@ -88,9 +89,9 @@ def fetchPointValuesFromField(dataDF, xyPoints, resType, interpMethod='bilinear'
         # read field
         field = IOf.readRaster(row[resType])
 
-        value, _ = gT.projectOnRaster(field, xyPoints, interp=interpMethod, inData="rasterData", outData="value")
-        dataDF.loc[index, ('pointValues_%s' % resType)] = value['value'][0]
+        value, _ = gT.projectOnRaster(
+            field, xyPoints, interp=interpMethod, inData="rasterData", outData="value"
+        )
+        dataDF.loc[index, ("pointValues_%s" % resType)] = value["value"][0]
 
     return dataDF
-
-