avaframe/com1DFA/DFAfunctionsCython.pyx
#!python
# cython: boundscheck=False, wraparound=False, cdivision=True
"""
function related to SPH calculations in com1DFA
to build: go to repository containing this file and run:
python setup.py build_ext --inplace
"""
# Load modules
import copy
import logging
import numpy as np
import cython
cimport numpy as np
from libc cimport math as math
# Local imports
import avaframe.com1DFA.DFAtools as DFAtls
cimport avaframe.com1DFA.DFAToolsCython as DFAtlsC
cimport avaframe.com1DFA.damCom1DFA as damCom1DFA
import avaframe.com1DFA.particleTools as particleTools
import avaframe.in3Utils.geoTrans as geoTrans
# create local logger
# change log level in calling module to DEBUG to see log messages
log = logging.getLogger(__name__)
def computeForceC(cfg, particles, fields, dem, int frictType):
""" compute forces acting on the particles (without the SPH component)
Cython implementation implementation
Parameters
----------
cfg: configparser
configuration for DFA simulation
particles : dict
particles dictionary at t
fields: dict
fields dictionary
dem : dict
dictionary with dem information
frictType: int
identifier for friction law to be used
Returns
-------
force : dict
force dictionary
particles : dict
"""
# read input parameters
cdef double enthRef = cfg.getfloat('enthRef')
cdef double muSamosAt = cfg.getfloat('musamosat')
cdef double tau0SamosAt = cfg.getfloat('tau0samosat')
cdef double Rs0SamosAt = cfg.getfloat('Rs0samosat')
cdef double kappaSamosAt = cfg.getfloat('kappasamosat')
cdef double BSamosAt = cfg.getfloat('Bsamosat')
cdef double RSamosAt = cfg.getfloat('Rsamosat')
cdef double muSamosAtSmall = cfg.getfloat('musamosatsmall')
cdef double tau0SamosAtSmall = cfg.getfloat('tau0samosatsmall')
cdef double Rs0SamosAtSmall = cfg.getfloat('Rs0samosatsmall')
cdef double kappaSamosAtSmall = cfg.getfloat('kappasamosatsmall')
cdef double BSamosAtSmall = cfg.getfloat('Bsamosatsmall')
cdef double RSamosAtSmall = cfg.getfloat('Rsamosatsmall')
cdef double muSamosAtMedium = cfg.getfloat('musamosatmedium')
cdef double tau0SamosAtMedium = cfg.getfloat('tau0samosatmedium')
cdef double Rs0SamosAtMedium = cfg.getfloat('Rs0samosatmedium')
cdef double kappaSamosAtMedium = cfg.getfloat('kappasamosatmedium')
cdef double BSamosAtMedium = cfg.getfloat('Bsamosatmedium')
cdef double RSamosAtMedium = cfg.getfloat('Rsamosatmedium')
cdef double entEroEnergy = cfg.getfloat('entEroEnergy')
cdef double entShearResistance = cfg.getfloat('entShearResistance')
cdef double entDefResistance = cfg.getfloat('entDefResistance')
cdef double rho = cfg.getfloat('rho')
cdef double rhoEnt = cfg.getfloat('rhoEnt')
cdef double hRes = cfg.getfloat('hRes')
cdef double gravAcc = cfg.getfloat('gravAcc')
cdef double xsiVoellmy = cfg.getfloat('xsivoellmy')
cdef double muVoellmy = cfg.getfloat('muvoellmy')
cdef double xsiVoellmyMinShear = cfg.getfloat('xsivoellmyminshear')
cdef double muVoellmyMinShear = cfg.getfloat('muvoellmyminshear')
cdef double tau0VoellmyMinShear = cfg.getfloat('tau0voellmyminshear')
cdef double muCoulomb = cfg.getfloat('mucoulomb')
cdef double muCoulombMinShear = cfg.getfloat('mucoulombminshear')
cdef double tau0CoulombMinShear = cfg.getfloat('tau0coulombminshear')
cdef double curvAccInFriction = cfg.getfloat('curvAccInFriction')
cdef double curvAccInTangent = cfg.getfloat('curvAccInTangent')
cdef int curvAccInGradient = cfg.getint('curvAccInGradient')
cdef double velMagMin = cfg.getfloat('velMagMin')
cdef double depMin = cfg.getfloat('depMin')
cdef int interpOption = cfg.getint('interpOption')
cdef int explicitFriction = cfg.getint('explicitFriction')
cdef int reprojMethod = cfg.getint('reprojMethodForce')
cdef int reprojectionIterations = cfg.getint('reprojectionIterations')
cdef double thresholdProjection = cfg.getfloat('thresholdProjection')
cdef double subgridMixingFactor = cfg.getfloat('subgridMixingFactor')
cdef int viscOption = cfg.getint('viscOption')
cdef double dt = particles['dt']
cdef double mu0 = cfg.getfloat('mu0wetsnow')
cdef double xsiWetSnow = cfg.getfloat('xsiwetsnow')
cdef int nPart = particles['nPart']
cdef double csz = dem['header']['cellsize']
cdef int nrows = dem['header']['nrows']
cdef int ncols = dem['header']['ncols']
cdef double[:, :] ZDEM = dem['rasterData']
cdef double[:, :] nxArray = dem['Nx']
cdef double[:, :] nyArray = dem['Ny']
cdef double[:, :] nzArray = dem['Nz']
cdef double[:, :] areaRatser = dem['areaRaster']
cdef np.ndarray[np.uint8_t, ndim = 1, cast=True] outOfDEM
outOfDEM = np.array(dem['outOfDEM'], dtype=bool)
# read particles and fields
cdef double[:] mass = particles['m']
cdef double[:] hArray = particles['h']
cdef double[:] xArray = particles['x']
cdef double[:] yArray = particles['y']
cdef double[:] zArray = particles['z']
cdef double[:] uxArray = particles['ux']
cdef double[:] uyArray = particles['uy']
cdef double[:] uzArray = particles['uz']
cdef long[:] ID = particles['ID']
cdef double[:] totalEnthalpyArray = particles['totalEnthalpy']
cdef double[:, :] VX = fields['Vx']
cdef double[:, :] VY = fields['Vy']
cdef double[:, :] VZ = fields['Vz']
cdef double[:, :] entrMassRaster = fields['entrMassRaster']
cdef double[:, :] entrEnthRaster = fields['entrEnthRaster']
cdef double[:, :] cResRaster = fields['cResRaster']
cdef int[:] indXDEM = particles['indXDEM']
cdef int[:] indYDEM = particles['indYDEM']
# initialize outputs
cdef double[:] forceX = np.zeros(nPart, dtype=np.float64)
cdef double[:] forceY = np.zeros(nPart, dtype=np.float64)
cdef double[:] forceZ = np.zeros(nPart, dtype=np.float64)
cdef double[:] forceFrict = np.zeros(nPart, dtype=np.float64)
cdef double[:] dM = np.zeros(nPart, dtype=np.float64)
cdef double[:] gEff = np.zeros(nPart, dtype=np.float64)
cdef double[:] curvAcc = np.zeros(nPart, dtype=np.float64)
# declare intermediate step variables
cdef int indCellX, indCellY
cdef double areaPart, areaCell, araEntrPart, cResCell, cResPart, uMag, uMagRes, m, dm, h, entrMassCell, entrEnthCell, dEnergyEntr, dis
cdef double x, y, z, xEnd, yEnd, zEnd, ux, uy, uz, uxDir, uyDir, uzDir, totalEnthalpy, enthalpy, dTotalEnthalpy
cdef double nx, ny, nz, nxEnd, nyEnd, nzEnd, nxAvg, nyAvg, nzAvg
cdef double gravAccNorm, accNormCurv, effAccNorm, gravAccTangX, gravAccTangY, gravAccTangZ, forceBotTang, sigmaB, tau
# variables for interpolation
cdef int Lx0, Ly0, LxEnd0, LyEnd0, iCell, iCellEnd
cdef double w[4]
cdef double wEnd[4]
cdef int k
force = {}
# loop on particles
for k in range(nPart):
m = mass[k]
x = xArray[k]
y = yArray[k]
z = zArray[k]
h = hArray[k]
if h < depMin:
h = depMin
ux = uxArray[k]
uy = uyArray[k]
uz = uzArray[k]
indCellX = indXDEM[k]
indCellY = indYDEM[k]
# deduce area
areaPart = m / (h * rho)
# get cell and weights
Lx0, Ly0, iCell, w[0], w[1], w[2], w[3] = DFAtlsC.getCellAndWeights(x, y, ncols, nrows, csz, interpOption)
# get normal at the particle location
nx, ny, nz = DFAtlsC.getVector(Lx0, Ly0, w[0], w[1], w[2], w[3], nxArray, nyArray, nzArray)
nx, ny, nz = DFAtlsC.normalize(nx, ny, nz)
if viscOption == 1:
# add artificial viscosity
ux, uy, uz = addArtificialViscosity(m, h, dt, rho, ux, uy, uz, subgridMixingFactor, Lx0, Ly0,
w[0], w[1], w[2], w[3], VX, VY, VZ, nx, ny, nz)
# get normal at the particle estimated end location
xEnd = x + dt * ux
yEnd = y + dt * uy
zEnd = z + dt * uz
# this point is not necessarily on the surface, project it on the surface
if reprojMethod == 0:
# Project vertically on the dem
iCellEnd = DFAtlsC.getCells(xEnd, yEnd, ncols, nrows, csz)
if iCellEnd >= 0 and outOfDEM[iCellEnd] == 0:
LxEnd0, LyEnd0, iCellEnd, wEnd[0], wEnd[1], wEnd[2], wEnd[3] = DFAtlsC.getCellAndWeights(xEnd, yEnd, ncols, nrows, csz, interpOption)
elif reprojMethod == 1:
# project trying to keep the travelled distance constant
xEnd, yEnd, zEnd, iCellEnd, LxEnd0, LyEnd0, wEnd[0], wEnd[1], wEnd[2], wEnd[3] = DFAtlsC.distConservProjectionIteratrive(
x, y, z, ZDEM, nxArray, nyArray, nzArray, xEnd, yEnd, zEnd, csz, ncols, nrows, interpOption,
reprojectionIterations, thresholdProjection)
elif reprojMethod == 2:
# project using samos method
xEnd, yEnd, iCellEnd, LxEnd0, LyEnd0, wEnd[0], wEnd[1], wEnd[2], wEnd[3] = DFAtlsC.samosProjectionIteratrive(
xEnd, yEnd, zEnd, ZDEM, nxArray, nyArray, nzArray, csz, ncols, nrows, interpOption, reprojectionIterations)
if iCellEnd < 0:
# if not on the DEM take x, y as end point
LxEnd0 = Lx0
LyEnd0 = Ly0
wEnd = w
elif outOfDEM[iCellEnd]:
# if in a noData area from the DEM take x, y as end point
LxEnd0 = Lx0
LyEnd0 = Ly0
wEnd = w
# get the normal at this location
nxEnd, nyEnd, nzEnd = DFAtlsC.getVector(LxEnd0, LyEnd0, wEnd[0], wEnd[1], wEnd[2], wEnd[3], nxArray, nyArray, nzArray)
nxEnd, nyEnd, nzEnd = DFAtlsC.normalize(nxEnd, nyEnd, nzEnd)
# get average of those normals
nxAvg = nx + nxEnd
nyAvg = ny + nyEnd
nzAvg = nz + nzEnd
nxAvg, nyAvg, nzAvg = DFAtlsC.normalize(nxAvg, nyAvg, nzAvg)
# acceleration due to curvature
accNormCurv = (ux*(nxEnd-nx) + uy*(nyEnd-ny) + uz*(nzEnd-nz)) / dt
# normal component of the acceleration of gravity
gravAccNorm = - gravAcc * nzAvg # use nzAvg because we want the average gNorm on the time step
# turn off curvature with the curvAccInFriction coefficient
effAccNorm = gravAccNorm + curvAccInFriction * accNormCurv
# save the normal component of the gravity (plus the curvature term if desiered)
# this is needed to compute the pressure gradient
# save the norm of the gEff, because we know in which direction it goes (minus the normal vector)
if curvAccInGradient == 1:
# use gravity plus curvature (if the effAccNorm > 0, material is detatched, then gEff = 0)
if(effAccNorm <= 0.0):
gEff[k] = -effAccNorm
else:
gEff[k] = 0
else:
# only use gravity
gEff[k] = -gravAccNorm
# body forces (tangential component of acceleration of gravity with the term due to curvature)
gravAccTangX = - (gravAccNorm + curvAccInTangent * accNormCurv) * nxAvg
gravAccTangY = - (gravAccNorm + curvAccInTangent * accNormCurv) * nyAvg
gravAccTangZ = - gravAcc - (gravAccNorm + curvAccInTangent * accNormCurv) * nzAvg
# adding gravity force contribution
forceX[k] = forceX[k] + gravAccTangX * m
forceY[k] = forceY[k] + gravAccTangY * m
forceZ[k] = forceZ[k] + gravAccTangZ * m
# Calculating bottom shear and normal stress
# get new velocity magnitude (leave 0 if uMag is 0)
# this is important because uMag is first used to compute tau
uMag = DFAtlsC.norm(ux, uy, uz)
if(-effAccNorm < 0.0):
# if fluid detatched
# log.info('fluid detatched for particle %s, effAccNorm %.16g' % (k, effAccNorm))
tau = 0.0
else:
# bottom normal stress sigmaB
sigmaB = - effAccNorm * rho * h
if frictType == 1:
# SamosAT friction type (bottom shear stress)
tau = DFAtlsC.SamosATfric(rho, tau0SamosAt, Rs0SamosAt, muSamosAt, kappaSamosAt, BSamosAt, RSamosAt, uMag, sigmaB, h)
elif frictType == 2:
# coulomb friction type (bottom shear stress)
tau = muCoulomb * sigmaB
elif frictType == 3:
# voellmy friction type
tau = muVoellmy * sigmaB + rho * uMag * uMag * gravAcc / xsiVoellmy
elif frictType == 4:
# add enthalpy dependent mu if wetSnow is activated
totalEnthalpy = totalEnthalpyArray[k]
enthalpy = totalEnthalpy - gravAcc * z - 0.5 * uMag * uMag
mu = mu0 * math.exp(-enthalpy / enthRef)
tau = mu * sigmaB + rho * uMag * uMag * gravAcc / xsiWetSnow
elif frictType == 5:
# SamosAT friction type (bottom shear stress) - for small ava calibration parameters
tau = DFAtlsC.SamosATfric(rho, tau0SamosAtSmall, Rs0SamosAtSmall, muSamosAtSmall, kappaSamosAtSmall, BSamosAtSmall, RSamosAtSmall, uMag, sigmaB, h)
elif frictType == 6:
# SamosAT friction type (bottom shear stress) - for medium ava calibration parameters
tau = DFAtlsC.SamosATfric(rho, tau0SamosAtMedium, Rs0SamosAtMedium, muSamosAtMedium, kappaSamosAtMedium, BSamosAtMedium, RSamosAtMedium, uMag, sigmaB, h)
elif frictType == 7:
# voellmy MinShear friction type
tau = muVoellmyMinShear * sigmaB + rho * uMag * uMag * gravAcc / xsiVoellmyMinShear + tau0VoellmyMinShear
elif frictType == 8:
# coulomb MinShear friction type
tau = muCoulombMinShear * sigmaB + tau0CoulombMinShear
else:
tau = 0.0
# adding bottom shear resistance contribution
# norm of the friction force >=0 (if 0 -> detatchment)
forceBotTang = - areaPart * tau
if explicitFriction == 1:
# explicit formulation
forceFrict[k] = forceFrict[k] - forceBotTang
elif explicitFriction == 0:
# make sure uMag is not 0
if uMag<velMagMin:
uMagRes = velMagMin
else:
uMagRes = uMag
forceFrict[k] = forceFrict[k] - forceBotTang/uMagRes
# compute entrained mass and if wetSnow case compute total enthalpy of particles with entrainment
entrMassCell = entrMassRaster[indCellY, indCellX]
if entrMassCell > 0.0:
# compute entrained mass
dm, areaEntrPart = computeEntMassAndForce(dt, entrMassCell, areaPart, uMag, tau, entEroEnergy, rhoEnt)
# speed loss due to energy loss due to entrained mass
dEnergyEntr = areaEntrPart * entShearResistance + dm * entDefResistance
if frictType == 4 and dm > 0.0:
# enthalpy change due to entrained mass
entrEnthCell = entrEnthRaster[indCellY, indCellX]
dTotalEnthalpy = (entrEnthCell + gravAcc * z) * dm
totalEnthalpy = (totalEnthalpy * m + dTotalEnthalpy)
# update velocity
ux = ux * m / (m + dm)
uy = uy * m / (m + dm)
uz = uz * m / (m + dm)
# update mass
m = m + dm
mass[k] = m
dM[k] = dm
# speed loss due to energy loss due to entrained mass
if dEnergyEntr > 0.0:
dis = 1.0 - dEnergyEntr / (0.5 * m * (uMag*uMag + velMagMin))
if dis < 0.0:
dis = 0.0
# update velocity
ux = ux * dis
uy = uy * dis
uz = uz * dis
# if wetSnow update total enthalpy according to new particle mass
if frictType == 4 and dm > 0.0:
# update specific enthalpy of particle
totalEnthalpyArray[k] = totalEnthalpy / m
# adding resistance force due to obstacles
cResCell = cResRaster[indCellY][indCellX]
cResPart = computeResForce(hRes, h, areaPart, rho, cResCell, uMag, explicitFriction)
forceFrict[k] = forceFrict[k] - cResPart
uxArray[k] = ux
uyArray[k] = uy
uzArray[k] = uz
# save results
force['dM'] = np.asarray(dM)
force['forceX'] = np.asarray(forceX)
force['forceY'] = np.asarray(forceY)
force['forceZ'] = np.asarray(forceZ)
force['forceFrict'] = np.asarray(forceFrict)
particles['gEff'] = np.asarray(gEff)
particles['curvAcc'] = np.asarray(curvAcc)
particles['ux'] = np.asarray(uxArray)
particles['uy'] = np.asarray(uyArray)
particles['uz'] = np.asarray(uzArray)
particles['m'] = np.asarray(mass)
particles['totalEnthalpy'] = np.asarray(totalEnthalpyArray)
# update mass available for entrainement
# TODO: this allows to entrain more mass then available...
for k in range(nPart):
indCellX = indXDEM[k]
indCellY = indYDEM[k]
entrMassCell = entrMassRaster[indCellY, indCellX]
areaCell = areaRatser[indCellY, indCellX]
dm = dM[k]
# update surface entrainment mass available
entrMassCell = entrMassCell - dm/areaCell
if entrMassCell < 0:
entrMassCell = 0
entrMassRaster[indCellY, indCellX] = entrMassCell
fields['entrMassRaster'] = np.asarray(entrMassRaster)
return particles, force, fields
cpdef (double, double) computeEntMassAndForce(double dt, double entrMassCell,
double areaPart, double uMag,
double tau, double entEroEnergy,
double rhoEnt):
""" compute force component due to entrained mass
Parameters
----------
dt: float
time step
entrMassCell : float
available mass for entrainement
areaPart : float
particle area
uMag : float
particle speed (velocity magnitude)
tau : float
bottom shear stress
Returns
-------
dm : float
entrained mass
areaEntrPart : float
Area for entrainement energy loss computation
entEroEnergy: float
erosion entrainment energy constant
rhoEnt: float
entrainement density
"""
cdef double width, ABotSwiped, areaEntrPart
# compute entrained mass
cdef double dm = 0
if entrMassCell > 0:
# either erosion or ploughing but not both
if(entEroEnergy > 0):
# erosion: erode according to shear and erosion energy
dm = areaPart * tau * uMag * dt / entEroEnergy
areaEntrPart = areaPart
# TODO why not this?
#areaEntrPart = math.sqrt(areaPart) * uMag * dt
else:
# ploughing in at avalanche front: erode full area weight
# mass available in the cell [kg/m²]
# width of the particle
width = math.sqrt(areaPart)
# bottom area covered by the particle during dt
ABotSwiped = width * uMag * dt
dm = entrMassCell * ABotSwiped
areaEntrPart = entrMassCell / rhoEnt
return dm, areaEntrPart
cpdef double computeResForce(double hRes, double h, double areaPart, double rho,
double cResCell, double uMag, int explicitFriction):
""" compute force component due to resistance
Parameters
----------
hRes: float
resistance height
h : float
particle flow thickness
areaPart : float
particle area
rho : float
snow density
cResCell : float
resisance coefficient of cell
uMag : float
particle speed (velocity magnitude)
explicitFriction: int
if 1 add resistance with an explicit method. Implicit otherwise
Returns
-------
cResPart : float
resistance component for particle
"""
cdef double hResEff = hRes
cdef double cRecResPart
if(h < hRes):
hResEff = h
if explicitFriction == 1:
# explicit formulation
cRecResPart = - rho * areaPart * hResEff * cResCell * uMag * uMag
elif explicitFriction == 0:
cRecResPart = - rho * areaPart * hResEff * cResCell * uMag
return cRecResPart
cdef (double, double, double) addArtificialViscosity(double m, double h, double dt, double rho,
double ux, double uy, double uz, double subgridMixingFactor,
int Lx0, int Ly0, double w0, double w1, double w2, double w3,
double[:, :] VX, double[:, :] VY, double[:, :] VZ,
double nx, double ny, double nz):
""" add artificial viscosity
Add the artificial viscosity in an implicit way and this before adding the other forces.
Parameters
----------
m : float
mass of the particle
h : float
flow thickness of the particle
dt : float
time step
rho : float
density
subgridMixingFactor : float
viscosity coefficient
Lx0: int
colomn of the nearest lower left cell
Ly0: int
row of the nearest lower left cell
w: float[4]
corresponding weights
location in the y location of desiered interpolation
VX: 2D numpy array
x component of the velocity vector field at the grid nodes
VY: 2D numpy array
y component of the velocity vector field at the grid nodes
VZ: 2D numpy array
z component of the velocity vector field at the grid nodes
nx: float
x component of the interpolated vector field at particle position
ny: float
y component of the interpolated vector field at particle position
nz: float
z component of the interpolated vector field at particle position
Returns
-------
ux: float
x component of the velocity uptated with the viscous force
uy: float
y component of the uptated with the viscous force
uz: float
z component of the uptated with the viscous force
"""
cdef double vMeanx, vMeany, vMeanz, vMeanNorm, dvX, dvY, dvZ
vMeanx, vMeany, vMeanz = DFAtlsC.getVector(Lx0, Ly0, w0, w1, w2, w3, VX, VY, VZ)
# compute normal component of the velocity
vMeanNorm = DFAtlsC.scalProd(vMeanx, vMeany, vMeanz, nx, ny, nz)
# remove normal component (make sure vMean is in the tangent plane)
vMeanx = vMeanx - vMeanNorm * nx
vMeany = vMeany - vMeanNorm * ny
vMeanz = vMeanz - vMeanNorm * nz
# compute particle to field velocity difference
dvX = vMeanx - ux
dvY = vMeany - uy
dvZ = vMeanz - uz
dvMag = DFAtlsC.norm(dvX, dvY, dvZ)
Alat = 2.0 * math.sqrt((m * h) / rho)
fDrag = (subgridMixingFactor * 0.5 * rho * dvMag * Alat * dt) / m
# update velocity with artificial viscosity - implicit method
ux = ux + fDrag * vMeanx
uy = uy + fDrag * vMeany
uz = uz + fDrag * vMeanz
ux = ux / (1.0 + fDrag)
uy = uy / (1.0 + fDrag)
uz = uz / (1.0 + fDrag)
return ux, uy, uz
def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0):
""" update particle position using euler forward scheme
Cython implementation
Parameters
----------
cfg: configparser
configuration for DFA simulation
particles : dict
particles dictionary at t
dem : dict
dictionary with dem information
force : dict
force dictionary
fields : dict
fields dictionary with flow thickness (needed if there is a dam)
typeStop: int
0 if standard stopping criterion, if 1 stopping criterion based on SPHforce - used for iniStep
Returns
-------
particles : dict
particles dictionary at t + dt
"""
# read input parameters
cdef double stopCrit = cfg.getfloat('stopCrit')
cdef double stopCritIni = cfg.getfloat('stopCritIni')
cdef double stopCritIniSmall = cfg.getfloat('stopCritIniSmall')
cdef double uFlowingThreshold = cfg.getfloat('uFlowingThreshold')
cdef double dt = particles['dt']
log.debug('dt used now is %f' % dt)
cdef double gravAcc = cfg.getfloat('gravAcc')
cdef double velMagMin = cfg.getfloat('velMagMin')
cdef double rho = cfg.getfloat('rho')
cdef int interpOption = cfg.getint('interpOption')
cdef int explicitFriction = cfg.getint('explicitFriction')
cdef int reprojMethod = cfg.getint('reprojMethodPosition')
cdef int reprojectionIterations = cfg.getint('reprojectionIterations')
cdef double thresholdProjection = cfg.getfloat('thresholdProjection')
cdef double centeredPosition = cfg.getfloat('centeredPosition')
cdef int snowSlide = cfg.getint('snowSlide')
cdef int dissDam = cfg.getint('dissDam')
cdef double csz = dem['header']['cellsize']
cdef int nrows = dem['header']['nrows']
cdef int ncols = dem['header']['ncols']
cdef int nPart = particles['nPart']
cdef double[:, :] nxArray = dem['Nx']
cdef double[:, :] nyArray = dem['Ny']
cdef double[:, :] nzArray = dem['Nz']
cdef double[:, :] ZDEM = dem['rasterData']
cdef np.ndarray[np.uint8_t, ndim = 1, cast=True] outOfDEM
outOfDEM = np.array(dem['outOfDEM'], dtype=bool)
# read particles and fields
cdef double[:] mass = particles['m']
cdef double[:] idFixed = particles['idFixed']
cdef double[:] sArray = particles['trajectoryLengthXY']
cdef double[:] sCorArray = particles['trajectoryLengthXYCor']
cdef double[:] lArray = particles['trajectoryLengthXYZ']
cdef double[:] xArray = particles['x']
cdef double[:] yArray = particles['y']
cdef double[:] zArray = particles['z']
cdef double[:] uxArray = particles['ux']
cdef double[:] uyArray = particles['uy']
cdef double[:] uzArray = particles['uz']
cdef double[:] velocityMagArray = particles['velocityMag']
cdef double[:] uAccArray = particles['uAcc']
cdef double[:] totalEnthalpyArray = particles['totalEnthalpy']
cdef double TotkinEne = particles['kineticEne']
cdef double TotpotEne = particles['potentialEne']
cdef double peakKinEne = particles['peakKinEne']
cdef double peakForceSPH = particles['peakForceSPH']
cdef double totKForceSPH = particles['forceSPHIni']
# read fields
cdef double[:] forceX = force['forceX']
cdef double[:] forceY = force['forceY']
cdef double[:] forceZ = force['forceZ']
cdef double[:] forceFrict = force['forceFrict']
cdef double[:] forceSPHX = force['forceSPHX']
cdef double[:] forceSPHY = force['forceSPHY']
cdef double[:] forceSPHZ = force['forceSPHZ']
cdef double[:] dM = force['dM']
# read dam
dam = dem['damLine']
cdef int flagDam = dam['dam']
cdef int restitutionCoefficient = dam['restitutionCoefficient']
cdef int nIterDam = dam['nIterDam']
cdef int nDamPoints = dam['nPoints']
cdef long[:] cellsCrossed = dam['cellsCrossed']
cdef double[:] xFootArray = dam['x']
cdef double[:] yFootArray = dam['y']
cdef double[:] zFootArray = dam['z']
cdef double[:] xCrownArray = dam['xCrown']
cdef double[:] yCrownArray = dam['yCrown']
cdef double[:] zCrownArray = dam['zCrown']
cdef double[:] xTangentArray = dam['xTangent']
cdef double[:] yTangentArray = dam['yTangent']
cdef double[:] zTangentArray = dam['zTangent']
# read fields
cdef double[:, :] FD = fields['FT']
# initialize outputs
cdef double TotkinEneNew = 0
cdef double TotpotEneNew = 0
cdef double totForceSPHNew = 0
cdef double[:] mNewArray = np.zeros(nPart, dtype=np.float64)
cdef double[:] xNewArray = np.zeros(nPart, dtype=np.float64)
cdef double[:] yNewArray = np.zeros(nPart, dtype=np.float64)
cdef double[:] zNewArray = np.zeros(nPart, dtype=np.float64)
cdef double[:] sNewArray = np.zeros(nPart, dtype=np.float64)
cdef double[:] sCorNewArray = np.zeros(nPart, dtype=np.float64)
cdef double[:] lNewArray = np.zeros(nPart, dtype=np.float64)
cdef double[:] uxArrayNew = np.zeros(nPart, dtype=np.float64)
cdef double[:] uyArrayNew = np.zeros(nPart, dtype=np.float64)
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)
# 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, massFlowing = 0, dissEm = 0
cdef int k, inter
cdef int nRemove = 0
# variables for interpolation
cdef int Lx0, Ly0, LxNew0, LyNew0, iCell, iCellNew
cdef double w[4]
cdef double wNew[4]
# loop on particles
for k in range(nPart):
m = mass[k]
x = xArray[k]
y = yArray[k]
z = zArray[k]
ux = uxArray[k]
uy = uyArray[k]
uz = uzArray[k]
s = sArray[k]
sCor = sCorArray[k]
l = lArray[k]
idfixed = idFixed[k]
# Force magnitude (without friction)
ForceDriveX = forceX[k] + forceSPHX[k]
ForceDriveY = forceY[k] + forceSPHY[k]
ForceDriveZ = forceZ[k] + forceSPHZ[k]
# velocity magnitude
uMag = DFAtlsC.norm(ux, uy, uz)
uMagt0 = DFAtlsC.norm(ux, uy, uz)
# procede to time integration
# operator splitting
# estimate new velocity due to driving force
uxNew = ux + ForceDriveX * dt / m
uyNew = uy + ForceDriveY * dt / m
uzNew = uz + ForceDriveZ * dt / m
# take friction force into account
if typeStop != 1:
uxNew, uyNew, uzNew, dtStop = account4FrictionForce(uxNew, uyNew, uzNew, m, dt, forceFrict[k], uMag, explicitFriction)
else:
dtStop = dt
# update mass (already done un computeForceC)
mNew = m
massEntrained = massEntrained + dM[k]
# update position
if centeredPosition:
xNew = x + dtStop * 0.5 * (ux + uxNew)
yNew = y + dtStop * 0.5 * (uy + uyNew)
zNew = z + dtStop * 0.5 * (uz + uzNew)
else:
xNew = x + dtStop * uxNew
yNew = y + dtStop * uyNew
zNew = z + dtStop * uzNew
# make sure particle is on the mesh (normal reprojection!!)
if reprojMethod == 0:
# Project vertically on the dem
iCellNew = DFAtlsC.getCells(xNew, yNew, ncols, nrows, csz)
if iCellNew >= 0 and outOfDEM[iCellNew] == 0:
Lx0, Ly0, iCellNew, wNew[0], wNew[1], wNew[2], wNew[3] = DFAtlsC.getCellAndWeights(xNew, yNew, ncols, nrows, csz, interpOption)
zNew = DFAtlsC.getScalar(Lx0, Ly0, wNew[0], wNew[1], wNew[2], wNew[3], ZDEM)
elif reprojMethod == 1:
# project trying to keep the travelled distance constant
xNew, yNew, zNew, iCellNew, LxNew0, LyNew0, wNew[0], wNew[1], wNew[2], wNew[3] = DFAtlsC.distConservProjectionIteratrive(
x, y, z, ZDEM, nxArray, nyArray, nzArray, xNew, yNew, zNew, csz, ncols, nrows, interpOption,
reprojectionIterations, thresholdProjection)
elif reprojMethod == 2:
# project using samos method
xNew, yNew, iCellNew, LxNew0, LyNew0, wNew[0], wNew[1], wNew[2], wNew[3] = DFAtlsC.samosProjectionIteratrive(
xNew, yNew, zNew, ZDEM, nxArray, nyArray, nzArray, csz, ncols, nrows, interpOption, reprojectionIterations)
if iCellNew >= 0:
zNew = DFAtlsC.getScalar(LxNew0, LyNew0, wNew[0], wNew[1], wNew[2], wNew[3], ZDEM)
if iCellNew < 0:
# if the particle is not on the DEM, memorize it and remove it at the next update
keepParticle[k] = 0
LxNew0 = 0
LyNew0 = 0
wNew = [0, 0, 0, 0]
nRemove = nRemove + 1
continue # this particle will be removed, skip the what is bellow and directly go to the next particle
elif outOfDEM[iCellNew]:
# if the particle is on the DEM but in a noData area,
# memorize it and remove it at the next update
keepParticle[k] = 0
nRemove = nRemove + 1
continue # this particle will be removed, skipp to the next particle
# get cell and weights at old position
Lx0, Ly0, iCell, w[0], w[1], w[2], w[3] = DFAtlsC.getCellAndWeights(x, y, ncols, nrows, csz, interpOption)
# get normal at the old particle location
nx, ny, nz = DFAtlsC.getVector(Lx0, Ly0, w[0], w[1], w[2], w[3], nxArray, nyArray, nzArray)
# get normal at the new particle location
nxNew, nyNew, nzNew = DFAtlsC.getVector(LxNew0, LyNew0, wNew[0], wNew[1], wNew[2], wNew[3], nxArray, nyArray, nzArray)
# get average normal between old and new position
nx, ny, nz = DFAtlsC.normalize(nx+nxNew, ny+nyNew, nz+nzNew)
# normalize new normal
nxNew, nyNew, nzNew = DFAtlsC.normalize(nxNew, nyNew, nzNew)
# compute the distance traveled by the particle
dl = DFAtlsC.norm((xNew-x), (yNew-y), (zNew-z))
lNew = l + dl
# compute the horizontal distance traveled by the particle
ds = DFAtlsC.norm((xNew-x), (yNew-y), 0)
sNew = s + ds
# compute the horizontal distance traveled by the particle (corrected with
# the angle difference between the slope and the normal)
sCorNew = sCor + nz*dl
# reproject velocity
uxNew, uyNew, uzNew, uMag = DFAtlsC.reprojectVelocity(uxNew, uyNew, uzNew, nxNew, nyNew, nzNew, velMagMin, 1)
# ############## Start add Dam interaction ##############################################
if flagDam and cellsCrossed[iCell] == 1:
# if there is an interaction with the dam, update position and velocity
inter, xNew, yNew, zNew, uxNew, uyNew, uzNew, txWall, tyWall, tzWall, dissEm = damCom1DFA.getWallInteraction(x, y, z,
xNew, yNew, zNew, uxNew, uyNew, uzNew, nDamPoints, xFootArray, yFootArray, zFootArray,
xCrownArray, yCrownArray, zCrownArray, xTangentArray, yTangentArray, zTangentArray,
ncols, nrows, csz, interpOption, restitutionCoefficient, nIterDam, nxArray, nyArray, nzArray, ZDEM, FD)
# if there was an interaction with the dam, reproject and take dEM into account
if inter == 1:
LxNew0, LyNew0, iCellNew, wNew[0], wNew[1], wNew[2], wNew[3] = DFAtlsC.getCellAndWeights(xNew, yNew, ncols, nrows, csz, interpOption)
if iCellNew >= 0:
zNew = DFAtlsC.getScalar(LxNew0, LyNew0, wNew[0], wNew[1], wNew[2], wNew[3], ZDEM)
# do we need to remove the particle?
if iCellNew < 0:
# if the particle is not on the DEM, memorize it and remove it at the next update
keepParticle[k] = 0
LxNew0 = 0
LyNew0 = 0
wNew = [0, 0, 0, 0]
nRemove = nRemove + 1
continue # this particle will be removed, skipp to the next particle
elif outOfDEM[iCellNew]:
# if the particle is on the DEM but in a noData area,
# memorize it and remove it at the next update
keepParticle[k] = 0
nRemove = nRemove + 1
continue # this particle will be removed, skipp to the next particle
# reduce velocity normal to footline for energy loss due to flowing over dam
# dEm>0 means the snow thickness exceeds the dam height???
# reduce velocity normal to footline for energy loss due to flowing over dam
# (but there is maybe nothing flowing over...)
if dissDam == 1 and dissEm < 0.0:
# First multiply dissEm by the gravAcc magnitude (because this was not done in the dam function)
dissEm = gravAcc*dissEm
fNx, fNy, fNz = DFAtlsC.crossProd(nxNew, nyNew, nzNew, txWall, tyWall, tzWall)
fNx, fNy, fNz = DFAtlsC.normalize(fNx, fNy, fNz)
dv = math.sqrt(-2.0 * dissEm)
uN = -DFAtlsC.scalProd(uxNew, uyNew, uzNew, fNx, fNy, fNz)
if uN < dv:
dv = uN
if dv > 0.0:
uxNew = uxNew + dv * fNx
uyNew = uyNew + dv * fNy
uzNew = uzNew + dv * fNz
# reproject velocity
uxNew, uyNew, uzNew, uMag = DFAtlsC.reprojectVelocity(uxNew, uyNew, uzNew, nxNew, nyNew, nzNew, velMagMin, 1)
# ############## End add Dam interaction ##############################################
# prepare for stopping criterion
if uMag > uFlowingThreshold:
# if velocity is bigger then threshold add to flowing mass
massFlowing = massFlowing + mNew
TotkinEneNew = TotkinEneNew + 0.5 * m * uMag * uMag
TotpotEneNew = TotpotEneNew + mNew * gravAcc * zNew
totForceSPHNew = totForceSPHNew + mNew * DFAtlsC.norm(forceSPHX[k], forceSPHY[k], forceSPHZ[k])
if idfixed == 1:
# idfixed = 1 if particles belong to 'fixed' boundary particles - so zero velocity and fixed position
xNewArray[k] = x
yNewArray[k] = y
zNewArray[k] = z
uxArrayNew[k] = ux
uyArrayNew[k] = uy
uzArrayNew[k] = uz
lNewArray[k] = l
sNewArray[k] = s
sCorNewArray[k] = sCor
mNewArray[k] = m
else:
# idfixed = 0 particles belong to the actual releae area
xNewArray[k] = xNew
yNewArray[k] = yNew
zNewArray[k] = zNew
uxArrayNew[k] = uxNew
uyArrayNew[k] = uyNew
uzArrayNew[k] = uzNew
lNewArray[k] = lNew
sNewArray[k] = sNew
sCorNewArray[k] = sCorNew
mNewArray[k] = mNew
# compute acceleration
uMagt1 = DFAtlsC.norm(uxNew, uyNew, uzNew)
uAcc = (uMagt1 - uMagt0) / dtStop
uAccArray[k] = uAcc
velocityMagArrayNew[k] = uMagt1
particles['ux'] = np.asarray(uxArrayNew)
particles['uy'] = np.asarray(uyArrayNew)
particles['uz'] = np.asarray(uzArrayNew)
particles['velocityMag'] = np.asarray(velocityMagArrayNew)
particles['uAcc'] = np.asarray(uAccArray)
particles['trajectoryLengthXYZ'] = np.asarray(lNewArray)
particles['trajectoryLengthXY'] = np.asarray(sNewArray)
particles['trajectoryLengthXYCor'] = np.asarray(sCorNewArray)
particles['m'] = np.asarray(mNewArray)
particles['mTot'] = np.sum(particles['m'])
particles['x'] = np.asarray(xNewArray)
particles['y'] = np.asarray(yNewArray)
particles['z'] = np.asarray(zNewArray)
particles['massEntrained'] = massEntrained
log.debug('Entrained DFA mass: %s kg', np.asarray(massEntrained))
particles['kineticEne'] = TotkinEneNew
particles['potentialEne'] = TotpotEneNew
if typeStop == 1:
# typeStop = 1 for initialisation step where particles are redistributed to reduce SPH force
# here stop criterion based on SPHForce
value = totForceSPHNew
peakValue = particles['peakForceSPH']
oldValue = particles['forceSPHIni']
particles['forceSPHIni'] = totForceSPHNew
if oldValue == 0.0:
stop = False
elif value < 1.:
stop = oldValue/value < stopCritIniSmall
else:
stop = value <= stopCritIni*peakValue
log.debug('SPHFORCE value %f and stop value %f' % (totForceSPHNew, stopCritIni*peakValue))
if peakForceSPH < totForceSPHNew:
particles['peakForceSPH'] = totForceSPHNew
else:
# avalanche computation stop criterion based on kinetic energy or massFlowing
if peakKinEne < TotkinEneNew:
particles['peakKinEne'] = TotkinEneNew
if particles['peakMassFlowing'] < massFlowing:
particles['peakMassFlowing'] = massFlowing
if cfg['stopCritType'] == 'massFlowing':
value = massFlowing
peakValue = particles['peakMassFlowing']
stop = (value <= stopCrit*peakValue) and peakValue > 0
elif cfg['stopCritType'] == 'kinEnergy':
value = TotkinEneNew
peakValue = particles['peakKinEne']
stop = value <= stopCrit*peakValue
else:
message = 'stopCritType %s is not defined' % cfg['stopCritType']
log.error(message)
raise AssertionError(message)
if stop:
particles['iterate'] = False
if typeStop == 1:
log.debug('stopping initial particle distribution')
else:
log.debug('stopping because of %s stopCriterion.' % (cfg['stopCritType']))
# 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)
return particles
cpdef (double, double, double, double) account4FrictionForce(double ux, double uy,
double uz, double m,
double dt, double forceFrict,
double uMag, int explicitFriction):
""" update velocity with friction force
Parameters
----------
ux: float
x velocity
uy: float
y velocity
uz: float
z velocity
m : float
particle area
dt : float
time step
forceFrict : float
friction force (actually a force for the explicit method, force/vel for implicit)
uMag : float
particle speed (velocity magnitude)
explicitFriction: int
if 1 add resistance with an explicit method. Implicit otherwise
Returns
-------
uxNew: float
x velocity
uyNew: float
y velocity
uzNew: float
z velocity
dtStop : float
time step (smaller then dt if the particle stops during this time step)
"""
cdef double xDir, yDir, zDir, uxNew, uyNew, uzNew, uMagNew, dtStop
if explicitFriction == 1:
# explicite method
uMagNew = DFAtlsC.norm(ux, uy, uz)
# will friction force stop the particle
if uMagNew<dt*forceFrict/m:
# stop the particle
uxNew = 0
uyNew = 0
uzNew = 0
# particle stops after
if uMag<=0:
dtStop = 0
else:
dtStop = m * uMagNew / (forceFrict)
else:
# add friction force in the opposite direction of the motion
xDir, yDir, zDir = DFAtlsC.normalize(ux, uy, uz)
uxNew = ux - xDir * forceFrict * dt / m
uyNew = uy - yDir * forceFrict * dt / m
uzNew = uz - zDir * forceFrict * dt / m
dtStop = dt
elif explicitFriction == 0:
# implicite method
uxNew = ux / (1.0 + dt * forceFrict / m)
uyNew = uy / (1.0 + dt * forceFrict / m)
uzNew = uz / (1.0 + dt * forceFrict / m)
dtStop = dt
return uxNew, uyNew, uzNew, dtStop
def updateFieldsC(cfg, particles, dem, fields):
""" update fields and particles flow thickness
Cython implementation
Parameters
----------
cfg: configparser
configuration for DFA simulation
particles : dict
particles dictionary
dem : dict
dictionary with dem information
fields : dict
fields dictionary
Returns
-------
particles : dict
particles dictionary
fields : dict
fields dictionary
"""
# read input parameters
cdef double rho = cfg.getfloat('rho')
cdef int interpOption = cfg.getint('interpOption')
header = dem['header']
cdef int nrows = header['nrows']
cdef int ncols = header['ncols']
cdef double xllc = 0
cdef double yllc = 0
cdef double csz = header['cellsize']
cdef int nPart = np.size(particles['x'])
cdef double[:, :] areaRaster = dem['areaRaster']
# read particles and fields
cdef double[:] mass = particles['m']
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 bint computeTA = fields['computeTA']
cdef bint computeKE = fields['computeKE']
cdef bint computeP = fields['computeP']
cdef double[:, :] PFV = fields['pfv']
cdef double[:, :] PP = fields['ppr']
cdef double[:, :] PFT = fields['pft']
cdef double[:, :] PTA = fields['pta']
cdef double[:, :] PKE = fields['pke']
# initialize outputs
cdef double[:, :] MassBilinear = 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))
cdef double[:, :] MomBilinearX = np.zeros((nrows, ncols))
cdef double[:, :] MomBilinearY = np.zeros((nrows, ncols))
cdef double[:, :] MomBilinearZ = np.zeros((nrows, ncols))
cdef double[:, :] VXBilinear = np.zeros((nrows, ncols))
cdef double[:, :] VYBilinear = np.zeros((nrows, ncols))
cdef double[:, :] VZBilinear = np.zeros((nrows, ncols))
cdef double[:, :] kineticEnergy = np.zeros((nrows, ncols))
cdef double[:, :] travelAngleField = np.zeros((nrows, ncols))
# declare intermediate step variables
cdef double[:] hBB = np.zeros((nPart))
cdef double m, h, x, y, z, s, ux, uy, uz, nx, ny, nz, hbb, hLim, areaPart, trajectoryAngle
cdef int k, i
cdef int indx, indy
cdef int ind1[4]
cdef int ind2[4]
ind1[:] = [0, 1, 0, 1]
ind2[:] = [0, 0, 1, 1]
# variables for interpolation
cdef int Lx0, Ly0, iCell
cdef double w[4]
cdef double mwi
for k in range(nPart):
x = xArray[k]
y = yArray[k]
ux = uxArray[k]
uy = uyArray[k]
uz = uzArray[k]
m = mass[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
Lx0, Ly0, iCell, w[0], w[1], w[2], w[3] = DFAtlsC.getCellAndWeights(x, y, ncols, nrows, csz, interpOption)
# for the travel angle we simply do a nearest interpolation
indx = <int>math.round(x / csz)
indy = <int>math.round(y / csz)
if computeTA:
trajectoryAngle = trajectoryAngleArray[k]
travelAngleField[indy, indx] = max(travelAngleField[indy, indx], trajectoryAngle)
# add the component of the points value to the 4 neighbour grid points
# TODO : check if giving the arrays [0 1 0 1].. is faster
for i in range(4):
indx = Lx0 + ind1[i]
indy = Ly0 + ind2[i]
mwi = m * w[i]
MassBilinear[indy, indx] = MassBilinear[indy, indx] + mwi
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 i in range(ncols):
for j in range(nrows):
m = MassBilinear[j, i]
if m > 0:
# TODO: here we devide by the area of the vertex, would it not make
# more sense to devide by the area of the cell in the previous loop?
FTBilinear[j, i] = m / (areaRaster[j, i] * rho)
VXBilinear[j, i] = MomBilinearX[j, i]/m
VYBilinear[j, i] = MomBilinearY[j, i]/m
VZBilinear[j, i] = MomBilinearZ[j, i]/m
VBilinear[j, i] = DFAtlsC.norm(VXBilinear[j, i], VYBilinear[j, i], VZBilinear[j, i])
if VBilinear[j, i] > PFV[j, i]:
PFV[j, i] = VBilinear[j, i]
if FTBilinear[j, i] > PFT[j, i]:
PFT[j, i] = FTBilinear[j, i]
if computeP:
PBilinear[j, i] = computePressure(VBilinear[j, i], rho)
if PBilinear[j, i] > PP[j, i]:
PP[j, i] = PBilinear[j, i]
if computeTA:
if travelAngleField[j, i] > PTA[j, i]:
PTA[j, i] = travelAngleField[j, i]
if computeKE:
# in J/cell (this is not normalized yet and depends on the cell size used for the computation)
kineticEnergy[j, i] = 0.5*m*VBilinear[j, i]*VBilinear[j, i]
if kineticEnergy[j, i] > PKE[j, i]:
PKE[j, i] = kineticEnergy[j, i]
fields['FM'] = np.asarray(MassBilinear)
fields['FV'] = np.asarray(VBilinear)
fields['Vx'] = np.asarray(VXBilinear)
fields['Vy'] = np.asarray(VYBilinear)
fields['Vz'] = np.asarray(VZBilinear)
fields['FT'] = np.asarray(FTBilinear)
fields['pfv'] = np.asarray(PFV)
fields['pft'] = np.asarray(PFT)
if computeP:
fields['ppr'] = np.asarray(PP)
fields['P'] = np.asarray(PBilinear)
if computeTA:
fields['TA'] = np.asarray(travelAngleField)
fields['pta'] = np.asarray(PTA)
if computeKE:
fields['pke'] = np.asarray(PKE)
for k in range(nPart):
x = xArray[k]
y = yArray[k]
Lx0, Ly0, iCell, w[0], w[1], w[2], w[3] = DFAtlsC.getCellAndWeights(x, y, ncols, nrows, csz, interpOption)
hbb = DFAtlsC.getScalar(Lx0, Ly0, w[0], w[1], w[2], w[3], FTBilinear)
hBB[k] = hbb
particles['h'] = np.asarray(hBB)
return particles, fields
cpdef double computePressure(double v, double rho):
"""Compute pressure using the p = rho*v² equation
Parameters
----------
v : float
velocity magnitude
rho : float
snow density
Returns
-------
p : float
pressure computed using p = rho*v²
"""
cdef double p
p = rho * v * v
return p
def computeTrajectoryAngleC(particles, zPartArray0):
"""Compute the travel angle associated to the particles
Parameters
----------
particles : dict
particles dictionary at t
zPartArray0 : dict
z coordinate of particles at t=0s
Returns
-------
particles : dict
particles dictionary updated with the travel angle
"""
cdef int[:] parentIDArray = particles['parentID'].astype('intc')
cdef int nPart = particles['nPart']
cdef double[:] zArray = particles['z']
cdef double[:] sArray = particles['trajectoryLengthXY']
cdef double[:] z0Array = zPartArray0
cdef double[:] gammaArray = np.zeros(nPart)
cdef int parentID, j
cdef double tanGamma, gamma, s, z, z0
# get particle location
# first compute travel angle for each particle
for k in range(nPart):
# get parent Id in order to get the first z position
parentID = parentIDArray[k]
# get z0
z0 = z0Array[parentID]
# compute tan of the travel angle
z = zArray[k]
s = sArray[k]
if s > 0:
tanGamma = (z0 - z) / s
else:
tanGamma = 0.0
# get the travel angle
gamma = math.atan(tanGamma) * 180.0 / math.pi
gammaArray[k] = gamma
particles['trajectoryAngle'] = np.asarray(gammaArray)
return particles
def getNeighborsC(particles, dem):
""" Locate particles on DEM and neighbour search grid
Ĺocate each particle in a grid cell (both DEM and neighbour search grid) and build the
indPartInCell and partInCell arrays for SPH neighbourSearch and
InCell, IndX and IndY arrays location on the DEM grid.
See issue #200 and documentation for details
Parameters
----------
particles : dict
dem : dict
dem dict with neighbour search grid header (information about neighbour search grid)
Returns
-------
particles : dict
updated particles dictionary with
neighbours info indPartInCell and partInCell arrays for neighbour search and
with InCell, IndX and IndY arrays location on the DEM grid
"""
# get DEM grid information
header = dem['header']
cdef int nColsDEM = header['ncols']
cdef int nRowsDEM = header['nrows']
cdef float cszDEM = header['cellsize']
# get neighbour search grid information
headerNeighbourGrid = dem['headerNeighbourGrid']
cdef int nColsNeighbourGrid = headerNeighbourGrid['ncols']
cdef int nRowsNeighbourGrid = headerNeighbourGrid['nrows']
cdef float cszNeighbourGrid = headerNeighbourGrid['cellsize']
# get particle location
cdef int nPart = particles['nPart']
cdef int k
cdef double[:] xArray = particles['x']
cdef double[:] yArray = particles['y']
# initialize outputs
cdef int nCellsNeighbourGrid = nColsNeighbourGrid*nRowsNeighbourGrid
cdef int[:] indPartInCell = np.zeros(nCellsNeighbourGrid + 1).astype('intc')
cdef int[:] indPartInCell2 = np.zeros(nCellsNeighbourGrid + 1).astype('intc')
cdef int[:] partInCell = np.zeros(nPart).astype('intc')
cdef int[:] indXDEM = np.zeros(nPart).astype('intc')
cdef int[:] indYDEM = np.zeros(nPart).astype('intc')
cdef int[:] inCellDEM = np.zeros(nPart).astype('intc')
# Count number of particles in each SPH grid cell
cdef int indx, indy, ic
for k in range(nPart):
indx = <int>math.round(xArray[k] / cszNeighbourGrid)
indy = <int>math.round(yArray[k] / cszNeighbourGrid)
# get index of cell containing the particle
ic = indx + nColsNeighbourGrid * indy
indPartInCell[ic+1] = indPartInCell[ic+1] + 1
for ic in range(nCellsNeighbourGrid):
indPartInCell[ic+1] = indPartInCell[ic] + indPartInCell[ic+1]
indPartInCell2[ic+1] = indPartInCell[ic+1]
# make the list of which particles are in which cell
for k in range(nPart):
indx = <int>math.round(xArray[k] / cszNeighbourGrid)
indy = <int>math.round(yArray[k] / cszNeighbourGrid)
ic = indx + nColsNeighbourGrid * indy
partInCell[indPartInCell2[ic+1]-1] = k
indPartInCell2[ic+1] = indPartInCell2[ic+1] - 1
indXDEM[k] = <int>math.round(xArray[k] / cszDEM)
indYDEM[k] = <int>math.round(yArray[k] / cszDEM)
# get index of cell containing the particle
inCellDEM[k] = indXDEM[k] + nColsDEM * indYDEM[k]
particles['inCellDEM'] = np.asarray(inCellDEM)
particles['indXDEM'] = np.asarray(indXDEM)
particles['indYDEM'] = np.asarray(indYDEM)
particles['indPartInCell'] = np.asarray(indPartInCell)
particles['partInCell'] = np.asarray(partInCell)
return particles
def computeCohesionForceC(cfg, particles, force):
""" compute elastic cohesion forces acting on the particles
this is computed when the snow slide option is activated (snowSlide = 1
)
Parameters
----------
cfg: configparser
configuration for DFA simulation
particles : dict
particles dictionary at time t
force : dict
force dictionary
Returns
-------
particles : dict
particles dictionary at t (updated bondDist due to breaking)
force : dict
force dictionary with bonding elastic force
"""
# read input parameters
cdef double dt = particles['dt']
cdef double cohesiveSurfaceTension = cfg.getfloat('cohesiveSurfaceTension')
cdef double cohesionMaxStrain = cfg.getfloat('cohesionMaxStrain')
cdef double minDistCohesion = cfg.getfloat('minDistCohesion')
cdef int nPart = particles['nPart']
# read particles and fields
cdef double[:] mass = particles['m']
cdef double[:] hArray = particles['h']
cdef double[:] xArray = particles['x']
cdef double[:] yArray = particles['y']
cdef double[:] zArray = particles['z']
cdef double[:] uxArray = particles['ux']
cdef double[:] uyArray = particles['uy']
cdef double[:] uzArray = particles['uz']
cdef int[:] bondStart = particles['bondStart'].astype('intc')
cdef double[:] bondDist = particles['bondDist']
cdef int[:] bondPart = particles['bondPart'].astype('intc')
cdef double[:] forceSPHX = force['forceSPHX']
cdef double[:] forceSPHY = force['forceSPHY']
cdef double[:] forceSPHZ = force['forceSPHZ']
# initialize outputs
cdef int k, l, ib
cdef double dist0, dist, vx, vy, vz, xk, yk, zk, xl, yl, zl
cdef double hk, hl, ux, uy, uz, dx, dy, dz
cdef double epsilon, mkl, Akl, contact, forceCohesion
# loop on particles
for k in range(nPart):
hk = hArray[k]
xk = xArray[k]
yk = yArray[k]
zk = zArray[k]
ux = uxArray[k]
uy = uyArray[k]
uz = uzArray[k]
# loop on all bonded particles
for ib in range(bondStart[k], bondStart[k + 1]):
# get initial distance between particles
dist0 = bondDist[ib]
# if bond already broken, just skip it
if (dist0 < 0.0):
continue
# get the bonded particle
l = bondPart[ib]
xl = xArray[l]
yl = yArray[l]
zl = zArray[l]
vx = uxArray[l]
vy = uyArray[l]
vz = uzArray[l]
dx = xl - xk + dt * (vx - ux)
dy = yl - yk + dt * (vy - uy)
dz = zl - zk + dt * (vz - uz)
dist = DFAtlsC.norm(dx, dy, dz)
# only if the distance is bigger than a threshold (avoid dividing by 0)
if (dist > minDistCohesion):
# strain
epsilon = dist / dist0 - 1.0
# allow breaking in both compression and extension
if (abs(epsilon) > cohesionMaxStrain):
# break bond
bondDist[ib] = -1
# if not broken add elastic force
else:
hl = hArray[l]
# compute average depth and mass
hkl = 0.5 * (hk + hl)
# sqrt(3) from 6-neighbors-assumption (Iwould add a 2 here... 2/sqrt(3))
Akl = hkl * dist / math.sqrt(3)
# cohesiveSurfaceTension used as elasticity modulus (Pa)
contact = Akl * cohesiveSurfaceTension
# dx / dist is the unit direction vector towards bonded particle l
# cohesion force (including the normalizing factor 1/dist)
forcecohesion = (1 / dist) * epsilon * contact
forceSPHX[k] = forceSPHX[k] + forcecohesion * dx
forceSPHY[k] = forceSPHY[k] + forcecohesion * dy
forceSPHZ[k] = forceSPHZ[k] + forcecohesion * dz
force['forceSPHX'] = np.asarray(forceSPHX)
force['forceSPHY'] = np.asarray(forceSPHY)
force['forceSPHZ'] = np.asarray(forceSPHZ)
particles['bondDist'] = np.asarray(bondDist)
return force, particles
def plotBondC(particles):
""" update edges for plot (bonds that still exist)
Cython implementation
Parameters
----------
particles : dict
particles dictionary
Returns
-------
edges : 2D numpy array
bonds that still exist
"""
# read input parameters
cdef int nPart = particles['nPart']
# read particles and fields
cdef int[:] bondStart = particles['bondStart'].astype('intc')
cdef double[:] bondDist = particles['bondDist']
cdef int[:] bondPart = particles['bondPart'].astype('intc')
cdef int nEdges = bondPart.shape[0]
cdef int[:, :] edges = np.zeros((nEdges, 2), dtype=np.int32)
# initialize outputs
cdef int k, l, ib, count = 0
cdef double dist0
# loop on particles
for k in range(nPart):
# loop on all bonded particles
for ib in range(bondStart[k], bondStart[k + 1]):
# get initial distance between particles
dist0 = bondDist[ib]
# if bond already broken, just skip it
if (dist0 > 0.0):
l = bondPart[ib]
edges[count, 0] = k
edges[count, 1] = l
count = count + 1
return edges
def initializeBondsC(particles, triangles):
""" Initialize bonds if snowSlide is activated
Parameters
----------
particles : dict
particles dictionary at t start
triangles : matplotlib object
result from tri.Triangulation(x, y)
Returns
-------
particles : dict
particles dictionary updated with the bonds
"""
# get all inputs
# particle information
cdef double[:] xArray = particles['x']
cdef double[:] yArray = particles['y']
cdef double[:] zArray = particles['z']
cdef int[:, :] tri = triangles.triangles
cdef int[:, :] edges = triangles.edges
cdef int nPart = xArray.shape[0]
cdef int nTri = tri.shape[0]
cdef int nEdges = edges.shape[0]
# initialize variables and outputs
cdef int[:] bondStart = np.zeros(nPart + 1, dtype=np.int32)
cdef int[:] bondStart2 = np.zeros(nPart + 1, dtype=np.int32)
cdef int[:] bondPart = np.zeros(2 * nEdges, dtype=np.int32)
cdef double[:] bondDist = np.zeros(2 * nEdges, dtype=np.float64)
cdef double dx, dy, dz
cdef int k, e, p1, p2
cdef double distanceIni
# count bonds for each particle
for e in range(nEdges):
p1 = edges[e, 0]
p2 = edges[e, 1]
bondStart[p1+1] = bondStart[p1+1] + 1
bondStart[p2+1] = bondStart[p2+1] + 1
# create cumulative sum
for k in range(nPart):
bondStart[k+1] = bondStart[k] + bondStart[k+1]
bondStart2[k+1] = bondStart[k+1]
# make the list of which particles are bonded with which particle
for e in range(nEdges):
p1 = edges[e, 0]
p2 = edges[e, 1]
dx = xArray[p1] - xArray[p2]
dy = yArray[p1] - yArray[p2]
dz = zArray[p1] - zArray[p2]
distanceIni = DFAtlsC.norm(dx, dy, dz)
bondPart[bondStart2[p1+1]-1] = p2
bondDist[bondStart2[p1+1]-1] = distanceIni
bondStart2[p1+1] = bondStart2[p1+1] - 1
bondPart[bondStart2[p2+1]-1] = p1
bondDist[bondStart2[p2+1]-1] = distanceIni
bondStart2[p2+1] = bondStart2[p2+1] - 1
particles['bondStart'] = np.asarray(bondStart)
particles['bondPart'] = np.asarray(bondPart)
particles['bondDist'] = np.asarray(bondDist)
return particles
def countRemovedBonds(particles, mask, nRemove):
""" count number of bonds to remove
Parameters
----------
particles : dict
particles dictionary at t start
mask : 1D numpy array
particles to keep
nRemove : int
number of particles removed
Returns
-------
nBondRemove : int
number of bonds removed
"""
# read input parameters
cdef int nPart = particles['nPart']
cdef int[:] bondStart = particles['bondStart'].astype('intc')
# read particles and fields
cdef int[:] keepParticle = mask.astype('intc')
cdef int nBondRemove = 0
cdef int k, ib
# count bonds for each particle
for k in range(nPart):
if keepParticle[k] == 0:
# loop on all bonded particles
for ib in range(bondStart[k], bondStart[k + 1]):
nBondRemove = nBondRemove + 2
return nBondRemove
def removedBonds(particles, mask, nRemove, nBondRemove):
""" Remove bonds of removed particles
Parameters
----------
particles : dict
particles dictionary at t start
mask : 1D numpy array
particles to keep
nRemove : int
number of particles removed
Returns
-------
particles : dict
particles dictionary according to removed particles
"""
# read input parameters
cdef int nPart = particles['nPart']
# read particles and fields
cdef int[:] bondStart = particles['bondStart'].astype('intc')
cdef double[:] bondDist = particles['bondDist']
cdef int[:] bondPart = particles['bondPart'].astype('intc')
cdef int[:] keepParticle = mask.astype('intc')
cdef int nEdges = bondPart.shape[0]
cdef int nPartRemoved = nRemove
cdef int[:] bondStartNew = np.zeros(nPart + 1, dtype=np.int32)
cdef int[:] bondPartNew = np.zeros(nEdges - nBondRemove, dtype=np.int32)
cdef double[:] bondDistNew = np.zeros(nEdges - nBondRemove, dtype=np.float64)
cdef int countBond = 0
cdef int countBondNew = 0
cdef int k, ib, l, nBond
# countBond bonds for each particle
for k in range(nPart):
# count the bonds of the particle
nBond = 0
if keepParticle[k] == 1:
# loop on all bonded particles to k
for ib in range(bondStart[k], bondStart[k + 1]):
l = bondPart[ib]
if keepParticle[l] == 1:
# we keep the particle l so add it to the arrays
bondPartNew[countBondNew] = l
bondDistNew[countBondNew] = bondDist[ib]
countBondNew = countBondNew + 1
countBond = countBond + 1
nBond = nBond + 1
else:
countBond = countBond + 1
bondStartNew[k + 1] = bondStartNew[k] + nBond
else:
# we do not keep the particle
bondStartNew[k + 1] = bondStartNew[k]
# loop on all bonded particles
for ib in range(bondStart[k], bondStart[k + 1]):
# increment the counter without doing any thing else
countBond = countBond + 1
particles['bondStart'] = np.asarray(bondStartNew)
particles['bondPart'] = np.asarray(bondPartNew)
particles['bondDist'] = np.asarray(bondDistNew)
return particles
def computeForceSPHC(cfg, particles, force, dem, int sphOption, gradient=0):
""" Prepare data for C computation of lateral forces (SPH component)
acting on the particles (SPH component)
Parameters
----------
cfg: configparser
configuration for DFA simulation
particles : dict
particles dictionary at t
force : dict
force dictionary
dem : dict
dictionary with dem information
sphOption: int
which sphOption should be use
Returns
-------
particles : dict
particles dictionary at t
force : dict
force dictionary
"""
# get SPH grid information
headerNeighbourGrid = dem['headerNeighbourGrid']
# get DEM header and normal field
headerNormalGrid = dem['header']
nxArray = dem['Nx']
nyArray = dem['Ny']
nzArray = dem['Nz']
forceSPHX, forceSPHY, forceSPHZ = computeGradC(cfg, particles, headerNeighbourGrid, headerNormalGrid, nxArray, nyArray, nzArray, gradient, sphOption)
forceSPHX = np.asarray(forceSPHX)
forceSPHY = np.asarray(forceSPHY)
forceSPHZ = np.asarray(forceSPHZ)
force['forceSPHX'] = forceSPHX
force['forceSPHY'] = forceSPHY
force['forceSPHZ'] = forceSPHZ
return particles, force
def computeGradC(cfg, particles, headerNeighbourGrid, headerNormalGrid, double[:, :] nxArray, double[:, :] nyArray,
double[:, :] nzArray, gradient, int SPHoption):
""" compute lateral forces acting on the particles (SPH component)
Cython implementation
Parameters
----------
cfg: configparser
configuration for DFA simulation
particles : dict
particles dictionary at t
headerNeighbourGrid : dict
neighbour search header dictionary (information about SPH grid)
headerNormalGrid : double
normal grid header dictionary (information about the DEM grid)
nxArray : 2D numpy array
x component of the normal vector of the DEM
nyArray : 2D numpy array
y component of the normal vector of the DEM
nzArray : 2D numpy array
z component of the normal vector of the DEM
gradient : int
Return the gradient (if 1) or the force associated (if 0, default)
Returns
-------
GHX : 1D numpy array
x component of the lateral force
GHY : 1D numpy array
y component of the lateral force
GHZ : 1D numpy array
z component of the lateral force
"""
# get all inputs
# configuration parameters
cdef double rho = cfg.getfloat('rho')
cdef double minRKern = cfg.getfloat('minRKern')
cdef double velMagMin = cfg.getfloat('velMagMin')
cdef double gravAcc = cfg.getfloat('gravAcc')
cdef int interpOption = cfg.getint('interpOption')
cdef int viscOption = cfg.getint('viscOption')
# grid normal raster information
cdef double cszNormal = headerNormalGrid['cellsize']
cdef int nRowsNormal = headerNormalGrid['nrows']
cdef int nColsNormal = headerNormalGrid['ncols']
# neighbour search grid information and neighbour information
cdef double cszNeighbourGrid = headerNeighbourGrid['cellsize']
cdef int nRowsNeighbourGrid = headerNeighbourGrid['nrows']
cdef int nColsNeighbourGrid = headerNeighbourGrid['ncols']
cdef int[:] indPartInCell = particles['indPartInCell']
cdef int[:] partInCell = particles['partInCell']
# SPH kernel
# use "spiky" kernel: w = (rKernel - r)**3 * 10/(pi*rKernel**5)
cdef double rKernel = cszNeighbourGrid
cdef double facKernel = 10.0 / (math.pi * rKernel * rKernel * rKernel * rKernel * rKernel)
cdef double dfacKernel = - 3.0 * facKernel
# particle information
cdef double[:] gEff = particles['gEff']
cdef double[:] mass = particles['m']
cdef double[:] hArray = particles['h']
cdef double[:] xArray = particles['x']
cdef double[:] yArray = particles['y']
cdef double[:] zArray = particles['z']
cdef double[:] uxArray = particles['ux']
cdef double[:] uyArray = particles['uy']
cdef double[:] uzArray = particles['uz']
cdef int N = xArray.shape[0]
# initialize variables and outputs
cdef double[:] GHX = np.zeros(N, dtype=np.float64)
cdef double[:] GHY = np.zeros(N, dtype=np.float64)
cdef double[:] GHZ = np.zeros(N, dtype=np.float64)
cdef double K1 = 1
cdef double K2 = 1
cdef double gravAcc3
cdef double gradhX, gradhY, gradhZ, uMag, nx, ny, nz, G1, G2, mdwdrr
cdef double g1, g2, g11, g12, g22, g33
cdef double x, y, z, ux, uy, uz, vx, vy, wx, wy, uxOrtho, uyOrtho, uzOrtho
cdef double dx, dy, dz, dux, duy, duz, dn, r, hr, dwdr, wKernel
cdef int Lx0, Ly0, iCell
cdef double w[4]
cdef int lInd, rInd
cdef int indx, indy
cdef int k, ic, n, p, l, imax, imin, iPstart, iPend
cdef int grad = gradient
# artificial viscosity parameters
cdef double dwdrr, area
cdef double pikl, flux
cdef double mk, ml, hk, hl, ck, cl, lambdakl
# loop on particles
for k in range(N):
gradhX = 0
gradhY = 0
gradhZ = 0
pikl = 0
G1 = 0
G2 = 0
m11 = 0
m12 = 0
m22 = 0
x = xArray[k]
y = yArray[k]
z = zArray[k]
ux = uxArray[k]
uy = uyArray[k]
uz = uzArray[k]
hk = hArray[k]
mk = mass[k]
# locate particle in SPH grid
indx = <int>math.round(x / cszNeighbourGrid)
indy = <int>math.round(y / cszNeighbourGrid)
if SPHoption > 1:
# get normal vector
Lx0, Ly0, iCell, w[0], w[1], w[2], w[3] = DFAtlsC.getCellAndWeights(x, y, nColsNormal, nRowsNormal, cszNormal, interpOption)
nx, ny, nz = DFAtlsC.getVector(Lx0, Ly0, w[0], w[1], w[2], w[3], nxArray, nyArray, nzArray)
nx, ny, nz = DFAtlsC.normalize(nx, ny, nz)
# projection of gravity on normal vector or use the effective gravity (gravity + curvature acceleration part)
# This was done I computeForce and is passed over here
gravAcc3 = gEff[k]
if SPHoption > 2:
uMag = DFAtlsC.norm(ux, uy, uz)
if uMag < velMagMin:
ux = 1
uy = 0
uz = -(1*nx + 0*ny) / nz
ux, uy, uz = DFAtlsC.normalize(ux, uy, uz)
K1 = 1
K2 = 1
else:
ux, uy, uz = DFAtlsC.normalize(ux, uy, uz)
uxOrtho, uyOrtho, uzOrtho = DFAtlsC.crossProd(nx, ny, nz, ux, uy, uz)
uxOrtho, uyOrtho, uzOrtho = DFAtlsC.normalize(uxOrtho, uyOrtho, uzOrtho)
g1 = nx/(nz)
g2 = ny/(nz)
# check if we are on the bottom ot top row!!!
lInd = -1
rInd = 2
if indy == 0:
lInd = 0
if indy == nRowsNeighbourGrid - 1:
rInd = 1
for n in range(lInd, rInd):
ic = (indx - 1) + nColsNeighbourGrid * (indy + n)
# make sure not to take particles from the other edge
imax = max(ic, nColsNeighbourGrid * (indy + n))
imin = min(ic+3, nColsNeighbourGrid * (indy + n + 1))
iPstart = indPartInCell[imax]
iPend = indPartInCell[imin]
# loop on all particles in neighbour boxes
for p in range(iPstart, iPend):
# index of particle in neighbour box
l = partInCell[p]
if k != l:
dx = xArray[l] - x
dy = yArray[l] - y
dz = zArray[l] - z
#----------------------------SPHOPTION = 1--------------------------------------
# SamosAT style (no reprojecion on the surface, dz = 0 and gz is used)
if SPHoption == 1:
dz = 0
# get norm of r = xk - xl
r = DFAtlsC.norm(dx, dy, dz)
if r < minRKern * rKernel:
# impose a minimum distance between particles
r = minRKern * rKernel
if r < rKernel:
hr = rKernel - r
dwdr = dfacKernel * hr * hr
mdwdrr = mass[l] * dwdr / r
gradhX = gradhX + mdwdrr*dx
gradhY = gradhY + mdwdrr*dy
gradhZ = gradhZ + mdwdrr*dz
gravAcc3 = gravAcc
#----------------------------SPHOPTION = 2--------------------------------------
# Compute the gradient in the cartesian coord system with reprojecion on the surface
# and dz != 0 and g3 are used
# It is here also possible to add the artificial viscosity comming from the Ata theory
elif SPHoption == 2:
# like option 1 but with dz!=0
# get norm of r = xk - xl
r = DFAtlsC.norm(dx, dy, dz)
if r < minRKern * rKernel:
# impose a minimum distance between particles
dx = minRKern * rKernel * dx
dy = minRKern * rKernel * dy
dz = minRKern * rKernel * dz
r = minRKern * rKernel
if r < rKernel:
hl = hArray[l]
hr = rKernel - r
dwdrr = dfacKernel * hr * hr / r
ml = mass[l]
area = ml/(rho*hl)
flux = gravAcc3*hl
if viscOption == 2:
# ATA artificial viscosity
dux = uxArray[l] - ux
duy = uyArray[l] - uy
duz = uzArray[l] - uz
ck = math.sqrt(gravAcc3*hk)
cl = math.sqrt(gravAcc3*hl)
lamdbakl = (ck+cl)/2
pikl = - lamdbakl * DFAtlsC.scalProd(dux, duy, duz, dx, dy, dz) / r
# SPH gradient computation - standard SPH formulation
gradhX = gradhX + (flux + pikl)*dwdrr*dx*area
gradhY = gradhY + (flux + pikl)*dwdrr*dy*area
gradhZ = gradhZ + (flux + pikl)*dwdrr*dz*area
#----------------------------SPHOPTION = 3--------------------------------------
# Compute the gradient in the local coord system (will allow us to add the earth pressure coef)
# and this time reprojecion on the surface, dz != 0 and g3 are used
if SPHoption == 3:
# get coordinates in local coord system
r1 = DFAtlsC.scalProd(dx, dy, dz, ux, uy, uz)
r2 = DFAtlsC.scalProd(dx, dy, dz, uxOrtho, uyOrtho, uzOrtho)
# impse r3=0 even if the particle is not exactly on the tengent plane
# get norm of r = xk - xl
r = DFAtlsC.norm(r1, r2, 0)
if r < minRKern * rKernel:
# impose a minimum distance between particles
r1 = minRKern * rKernel * r1
r2 = minRKern * rKernel * r2
r = minRKern * rKernel
if r < rKernel:
hr = rKernel - r
dwdr = dfacKernel * hr * hr
ml = mass[l]
mdwdrr = ml * dwdr / r
G1 = mdwdrr * K1*r1
G2 = mdwdrr * K2*r2
gradhX = gradhX + ux*G1 + uxOrtho*G2
gradhY = gradhY + uy*G1 + uyOrtho*G2
gradhZ = gradhZ + (- g1*(ux*G1 + uxOrtho*G2) - g2*(uy*G1 + uyOrtho*G2))
if grad == 1:
if SPHoption == 2:
GHX[k] = GHX[k] + gradhX / gravAcc3
GHY[k] = GHY[k] + gradhY / gravAcc3
GHZ[k] = GHZ[k] + gradhZ / gravAcc3
else:
GHX[k] = GHX[k] - gradhX / rho
GHY[k] = GHY[k] - gradhY / rho
GHZ[k] = GHZ[k] - gradhZ / rho
elif grad == 0:
if SPHoption == 2:
# grad here is g*grad(h)
GHX[k] = GHX[k] + gradhX * mk
GHY[k] = GHY[k] + gradhY * mk
GHZ[k] = GHZ[k] + gradhZ * mk
else :
GHX[k] = GHX[k] + gradhX*gravAcc3 / rho* mk
GHY[k] = GHY[k] + gradhY*gravAcc3 / rho* mk
GHZ[k] = GHZ[k] + gradhZ*gravAcc3 / rho* mk
return GHX, GHY, GHZ
def computeIniMovement(cfg, particles, dem, dT, fields):
""" add artifical viscosity effect on velocity
Parameters
------------
cfg: configparser
configuration for DFA simulation
particles : dict
particles dictionary at t
dem : dict
dictionary with dem information
dT : float
time step
fields: dict
fields dictionary
Returns
--------
particles: dict
updated particle dictionary
force: dict
force dictionary
"""
cdef int interpOption = cfg.getint('interpOption')
cdef double rho = cfg.getfloat('rho')
cdef double subgridMixingFactor = cfg.getfloat('subgridMixingFactorIni')
cdef double gravAcc = cfg.getfloat('gravAcc')
cdef double dt = dT
cdef int nPart = particles['nPart']
cdef double csz = dem['header']['cellsize']
cdef int nrows = dem['header']['nrows']
cdef int ncols = dem['header']['ncols']
cdef double[:, :] ZDEM = dem['rasterData']
cdef double[:, :] nxArray = dem['Nx']
cdef double[:, :] nyArray = dem['Ny']
cdef double[:, :] nzArray = dem['Nz']
cdef double[:, :] areaRatser = dem['areaRaster']
cdef double[:] mass = particles['m']
cdef double[:] hArray = particles['h']
cdef double[:] xArray = particles['x']
cdef double[:] yArray = particles['y']
cdef double[:] zArray = particles['z']
cdef double[:] uxArray = particles['ux']
cdef double[:] uyArray = particles['uy']
cdef double[:] uzArray = particles['uz']
cdef int[:] indXDEM = particles['indXDEM']
cdef int[:] indYDEM = particles['indYDEM']
cdef double[:, :] VX = fields['Vx']
cdef double[:, :] VY = fields['Vy']
cdef double[:, :] VZ = fields['Vz']
cdef double[:] forceX = np.zeros(nPart, dtype=np.float64)
cdef double[:] forceY = np.zeros(nPart, dtype=np.float64)
cdef double[:] forceZ = np.zeros(nPart, dtype=np.float64)
cdef double[:] forceFrict = np.zeros(nPart, dtype=np.float64)
cdef double[:] dM = np.zeros(nPart, dtype=np.float64)
cdef double[:] curvAcc = np.zeros(nPart, dtype=np.float64)
cdef double[:] gEff = np.zeros(nPart, dtype=np.float64)
cdef int indCellX, indCellY
cdef double areaPart, uMag, m, h
cdef double vMeanx, vMeany, vMeanz, vMeanNorm, dvX, dvY, dvZ
cdef double x, y, z, xEnd, yEnd, zEnd, ux, uy, uz, uxDir, uyDir, uzDir
cdef double nx, ny, nz, nxEnd, nyEnd, nzEnd, nxAvg, nyAvg, nzAvg
# variables for interpolation
cdef int Lx0, Ly0, LxEnd0, LyEnd0, iCell, iCellEnd
cdef double w[4]
cdef double wEnd[4]
cdef int k
force = {}
for k in range(nPart):
m = mass[k]
x = xArray[k]
y = yArray[k]
z = zArray[k]
h = hArray[k]
ux = uxArray[k]
uy = uyArray[k]
uz = uzArray[k]
indCellX = indXDEM[k]
indCellY = indYDEM[k]
# deduce area
areaPart = m / (h * rho)
# get cell and weights
Lx0, Ly0, iCell, w[0], w[1], w[2], w[3] = DFAtlsC.getCellAndWeights(x, y, ncols, nrows, csz, interpOption)
# get normal at the particle location
nx, ny, nz = DFAtlsC.getVector(Lx0, Ly0, w[0], w[1], w[2], w[3], nxArray, nyArray, nzArray)
nx, ny, nz = DFAtlsC.normalize(nx, ny, nz)
# add artificial viscosity
ux, uy, uz = addArtificialViscosity(m, h, dt, rho, ux, uy, uz, subgridMixingFactor, Lx0, Ly0,
w[0], w[1], w[2], w[3], VX, VY, VZ, nx, ny, nz)
# update velocity
uxArray[k] = ux
uyArray[k] = uy
uzArray[k] = uz
gEff[k] = gravAcc * nz
# save results
force['dM'] = np.asarray(dM)
force['forceX'] = np.asarray(forceX)
force['forceY'] = np.asarray(forceY)
force['forceZ'] = np.asarray(forceZ)
force['forceFrict'] = np.asarray(forceFrict)
particles['gEff'] = np.asarray(gEff)
particles['curvAcc'] = np.asarray(curvAcc)
particles['ux'] = np.asarray(uxArray)
particles['uy'] = np.asarray(uyArray)
particles['uz'] = np.asarray(uzArray)
particles['m'] = np.asarray(mass)
return particles, force