diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 78a389c8d..34fd07549 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -391,9 +391,9 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType particles['uz'] = np.asarray(uzArray) particles['m'] = np.asarray(mass) particles['dmDet'] = np.asarray(dMDet) + particles['dmEnt'] = np.asarray(dM) particles['totalEnthalpy'] = np.asarray(totalEnthalpyArray) - # update mass available for entrainement # TODO: this allows to entrain more mass then available... for k in range(nPart): @@ -466,6 +466,7 @@ cpdef (double, double) computeEntMassAndForce(double dt, double entrMassCell, return dm, areaEntrPart + cpdef double computeDetMass(double dt, double detCell, double areaPart, double uMag): """ compute detrained mass @@ -657,6 +658,10 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): cdef double thresholdProjection = cfg.getfloat('thresholdProjection') cdef double centeredPosition = cfg.getfloat('centeredPosition') cdef int snowSlide = cfg.getint('snowSlide') + cdef int delStoppedParticles = cfg.getint('delStoppedParticles') + cdef int adaptSfcStopped = cfg.getint('adaptSfcStopped') + cdef int adaptSfcDetrainment = cfg.getint('adaptSfcDetrainment') + cdef int adaptSfcEntrainment = cfg.getint('adaptSfcEntrainment') cdef int dissDam = cfg.getint('dissDam') cdef double csz = dem['header']['cellsize'] cdef int nrows = dem['header']['nrows'] @@ -688,6 +693,14 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): cdef double peakKinEne = particles['peakKinEne'] cdef double peakForceSPH = particles['peakForceSPH'] cdef double totKForceSPH = particles['forceSPHIni'] + cdef long long[:] ID = particles['ID'] + # Stopped particles + cdef double[:] xStoppedArray = np.empty(0, dtype=np.float64) + cdef double[:] yStoppedArray = np.empty(0, dtype=np.float64) + cdef double[:] mStoppedArray = np.empty(0, dtype=np.float64) + cdef double[:] idStoppedArray = np.empty(0, dtype=np.float64) + cdef double[:] uMagStoppedArray = np.empty(0, dtype=np.float64) + cdef long[:] indRemoveParticle # read fields cdef double[:] forceX = force['forceX'] cdef double[:] forceY = force['forceY'] @@ -732,14 +745,16 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): cdef double[:] uzArrayNew = np.zeros(nPart, dtype=np.float64) cdef double[:] velocityMagArrayNew = np.zeros(nPart, dtype=np.float64) cdef int[:] keepParticle = np.ones(nPart, dtype=np.int32) + cdef int[:] notStopParticle = np.ones(nPart, dtype=np.int32) # declare intermediate step variables cdef double m, h, x, y, z, sCor, s, l, ux, uy, uz, nx, ny, nz, dtStop, idfixed cdef double mNew, xNew, yNew, zNew, uxNew, uyNew, uzNew, txWall, tyWall, tzWall, totalEnthalpy, totalEnthalpyNew cdef double sCorNew, sNew, lNew, ds, dl, uN, uMag, uMagNew, fNx, fNy, fNz, dv, uMagt0, uMagt1 cdef double ForceDriveX, ForceDriveY, ForceDriveZ - cdef double massEntrained = 0, massDetrained = 0, massFlowing = 0, dissEm = 0 + cdef double massEntrained = 0, massDetrained = 0, massStopped = 0, massFlowing = 0, dissEm = 0 cdef int k, inter cdef int nRemove = 0 + cdef int nStop = 0 # variables for interpolation cdef int Lx0, Ly0, LxNew0, LyNew0, iCell, iCellNew cdef double w[4] @@ -780,10 +795,26 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): else: dtStop = dt + uMagNew = DFAtlsC.norm(uxNew, uyNew, uzNew) + # update mass (already done un computeForceC) mNew = m massEntrained = massEntrained + dM[k] massDetrained = massDetrained + dMDet[k] + + # Stopped particles with velocity = 0 or mass = 0 + if delStoppedParticles == 1: + if uMagNew == 0 or mNew == 0: + xStoppedArray = np.append(xStoppedArray, xArray[k]) + yStoppedArray = np.append(yStoppedArray, yArray[k]) + mStoppedArray = np.append(mStoppedArray, mass[k]) + idStoppedArray = np.append(idStoppedArray, ID[k]) + uMagStoppedArray = np.append(uMagStoppedArray, uMagNew) + massStopped = massStopped + mass[k] + notStopParticle[k] = 0 # particle is deleted + nStop = nStop + 1 + continue + # update position if centeredPosition: xNew = x + dtStop * 0.5 * (ux + uxNew) @@ -960,6 +991,11 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): particles['kineticEne'] = TotkinEneNew particles['potentialEne'] = TotpotEneNew + particles['stoppedParticles']['x'] = np.asarray(xStoppedArray) + particles['stoppedParticles']['y'] = np.asarray(yStoppedArray) + particles['stoppedParticles']['m'] = np.asarray(mStoppedArray) + particles['stoppedParticles']['ID'] = np.asarray(idStoppedArray) + if typeStop == 1: # typeStop = 1 for initialisation step where particles are redistributed to reduce SPH force # here stop criterion based on SPHForce @@ -998,16 +1034,33 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): if stop: particles['iterate'] = False + particles['stoppedParticles'] = particles + massStopped = massStopped + massFlowing + if typeStop == 1: log.debug('stopping initial particle distribution') else: log.debug('stopping because of %s stopCriterion.' % (cfg['stopCritType'])) + particles['massStopped'] = - massStopped + # remove particles that are not located on the mesh any more if nRemove > 0: mask = np.array(np.asarray(keepParticle), dtype=bool) particles = particleTools.removePart(particles, mask, nRemove, 'because they exited the domain', snowSlide=snowSlide) + # remove particles that have mass = 0 or velocity = 0 + if nStop > 0 and delStoppedParticles == 1: + indRemoveParticle = np.array([], dtype=np.int64) + for k in range(len(keepParticle)): + # consider particles that are removed because they exit the domain + if keepParticle[k] == 0: + indRemoveParticle = np.append(indRemoveParticle, int(k)) + if notStopParticle[k] == 0: + nStop = nStop - 1 + notStopParticle = np.delete(notStopParticle, indRemoveParticle) + mask = np.array(np.asarray(notStopParticle), dtype=bool) + particles = particleTools.removePart(particles, mask, nStop, 'because their mass or velocity is zero', snowSlide=snowSlide) return particles @@ -1105,6 +1158,7 @@ def updateFieldsC(cfg, particles, dem, fields): """ # read input parameters cdef double rho = cfg.getfloat('rho') + cdef double rhoEnt = cfg.getfloat('rhoEnt') cdef int interpOption = cfg.getint('interpOption') header = dem['header'] cdef int nrows = header['nrows'] @@ -1117,12 +1171,17 @@ def updateFieldsC(cfg, particles, dem, fields): # read particles and fields cdef double[:] mass = particles['m'] cdef double[:] massDet = particles['dmDet'] + cdef double[:] massEnt = particles['dmEnt'] 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[:] trajectoryAngleArray = particles['trajectoryAngle'] + cdef double[:] mStoppedArray = particles['stoppedParticles']['m'] + cdef double[:] xStoppedArray = particles['stoppedParticles']['x'] + cdef double[:] yStoppedArray = particles['stoppedParticles']['y'] + cdef bint computeTA = fields['computeTA'] cdef bint computeKE = fields['computeKE'] cdef bint computeP = fields['computeP'] @@ -1135,6 +1194,8 @@ def updateFieldsC(cfg, particles, dem, fields): # initialize outputs cdef double[:, :] MassBilinear = np.zeros((nrows, ncols)) cdef double[:, :] MassDetBilinear = np.zeros((nrows, ncols)) + cdef double[:, :] MassStopBilinear = np.zeros((nrows, ncols)) + cdef double[:, :] MassEntBilinear = np.zeros((nrows, ncols)) cdef double[:, :] VBilinear = np.zeros((nrows, ncols)) cdef double[:, :] PBilinear = np.zeros((nrows, ncols)) cdef double[:, :] FTBilinear = np.zeros((nrows, ncols)) @@ -1146,10 +1207,14 @@ def updateFieldsC(cfg, particles, dem, fields): cdef double[:, :] VZBilinear = np.zeros((nrows, ncols)) cdef double[:, :] kineticEnergy = np.zeros((nrows, ncols)) cdef double[:, :] travelAngleField = np.zeros((nrows, ncols)) + cdef double[:, :] FTDetBilinear = np.zeros((nrows, ncols)) + cdef double[:, :] FTStopBilinear = np.zeros((nrows, ncols)) + cdef double[:, :] FTEntBilinear = np.zeros((nrows, ncols)) # declare intermediate step variables cdef double[:] hBB = np.zeros((nPart)) - cdef double m, dm, h, x, y, z, s, ux, uy, uz, nx, ny, nz, hbb, hLim, areaPart, trajectoryAngle - cdef int k, i + 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 cdef int indx, indy cdef int ind1[4] cdef int ind2[4] @@ -1158,7 +1223,7 @@ def updateFieldsC(cfg, particles, dem, fields): # variables for interpolation cdef int Lx0, Ly0, iCell cdef double w[4] - cdef double mwi, dmwi + cdef double mwi, dmDetWi, dmEntWi for k in range(nPart): x = xArray[k] @@ -1167,7 +1232,8 @@ def updateFieldsC(cfg, particles, dem, fields): uy = uyArray[k] uz = uzArray[k] m = mass[k] - dm = massDet[k] + dmDet = massDet[k] + dmEnt = massEnt[k] # find coordinates in normalized ref (origin (0,0) and cellsize 1) # find coordinates of the 4 nearest cornes on the raster # prepare for bilinear interpolation @@ -1184,18 +1250,35 @@ def updateFieldsC(cfg, particles, dem, fields): indx = Lx0 + ind1[i] indy = Ly0 + ind2[i] mwi = m * w[i] - dmwi = dm * w[i] + dmDetWi = dmDet * w[i] + dmEntWi = dmEnt * w[i] MassBilinear[indy, indx] = MassBilinear[indy, indx] + mwi - MassDetBilinear[indy, indx] = MassDetBilinear[indy, indx] + dmwi # PS TODO: sinnvoll??? + MassDetBilinear[indy, indx] = MassDetBilinear[indy, indx] + dmDetWi + MassEntBilinear[indy, indx] = MassEntBilinear[indy, indx] + dmEntWi MomBilinearX[indy, indx] = MomBilinearX[indy, indx] + mwi * ux MomBilinearY[indy, indx] = MomBilinearY[indy, indx] + mwi * uy MomBilinearZ[indy, indx] = MomBilinearZ[indy, indx] + mwi * uz + for l in range(len(mStoppedArray)): + xStop = xStoppedArray[l] + yStop = yStoppedArray[l] + mStop = - mStoppedArray[l] # Stopped mass is negativ for the flow + + Lx0, Ly0, iCell, w[0], w[1], w[2], w[3] = DFAtlsC.getCellAndWeights(xStop, yStop, ncols, nrows, csz, interpOption) + # for the travel angle we simply do a nearest interpolation + indx = math.round(xStop / csz) + indy = math.round(yStop / csz) + + for i in range(4): + indx = Lx0 + ind1[i] + indy = Ly0 + ind2[i] + mwi = mStop * w[i] + MassStopBilinear[indy, indx] = MassStopBilinear[indy, indx] + mwi + for i in range(ncols): for j in range(nrows): m = MassBilinear[j, i] - #PS: TODO: sinnvoll?? DMDet[j, i] = DMDet[j, i] + MassDetBilinear[j, i] if m > 0: @@ -1223,6 +1306,11 @@ def updateFieldsC(cfg, particles, dem, fields): if kineticEnergy[j, i] > PKE[j, i]: PKE[j, i] = kineticEnergy[j, i] + # thickness change due to detrainment, stopping and entrainment + 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) + fields['FM'] = np.asarray(MassBilinear) fields['FV'] = np.asarray(VBilinear) fields['Vx'] = np.asarray(VXBilinear) @@ -1232,6 +1320,10 @@ def updateFieldsC(cfg, particles, dem, fields): fields['pfv'] = np.asarray(PFV) fields['pft'] = np.asarray(PFT) fields['dmDet'] = np.asarray(DMDet) + fields['FTStop'] = np.asarray(FTStopBilinear) + fields['FTDet'] = np.asarray(FTDetBilinear) + fields['FTEnt'] = np.asarray(FTEntBilinear) + if computeP: fields['ppr'] = np.asarray(PP) fields['P'] = np.asarray(PBilinear) diff --git a/avaframe/com1DFA/checkCfg.py b/avaframe/com1DFA/checkCfg.py index 13bf2dc81..6059a7228 100644 --- a/avaframe/com1DFA/checkCfg.py +++ b/avaframe/com1DFA/checkCfg.py @@ -1,5 +1,5 @@ """ - Check if the .ini file is consistent for com1DFA computations +Check if the .ini file is consistent for com1DFA computations """ import logging @@ -10,8 +10,8 @@ def checkCfgConsistency(cfg): - """ Check the provided configuration for necessary consistency/relation between - parameters for com1DFA. + """Check the provided configuration for necessary consistency/relation between + parameters for com1DFA and for plausibility. Parameters ----------- @@ -23,24 +23,38 @@ def checkCfgConsistency(cfg): """ # Check if Ata Parameters are consistent - sphOption = float(cfg['GENERAL']['sphOption']) - viscOption = float(cfg['GENERAL']['viscOption']) + sphOption = float(cfg["GENERAL"]["sphOption"]) + viscOption = float(cfg["GENERAL"]["viscOption"]) if viscOption == 2: if sphOption != 2: # raise an error - message = ('If viscOption is set to 2 (Ata viscosity), sphOption = 2 is needed ' - '(or implement the Ata viscosity for other sphOption values)') + message = ( + "If viscOption is set to 2 (Ata viscosity), sphOption = 2 is needed " + "(or implement the Ata viscosity for other sphOption values)" + ) log.error(message) raise AssertionError(message) else: - log.debug('Ata parameters are consistent') + log.debug("Ata parameters are consistent") + + num = float(cfg["GENERAL"]["methodMeshNormal"]) + if num not in [1, 4, 6, 8]: + # raise an error + message = ( + "The methodMeshNormal set in the ini file does NOT have a plausible value. " + "Please set the parameter either to 1, 4, 6 or 8." + ) + log.error(message) + raise AssertionError(message) + else: + log.debug("MethodMeshNormal value is plausible") return True def checkCellSizeKernelRadius(cfg): - """ Check if sphKernelRadius is set to match the meshCellSize if so adapt sphKernelRadius value + """Check if sphKernelRadius is set to match the meshCellSize if so adapt sphKernelRadius value Parameters ----------- @@ -55,72 +69,105 @@ def checkCellSizeKernelRadius(cfg): """ # Check if Ata Parameters are consistent - sphKernelRadius = cfg['GENERAL']['sphKernelRadius'] - meshCellSize = cfg['GENERAL']['meshCellSize'] + sphKernelRadius = cfg["GENERAL"]["sphKernelRadius"] + meshCellSize = cfg["GENERAL"]["meshCellSize"] - if sphKernelRadius == 'meshCellSize': - cfg['GENERAL']['sphKernelRadius'] = cfg['GENERAL']['meshCellSize'] - log.info('sphKernelRadius is set to match meshCellSize of %s meters' % meshCellSize) + if sphKernelRadius == "meshCellSize": + cfg["GENERAL"]["sphKernelRadius"] = cfg["GENERAL"]["meshCellSize"] + log.info("sphKernelRadius is set to match meshCellSize of %s meters" % meshCellSize) return cfg -def checkCfgFrictionModel(cfg, inputSimFiles, relVolume=''): - """ check which friction model is chosen and if friction model parameters are of valid type - if samosATAuto - check if +def checkCfgFrictionModel(cfg, inputSimFiles, relVolume=""): + """check which friction model is chosen and if friction model parameters are of valid type + if samosATAuto - check if - relVolume < volClassSmall - set frictModel=samosATSmall - volClassSmall <= relVolume < volClassMedium - set frictModel= samosAtMedium - relVolume >= volClassMedium - set frictModel=samosAT + relVolume < volClassSmall - set frictModel=samosATSmall + volClassSmall <= relVolume < volClassMedium - set frictModel= samosAtMedium + relVolume >= volClassMedium - set frictModel=samosAT - Parameters - ------------- - cfg: dict - configuration settings - inputSimFiles: dict - info about available input files for simulation - relVolume: float - The release volume; optional - only needed if samosATAuto is set + Parameters + ------------- + cfg: dict + configuration settings + inputSimFiles: dict + info about available input files for simulation + relVolume: float + The release volume; optional - only needed if samosATAuto is set - Returns - -------- - cfg: dict - upated configuration settings + Returns + -------- + cfg: dict + upated configuration settings """ - frictParameters = ['musamosat', 'tau0samosat', 'Rs0samosat', 'kappasamosat', 'Rsamosat', - 'Bsamosat', - 'muvoellmy', 'xsivoellmy', - 'mucoulomb', - 'mu0wetsnow', 'xsiwetsnow', - 'musamosatsmall', 'tau0samosatsmall', 'Rs0samosatsmall', 'kappasamosatsmall', 'Rsamosatsmall', - 'Bsamosatsmall', - 'musamosatmedium', 'tau0samosatmedium', 'Rs0samosatmedium', 'kappasamosatmedium', 'Rsamosatmedium', - 'Bsamosatmedium', - 'mucoulombminshear', 'tau0coulombminshear', - 'muvoellmyminshear', 'xsivoellmyminshear'] + frictParameters = [ + "musamosat", + "tau0samosat", + "Rs0samosat", + "kappasamosat", + "Rsamosat", + "Bsamosat", + "muvoellmy", + "xsivoellmy", + "mucoulomb", + "mu0wetsnow", + "xsiwetsnow", + "musamosatsmall", + "tau0samosatsmall", + "Rs0samosatsmall", + "kappasamosatsmall", + "Rsamosatsmall", + "Bsamosatsmall", + "musamosatmedium", + "tau0samosatmedium", + "Rs0samosatmedium", + "kappasamosatmedium", + "Rsamosatmedium", + "Bsamosatmedium", + "mucoulombminshear", + "tau0coulombminshear", + "muvoellmyminshear", + "xsivoellmyminshear", + ] # if samosATAuto check for release volume and volume class to determine which parameter setup - if cfg['GENERAL']['frictModel'].lower() == 'samosatauto': - if relVolume < float(cfg['GENERAL']['volClassSmall']): - cfg['GENERAL']['frictModel'] = 'samosATSmall' - elif float(cfg['GENERAL']['volClassSmall']) <= relVolume < float(cfg['GENERAL']['volClassMedium']): - cfg['GENERAL']['frictModel'] = 'samosATMedium' - elif relVolume > float(cfg['GENERAL']['volClassMedium']): - cfg['GENERAL']['frictModel'] = 'samosAT' - log.info('samosATAuto - %.2f meter grid based release volume is %.2f and hence friction model: %s is chosen' % - (float(cfg['GENERAL']['meshCellSize']), relVolume, cfg['GENERAL']['frictModel'])) + if cfg["GENERAL"]["frictModel"].lower() == "samosatauto": + if relVolume < float(cfg["GENERAL"]["volClassSmall"]): + cfg["GENERAL"]["frictModel"] = "samosATSmall" + elif float(cfg["GENERAL"]["volClassSmall"]) <= relVolume < float(cfg["GENERAL"]["volClassMedium"]): + cfg["GENERAL"]["frictModel"] = "samosATMedium" + elif relVolume > float(cfg["GENERAL"]["volClassMedium"]): + cfg["GENERAL"]["frictModel"] = "samosAT" + log.info( + "samosATAuto - %.2f meter grid based release volume is %.2f and hence friction model: %s is chosen" + % ( + float(cfg["GENERAL"]["meshCellSize"]), + relVolume, + cfg["GENERAL"]["frictModel"], + ) + ) # fetch friction model - frictModel = cfg['GENERAL']['frictModel'] + frictModel = cfg["GENERAL"]["frictModel"] # remove friction parameter values that are not used - if frictModel.lower() == 'samosat': - frictParameterIn = [frictModel.lower() in frictP and 'small' not in frictP and 'medium' not in frictP for frictP in frictParameters] + if frictModel.lower() == "samosat": + frictParameterIn = [ + frictModel.lower() in frictP and "small" not in frictP and "medium" not in frictP + for frictP in frictParameters + ] else: - if frictModel.lower() == 'coulomb': - frictParameterIn = [frictModel.lower() in frictP and 'coulombminshear' not in frictP for frictP in frictParameters] - elif frictModel.lower() == 'voellmy': - frictParameterIn = [frictModel.lower() in frictP and 'voellmyminshear' not in frictP for frictP in frictParameters] + if frictModel.lower() == "coulomb": + frictParameterIn = [ + frictModel.lower() in frictP and "coulombminshear" not in frictP + for frictP in frictParameters + ] + elif frictModel.lower() == "voellmy": + frictParameterIn = [ + frictModel.lower() in frictP and "voellmyminshear" not in frictP + for frictP in frictParameters + ] else: frictParameterIn = [frictModel.lower() in frictP for frictP in frictParameters] @@ -128,32 +175,47 @@ def checkCfgFrictionModel(cfg, inputSimFiles, relVolume=''): if frictParameterIn[indexP]: noF = False try: - float(cfg['GENERAL'][frictP]) + float(cfg["GENERAL"][frictP]) except ValueError: noF = True if noF: - message = 'Friction model used %s, but %s is not of valid float type' % (frictModel, cfg['GENERAL'][frictP]) + message = "Friction model used %s, but %s is not of valid float type" % ( + frictModel, + cfg["GENERAL"][frictP], + ) log.error(message) raise ValueError(message) - elif np.isnan(float(cfg['GENERAL'][frictP])): - message = 'Friction model used %s, but %s is nan - not valid' % (frictModel, frictP) + elif np.isnan(float(cfg["GENERAL"][frictP])): + message = "Friction model used %s, but %s is nan - not valid" % ( + frictModel, + frictP, + ) log.error(message) raise ValueError(message) else: - log.info('Friction model parameter used: %s with value %s' % (frictP, cfg['GENERAL'][frictP])) + log.info( + "Friction model parameter used: %s with value %s" % (frictP, cfg["GENERAL"][frictP]) + ) # if spatialVoellmy check if mu and xi file are provided - if frictModel.lower() == 'spatialvoellmy': - if inputSimFiles['muFile'] is None: - message = 'No file for mu found' + if frictModel.lower() == "spatialvoellmy": + if inputSimFiles["muFile"] is None: + message = "No file for mu found" log.error(message) raise FileNotFoundError(message) - elif inputSimFiles['xiFile'] is None: - message = 'No file for xi found' + elif inputSimFiles["xiFile"] is None: + message = "No file for xi found" log.error(message) raise FileNotFoundError(message) else: - log.info('Mu field initialized from: %s and xi field from: %s and read from: %s, %s' % - (inputSimFiles['muFile'].name, inputSimFiles['xiFile'].name, cfg['INPUT']['muFile'], cfg['INPUT']['xiFile'])) + log.info( + "Mu field initialized from: %s and xi field from: %s and read from: %s, %s" + % ( + inputSimFiles["muFile"].name, + inputSimFiles["xiFile"].name, + cfg["INPUT"]["muFile"], + cfg["INPUT"]["xiFile"], + ) + ) return cfg diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index 8750f7822..f7f6aa3fd 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -151,6 +151,7 @@ def com1DFAMain(cfgMain, cfgInfo=""): for key in simDict: log.info("Simulation: %s" % key) exportFlag = simDict[key]["cfgSim"]["EXPORTS"].getboolean("exportData") + adaptDemPlot = simDict[key]["cfgSim"]["GENERAL"].getboolean("adaptDemPlot") # initialize reportDict list reportDictList = list() @@ -189,6 +190,7 @@ def com1DFAMain(cfgMain, cfgInfo=""): log.info("Overall (parallel) com1DFA computation took: %s s " % timeNeeded) log.info("--- ENDING (potential) PARALLEL PART ----") + dem = com1DFATools.chooseDemPlot(dem, adaptedDemBackground=adaptDemPlot) # postprocessing: writing report, creating plots dem, plotDict, reportDictList, simDFNew = com1DFAPostprocess( simDF, @@ -219,7 +221,7 @@ def com1DFACoreTask(simDict, inputSimFiles, avalancheDir, outDir, cuSim): # load configuration object for current sim cfg = simDict[cuSim]["cfgSim"] - # check configuraton for consistency + # check configuraton for consistency and plausibility checkCfg.checkCfgConsistency(cfg) # fetch simHash for current sim @@ -317,6 +319,7 @@ def com1DFAPostprocess(simDF, tCPUDF, simDFExisting, cfgMain, dem, reportDictLis # create plots and report reportDir = pathlib.Path(avalancheDir, "Outputs", modName, "reports") fU.makeADir(reportDir) + # Generate plots for all peakFiles if exportData: plotDict = oP.plotAllPeakFields(avalancheDir, cfgMain["FLAGS"], modName, demData=dem) @@ -895,6 +898,7 @@ def initializeMesh(cfg, demOri, num): # set origin to 0, 0 for computations, store original origin dem = setDEMoriginToZero(demOri) dem["originalHeader"] = demOri["header"].copy() + dem["originalRasterData"] = dem["rasterData"].copy() # read dem header headerDEM = dem["header"] @@ -1054,11 +1058,11 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName): releaseLine["header"] = dem["originalHeader"] inputSimLines["releaseLine"]["header"] = dem["originalHeader"] # export release area raster to file - if cfg['EXPORTS'].getboolean('exportRasters'): - outDir = pathlib.Path(cfgGen['avalancheDir'], 'Outputs', 'internalRasters') + if cfg["EXPORTS"].getboolean("exportRasters"): + outDir = pathlib.Path(cfgGen["avalancheDir"], "Outputs", "internalRasters") fU.makeADir(outDir) - IOf.writeResultToRaster(dem["originalHeader"], relRaster, (outDir / 'releaseRaster'), flip=True) - log.info('Release area raster derived from shp file saved to %s' % str(outDir / 'releaseRaster')) + IOf.writeResultToRaster(dem["originalHeader"], relRaster, (outDir / "releaseRaster"), flip=True) + log.info("Release area raster derived from shp file saved to %s" % str(outDir / "releaseRaster")) particles = initializeParticles( cfgGen, releaseLine, @@ -1358,6 +1362,7 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel particles["yllcenter"] = dem["originalHeader"]["yllcenter"] particles["nExitedParticles"] = 0.0 particles["dmDet"] = np.zeros(np.shape(hPartArray)) + particles["dmEnt"] = np.zeros(np.shape(hPartArray)) # remove particles that might lay outside of the release polygon if not cfg.getboolean("iniStep") and not cfg.getboolean("initialiseParticlesFromFile"): @@ -1403,6 +1408,14 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel if debugPlot: debPlot.plotPartIni(particles, dem) + # space for deposited particles TODO: extra dict for them? + particles["stoppedParticles"] = { + "x": np.empty(0), + "y": np.empty(0), + "m": np.empty(0), + "ID": np.empty(0), + "velocityMag": np.empty(0), + } return particles @@ -1476,6 +1489,12 @@ def initializeFields(cfg, dem, particles, releaseLine): 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)) # 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): @@ -1664,10 +1683,17 @@ def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPoin log.info("Entrainment area features: %s" % (entLine["Name"])) entLine = geoTrans.prepareArea(entLine, dem, thresholdPointInPoly, thList=entLine["thickness"]) entrMassRaster = entLine["rasterData"] - if cfg['EXPORTS'].getboolean('exportRasters'): - outDir = pathlib.Path(cfg['GENERAL']['avalancheDir'], 'Outputs', 'internalRasters') - IOf.writeResultToRaster(dem["originalHeader"], entrMassRaster, (outDir / 'entrainmentRaster'), flip=True) - log.info('Release area raster derived from shp file saved to %s' % str(outDir / 'entrainmentRaster')) + if cfg["EXPORTS"].getboolean("exportRasters"): + outDir = pathlib.Path(cfg["GENERAL"]["avalancheDir"], "Outputs", "internalRasters") + IOf.writeResultToRaster( + dem["originalHeader"], + entrMassRaster, + (outDir / "entrainmentRaster"), + flip=True, + ) + log.info( + "Release area raster derived from shp file saved to %s" % str(outDir / "entrainmentRaster") + ) # ToDo: not used in samos but implemented # tempRaster = cfg['GENERAL'].getfloat('entTempRef') + (dem['rasterData'] - cfg['GENERAL'].getfloat('entMinZ')) # * cfg['GENERAL'].getfloat('entTempGrad') @@ -1675,7 +1701,7 @@ def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPoin # tempRaster*cfg['GENERAL'].getfloat('cpWtr')/cfg['GENERAL'].getfloat('hFusion')) entrEnthRaster = np.where( entrMassRaster > 0, - cfg['GENERAL'].getfloat("entTempRef") * cfg['GENERAL'].getfloat("cpIce"), + cfg["GENERAL"].getfloat("entTempRef") * cfg["GENERAL"].getfloat("cpIce"), 0, ) reportAreaInfo["entrainment"] = "Yes" @@ -1684,7 +1710,7 @@ def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPoin entrEnthRaster = np.zeros((nrows, ncols)) reportAreaInfo["entrainment"] = "No" - entrMassRaster = entrMassRaster * cfg['GENERAL'].getfloat("rhoEnt") + entrMassRaster = entrMassRaster * cfg["GENERAL"].getfloat("rhoEnt") return entrMassRaster, entrEnthRaster, reportAreaInfo @@ -1847,6 +1873,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si timeM = [] massEntrained = [] massDetrained = [] + massStopped = [] massTotal = [] pfvTimeMax = [] @@ -1914,7 +1941,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si startTime = time.time() log.debug("Computing time step t = %f s, dt = %f s" % (t, dt)) # Perform computations - particles, fields, zPartArray0, tCPU = computeEulerTimeStep( + particles, fields, zPartArray0, tCPU, dem = computeEulerTimeStep( cfgGen, particles, fields, zPartArray0, dem, tCPU, frictType, resistanceType ) # set max values of fields to dataframe @@ -1930,6 +1957,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si # write mass balance info massEntrained.append(particles["massEntrained"]) massDetrained.append(particles["massDetrained"]) + massStopped.append(particles["massStopped"]) massTotal.append(particles["mTot"]) timeM.append(t) pfvTimeMax.append(np.nanmax(fields["FV"])) @@ -2004,7 +2032,6 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si nIter0 = nIter0 + 1 tCPUtimeLoop = time.time() - startTime tCPU["timeLoop"] = tCPU["timeLoop"] + tCPUtimeLoop - tCPU["nIter"] = nIter log.info("Ending computation at time t = %f s", t - dt) log.debug("Saving results for time step t = %f s", t - dt) @@ -2041,6 +2068,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si infoDict = { "massEntrained": massEntrained, "massDetrained": massDetrained, + "massStopped": massStopped, "timeStep": timeM, "massTotal": massTotal, "tCPU": tCPU, @@ -2278,6 +2306,7 @@ def writeMBFile(infoDict, avaDir, logName): t = infoDict["timeStep"] massEntrained = infoDict["massEntrained"] massDetrained = infoDict["massDetrained"] + massStopped = infoDict["massStopped"] massTotal = infoDict["massTotal"] massDetrainedTotal = np.zeros(len(massDetrained)) for m in range(1, len(massDetrained)): @@ -2290,16 +2319,17 @@ def writeMBFile(infoDict, avaDir, logName): massDir = pathlib.Path(avaDir, "Outputs", "com1DFA") fU.makeADir(massDir) with open(massDir / ("mass_%s.txt" % logName), "w") as mFile: - mFile.write("time, current, entrained, detrained\n") + mFile.write("time, current, entrained, detrained, detrainedTotal, stopped\n") for m in range(len(t)): mFile.write( - "%.02f, %.06f, %.06f, %.06f, %.06f\n" + "%.02f, %.06f, %.06f, %.06f, %.06f, %.06f\n" % ( t[m], massTotal[m], massEntrained[m], massDetrained[m], massDetrainedTotal[m], + massStopped[m], ) ) @@ -2334,6 +2364,8 @@ def computeEulerTimeStep(cfg, particles, fields, zPartArray0, dem, tCPU, frictTy fields dictionary at t + dt tCPU : dict computation time dictionary + dem: dict + dictionary with dem information including the adapted DEM """ # update cRes and detK rasters according to thresholds of FV and FT @@ -2418,10 +2450,21 @@ def computeEulerTimeStep(cfg, particles, fields, zPartArray0, dem, tCPU, frictTy if fields["computeTA"]: particles = DFAfunC.computeTrajectoryAngleC(particles, zPartArray0) particles, fields = DFAfunC.updateFieldsC(cfg, particles, dem, fields) + + # adapt DEM considering erosion and deposition + # only adapt DEM when in one grid cell the changing height > threshold + thresholdAdaptSfc = cfg.getfloat("thresholdAdaptSfc") + adaptStop = cfg.getboolean("adaptSfcStopped") and np.any(abs(fields["FTStop"]) > thresholdAdaptSfc) + adaptDet = cfg.getboolean("adaptSfcDetrainment") and np.any(abs(fields["FTDet"]) > thresholdAdaptSfc) + adaptEnt = cfg.getboolean("adaptSfcEntrainment") and np.any(abs(fields["FTEnt"]) > thresholdAdaptSfc) + if particles["t"] > 0: + if adaptStop or adaptDet or adaptEnt: + dem, fields = adaptDEM(dem, fields, cfg) + tCPUField = time.time() - startTime tCPU["timeField"] = tCPU["timeField"] + tCPUField - return particles, fields, zPartArray0, tCPU + return particles, fields, zPartArray0, tCPU, dem def releaseSecRelArea(cfg, particles, fields, dem, zPartArray0): @@ -2707,6 +2750,8 @@ def exportFields( if resType == "FTDet": dmDet = fields["dmDet"] resField = dmDet / (cfg["GENERAL"].getfloat("rho") * dem["areaRaster"]) + elif TSave == "final" and resType == "sfcChange": + resField = fields["sfcChangeTotal"] else: resField = fields[resType] if resType == "ppr": @@ -3260,3 +3305,61 @@ def saveContToPickle(contourDictXY, outDir, cuSimName): fi = open(outDir / ("contourDictXY_%s.pickle" % (cuSimName)), "wb") pickle.dump(contourDictXY, fi) fi.close() + + +def adaptDEM(dem, fields, cfg): + """adapt topography in respect to erosion and deposition + + Parameters + dem: dict + dictionary with info on DEM data + fields : dict + fields dictionary + cfg: dict + configuration settings + + Returns + --------- + dem: dict + dictionary with info on DEM data containing adapted topography + fields : dict + fields dictionary containing adapted DEM + """ + + ZDEM = dem["rasterData"].copy() + FTDet = fields["FTDet"] + FTStop = fields["FTStop"] + FTEnt = fields["FTEnt"] + sfcChangeTotal = fields["sfcChangeTotal"] + sfcChange = np.zeros_like(FTDet) + ZDEMadapt = ZDEM + + _, _, NzNormed = DFAtls.normalize(dem["Nx"].copy(), dem["Ny"].copy(), dem["Nz"].copy()) + + if cfg.getboolean("adaptSfcStopped"): + # compute thickness to depth + depthStop = FTStop / NzNormed + ZDEMadapt += depthStop + sfcChange += depthStop + if cfg.getboolean("adaptSfcDetrainment"): + # compute thickness to depth + depthDet = FTDet / NzNormed + ZDEMadapt += depthDet + sfcChange += depthDet + if cfg.getboolean("adaptSfcEntrainment"): + # compute thickness to depth + depthEnt = FTEnt / NzNormed + ZDEMadapt += depthEnt + sfcChange += depthEnt + + dem["rasterData"] = ZDEMadapt + fields["demAdapted"] = ZDEMadapt + + num = cfg.getfloat("methodMeshNormal") + dem = geoTrans.getNormalMesh(dem, num=num) + dem = DFAtls.getAreaMesh(dem, num) + + fields["sfcChange"] = sfcChange + fields["sfcChangeTotal"] = sfcChangeTotal + sfcChange + + return dem, fields diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index 366035dcc..d1dceab93 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/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, particles) - separated by | +# 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 # 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) @@ -385,6 +385,20 @@ entEroEnergy = 5000 entShearResistance = 0 entDefResistance = 0 +#++++++++++++ Deposition, Erosion and adaptive surface +# delete stopped particles (velocity = 0, mass = 0 or stopping criterion) but save in a additional dictionary +# activate with 1 +delStoppedParticles = 0 +# adapt the topography every time step +# activate with 1 +adaptSfcStopped = 0 +adaptSfcDetrainment = 0 +adaptSfcEntrainment = 0 +# 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 +adaptDemPlot = False + #+++++++++++++++++ Dam Parameters # the dam foot print is given in Inputs/DAM as a shape file line # the dam is located on the left side of the dam line diff --git a/avaframe/com1DFA/com1DFATools.py b/avaframe/com1DFA/com1DFATools.py index 125deac61..9eef87aac 100644 --- a/avaframe/com1DFA/com1DFATools.py +++ b/avaframe/com1DFA/com1DFATools.py @@ -459,3 +459,29 @@ def updateResCoeffFields(fields, cfg, t, dem): fields["detRaster"] = detRaster return fields + + +def chooseDemPlot(dem, adaptedDemBackground=False): + """ + choose the DEM (rasterData) that is used for report plots + + Parameters + ----------- + adaptedDemBackground: bool + if True the adapted DEM is used as background in the plots + + Returns + ------- + demPlot: dict + contains the rasterdata that should be used in report plots as background + """ + + demPlot = dem.copy() + # choose if the input DEM or the adapted DEM is used as background + if adaptedDemBackground: + demPlot["rasterData"] = dem["rasterData"] + log.info("The resulting adapted DEM is used in the plots.") + else: + demPlot["rasterData"] = dem["originalRasterData"] + log.info("The original DEM is used in the plots.") + return demPlot diff --git a/avaframe/com1DFA/deriveParameterSet.py b/avaframe/com1DFA/deriveParameterSet.py index 84a8bd27d..372810384 100644 --- a/avaframe/com1DFA/deriveParameterSet.py +++ b/avaframe/com1DFA/deriveParameterSet.py @@ -151,6 +151,8 @@ def checkResType(fullCfg, section, key, value): "particles", "dmDet", "FTDet", + "sfcChange", + "demAdapted", ] message = "The parameter % s is not a valid resType. It will not be saved" newResType = [] diff --git a/avaframe/com1DFA/particleTools.py b/avaframe/com1DFA/particleTools.py index 4c61dbad7..f4d042631 100644 --- a/avaframe/com1DFA/particleTools.py +++ b/avaframe/com1DFA/particleTools.py @@ -1,5 +1,5 @@ """ - Tools for handling particles, splitting, merging and tracking. +Tools for handling particles, splitting, merging and tracking. """ # Load modules @@ -27,40 +27,40 @@ def initialiseParticlesFromFile(cfg, avaDir, releaseScenario): # TODO: this is for development purposes, change or remove in the future # If initialisation from file - if cfg['particleFile']: - inDirPart = pathlib.Path(cfg['particleFile']) + if cfg["particleFile"]: + inDirPart = pathlib.Path(cfg["particleFile"]) else: - inDirPart = pathlib.Path(avaDir, 'Outputs', 'com1DFAOrig') + inDirPart = pathlib.Path(avaDir, "Outputs", "com1DFAOrig") - searchDir = inDirPart / 'particles' - inDirPart = list(searchDir.glob( - ('*' + releaseScenario + '_' + '*' + cfg['simTypeActual'] + '*'))) + searchDir = inDirPart / "particles" + inDirPart = list(searchDir.glob(("*" + releaseScenario + "_" + "*" + cfg["simTypeActual"] + "*"))) if inDirPart == []: - messagePart = ('Initialise particles from file - no particles file found for releaseScenario: %s and simType: %s' % - (releaseScenario, cfg['simTypeActual'])) + messagePart = ( + "Initialise particles from file - no particles file found for releaseScenario: %s and simType: %s" + % (releaseScenario, cfg["simTypeActual"]) + ) log.error(messagePart) raise FileNotFoundError(messagePart) elif len(inDirPart) > 1: - log.warning('More than one file found for Initialise particle from file: took %s' % inDirPart[0]) + log.warning("More than one file found for Initialise particle from file: took %s" % inDirPart[0]) inDirPart = inDirPart[0] else: inDirPart = inDirPart[0] - log.info('Initial particle distribution read from file!! %s' % - (inDirPart)) + log.info("Initial particle distribution read from file!! %s" % (inDirPart)) Particles, timeStepInfo = readPartFromPickle(inDirPart) particles = Particles[0] - xPartArray = particles['x'] + xPartArray = particles["x"] hPartArray = np.ones(len(xPartArray)) - particles['nPart'] = len(xPartArray) - particles['trajectoryLengthXY'] = np.zeros(np.shape(xPartArray)) - particles['trajectoryLengthXYZ'] = np.zeros(np.shape(xPartArray)) - particles['idFixed'] = np.zeros(np.shape(xPartArray)) + particles["nPart"] = len(xPartArray) + particles["trajectoryLengthXY"] = np.zeros(np.shape(xPartArray)) + particles["trajectoryLengthXYZ"] = np.zeros(np.shape(xPartArray)) + particles["idFixed"] = np.zeros(np.shape(xPartArray)) return particles, hPartArray def placeParticles(hCell, aCell, indx, indy, csz, massPerPart, nPPK, rng, cfg, ratioArea): - """ Create particles in given cell + """Create particles in given cell Compute number of particles to create in a given cell. Place particles in cell according to the chosen pattern (random semirandom @@ -102,14 +102,14 @@ def placeParticles(hCell, aCell, indx, indy, csz, massPerPart, nPPK, rng, cfg, r particle area """ - rho = cfg.getfloat('rho') - thresholdMassSplit = cfg.getfloat('thresholdMassSplit') - initPartDistType = cfg['initPartDistType'].lower() - massPerParticleDeterminationMethod = cfg['massPerParticleDeterminationMethod'] + rho = cfg.getfloat("rho") + thresholdMassSplit = cfg.getfloat("thresholdMassSplit") + initPartDistType = cfg["initPartDistType"].lower() + massPerParticleDeterminationMethod = cfg["massPerParticleDeterminationMethod"] volCell = aCell * hCell massCell = volCell * rho - if initPartDistType == 'random': - if massPerParticleDeterminationMethod == 'MPPKR': + if initPartDistType == "random": + if massPerParticleDeterminationMethod == "MPPKR": # impose a number of particles within a kernel radius so impose number of particles in a cell nFloat = nPPK * aCell / (math.pi * csz**2) else: @@ -125,16 +125,16 @@ def placeParticles(hCell, aCell, indx, indy, csz, massPerPart, nPPK, rng, cfg, r # make sure we do not violate the (massCell / n) < thresholdMassSplit x massPerPart rule if (massCell / n) / massPerPart >= thresholdMassSplit: n = n + 1 - elif initPartDistType == 'semirandom' or initPartDistType == 'uniform': - n1 = (np.ceil(np.sqrt(massCell / massPerPart))).astype('int') - n = n1*n1 - d = csz/n1 - pos = np.linspace(0., csz-d, n1) + d/2. + elif initPartDistType == "semirandom" or initPartDistType == "uniform": + n1 = (np.ceil(np.sqrt(massCell / massPerPart))).astype("int") + n = n1 * n1 + d = csz / n1 + pos = np.linspace(0.0, csz - d, n1) + d / 2.0 x, y = np.meshgrid(pos, pos) x = x.flatten() y = y.flatten() - elif initPartDistType == 'triangular': - if massPerParticleDeterminationMethod == 'MPPKR': + elif initPartDistType == "triangular": + if massPerParticleDeterminationMethod == "MPPKR": # impose a number of particles within a kernel radius so impose number of particles in a cell # (regardless of the slope) nPPC = nPPK / math.pi @@ -143,18 +143,20 @@ def placeParticles(hCell, aCell, indx, indy, csz, massPerPart, nPPK, rng, cfg, r # ToDo: this only works if the release thickness is constant in a release area!!! nPPC = hCell * (csz**2 / ratioArea) * rho / massPerPart n = int(np.floor(nPPC)) - indx = indx - 1/2 - indy = indy - 1/2 + indx = indx - 1 / 2 + indy = indy - 1 / 2 # compute triangles properties Aparticle = csz**2 / n - sTri = math.sqrt(Aparticle/(math.sqrt(3)/2)) - hTri = sTri * math.sqrt(3)/2 - jMin = int(indy * csz/hTri) + sTri = math.sqrt(Aparticle / (math.sqrt(3) / 2)) + hTri = sTri * math.sqrt(3) / 2 + jMin = int(indy * csz / hTri) if jMin * hTri < indy * csz: jMin = jMin + 1 else: - log.warning('Chosen value for initial particle distribution type not available: %s uniform is used instead' % - initPartDistType) + log.warning( + "Chosen value for initial particle distribution type not available: %s uniform is used instead" + % initPartDistType + ) mPart = massCell / n aPart = aCell / n @@ -162,41 +164,41 @@ def placeParticles(hCell, aCell, indx, indy, csz, massPerPart, nPPK, rng, cfg, r # TODO make this an independent function ####################### # start ############### - if initPartDistType == 'random': + if initPartDistType == "random": # place particles randomly in the cell xPart = csz * (rng.random(n) - 0.5 + indx) yPart = csz * (rng.random(n) - 0.5 + indy) - elif initPartDistType == 'semirandom' or initPartDistType == 'uniform': + elif initPartDistType == "semirandom" or initPartDistType == "uniform": # place particles equaly distributed - xPart = csz * (- 0.5 + indx) + x - yPart = csz * (- 0.5 + indy) + y - if initPartDistType == 'semirandom': + xPart = csz * (-0.5 + indx) + x + yPart = csz * (-0.5 + indy) + y + if initPartDistType == "semirandom": # place particles equaly distributed with a small variation xPart = xPart + (rng.random(n) - 0.5) * d yPart = yPart + (rng.random(n) - 0.5) * d - elif initPartDistType == 'triangular': + elif initPartDistType == "triangular": xPart = np.empty(0) yPart = np.empty(0) n = 0 j = jMin - iTemp = int(indx * csz/sTri) + iTemp = int(indx * csz / sTri) if iTemp * sTri < indx * csz: iTemp = iTemp + 1 - while j * hTri < (indy+1) * csz: - i = iTemp - 1/2 * j%2 - while i * sTri < (indx+1) * csz: + while j * hTri < (indy + 1) * csz: + i = iTemp - 1 / 2 * j % 2 + while i * sTri < (indx + 1) * csz: if i * sTri >= indx * csz and j * hTri >= indy * csz: - xPart = np.append(xPart, i*sTri) - yPart = np.append(yPart, j*hTri) - n = n+1 - i = i+1 - j = j+1 + xPart = np.append(xPart, i * sTri) + yPart = np.append(yPart, j * hTri) + n = n + 1 + i = i + 1 + j = j + 1 return xPart, yPart, mPart, n, aPart -def removePart(particles, mask, nRemove, reasonString='', snowSlide=0): - """ remove given particles +def removePart(particles, mask, nRemove, reasonString="", snowSlide=0): + """remove given particles Parameters ---------- @@ -214,31 +216,33 @@ def removePart(particles, mask, nRemove, reasonString='', snowSlide=0): particles : dict particles dictionary """ - if reasonString != '': - log.debug('removed %s particles %s' % (nRemove, reasonString)) - if reasonString == 'because they exited the domain': - particles['nExitedParticles'] = particles['nExitedParticles'] + nRemove - nPart = particles['nPart'] + if reasonString != "": + log.debug("removed %s particles %s" % (nRemove, reasonString)) + if reasonString == "because they exited the domain": + particles["nExitedParticles"] = particles["nExitedParticles"] + nRemove + nPart = particles["nPart"] if snowSlide == 1: # if snowSlide is activated, we need to remove the particles as well as the bonds accordingly # we do this first befor nPart changes nBondRemove = DFAfunC.countRemovedBonds(particles, mask, nRemove) particles = DFAfunC.removedBonds(particles, mask, nRemove, nBondRemove) for key in particles: - if key == 'nPart': - particles['nPart'] = particles['nPart'] - nRemove + if key == "nPart": + particles["nPart"] = particles["nPart"] - nRemove + elif key == "stoppedParticles": + continue # for all keys in particles that are arrays of size nPart do: elif type(particles[key]).__module__ == np.__name__: if np.size(particles[key]) == nPart: particles[key] = np.extract(mask, particles[key]) - particles['mTot'] = np.sum(particles['m']) + particles["mTot"] = np.sum(particles["m"]) return particles def addParticles(particles, nAdd, ind, mNew, xNew, yNew, zNew): - """ add particles + """add particles Parameters ---------- @@ -263,38 +267,38 @@ def addParticles(particles, nAdd, ind, mNew, xNew, yNew, zNew): particles dictionary with modified particle and new ones """ # get old values - nPart = particles['nPart'] - nID = particles['nID'] + nPart = particles["nPart"] + nID = particles["nID"] # update total number of particles and number of IDs used so far - particles['nPart'] = particles['nPart'] + nAdd - particles['nID'] = nID + nAdd + particles["nPart"] = particles["nPart"] + nAdd + particles["nID"] = nID + nAdd # log.info('Spliting particle %s in %s' % (ind, nAdd)) for key in particles: # update splitted particle mass # first update the old particle - particles['m'][ind] = mNew - particles['x'][ind] = xNew[0] - particles['y'][ind] = yNew[0] - particles['z'][ind] = zNew[0] + particles["m"][ind] = mNew + particles["x"][ind] = xNew[0] + particles["y"][ind] = yNew[0] + particles["z"][ind] = zNew[0] # add new particles at the end of the arrays # for all keys in particles that are arrays of size nPart do: if type(particles[key]).__module__ == np.__name__: # create unique ID for the new particles - if key == 'ID': - particles['ID'] = np.append(particles['ID'], np.arange(nID, nID + nAdd, 1)) - elif key == 'x': + if key == "ID": + particles["ID"] = np.append(particles["ID"], np.arange(nID, nID + nAdd, 1)) + elif key == "x": particles[key] = np.append(particles[key], xNew[1:]) - elif key == 'y': + elif key == "y": particles[key] = np.append(particles[key], yNew[1:]) - elif key == 'z': + elif key == "z": particles[key] = np.append(particles[key], zNew[1:]) - elif key == 'bondStart': + elif key == "bondStart": # no bonds for added particles: - nBondsParts = np.size(particles['bondPart']) - particles[key] = np.append(particles[key], nBondsParts*np.ones((nAdd))) + nBondsParts = np.size(particles["bondPart"]) + particles[key] = np.append(particles[key], nBondsParts * np.ones((nAdd))) # set the parent properties to new particles due to splitting elif np.size(particles[key]) == nPart: - particles[key] = np.append(particles[key], particles[key][ind]*np.ones((nAdd))) + particles[key] = np.append(particles[key], particles[key][ind] * np.ones((nAdd))) # ToDo: maybe also update the h smartly return particles @@ -319,13 +323,13 @@ def splitPartMass(particles, cfg): particles dictionary """ - rho = cfg.getfloat('rho') - thresholdMassSplit = cfg.getfloat('thresholdMassSplit') - distSplitPart = cfg.getfloat('distSplitPart') - massPerPart = particles['massPerPart'] - mPart = particles['m'] + rho = cfg.getfloat("rho") + thresholdMassSplit = cfg.getfloat("thresholdMassSplit") + distSplitPart = cfg.getfloat("distSplitPart") + massPerPart = particles["massPerPart"] + mPart = particles["m"] # decide which particles to split - nSplit = np.ceil(mPart/(massPerPart*thresholdMassSplit)) + nSplit = np.ceil(mPart / (massPerPart * thresholdMassSplit)) Ind = np.where(nSplit > 1)[0] # loop on particles to split for ind in Ind: @@ -336,7 +340,7 @@ def splitPartMass(particles, cfg): # add new particles particles = addParticles(particles, nAdd, ind, mNew, xNew, yNew, zNew) - particles['mTot'] = np.sum(particles['m']) + particles["mTot"] = np.sum(particles["m"]) return particles @@ -362,42 +366,42 @@ def splitPartArea(particles, cfg, dem): """ # get cfg info - rho = cfg.getfloat('rho') - sphKernelRadius = cfg.getfloat('sphKernelRadius') - cMinNPPK = cfg.getfloat('cMinNPPK') - cMinMass = cfg.getfloat('cMinMass') - nSplit = cfg.getint('nSplit') + rho = cfg.getfloat("rho") + sphKernelRadius = cfg.getfloat("sphKernelRadius") + cMinNPPK = cfg.getfloat("cMinNPPK") + cMinMass = cfg.getfloat("cMinMass") + nSplit = cfg.getint("nSplit") # get dem info - csz = dem['header']['cellsize'] - Nx = dem['Nx'] - Ny = dem['Ny'] - Nz = dem['Nz'] + csz = dem["header"]["cellsize"] + Nx = dem["Nx"] + Ny = dem["Ny"] + Nz = dem["Nz"] # get the threshold area over which we split the particle - massPerPart = particles['massPerPart'] - nPPK = particles['nPPK'] + massPerPart = particles["massPerPart"] + nPPK = particles["nPPK"] aMax = math.pi * sphKernelRadius**2 / (cMinNPPK * nPPK) mMin = massPerPart * cMinMass # get particle area - mPart = particles['m'] - hPart = particles['h'] - aPart = mPart/(rho*hPart) + mPart = particles["m"] + hPart = particles["h"] + aPart = mPart / (rho * hPart) # find particles to split - tooBig = np.where((aPart > aMax) & (mPart/nSplit > mMin))[0] + tooBig = np.where((aPart > aMax) & (mPart / nSplit > mMin))[0] # count new particles nNewPart = 0 # loop on particles to split for ind in tooBig: # compute new mass and particles to add mNew = mPart[ind] / nSplit - nAdd = nSplit-1 + nAdd = nSplit - 1 nNewPart = nNewPart + nAdd # get position of new particles xNew, yNew, zNew = getSplitPartPosition(cfg, particles, aPart, Nx, Ny, Nz, csz, nSplit, ind) # add new particles particles = addParticles(particles, nAdd, ind, mNew, xNew, yNew, zNew) - log.debug('Added %s because of splitting' % (nNewPart)) + log.debug("Added %s because of splitting" % (nNewPart)) - particles['mTot'] = np.sum(particles['m']) + particles["mTot"] = np.sum(particles["m"]) return particles @@ -424,28 +428,30 @@ def mergePartArea(particles, cfg, dem): """ # get cfg info - rho = cfg.getfloat('rho') - sphKernelRadius = cfg.getfloat('sphKernelRadius') - cMaxNPPK = cfg.getfloat('cMaxNPPK') + rho = cfg.getfloat("rho") + sphKernelRadius = cfg.getfloat("sphKernelRadius") + cMaxNPPK = cfg.getfloat("cMaxNPPK") # get the threshold area under which we merge the particle - nPPK = particles['nPPK'] + nPPK = particles["nPPK"] aMin = math.pi * sphKernelRadius**2 / (cMaxNPPK * nPPK) # get particle area - mPart = particles['m'] - hPart = particles['h'] - xPart = particles['x'] - yPart = particles['y'] - zPart = particles['z'] - aPart = mPart/(rho*hPart) + mPart = particles["m"] + hPart = particles["h"] + xPart = particles["x"] + yPart = particles["y"] + zPart = particles["z"] + aPart = mPart / (rho * hPart) # find particles to merge tooSmall = np.where(aPart < aMin)[0] - keepParticle = np.ones((particles['nPart']), dtype=bool) + keepParticle = np.ones((particles["nPart"]), dtype=bool) nRemoved = 0 # loop on particles to merge for ind in tooSmall: if keepParticle[ind]: # find nearest particle - rMerge, neighbourInd = getClosestNeighbour(xPart, yPart, zPart, ind, sphKernelRadius, keepParticle) + rMerge, neighbourInd = getClosestNeighbour( + xPart, yPart, zPart, ind, sphKernelRadius, keepParticle + ) # only merge a particle if it is closer thant the kernel radius if rMerge < sphKernelRadius: # remove neighbourInd from tooSmall if possible @@ -454,18 +460,19 @@ def mergePartArea(particles, cfg, dem): # compute new mass and particles to add mNew = mPart[ind] + mPart[neighbourInd] # compute mass averaged values - for key in ['x', 'y', 'z', 'ux', 'uy', 'uz']: - particles[key][ind] = (mPart[ind]*particles[key][ind] + - mPart[neighbourInd]*particles[key][neighbourInd]) / mNew - particles['m'][ind] = mNew + for key in ["x", "y", "z", "ux", "uy", "uz"]: + particles[key][ind] = ( + mPart[ind] * particles[key][ind] + mPart[neighbourInd] * particles[key][neighbourInd] + ) / mNew + particles["m"][ind] = mNew # ToDo: mabe also update h - particles = removePart(particles, keepParticle, nRemoved, reasonString='') # 'because of colocation') + particles = removePart(particles, keepParticle, nRemoved, reasonString="") # 'because of colocation') return particles def getClosestNeighbour(xPartArray, yPartArray, zPartArray, ind, sphKernelRadius, keepParticle): - """ find closest neighbour + """find closest neighbour Parameters ---------- @@ -489,11 +496,15 @@ def getClosestNeighbour(xPartArray, yPartArray, zPartArray, ind, sphKernelRadius neighbourInd : int index of closest neighbour """ - r = DFAtls.norm(xPartArray-xPartArray[ind], yPartArray-yPartArray[ind], zPartArray-zPartArray[ind]) + r = DFAtls.norm( + xPartArray - xPartArray[ind], + yPartArray - yPartArray[ind], + zPartArray - zPartArray[ind], + ) # make sure you don't find the particle itself - r[ind] = 2*sphKernelRadius + r[ind] = 2 * sphKernelRadius # make sure you don't find a particle that has already been merged - r[~keepParticle] = 2*sphKernelRadius + r[~keepParticle] = 2 * sphKernelRadius # find nearest particle neighbourInd = np.argmin(r) rMerge = r[neighbourInd] @@ -517,18 +528,18 @@ def mergeParticleDict(particles1, particles2): """ particles = {} - nPart1 = particles1['nPart'] + nPart1 = particles1["nPart"] # loop on the keys from particles1 dicionary for key in particles1: # deal with specific cases # nPart: just sum them up - if key == 'nPart': - particles['nPart'] = particles1['nPart'] + particles2['nPart'] + if key == "nPart": + particles["nPart"] = particles1["nPart"] + particles2["nPart"] # massPerPart, should stay unchanged. If ever they are different take # the minimum # ToDo: are we sure we want the minimum? - elif key == 'massPerPart': - particles['massPerPart'] = min(particles1['massPerPart'], particles2['massPerPart']) + elif key == "massPerPart": + particles["massPerPart"] = min(particles1["massPerPart"], particles2["massPerPart"]) # now if the value is a numpy array and this key is also in particles2 elif (key in particles2) and (type(particles1[key]).__module__ == np.__name__): # deal with the specific cases: @@ -537,8 +548,8 @@ def mergeParticleDict(particles1, particles2): # here when we merge the 2 arrays we make sure to shift the value # of particles2 so that the ID stays a unique identifier and # that the parentID is consistent with this shift. - if (key == 'ID') or (key == 'parentID'): - particles[key] = np.append(particles1[key], (particles2[key] + particles1['nID'])) + if (key == "ID") or (key == "parentID"): + particles[key] = np.append(particles1[key], (particles2[key] + particles1["nID"])) # general case where the key value is an array with as many elements # as particles elif np.size(particles1[key]) == nPart1: @@ -591,19 +602,21 @@ def getSplitPartPosition(cfg, particles, aPart, Nx, Ny, Nz, csz, nSplit, ind): zNew : numpy array z components of the splitted particles """ - rng = np.random.default_rng(int(cfg['seed'])) - x = particles['x'] - y = particles['y'] - z = particles['z'] - ux = particles['ux'] - uy = particles['uy'] - uz = particles['uz'] + rng = np.random.default_rng(int(cfg["seed"])) + x = particles["x"] + y = particles["y"] + z = particles["z"] + ux = particles["ux"] + uy = particles["uy"] + uz = particles["uz"] rNew = np.sqrt(aPart[ind] / (math.pi * nSplit)) - alpha = 2*math.pi*(np.arange(nSplit)/nSplit + rng.random(1)) - cos = rNew*np.cos(alpha) - sin = rNew*np.sin(alpha) + alpha = 2 * math.pi * (np.arange(nSplit) / nSplit + rng.random(1)) + cos = rNew * np.cos(alpha) + sin = rNew * np.sin(alpha) nx, ny, nz = geoTrans.getNormalArray(np.array([x[ind]]), np.array([y[ind]]), Nx, Ny, Nz, csz) - e1x, e1y, e1z, e2x, e2y, e2z = getTangenVectors(nx, ny, nz, np.array([ux[ind]]), np.array([uy[ind]]), np.array([uz[ind]])) + e1x, e1y, e1z, e2x, e2y, e2z = getTangenVectors( + nx, ny, nz, np.array([ux[ind]]), np.array([uy[ind]]), np.array([uz[ind]]) + ) xNew = x[ind] + cos * e1x + sin * e2x yNew = y[ind] + cos * e1y + sin * e2y zNew = z[ind] + cos * e1z + sin * e2z @@ -634,24 +647,24 @@ def getSplitPartPositionSimple(particles, rho, distSplitPart, ind): zNew : numpy array z components of the splitted particles """ - mPart = particles['m'][ind] - hPart = particles['h'][ind] - xPart = particles['x'][ind] - yPart = particles['y'][ind] - zPart = particles['z'][ind] - uxPart = particles['ux'][ind] - uyPart = particles['uy'][ind] - uzPart = particles['uz'][ind] + mPart = particles["m"][ind] + hPart = particles["h"][ind] + xPart = particles["x"][ind] + yPart = particles["y"][ind] + zPart = particles["z"][ind] + uxPart = particles["ux"][ind] + uyPart = particles["uy"][ind] + uzPart = particles["uz"][ind] # get the area of the particle as well as the distance expected between the old and new particle # note that if we did not update the particles FT, we use here the h from the previous time step - aPart = mPart/(rho*hPart) - rNew = distSplitPart * np.sqrt(aPart/math.pi) - cos = rNew*np.array([-1, 1]) + aPart = mPart / (rho * hPart) + rNew = distSplitPart * np.sqrt(aPart / math.pi) + cos = rNew * np.array([-1, 1]) # compute velocity mag to get the direction of the flow (e_1) uMag = DFAtls.norm(uxPart, uyPart, uzPart) - xNew = xPart + cos * uxPart/uMag - yNew = yPart + cos * uyPart/uMag - zNew = zPart + cos * uzPart/uMag + xNew = xPart + cos * uxPart / uMag + yNew = yPart + cos * uyPart / uMag + zNew = zPart + cos * uzPart / uMag # toDo: do we need to reproject the particles on the dem? return xNew, yNew, zNew @@ -703,7 +716,7 @@ def getTangenVectors(nx, ny, nz, ux, uy, uz): # if vector u is zero use the tangent vector in x direction for e1 e1x = np.array([1]) e1y = np.array([0]) - e1z = -nx/nz + e1z = -nx / nz e1x, e1y, e1z = DFAtls.normalize(e1x, e1y, e1z) # compute the othe tengent vector e2x, e2y, e2z = DFAtls.crossProd(nx, ny, nz, e1x, e1y, e1z) @@ -712,7 +725,7 @@ def getTangenVectors(nx, ny, nz, ux, uy, uz): def findParticles2Track(particles, center, radius): - '''Find particles within a circle arround a given point + """Find particles within a circle arround a given point Parameters ---------- @@ -732,25 +745,25 @@ def findParticles2Track(particles, center, radius): array with Parent ID of particles to track track: boolean False if no particles are tracked - ''' + """ track = True - x = particles['x'] - y = particles['y'] - xc = center['x'] - yc = center['y'] - r = DFAtls.norm(x-xc, y-yc, 0) + x = particles["x"] + y = particles["y"] + xc = center["x"] + yc = center["y"] + r = DFAtls.norm(x - xc, y - yc, 0) index = np.where(r <= radius) - particles2Track = particles['parentID'][index] - log.info('Tracking %d particles' % len(index[0])) + particles2Track = particles["parentID"][index] + log.info("Tracking %d particles" % len(index[0])) if len(index[0]) < 1: - log.warning('Found No particles to track ') + log.warning("Found No particles to track ") track = False return particles2Track, track def getTrackedParticles(particlesList, particles2Track): - '''Track particles along time given the parentID of the particles to track + """Track particles along time given the parentID of the particles to track Parameters ---------- @@ -767,21 +780,21 @@ def getTrackedParticles(particlesList, particles2Track): are tracked) nPartTracked : int total number of tracked particles - ''' + """ nPartTracked = np.size(particles2Track) # add trackedParticles array to the particles dictionary for every saved time step for particles in particlesList: # find index of particles to track - index = [ind for ind, parent in enumerate(particles['parentID']) if parent in particles2Track] + index = [ind for ind, parent in enumerate(particles["parentID"]) if parent in particles2Track] nPartTracked = max(nPartTracked, len(index)) - trackedParticles = np.zeros(particles['nPart']) + trackedParticles = np.zeros(particles["nPart"]) trackedParticles[index] = 1 - particles['trackedParticles'] = trackedParticles + particles["trackedParticles"] = trackedParticles return particlesList, nPartTracked def getTrackedParticlesProperties(particlesList, nPartTracked, properties): - '''Get the desired properties for the tracked particles + """Get the desired properties for the tracked particles Parameters ---------- @@ -800,11 +813,11 @@ def getTrackedParticlesProperties(particlesList, nPartTracked, properties): properties = ['x', 'y'], the dictionary will have the keys 't', 'x' and 'y'. trackedPartProp['x'] will be a 2D numpy array, each line corresponds to the 'x' time series of a tracked particle) - ''' + """ # buid time series for desired properties of tracked particles nTimeSteps = len(particlesList) trackedPartProp = {} - trackedPartProp['t'] = np.zeros(nTimeSteps) + trackedPartProp["t"] = np.zeros(nTimeSteps) newProperties = [] # initialize @@ -813,15 +826,15 @@ def getTrackedParticlesProperties(particlesList, nPartTracked, properties): trackedPartProp[key] = np.zeros((nTimeSteps, nPartTracked)) newProperties.append(key) else: - log.warning('%s is not a particle property' % key) + log.warning("%s is not a particle property" % key) # extract wanted properties and build the time series trackedPartID = [] for particles, nTime in zip(particlesList, range(nTimeSteps)): - trackedParticles = particles['trackedParticles'] - trackedPartProp['t'][nTime] = particles['t'] + trackedParticles = particles["trackedParticles"] + trackedPartProp["t"][nTime] = particles["t"] index = np.where(trackedParticles == 1) - for ind, id in zip(index[0], particles['ID'][index]): + for ind, id in zip(index[0], particles["ID"][index]): if id not in trackedPartID: trackedPartID.append(id) indCol = trackedPartID.index(id) @@ -831,30 +844,30 @@ def getTrackedParticlesProperties(particlesList, nPartTracked, properties): return trackedPartProp -def readPartFromPickle(inDir, simName='', flagAvaDir=False, comModule='com1DFA'): - """ Read pickles within a directory and return List of dicionaries read from pickle - - Parameters - ----------- - inDir: str - path to input directory - simName : str - simulation name - flagAvaDir: bool - if True inDir corresponds to an avalanche directory and pickles are - read from avaDir/Outputs/com1DFA/particles - comModule: str - module that computed the particles +def readPartFromPickle(inDir, simName="", flagAvaDir=False, comModule="com1DFA"): + """Read pickles within a directory and return List of dicionaries read from pickle + + Parameters + ----------- + inDir: str + path to input directory + simName : str + simulation name + flagAvaDir: bool + if True inDir corresponds to an avalanche directory and pickles are + read from avaDir/Outputs/com1DFA/particles + comModule: str + module that computed the particles """ if flagAvaDir: - inDir = pathlib.Path(inDir, 'Outputs', comModule, 'particles') + inDir = pathlib.Path(inDir, "Outputs", comModule, "particles") # search for all pickles within directory if simName: - name = '*' + simName + '*.pickle' + name = "*" + simName + "*.pickle" else: - name = '*.pickle' + name = "*.pickle" PartDicts = sorted(list(inDir.glob(name))) # initialise list of particle dictionaries @@ -863,42 +876,42 @@ def readPartFromPickle(inDir, simName='', flagAvaDir=False, comModule='com1DFA') for particles in PartDicts: particles = pickle.load(open(particles, "rb")) Particles.append(particles) - timeStepInfo.append(particles['t']) + timeStepInfo.append(particles["t"]) return Particles, timeStepInfo def savePartToCsv(particleProperties, dictList, outDir, countParticleCsv=None): - """ Save each particle dictionary from a list to a csv file; - works also for one dictionary instead of list - - Parameters - --------- - particleProperties: str - all particle properties that shall be saved to csv file - (e.g.: m, velocityMagnitude, ux,..) - dictList: list or dict - list of dictionaries or single dictionary - outDir: str - path to output directory; particlesCSV will be created in this - outDir - countParticleCsv: int - number of particlesDict to be saved according to time step - ranging from 0...n where n is the total number - of saved particleDicts - if list of dicts is saved this parameter is ignored + """Save each particle dictionary from a list to a csv file; + works also for one dictionary instead of list + + Parameters + --------- + particleProperties: str + all particle properties that shall be saved to csv file + (e.g.: m, velocityMagnitude, ux,..) + dictList: list or dict + list of dictionaries or single dictionary + outDir: str + path to output directory; particlesCSV will be created in this + outDir + countParticleCsv: int + number of particlesDict to be saved according to time step - ranging from 0...n where n is the total number + of saved particleDicts - if list of dicts is saved this parameter is ignored """ # set output directory - outDir = outDir / 'particlesCSV' + outDir = outDir / "particlesCSV" fU.makeADir(outDir) # read particle properties to be saved - particleProperties = particleProperties.split('|') + particleProperties = particleProperties.split("|") # write particles locations and properties to csv file nParticles = len(dictList) if nParticles == 1: if countParticleCsv is None: - message = 'Indicator for step in timeseries for particlesDict to be saved is required, set countPartCsv' + message = "Indicator for step in timeseries for particlesDict to be saved is required, set countPartCsv" log.error(message) raise AssertionError(message) else: @@ -910,57 +923,57 @@ def savePartToCsv(particleProperties, dictList, outDir, countParticleCsv=None): particles = dictList[0] else: particles = dictList[count] - simName = particles['simName'] + simName = particles["simName"] csvData = {} - csvData['X'] = particles['x'] + particles['xllcenter'] - csvData['Y'] = particles['y'] + particles['yllcenter'] - csvData['Z'] = particles['z'] + csvData["X"] = particles["x"] + particles["xllcenter"] + csvData["Y"] = particles["y"] + particles["yllcenter"] + csvData["Z"] = particles["z"] for partProp in particleProperties: - if partProp == 'velocityMagnitude': - ux = particles['ux'] - uy = particles['uy'] - uz = particles['uz'] + if partProp == "velocityMagnitude": + ux = particles["ux"] + uy = particles["uy"] + uz = particles["uz"] csvData[partProp] = DFAtls.norm(ux, uy, uz) else: if partProp in particles: csvData[partProp] = particles[partProp] else: - log.warning('Particle property %s does not exist' % (partProp)) - csvData['time'] = particles['t'] + log.warning("Particle property %s does not exist" % (partProp)) + csvData["time"] = particles["t"] # create pandas dataFrame and save to csv - outFile = outDir / ('particles%s.csv.%04d' % (simName, count)) + outFile = outDir / ("particles%s.csv.%04d" % (simName, count)) particlesData = pd.DataFrame(data=csvData) particlesData.to_csv(outFile, index=False) count = count + 1 def reshapeParticlesDicts(particlesList, propertyList): - """ reshape particlesList from one dict per time step with all particle properties for each particle, - to one dict with an array of property values for all time steps for each particle - shape: nx x ny; nx time steps, ny number of particles - - Parameters - ----------- - particlesList: list - list of particle dicts, one dict per time step - propertyList: list - list of property names that shall be reshaped and saved to particlesTimeArrays - - Returns - -------- - particlesTimeArrays: dict - dict with time series of properties of particles - key: property, item: timeSteps x particlesID array of property values + """reshape particlesList from one dict per time step with all particle properties for each particle, + to one dict with an array of property values for all time steps for each particle + shape: nx x ny; nx time steps, ny number of particles + + Parameters + ----------- + particlesList: list + list of particle dicts, one dict per time step + propertyList: list + list of property names that shall be reshaped and saved to particlesTimeArrays + + Returns + -------- + particlesTimeArrays: dict + dict with time series of properties of particles + key: property, item: timeSteps x particlesID array of property values """ # create particlesTimeArrays particlesTimeArrays = {} - idParticles = [p['nPart'] for p in particlesList] + idParticles = [p["nPart"] for p in particlesList] if len(list(set(idParticles))) > 1: - message = 'Number of particles changed throughout simulation' + message = "Number of particles changed throughout simulation" log.error(message) raise AssertionError(message) @@ -969,22 +982,22 @@ def reshapeParticlesDicts(particlesList, propertyList): particlesTimeArrays[props] = np.zeros(len(particlesList)) particlesTimeArrays[props][:] = np.asarray([p[props] for p in particlesList]) else: - particlesTimeArrays[props] = np.zeros((len(particlesList), particlesList[0]['nPart'])) - for idx in particlesList[0]['ID']: - particlesTimeArrays[props][:,idx] = np.asarray([p[props][idx] for p in particlesList]) + particlesTimeArrays[props] = np.zeros((len(particlesList), particlesList[0]["nPart"])) + for idx in particlesList[0]["ID"]: + particlesTimeArrays[props][:, idx] = np.asarray([p[props][idx] for p in particlesList]) return particlesTimeArrays def savePartDictToPickle(partDict, fName): - """ save a single dict to a pickle - - Parameters - ----------- - partDict: dict - dict with info - fName: pathlib path - path to saving location + """save a single dict to a pickle + + Parameters + ----------- + partDict: dict + dict with info + fName: pathlib path + path to saving location """ fi = open(fName, "wb") diff --git a/avaframe/out3Plot/plotUtils.py b/avaframe/out3Plot/plotUtils.py index 7673ba8e3..6958cea5d 100644 --- a/avaframe/out3Plot/plotUtils.py +++ b/avaframe/out3Plot/plotUtils.py @@ -202,6 +202,19 @@ colorsProb = ["#FEF1F1", "#B2AB96", "#5B8BA3", "#2D5393", "#1A0C64"] cmapProbmap = copy.copy(cmapCrameri.lapaz.reversed()) +levSfcC = list(fU.splitIniValueToArraySteps(cfgPlotUtils["surfaceChangeLevels"])) +colorsSfcC = [ + "#fefeb2", + "#d2d184", + "#bab98d", + "#a0a598", + "#6f878d", + "#386982", + "#05598c", +] +cmapSfcChange = copy.copy(cmapCrameri.nuuk.reversed()) + + ############################################### # Set colormaps to use ############################################### @@ -222,6 +235,8 @@ cmapEnergy = {"cmap": cmapE, "colors": colorsE, "levels": levE} +cmapSfcChange = {"cmap": cmapSfcChange, "colors": colorsSfcC, "levels": levSfcC} + # 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 @@ -230,6 +245,7 @@ except AttributeError: cmapZdelta = {"cmap": copy.copy(cmapCrameri.lapaz), "colors": [], "levels": []} +cmapDmDet = copy.copy(cmapCrameri.lapaz.reversed()) colorMaps = { "ppr": cmapPres, @@ -249,8 +265,10 @@ "TA": cmapTravelAngle, "pke": cmapEnergy, "zdelta": cmapZdelta, - "dmDet": cmapProb, "FTDet": cmapThickness, + "dmDet": cmapDmDet, + "demAdapted": cmapGreys, + "sfcChange": cmapSfcChange, } cmapDEM = cmapGreys diff --git a/avaframe/out3Plot/plotUtilsCfg.ini b/avaframe/out3Plot/plotUtilsCfg.ini index 6cc8cd51b..8371e55a0 100644 --- a/avaframe/out3Plot/plotUtilsCfg.ini +++ b/avaframe/out3Plot/plotUtilsCfg.ini @@ -84,6 +84,8 @@ unitTA = ° unitM = kg unitzdelta = m unitdmDet = kg +unitsfcChange = m +unitdemAdapted = m # threshold levels elevMaxppr = 100 elevMaxpft = 1 @@ -108,6 +110,8 @@ energyColorLevels = 5|10|100|500|1000 travelAngleColorLevels = 28|29|30|31|32|33|34 # for proba, provide 5 levels probaColorLevels = 0|0.25|0.50|0.75|1. +# for thicknessChange +surfaceChangeLevels = -1|-0.3|-0.1|0|0.1|0.3|1 # 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 diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py index 210c7f2ec..047935c96 100644 --- a/avaframe/tests/test_DFAfunctionsCython.py +++ b/avaframe/tests/test_DFAfunctionsCython.py @@ -211,6 +211,10 @@ def test_updatePositionC(): "dissDam": "1", "snowSlide": "1", "wetSnow": "1", + "delStoppedParticles": "0", + "adaptSfcStopped": "0", + "adaptSfcDetrainment": "0", + "adaptSfcEntrainment": "0", } particles = { @@ -236,6 +240,9 @@ def test_updatePositionC(): "iterate": True, "totalEnthalpy": np.asarray([0.0, 0.0, 0.0]), "velocityMag": np.asarray([1.0, 1.0, 1.0]), + "h": np.asarray([1.0, 1.0, 1.0]), + "ID": np.array([0, 1, 2]), + "stoppedParticles": {}, } particles["potentialEne"] = np.sum(9.81 * particles["z"] * particles["m"]) @@ -344,6 +351,9 @@ def test_updatePositionC(): "peakMassFlowing": 0.0, "iterate": True, "totalEnthalpy": np.asarray([0.0, 0.0, 0.0]), + "h": np.asarray([1.0, 1.0, 1.0]), + "ID": np.array([0, 1, 2]), + "stoppedParticles": {}, } particles["potentialEne"] = np.sum(9.81 * particles["z"] * particles["m"]) @@ -387,6 +397,9 @@ def test_updatePositionC(): "peakMassFlowing": 0.0, "iterate": True, "totalEnthalpy": np.asarray([0.0, 0.0, 0.0]), + "h": np.asarray([1.0, 1.0, 1.0]), + "ID": np.array([0, 1, 2]), + "stoppedParticles": {}, } particles["potentialEne"] = np.sum(9.81 * particles["z"] * particles["m"]) typeStop = 1 @@ -450,6 +463,9 @@ def test_updatePositionC(): "peakMassFlowing": 0.0, "iterate": True, "totalEnthalpy": np.asarray([0.0, 0.0, 0.0]), + "h": np.asarray([1.0, 1.0, 1.0]), + "ID": np.array([0, 1, 2]), + "stoppedParticles": {}, } particles["potentialEne"] = np.sum(9.81 * particles["z"] * particles["m"]) typeStop = 1 @@ -480,6 +496,62 @@ def test_updatePositionC(): assert (potEneNew - 1.0e-4) < particles["potentialEne"] < (potEneNew + 1.0e-4) assert particles["iterate"] == True + particles = { + "dt": 1.0, + "m": np.asarray([10.0, 10.0, 10.0]), + "idFixed": np.asarray([0.0, 0.0, 0.0]), + "trajectoryLengthXY": np.asarray([0.0, 0.0, 0.0]), + "trajectoryLengthXYCor": np.asarray([0.0, 0.0, 0.0]), + "trajectoryLengthXYZ": np.asarray([0.0, 0.0, 0.0]), + "x": np.asarray([0.0, 1.0, 2.0]), + "y": np.asarray([2.0, 3.0, 4.0]), + "z": np.asarray([1.0, 1.0, 1.0]), + "ux": np.asarray([0.0, 1.0, 1.0]), + "uy": np.asarray([0.0, 1.0, 1.0]), + "uz": np.asarray([0.0, 0.0, 0.0]), + "uAcc": np.asarray([0.0, 0.0, 0.0]), + "kineticEne": 0.0, + "peakKinEne": 0.0, + "peakForceSPH": 0.0, + "forceSPHIni": 0.0, + "nPart": 3, + "peakMassFlowing": 0.0, + "iterate": True, + "totalEnthalpy": np.asarray([0.0, 0.0, 0.0]), + "velocityMag": np.asarray([0.0, 1.0, 1.0]), + "h": np.asarray([1.0, 1.0, 1.0]), + "ID": np.array([0, 1, 2]), + "stoppedParticles": {}, + } + particles["potentialEne"] = np.sum(9.81 * particles["z"] * particles["m"]) + + force = { + "forceZ": np.asarray([0.0, 0.0, 0.0]), + "forceFrict": np.asarray([10.0, 10.0, 10.0]), + "dM": np.asarray([0.0, 0.0, 0.0]), + "dMDet": np.asarray([0.0, 0.0, 0.0]), + "forceX": np.asarray([0.0, 50.0, 50.0]), + "forceY": np.asarray([0.0, 50.0, 50.0]), + "forceSPHX": np.asarray([0.0, 50.0, 50.0]), + "forceSPHY": np.asarray([0.0, 50.0, 50.0]), + "forceSPHZ": np.asarray([0.0, 0.0, 0.0]), + } + + cfg["GENERAL"]["delStoppedParticles"] = "1" + cfg["GENERAL"]["explicitFriction"] = "1" + cfg["GENERAL"]["snowSlide"] = "0" + typeStop = 1 + + particles = DFAfunC.updatePositionC(cfg["GENERAL"], particles, dem, force, fields, typeStop=typeStop) + assert np.array_equal(particles["stoppedParticles"]["m"], np.array([10.0])) + assert particles["stoppedParticles"]["x"] == 0.0 + assert particles["stoppedParticles"]["y"] == 2.0 + assert particles["stoppedParticles"]["ID"] == 0 + assert np.array_equal(particles["ID"], np.array([1, 2])) + assert particles["nPart"] == 2 + assert np.all(particles["velocityMag"] > 0) + assert particles["massStopped"] == -np.sum(particles["stoppedParticles"]["m"]) + def test_computeTrajectoryAngle(): # first compute travel angle for each particle diff --git a/avaframe/tests/test_checkCfg.py b/avaframe/tests/test_checkCfg.py index be46e49c1..36ec4a6ac 100644 --- a/avaframe/tests/test_checkCfg.py +++ b/avaframe/tests/test_checkCfg.py @@ -1,5 +1,5 @@ """ - Pytest for module com1DFA +Pytest for module com1DFA """ # Load modules @@ -14,158 +14,263 @@ def test_checkCfgConsistency(): - """ test check if sphOption and viscOption settings are consistent """ + """test check if sphOption and viscOption settings are consistent""" # setup requuired input data cfg = configparser.ConfigParser() - cfg['GENERAL'] = {'viscOption': 2, 'sphOption': 2} + cfg["GENERAL"] = {"viscOption": 2, "sphOption": 2, "methodMeshNormal": 1} # call function to be tested testFlag = checkCfg.checkCfgConsistency(cfg) assert testFlag == True - cfg['GENERAL'] = {'viscOption': 2, 'sphOption': 1} + cfg["GENERAL"] = {"viscOption": 2, "sphOption": 1, "methodMeshNormal": 1} # call function to be tested with pytest.raises(AssertionError) as e: assert checkCfg.checkCfgConsistency(cfg) - assert 'If viscOption is set to 2' in str(e.value) + assert "If viscOption is set to 2" in str(e.value) + + cfg["GENERAL"] = {"viscOption": 2, "sphOption": 2, "methodMeshNormal": 2} + + # call function to be tested + with pytest.raises(AssertionError) as e: + assert checkCfg.checkCfgConsistency(cfg) + assert "The methodMeshNormal set in the ini file does NOT have" in str(e.value) def test_checkCellSizeKernelRadius(): - """ test check if sphOption and viscOption settings are consistent """ + """test check if sphOption and viscOption settings are consistent""" # setup requuired input data - cfg = {'GENERAL': {'sphKernelRadius': '5', 'meshCellSize': '10'}} + cfg = {"GENERAL": {"sphKernelRadius": "5", "meshCellSize": "10"}} # call function to be tested cfg = checkCfg.checkCellSizeKernelRadius(cfg) - assert cfg['GENERAL']['sphKernelRadius'] == '5' + assert cfg["GENERAL"]["sphKernelRadius"] == "5" # setup requuired input data - cfg = {'GENERAL': {'sphKernelRadius': 'meshCellSize', 'meshCellSize': '10'}} + cfg = {"GENERAL": {"sphKernelRadius": "meshCellSize", "meshCellSize": "10"}} # call function to be tested cfg = checkCfg.checkCellSizeKernelRadius(cfg) - assert cfg['GENERAL']['sphKernelRadius'] == '10' + assert cfg["GENERAL"]["sphKernelRadius"] == "10" def test_checkCfgFrictionModel(): - """ test if friction model parameters are set to nan if not used """ - - cfg = {'GENERAL': {'frictModel': 'samosAT', 'musamosat': '0.155', 'tau0samosat': '0', 'Rs0samosat': '0.222', - 'kappasamosat': '0.43', 'Rsamosat': '0.05', 'Bsamosat': '4.13', - 'muvoellmy': '4000.', 'xsivoellmy': '4000.'}} + """test if friction model parameters are set to nan if not used""" + + cfg = { + "GENERAL": { + "frictModel": "samosAT", + "musamosat": "0.155", + "tau0samosat": "0", + "Rs0samosat": "0.222", + "kappasamosat": "0.43", + "Rsamosat": "0.05", + "Bsamosat": "4.13", + "muvoellmy": "4000.", + "xsivoellmy": "4000.", + } + } inputSimFiles = {} cfg = checkCfg.checkCfgFrictionModel(cfg, inputSimFiles) - assert cfg['GENERAL']['frictModel'] == 'samosAT' - assert float(cfg['GENERAL']['musamosat']) == 0.155 - assert float(cfg['GENERAL']['tau0samosat']) == 0. - assert float(cfg['GENERAL']['kappasamosat']) == 0.43 - assert float(cfg['GENERAL']['Rsamosat']) == 0.05 - assert float(cfg['GENERAL']['Bsamosat']) == 4.13 - assert float(cfg['GENERAL']['muvoellmy']) == 4000. - assert float(cfg['GENERAL']['xsivoellmy']) == 4000. - - - cfg = {'GENERAL': {'frictModel': 'samosAT', 'musamosat': '0.155', 'tau0samosat': '0', 'Rs0samosat': '0.222', - 'kappasamosat': '0.43', 'Rsamosat': '0.05', 'Bsamosat': '9.13', - 'muvoellmy': '4000.', 'xsivoellmy': 'nan'}} + assert cfg["GENERAL"]["frictModel"] == "samosAT" + assert float(cfg["GENERAL"]["musamosat"]) == 0.155 + assert float(cfg["GENERAL"]["tau0samosat"]) == 0.0 + assert float(cfg["GENERAL"]["kappasamosat"]) == 0.43 + assert float(cfg["GENERAL"]["Rsamosat"]) == 0.05 + assert float(cfg["GENERAL"]["Bsamosat"]) == 4.13 + assert float(cfg["GENERAL"]["muvoellmy"]) == 4000.0 + assert float(cfg["GENERAL"]["xsivoellmy"]) == 4000.0 + + cfg = { + "GENERAL": { + "frictModel": "samosAT", + "musamosat": "0.155", + "tau0samosat": "0", + "Rs0samosat": "0.222", + "kappasamosat": "0.43", + "Rsamosat": "0.05", + "Bsamosat": "9.13", + "muvoellmy": "4000.", + "xsivoellmy": "nan", + } + } cfg = checkCfg.checkCfgFrictionModel(cfg, inputSimFiles) - assert cfg['GENERAL']['frictModel'] == 'samosAT' - assert float(cfg['GENERAL']['musamosat']) == 0.155 - assert float(cfg['GENERAL']['tau0samosat']) == 0 - assert float(cfg['GENERAL']['kappasamosat']) == 0.43 - assert float(cfg['GENERAL']['Rsamosat']) == 0.05 - assert float(cfg['GENERAL']['Bsamosat']) == 9.13 - assert float(cfg['GENERAL']['muvoellmy']) == 4000. - assert np.isnan(float(cfg['GENERAL']['xsivoellmy'])) - - cfg = {'GENERAL': {'frictModel': 'samosAT', 'musamosat': '0.155', 'tau0samosat': '0', 'Rs0samosat': '0.222', - 'kappasamosat': '0.43', 'Rsamosat': '0.05', 'Bsamosat': 'nan', - 'muvoellmy': '4000.', 'xsivoellmy': 'nan'}} + assert cfg["GENERAL"]["frictModel"] == "samosAT" + assert float(cfg["GENERAL"]["musamosat"]) == 0.155 + assert float(cfg["GENERAL"]["tau0samosat"]) == 0 + assert float(cfg["GENERAL"]["kappasamosat"]) == 0.43 + assert float(cfg["GENERAL"]["Rsamosat"]) == 0.05 + assert float(cfg["GENERAL"]["Bsamosat"]) == 9.13 + assert float(cfg["GENERAL"]["muvoellmy"]) == 4000.0 + assert np.isnan(float(cfg["GENERAL"]["xsivoellmy"])) + + cfg = { + "GENERAL": { + "frictModel": "samosAT", + "musamosat": "0.155", + "tau0samosat": "0", + "Rs0samosat": "0.222", + "kappasamosat": "0.43", + "Rsamosat": "0.05", + "Bsamosat": "nan", + "muvoellmy": "4000.", + "xsivoellmy": "nan", + } + } # call function to be tested with pytest.raises(ValueError) as e: assert checkCfg.checkCfgFrictionModel(cfg, inputSimFiles) - assert 'Friction model used samosAT, but Bsamosat is nan - not valid' in str(e.value) - - cfg = {'GENERAL': {'frictModel': 'samosAT', 'musamosat': '0.155', 'tau0samosat': '0', 'Rs0samosat': '0.222', - 'kappasamosat': '0.43', 'Rsamosat': '0.05', 'Bsamosat': 'test', - 'muvoellmy': '4000.', 'xsivoellmy': 'nan'}} + assert "Friction model used samosAT, but Bsamosat is nan - not valid" in str(e.value) + + cfg = { + "GENERAL": { + "frictModel": "samosAT", + "musamosat": "0.155", + "tau0samosat": "0", + "Rs0samosat": "0.222", + "kappasamosat": "0.43", + "Rsamosat": "0.05", + "Bsamosat": "test", + "muvoellmy": "4000.", + "xsivoellmy": "nan", + } + } # call function to be tested with pytest.raises(ValueError) as e: assert checkCfg.checkCfgFrictionModel(cfg, inputSimFiles) - assert 'Friction model used samosAT, but test is not of valid' in str(e.value) - - - cfg = {'GENERAL': {'frictModel': 'samosATAuto', 'musamosat': '0.155', 'tau0samosat': '0', 'Rs0samosat': '0.222', - 'kappasamosat': '0.43', 'Rsamosat': '0.05', 'Bsamosat': '9.13', - 'muvoellmy': '4000.', 'xsivoellmy': 'nan', 'musamosatsmall': '0.22', 'tau0samosatsmall': '0', - 'Rs0samosatsmall': '0.222', 'kappasamosatsmall': '0.43', 'Rsamosatsmall': '0.05', - 'Bsamosatsmall': '4.13', 'musamosatmedium': '0.17', 'tau0samosatmedium': '0', - 'Rs0samosatmedium': '0.222', 'kappasamosatmedium': '0.43', 'Rsamosatmedium': '0.05', - 'Bsamosatmedium': '4.13', 'volClassSmall': '25000.', 'volClassMedium': '60000.', - 'meshCellSize': '5'}} - - cfg = checkCfg.checkCfgFrictionModel(cfg, inputSimFiles, relVolume=24999.) - - assert cfg['GENERAL']['frictModel'] == 'samosATSmall' - assert float(cfg['GENERAL']['musamosatsmall']) == 0.22 - assert float(cfg['GENERAL']['tau0samosatsmall']) == 0 - assert float(cfg['GENERAL']['Rs0samosatsmall']) == 0.222 - assert float(cfg['GENERAL']['kappasamosatsmall']) == 0.43 - assert float(cfg['GENERAL']['Rsamosatsmall']) == 0.05 - assert float(cfg['GENERAL']['Bsamosatsmall']) == 4.13 - assert float(cfg['GENERAL']['muvoellmy']) == 4000. - assert np.isnan(float(cfg['GENERAL']['xsivoellmy'])) - - cfg = {'GENERAL': {'frictModel': 'samosATAuto', 'musamosat': '0.155', 'tau0samosat': '0', 'Rs0samosat': '0.222', - 'kappasamosat': '0.43', 'Rsamosat': '0.05', 'Bsamosat': '9.13', - 'muvoellmy': '4000.', 'xsivoellmy': 'nan', 'musamosatsmall': '0.22', 'tau0samosatsmall': '0', - 'Rs0samosatsmall': '0.222', 'kappasamosatsmall': '0.43', 'Rsamosatsmall': '0.05', - 'Bsamosatsmall': '4.13', 'musamosatmedium': '0.17', 'tau0samosatmedium': '0.9', - 'Rs0samosatmedium': '1.222', 'kappasamosatmedium': '2.43', 'Rsamosatmedium': '0.75', - 'Bsamosatmedium': '4.23', 'volClassSmall': '25000.', 'volClassMedium': '60000.', - 'meshCellSize': '5'}} - - cfg = checkCfg.checkCfgFrictionModel(cfg, inputSimFiles, relVolume=26999.) - - assert cfg['GENERAL']['frictModel'] == 'samosATMedium' - assert float(cfg['GENERAL']['musamosatmedium']) == 0.17 - assert float(cfg['GENERAL']['tau0samosatmedium']) == 0.9 - assert float(cfg['GENERAL']['Rs0samosatmedium']) == 1.222 - assert float(cfg['GENERAL']['kappasamosatmedium']) == 2.43 - assert float(cfg['GENERAL']['Rsamosatmedium']) == 0.75 - assert float(cfg['GENERAL']['Bsamosatmedium']) == 4.23 - assert float(cfg['GENERAL']['muvoellmy']) == 4000. - assert np.isnan(float(cfg['GENERAL']['xsivoellmy'])) - - cfg = {'GENERAL': {'frictModel': 'samosATAuto', 'musamosat': '0.155', 'tau0samosat': '0.8', 'Rs0samosat': '0.2227', - 'kappasamosat': '1.43', 'Rsamosat': '1.05', 'Bsamosat': '9.13', - 'muvoellmy': '4000.', 'xsivoellmy': 'nan', 'musamosatsmall': '0.22', 'tau0samosatsmall': '0', - 'Rs0samosatsmall': '0.222', 'kappasamosatsmall': '0.43', 'Rsamosatsmall': '0.05', - 'Bsamosatsmall': '4.13', 'musamosatmedium': '0.17', 'tau0samosatmedium': '0', - 'Rs0samosatmedium': '0.222', 'kappasamosatmedium': '0.43', 'Rsamosatmedium': '0.05', - 'Bsamosatmedium': '4.13', 'volClassSmall': '25000.', 'volClassMedium': '60000.', - 'meshCellSize': '5'}} - - cfg = checkCfg.checkCfgFrictionModel(cfg, inputSimFiles, relVolume=74999.) - - assert cfg['GENERAL']['frictModel'] == 'samosAT' - assert float(cfg['GENERAL']['musamosat']) == 0.155 - assert float(cfg['GENERAL']['tau0samosat']) == 0.8 - assert float(cfg['GENERAL']['Rs0samosat']) == 0.2227 - assert float(cfg['GENERAL']['kappasamosat']) == 1.43 - assert float(cfg['GENERAL']['Rsamosat']) == 1.05 - assert float(cfg['GENERAL']['Bsamosat']) == 9.13 - assert float(cfg['GENERAL']['muvoellmy']) == 4000. - assert np.isnan(float(cfg['GENERAL']['xsivoellmy'])) + assert "Friction model used samosAT, but test is not of valid" in str(e.value) + + cfg = { + "GENERAL": { + "frictModel": "samosATAuto", + "musamosat": "0.155", + "tau0samosat": "0", + "Rs0samosat": "0.222", + "kappasamosat": "0.43", + "Rsamosat": "0.05", + "Bsamosat": "9.13", + "muvoellmy": "4000.", + "xsivoellmy": "nan", + "musamosatsmall": "0.22", + "tau0samosatsmall": "0", + "Rs0samosatsmall": "0.222", + "kappasamosatsmall": "0.43", + "Rsamosatsmall": "0.05", + "Bsamosatsmall": "4.13", + "musamosatmedium": "0.17", + "tau0samosatmedium": "0", + "Rs0samosatmedium": "0.222", + "kappasamosatmedium": "0.43", + "Rsamosatmedium": "0.05", + "Bsamosatmedium": "4.13", + "volClassSmall": "25000.", + "volClassMedium": "60000.", + "meshCellSize": "5", + } + } + + cfg = checkCfg.checkCfgFrictionModel(cfg, inputSimFiles, relVolume=24999.0) + + assert cfg["GENERAL"]["frictModel"] == "samosATSmall" + assert float(cfg["GENERAL"]["musamosatsmall"]) == 0.22 + assert float(cfg["GENERAL"]["tau0samosatsmall"]) == 0 + assert float(cfg["GENERAL"]["Rs0samosatsmall"]) == 0.222 + assert float(cfg["GENERAL"]["kappasamosatsmall"]) == 0.43 + assert float(cfg["GENERAL"]["Rsamosatsmall"]) == 0.05 + assert float(cfg["GENERAL"]["Bsamosatsmall"]) == 4.13 + assert float(cfg["GENERAL"]["muvoellmy"]) == 4000.0 + assert np.isnan(float(cfg["GENERAL"]["xsivoellmy"])) + + cfg = { + "GENERAL": { + "frictModel": "samosATAuto", + "musamosat": "0.155", + "tau0samosat": "0", + "Rs0samosat": "0.222", + "kappasamosat": "0.43", + "Rsamosat": "0.05", + "Bsamosat": "9.13", + "muvoellmy": "4000.", + "xsivoellmy": "nan", + "musamosatsmall": "0.22", + "tau0samosatsmall": "0", + "Rs0samosatsmall": "0.222", + "kappasamosatsmall": "0.43", + "Rsamosatsmall": "0.05", + "Bsamosatsmall": "4.13", + "musamosatmedium": "0.17", + "tau0samosatmedium": "0.9", + "Rs0samosatmedium": "1.222", + "kappasamosatmedium": "2.43", + "Rsamosatmedium": "0.75", + "Bsamosatmedium": "4.23", + "volClassSmall": "25000.", + "volClassMedium": "60000.", + "meshCellSize": "5", + } + } + + cfg = checkCfg.checkCfgFrictionModel(cfg, inputSimFiles, relVolume=26999.0) + + assert cfg["GENERAL"]["frictModel"] == "samosATMedium" + assert float(cfg["GENERAL"]["musamosatmedium"]) == 0.17 + assert float(cfg["GENERAL"]["tau0samosatmedium"]) == 0.9 + assert float(cfg["GENERAL"]["Rs0samosatmedium"]) == 1.222 + assert float(cfg["GENERAL"]["kappasamosatmedium"]) == 2.43 + assert float(cfg["GENERAL"]["Rsamosatmedium"]) == 0.75 + assert float(cfg["GENERAL"]["Bsamosatmedium"]) == 4.23 + assert float(cfg["GENERAL"]["muvoellmy"]) == 4000.0 + assert np.isnan(float(cfg["GENERAL"]["xsivoellmy"])) + + cfg = { + "GENERAL": { + "frictModel": "samosATAuto", + "musamosat": "0.155", + "tau0samosat": "0.8", + "Rs0samosat": "0.2227", + "kappasamosat": "1.43", + "Rsamosat": "1.05", + "Bsamosat": "9.13", + "muvoellmy": "4000.", + "xsivoellmy": "nan", + "musamosatsmall": "0.22", + "tau0samosatsmall": "0", + "Rs0samosatsmall": "0.222", + "kappasamosatsmall": "0.43", + "Rsamosatsmall": "0.05", + "Bsamosatsmall": "4.13", + "musamosatmedium": "0.17", + "tau0samosatmedium": "0", + "Rs0samosatmedium": "0.222", + "kappasamosatmedium": "0.43", + "Rsamosatmedium": "0.05", + "Bsamosatmedium": "4.13", + "volClassSmall": "25000.", + "volClassMedium": "60000.", + "meshCellSize": "5", + } + } + + cfg = checkCfg.checkCfgFrictionModel(cfg, inputSimFiles, relVolume=74999.0) + + assert cfg["GENERAL"]["frictModel"] == "samosAT" + assert float(cfg["GENERAL"]["musamosat"]) == 0.155 + assert float(cfg["GENERAL"]["tau0samosat"]) == 0.8 + assert float(cfg["GENERAL"]["Rs0samosat"]) == 0.2227 + assert float(cfg["GENERAL"]["kappasamosat"]) == 1.43 + assert float(cfg["GENERAL"]["Rsamosat"]) == 1.05 + assert float(cfg["GENERAL"]["Bsamosat"]) == 9.13 + assert float(cfg["GENERAL"]["muvoellmy"]) == 4000.0 + assert np.isnan(float(cfg["GENERAL"]["xsivoellmy"])) diff --git a/avaframe/tests/test_com1DFA.py b/avaframe/tests/test_com1DFA.py index 77be236d2..369f5b467 100644 --- a/avaframe/tests/test_com1DFA.py +++ b/avaframe/tests/test_com1DFA.py @@ -21,6 +21,8 @@ from avaframe.com1DFA import com1DFA from avaframe.in2Trans.rasterUtils import transformFromASCHeader from avaframe.in3Utils import cfgUtils +import avaframe.in3Utils.geoTrans as geoTrans +import avaframe.com1DFA.DFAtools as DFAtls def test_prepareInputData(tmp_path): @@ -700,7 +702,7 @@ def test_initializeMassEnt(): "cpIce": "2050.", "TIni": "-10.", } - cfg['EXPORTS'] = {'exportRasters': False} + cfg["EXPORTS"] = {"exportRasters": False} simTypeActual = "entres" dirName = pathlib.Path(__file__).parents[0] @@ -1340,6 +1342,7 @@ def test_initializeParticles(): "entTempRef": "-10.", "cpIce": "2050.", "TIni": "-10.", + "rhoEnt": "200.", } demHeader = {} demHeader["cellsize"] = 1 @@ -1421,6 +1424,8 @@ def test_initializeParticles(): "velocityMag", "nExitedParticles", "tPlot", + "dmEnt", + "stoppedParticles", ] # call function to be tested @@ -1543,8 +1548,10 @@ def test_writeMBFile(tmp_path): infoDict = {"timeStep": np.asarray([0, 1, 2, 3, 4])} infoDict["massEntrained"] = np.asarray([0, 0, 10, 20, 30]) infoDict["massDetrained"] = np.asarray([0, 0, 0, 0, 0]) + infoDict["massDetrainedTotal"] = np.asarray([0, 0, 0, 0, 0]) infoDict["massTotal"] = np.asarray([60.0, 60.0, 70.0, 90.0, 120.0]) infoDict["pfvTimeMax"] = np.asarray([0, 0, 0, 0, 0]) + infoDict["massStopped"] = np.asarray([0.0, 0.0, 0.0, 0.0, 0.0]) avaName = "data/avaTest" avaDir = pathlib.Path(tmp_path, avaName) logName = "simTestName" @@ -1560,13 +1567,17 @@ def test_writeMBFile(tmp_path): assert np.array_equal(mbInfo[:, 0], infoDict["timeStep"]) assert np.array_equal(mbInfo[:, 2], infoDict["massEntrained"]) assert np.array_equal(mbInfo[:, 3], infoDict["massDetrained"]) + assert np.array_equal(mbInfo[:, 4], infoDict["massDetrainedTotal"]) assert np.array_equal(mbInfo[:, 1], infoDict["massTotal"]) + assert np.array_equal(mbInfo[:, 5], infoDict["massStopped"]) assert mbInfo.shape[0] == 5 - assert mbInfo.shape[1] == 5 + assert mbInfo.shape[1] == 6 infoDict["massEntrained"] = np.asarray([0, 0, 0, 0, 0]) infoDict["massDetrained"] = np.asarray([0, 10, 0, 30, 0]) infoDict["massTotal"] = np.asarray([60.0, 50.0, 50.0, 20.0, 20.0]) + infoDict["massStopped"] = np.asarray([10.0, 10.0, 10.0, 50.0, 0.0]) + infoDict["massDetrainedTotal"] = np.asarray([0, 10, 10, 40, 40]) com1DFA.writeMBFile(infoDict, avaDir, logName) mbFilePath = avaDir / "Outputs" / "com1DFA" / "mass_simTestName.txt" @@ -1575,13 +1586,17 @@ def test_writeMBFile(tmp_path): assert np.array_equal(mbInfo[:, 0], infoDict["timeStep"]) assert np.array_equal(mbInfo[:, 2], infoDict["massEntrained"]) assert np.array_equal(mbInfo[:, 3], infoDict["massDetrained"]) + assert np.array_equal(mbInfo[:, 4], infoDict["massDetrainedTotal"]) assert np.array_equal(mbInfo[:, 1], infoDict["massTotal"]) + assert np.array_equal(mbInfo[:, 5], infoDict["massStopped"]) assert mbInfo.shape[0] == 5 - assert mbInfo.shape[1] == 5 + assert mbInfo.shape[1] == 6 infoDict["massEntrained"] = np.asarray([0, 20, 0, 0, 10]) infoDict["massDetrained"] = np.asarray([0, 10, 0, 30, 0]) + infoDict["massDetrainedTotal"] = np.asarray([0, 10, 10, 40, 40]) infoDict["massTotal"] = np.asarray([60.0, 70.0, 70.0, 40.0, 50.0]) + infoDict["massStopped"] = np.asarray([0, 10, 0, 30, 0]) com1DFA.writeMBFile(infoDict, avaDir, logName) mbFilePath = avaDir / "Outputs" / "com1DFA" / "mass_simTestName.txt" @@ -1590,9 +1605,11 @@ def test_writeMBFile(tmp_path): assert np.array_equal(mbInfo[:, 0], infoDict["timeStep"]) assert np.array_equal(mbInfo[:, 2], infoDict["massEntrained"]) assert np.array_equal(mbInfo[:, 3], infoDict["massDetrained"]) + assert np.array_equal(mbInfo[:, 4], infoDict["massDetrainedTotal"]) assert np.array_equal(mbInfo[:, 1], infoDict["massTotal"]) + assert np.array_equal(mbInfo[:, 5], infoDict["massStopped"]) assert mbInfo.shape[0] == 5 - assert mbInfo.shape[1] == 5 + assert mbInfo.shape[1] == 6 def test_savePartToPickle(tmp_path): @@ -1801,11 +1818,18 @@ def test_initializeFields(): "uz": np.asarray([0.0, 0.0, 0.0]), "m": np.asarray([10.0, 10.0, 10.0]), "dmDet": np.asarray([0.0, 0.0, 0.0]), + "dmEnt": np.asarray([0.0, 0.0, 0.0]), "trajectoryAngle": np.asarray([0.0, 0.0, 0.0]), + "stoppedParticles": {"m": np.empty(0), "x": np.empty(0), "y": np.empty(0)}, } cfg = configparser.ConfigParser() cfg["REPORT"] = {"plotFields": "ppr|pft|pfv"} - cfg["GENERAL"] = {"rho": "200.", "interpOption": "2", "resType": "ppr|pft|pfv"} + cfg["GENERAL"] = { + "rho": "200.", + "interpOption": "2", + "resType": "ppr|pft|pfv", + "rhoEnt": 200, + } dem["originalHeader"] = dem["header"] dem["header"]["xllcenter"] = 0.0 @@ -1820,7 +1844,7 @@ def test_initializeFields(): # print("compute TA", fields["computeTA"]) # print("compute P", fields["computeP"]) - assert len(fields) == 17 + assert len(fields) == 23 assert fields["computeTA"] is False assert fields["computeKE"] is False assert fields["computeP"] @@ -1837,12 +1861,22 @@ def test_initializeFields(): assert np.sum(fields["FT"]) != 0.0 assert np.sum(fields["FM"]) != 0.0 assert np.sum(fields["dmDet"]) == 0.0 + assert np.sum(fields["sfcChange"]) == 0.0 + assert np.sum(fields["demAdapted"]) == 0.0 + assert np.sum(fields["FTDet"]) == 0.0 + assert np.sum(fields["FTStop"]) == 0.0 + assert np.sum(fields["FTEnt"]) == 0.0 cfg["REPORT"] = {"plotFields": "pft|pfv"} - cfg["GENERAL"] = {"resType": "pke|pta|pft|pfv", "rho": "200.", "interpOption": "2"} + cfg["GENERAL"] = { + "resType": "pke|pta|pft|pfv", + "rho": "200.", + "interpOption": "2", + "rhoEnt": 200, + } # call function to be tested particles, fields = com1DFA.initializeFields(cfg, dem, particles, "") - assert len(fields) == 17 + assert len(fields) == 23 assert fields["computeTA"] assert fields["computeKE"] assert fields["computeP"] is False @@ -2123,7 +2157,7 @@ def test_initializeSimulation(tmp_path): "TIni": "-10.", "ResistanceModel": "cRes", } - cfg['EXPORTS'] = {'exportRasters': "False"} + cfg["EXPORTS"] = {"exportRasters": "False"} # setup dem input demHeader = {} demHeader["xllcenter"] = 1.0 @@ -2293,7 +2327,7 @@ def test_initializeSimulation(tmp_path): "restitutionCoefficient": 1, "nIterDam": 1, } - cfg['EXPORTS'] = {'exportRasters': "False"} + cfg["EXPORTS"] = {"exportRasters": "False"} releaseLine = { "x": np.asarray([6.9, 8.5, 8.5, 6.9, 6.9]), "y": np.asarray([7.9, 7.9, 9.5, 9.5, 7.9]), @@ -2415,6 +2449,9 @@ def test_runCom1DFA(tmp_path, caplog): "dmDet", "massDetrained", "tPlot", + "dmEnt", + "massStopped", + "stoppedParticles", ] # read one particles dictionary @@ -2598,3 +2635,167 @@ def test_fetchRelVolume(tmp_path): relVolume = com1DFA.fetchRelVolume(rel1, cfg, demPath, None) assert relVolume == 38.0 + + +def test_adaptDEM(): + cfg = configparser.ConfigParser() + cfg["GENERAL"] = { + "methodMeshNormal": 1, + "adaptSfcStopped": 0, + "adaptSfcDetrainment": 0, + "adaptSfcEntrainment": 0, + } + + header = { + "nrows": 5, + "ncols": 5, + "cellsize": 5, + } + + data = np.array( + [ + [1.0, 2.0, 3.0, 4.0, 5.0], + [1.0, 2.0, 3.0, 4.0, 5.0], + [1.0, 2.0, 3.0, 4.0, 5.0], + [1.0, 2.0, 3.0, 4.0, 5.0], + [1.0, 2.0, 3.0, 4.0, 5.0], + ] + ) + dem = { + "header": header, + "rasterData": data, + "originalHeader": {}, + "originalRasterData": data, + "headerNeighbourGrid": {}, + "damLine": {}, + } + + fields = { + "FTDet": np.zeros_like(data), + "FTStop": np.zeros_like(data), + "FTEnt": np.zeros_like(data), + "demAdapted": data, + "sfcChangeTotal": np.zeros_like(data), + "sfcChange": np.zeros_like(data), + } + + dem = geoTrans.getNormalMesh(dem, num=cfg["GENERAL"].getfloat("methodMeshNormal")) + dem = DFAtls.getAreaMesh(dem, cfg["GENERAL"].getfloat("methodMeshNormal")) + + _, _, NzNormed = DFAtls.normalize(dem["Nx"].copy(), dem["Ny"].copy(), dem["Nz"].copy()) + + demInput = dem.copy() + fieldsInput = fields.copy() + + demAdapted, fieldsAdapted = com1DFA.adaptDEM(demInput, fieldsInput, cfg["GENERAL"]) + for key in demAdapted.keys(): + assert np.all(demAdapted[key] == dem[key]) + for key in fieldsAdapted.keys(): + assert np.all(fieldsAdapted[key] == fields[key]) + + cfg = configparser.ConfigParser() + cfg["GENERAL"] = { + "methodMeshNormal": 1, + "adaptSfcStopped": 1, + "adaptSfcDetrainment": 1, + "adaptSfcEntrainment": 1, + } + + # all rasters for depth changes are zero + demAdapted, fieldsAdapted = com1DFA.adaptDEM(demInput, fieldsInput, cfg["GENERAL"]) + for key in demAdapted.keys(): + assert np.all(demAdapted[key] == dem[key]) + for key in fieldsAdapted.keys(): + assert np.all(fieldsAdapted[key] == fields[key]) + + fields["FTDet"] += 1 + fieldsInput = fields.copy() + demInput = dem.copy() + + demAdapted, fieldsAdapted = com1DFA.adaptDEM(demInput, fieldsInput, cfg["GENERAL"]) + + for key in demAdapted.keys(): + if key == "rasterData": + assert np.all(demAdapted[key] == dem[key] + 1 / NzNormed) + elif key == "header": + assert np.all(demAdapted[key] == dem[key]) + + for key in fieldsAdapted.keys(): + if key == "demAdapted": + assert np.all(fieldsAdapted[key] == fields[key] + 1 / NzNormed) + elif key == "sfcChange": + assert np.all(fieldsAdapted[key] == fields["FTDet"] / NzNormed) + elif key == "sfcChangeTotal": + assert np.all(fieldsAdapted[key] == fields["FTDet"] / NzNormed) + else: + assert np.all(fieldsAdapted[key] == fields[key]) + + cfg = configparser.ConfigParser() + cfg["GENERAL"] = { + "methodMeshNormal": 1, + "adaptSfcStopped": 1, + "adaptSfcDetrainment": 0, + "adaptSfcEntrainment": 1, + } + + fieldsInput = fields.copy() + demInput = dem.copy() + demAdapted, fieldsAdapted = com1DFA.adaptDEM(demInput, fieldsInput, cfg["GENERAL"]) + + for key in demAdapted.keys(): + assert np.all(demAdapted[key] == dem[key]) + for key in fieldsAdapted.keys(): + assert np.all(fieldsAdapted[key] == fields[key]) + + cfg = configparser.ConfigParser() + cfg["GENERAL"] = { + "methodMeshNormal": 1, + "adaptSfcStopped": 1, + "adaptSfcDetrainment": 1, + "adaptSfcEntrainment": 1, + } + + fields["FTEnt"] -= 1 + fieldsInput = fields.copy() + demInput = dem.copy() + + demAdapted, fieldsAdapted = com1DFA.adaptDEM(demInput, fieldsInput, cfg["GENERAL"]) + + for key in demAdapted.keys(): + assert np.all(demAdapted[key] == dem[key]) + for key in fieldsAdapted.keys(): + assert np.all(fieldsAdapted[key] == fields[key]) + + fields["FTEnt"] = np.zeros_like(fields["FTDet"]) + fields["FTDet"] = np.array( + [ + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [1, 1, 1, 1, 1], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + ], + dtype=float, + ) + fieldsInput = fields.copy() + demInput = dem.copy() + demAdapted, fieldsAdapted = com1DFA.adaptDEM(demInput, fieldsInput, cfg["GENERAL"]) + + assert np.all( + demAdapted["rasterData"] + == np.array( + [ + [1.0, 2.0, 3.0, 4.0, 5.0], + [1.0, 2.0, 3.0, 4.0, 5.0], + [1.0, 2.0, 3.0, 4.0, 5.0] + 1 / NzNormed[2], + [1.0, 2.0, 3.0, 4.0, 5.0], + [1.0, 2.0, 3.0, 4.0, 5.0], + ] + ) + ) + assert np.any(demAdapted["Nx"] != dem["Nx"]) + assert np.any(demAdapted["Ny"] != dem["Ny"]) + assert np.all(demAdapted["Nz"] == dem["Nz"]) + assert np.any(dem["areaRaster"] != demAdapted["areaRaster"]) + assert np.all(fieldsAdapted["sfcChange"] == fields["FTDet"] / NzNormed) + assert np.all(fieldsAdapted["sfcChangeTotal"] == fields["FTDet"] / NzNormed) diff --git a/docs/com1DFAAlgorithm.rst b/docs/com1DFAAlgorithm.rst index e13ecf1b2..405a7102c 100644 --- a/docs/com1DFAAlgorithm.rst +++ b/docs/com1DFAAlgorithm.rst @@ -166,6 +166,10 @@ Other particles properties are also initialized here: - ``dmDet`` - detrained mass of particle [kg] + - ``dmEnt`` - entrained mass of particle [kg] + + - ``stoppedParticles`` - dictionary with particles (containing x-, y-coordinates, mass and ID) that are stopped (mass or velocity is zero) and deleted from the particles in each time step + For more details, see :py:func:`com1DFA.com1DFA.initializeParticles`. Go back to :ref:`com1DFAAlgorithm:Algorithm graph` diff --git a/docs/moduleCom1DFA.rst b/docs/moduleCom1DFA.rst index 2e95d63e0..e76bedd31 100644 --- a/docs/moduleCom1DFA.rst +++ b/docs/moduleCom1DFA.rst @@ -295,6 +295,8 @@ The result types that can be chosen to be exported are (all correspond to fields * TA - travel angle * dmDet - detrained mass * 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 * particles (:ref:`com1DFAAlgorithm:Particle properties`) Have a look at the designated subsection Output in ``com1DFA/com1DFACfg.ini``.