##~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~##
## ##
## This file forms part of the Badlands surface processes modelling application. ##
## ##
## For full license and copyright information, please refer to the LICENSE.md file ##
## located at the project root, or contact the authors. ##
## ##
##~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~##
"""
This module defines the **stratigraphic layers** based on a **regular mesh**.
"""
import os
import glob
import time
import h5py
import numpy
from scipy import interpolate
from scipy.spatial import cKDTree
from scipy.interpolate import RegularGridInterpolator
if "READTHEDOCS" not in os.environ:
from badlands import flowalgo
[docs]class strataMesh:
"""
This class builds stratigraphic layer on each depositional point of the regular mesh.
Args:
sdx: discretisation value [m]
bbX: extent of stratal regular grid along X
bbY: extent of stratal regular grid along Y
layNb: total number of stratigraphic layers
xyTIN: numpy float-type array containing the coordinates for each nodes in the TIN (in m)
folder: name of the output folder.
h5file: first part of the hdf5 file name.
poro0: initial porosity.
poroC: porosity evolution coefficient.
cumdiff: numpy array containing cumulative erosion/deposition from previous simulation.
rfolder: restart folder.
rstep: restart step.
"""
def __init__(
self,
sdx,
bbX,
bbY,
layNb,
xyTIN,
folder,
h5file,
poro0,
poroC,
cumdiff=0,
rfolder=None,
rstep=0,
):
self.ids = None
self.ptsNb = None
self.oldload = None
self.tree = None
self.folder = folder
self.h5file = h5file + ".time"
self.step = 0
self.upper = None
self.lower = None
self.poro0 = poro0
self.poroC = poroC
self.xyTIN = xyTIN
# User defined parameter
self.dx = sdx
# Create stratal regular grid
self.nx = int(round((bbX[1] - bbX[0]) / sdx - 0.5) + 1)
self.ny = int(round((bbY[1] - bbY[0]) / sdx - 0.5) + 1)
self.xgrid = numpy.linspace(bbX[0], bbX[1], num=self.nx)
self.ygrid = numpy.linspace(bbY[0], bbY[1], num=self.ny)
xi, yi = numpy.meshgrid(self.xgrid, self.ygrid)
self.xyi = numpy.dstack([xi.flatten(), yi.flatten()])[0]
# Partition mesh
self.buildPartition(bbX, bbY)
self.ptsNb = len(self.ids)
if rstep > 0:
if os.path.exists(rfolder):
folder = rfolder + "/h5/"
fileCPU = "sed.time%s.hdf5" % (rstep)
restartncpus = len(glob.glob1(folder, fileCPU))
if restartncpus == 0:
raise ValueError(
"The requested time step for the restart simulation cannot be found in the restart folder."
)
else:
raise ValueError(
"The restart folder is missing or the given path is incorrect."
)
if restartncpus != 1:
raise ValueError(
"When using the stratal model you need to run the restart simulation with the same number of processors as the previous one."
)
df = h5py.File("%s/h5/sed.time%s.hdf5" % (rfolder, rstep), "r")
layDepth = numpy.array((df["/layDepth"]))
layElev = numpy.array((df["/layElev"]))
layThick = numpy.array((df["/layThick"]))
layPoro = numpy.array((df["/layPoro"]))
rstlays = layDepth.shape[1]
layNb += rstlays
self.step = rstlays
# Define global stratigraphic dataset
self.stratIn = numpy.zeros([self.ptsNb], dtype=int)
self.stratElev = numpy.zeros([self.ptsNb, layNb])
self.stratThick = numpy.zeros([self.ptsNb, layNb])
self.stratPoro = numpy.zeros([self.ptsNb, layNb])
self.stratDepth = numpy.zeros([self.ptsNb, layNb])
self.stratPoro.fill(self.poro0)
if rstep > 0:
self.stratDepth[:, :rstlays] = layDepth
self.stratElev[:, :rstlays] = layElev
self.stratThick[:, :rstlays] = layThick
self.stratPoro[:, :rstlays] = layPoro
# Define TIN grid kdtree for interpolation
self.tree = cKDTree(xyTIN)
tindx = xyTIN[1, 0] - xyTIN[0, 0]
self.searchpts = max(int(sdx * sdx / (tindx * tindx)), 4)
if rstep > 0:
distances, indices = self.tree.query(self.xyi, k=self.searchpts)
weights = 1.0 / distances ** 2
if len(cumdiff[indices].shape) == 3:
cum_vals = cumdiff[indices][:, :, 0]
else:
cum_vals = cumdiff[indices]
fcum = numpy.average(cum_vals, weights=weights, axis=1)
onIDs = numpy.where(distances[:, 0] == 0)[0]
if len(onIDs) > 0:
fcum[onIDs] = cumdiff[indices[onIDs, 0]]
self.oldload = cumdiff
return
[docs] def update_TIN(self, xyTIN):
"""
Update TIN mesh after 3D displacements.
Args:
xyTIN: numpy float-type array containing the coordinates for each nodes in the TIN (in m)
"""
# Update TIN grid kdtree for interpolation
self.tree = cKDTree(xyTIN)
self.xyTIN = xyTIN
return
[docs] def move_mesh(self, dispX, dispY, cumdiff, verbose=False):
"""
Update stratal mesh after 3D displacements.
Args:
dispX: numpy float-type array containing X-displacement for each nodes in the stratal mesh
dispY: numpy float-type array containing Y-displacement for each nodes in the stratal mesh
cumdiff: numpy float-type array containing the cumulative erosion/deposition of the nodes in the TIN
verbose : (bool) when :code:`True`, output additional debug information (default: :code:`False`).
"""
# Move coordinates
walltime = time.process_time()
st_time = walltime
moveXY = numpy.zeros([self.xyi.shape[0], 2])
moveXY[:, 0] = self.xyi[:, 0] + dispX
moveXY[:, 1] = self.xyi[:, 1] + dispY
# Define point ids in local partition which needs to be send to neighbourhood
l0 = self.nx
l1 = l0 + self.nx
u1 = self.ptsNb - self.nx
u0 = u1 - self.nx
shape = (l1 - l0, self.step + 1)
if verbose:
print(" - move stratal mesh ", time.process_time() - walltime)
walltime = time.process_time()
deformXY = moveXY
deformThick = self.stratThick[:, : self.step + 1]
deformPoro = self.stratPoro[:, : self.step + 1]
deformElev = self.stratElev[:, : self.step + 1]
if verbose:
print(
" - create deformed stratal mesh arrays ",
time.process_time() - walltime,
)
# Build the kd-tree
walltime = time.process_time()
deformtree = cKDTree(deformXY)
if verbose:
print(
" - create deformed stratal mesh kd-tree ",
time.process_time() - walltime,
)
walltime = time.process_time()
distances, indices = deformtree.query(self.xyi, k=4)
if verbose:
print(" - query stratal mesh kd-tree ", time.process_time() - walltime)
# Compute inverse weighting distance
walltime = time.process_time()
with numpy.errstate(divide="ignore"):
w = 1.0 / distances ** 2
w3D = w.reshape((len(self.xyi), 4, 1))
weights = numpy.tile(w3D, (1, 1, self.step + 1))
# Perform interpolation
tmpIDs = numpy.where(distances[:, 0] == 0)[0]
if len(tmpIDs) > 0:
self.stratThick[tmpIDs, : self.step + 1] = deformThick[
indices[tmpIDs, 0], : self.step + 1
]
self.stratPoro[tmpIDs, : self.step + 1] = deformPoro[
indices[tmpIDs, 0], : self.step + 1
]
self.stratElev[tmpIDs, : self.step + 1] = deformElev[
indices[tmpIDs, 0], : self.step + 1
]
tmpID = numpy.where(distances[:, 0] > 0)[0]
if len(tmpID) > 0:
self.stratThick[tmpID, : self.step + 1] = numpy.average(
deformThick[indices[tmpID, :], :], weights=weights[tmpID, :], axis=1
)
self.stratPoro[tmpID, : self.step + 1] = numpy.average(
deformPoro[indices[tmpID, :], :], weights=weights[tmpID, :], axis=1
)
self.stratElev[tmpID, : self.step + 1] = numpy.average(
deformElev[indices[tmpID, :], :], weights=weights[tmpID, :], axis=1
)
else:
self.stratThick[:, : self.step + 1] = numpy.average(
deformThick[indices, :], weights=weights, axis=1
)
self.stratPoro[:, : self.step + 1] = numpy.average(
deformPoro[indices, :], weights=weights, axis=1
)
self.stratElev[:, : self.step + 1] = numpy.average(
deformElev[indices, :], weights=weights, axis=1
)
# Reset depostion flag
self.stratIn.fill(0)
tmpID = numpy.where(
numpy.amax(self.stratThick[:, : self.step + 1], axis=1) > 0
)[0]
self.stratIn[tmpID] = 1
if verbose:
print(
" - perform stratal mesh interpolation ", time.process_time() - walltime
)
if verbose:
print(" - moving stratal mesh function ", time.process_time() - st_time)
self.oldload = numpy.copy(cumdiff)
return
[docs] def buildStrata(self, elev, cumdiff, sea, boundsPt, write=0, outstep=0):
"""
Build the stratigraphic layer on the regular grid.
Args:
elev: numpy float-type array containing the elevation of the nodes in the TIN
cumdiff: numpy float-type array containing the cumulative erosion/deposition of the nodes in the TIN
sea: sea level elevation
boundPts: number of nodes on the edges of the TIN surface.
write: flag for output generation
outstep: sStep for output generation
Returns:
- sub_poro - numpy array containing the subsidence induced by porosity change.
"""
selev = numpy.zeros(len(self.xyi))
distances, indices = self.tree.query(self.xyi, k=self.searchpts)
with numpy.errstate(divide="ignore"):
weights = 1.0 / distances ** 2
if self.oldload is not None:
load_diff = cumdiff - self.oldload
else:
load_diff = cumdiff
if len(elev[indices].shape) == 3:
elev_vals = elev[indices][:, :, 0]
cum_vals = load_diff[indices][:, :, 0]
else:
elev_vals = elev[indices]
cum_vals = load_diff[indices]
felev = numpy.average(elev_vals, weights=weights, axis=1)
fcum = numpy.average(cum_vals, weights=weights, axis=1)
onIDs = numpy.where(distances[:, 0] == 0)[0]
if len(onIDs) > 0:
felev[onIDs] = elev[indices[onIDs, 0]]
fcum[onIDs] = load_diff[indices[onIDs, 0]]
self.oldload = numpy.copy(cumdiff)
selev = felev
# Update stratal elevation
self.stratElev[self.ids, self.step] = selev[self.ids] - sea
localCum = fcum[self.ids]
# Update stratal erosion
eroIDs = numpy.where(localCum < 0.0)[0]
self.eroLayer(self.ids[eroIDs], localCum)
# Update stratal deposition
depIDs = numpy.where(localCum > 0.0)[0]
subs = self.depoLayer(self.ids[depIDs], localCum)
subsi = numpy.reshape(subs, (len(self.ygrid), len(self.xgrid)))
subs_values = RegularGridInterpolator((self.ygrid, self.xgrid), subsi)
sub_poro = numpy.zeros(len(self.xyTIN[:, 0]))
sub_poro[boundsPt:] = subs_values(
(self.xyTIN[boundsPt:, 1], self.xyTIN[boundsPt:, 0])
)
sub_poro[sub_poro > 0.0] = 0.0
self.oldload += sub_poro
if write > 0:
self.layerMesh(selev[self.ids] + subs[self.ids])
self.write_hdf5_stratal(outstep - 1)
self.step += 1
return sub_poro
[docs] def buildPartition(self, bbX, bbY):
"""
Define a partition for the stratal mesh.
Args:
bbX: box boundary for X-axis
bbY: box boundary for Y-axis
"""
size = 1
# extent of X partition
Yst = numpy.zeros(1)
Yed = numpy.zeros(1)
partYID = numpy.zeros((1, 2))
nbY = int((self.ny - 1) / 1)
p = 0
Yst[p] = bbY[0]
Yed[p] = Yst[p] + nbY * self.dx
partYID[p, 0] = 0
partYID[p, 1] = (nbY + 1) * self.nx
Yed[size - 1] = bbY[1]
partYID[size - 1, 1] = self.ny * self.nx
# Get send/receive data ids for each processors
self.upper = numpy.zeros((1, 2))
self.lower = numpy.zeros((1, 2))
self.upper[0, 0] = partYID[0, 1]
self.upper[0, 1] = partYID[0, 1] + self.nx
self.lower[0, 0] = partYID[0, 0]
self.lower[0, 1] = partYID[0, 0] - self.nx
# Define partitions ID globally
Xst = numpy.zeros(1)
Xed = numpy.zeros(1)
Xst += bbX[0]
Xed += bbX[1]
# Loop over node coordinates and find if they belong to local partition
# Note: used a Cython/Fython class to increase search loop performance... in utils
partID = flowalgo.overlap(
self.xyi[:, 0], self.xyi[:, 1], Xst[0], Yst[0], Xed[0], Yed[0]
)
# Extract local domain nodes global ID
self.ids = numpy.where(partID > -1)[0]
return
[docs] def depoLayer(self, ids, depo):
"""
Add deposit to current stratigraphic layer.
Args:
ids: index of points subject to deposition
depo: value of the deposition for the given point (m)
Returns:
- subs - numpy array containing the subsidence induced by porosity change.
"""
# Initialise node deposition flag
tmpIDs = numpy.where(self.stratIn[ids] == 0)
self.stratIn[ids[tmpIDs]] = 1
# Add deposit to the considered layer time
self.stratThick[ids, self.step] += depo[ids]
# Define porosity values
subs = numpy.zeros(self.ptsNb)
if self.poro0 > 0:
self.stratPoro[ids, self.step] = self.poro0
cumThick = numpy.cumsum(self.stratThick[ids, self.step :: -1], axis=1)[
:, ::-1
]
poro = self.poro0 * numpy.exp(-self.poroC * cumThick / 1000.0)
poro[poro < 0.1] = 0.1
# tmpid = numpy.where(self.stratPoro[ids,:self.step+1]<poro)[0]
tmp1, tmp2 = numpy.where(
numpy.logical_and(
self.stratPoro[ids, : self.step + 1] < poro,
self.stratThick[ids, : self.step + 1] > 0.0,
)
)
if len(tmp1) > 0:
poro[tmp1, tmp2] = self.stratPoro[ids[tmp1], tmp2]
nh = (
self.stratThick[ids, : self.step + 1]
* (1.0 - self.stratPoro[ids, : self.step + 1])
/ (1.0 - poro)
)
nh = numpy.minimum(self.stratThick[ids, : self.step + 1], nh)
nh[nh < 0.0] = 0.0
# Update layer thickness
self.stratThick[ids, : self.step + 1] = nh
# Update porosity
self.stratPoro[ids, : self.step + 1] = poro
# Subsidence due to porosity change
cumThick2 = numpy.cumsum(self.stratThick[ids, self.step :: -1], axis=1)[
:, ::-1
]
subs[ids] = (cumThick2 - cumThick)[:, 0]
# subs[ids] = numpy.sum(nh-self.stratThick[ids,:self.step+1],axis=1)
subs[subs > 0.0] = 0.0
return subs
[docs] def eroLayer(self, nids, erosion):
"""
Erode top stratigraphic layers.
Args:
nids: index of points subject to erosion
erosion: value of the erosion for the given points (m)
"""
# Perform erosion on nodes containing stratigraphic layers
tmpIDs = numpy.where(self.stratIn[nids] == 1)[0]
if len(tmpIDs) == 0:
return
# Update node indices and associated erosion values
ids = nids[tmpIDs]
ero = -erosion[ids]
# Compute cumulative stratal thicknesses
cumThick = numpy.cumsum(self.stratThick[ids, self.step :: -1], axis=1)[:, ::-1]
# Find nodes with no remaining stratigraphic thicknesses
tmpIDs = numpy.where(ero >= cumThick[:, 0])[0]
self.stratIn[ids[tmpIDs]] = 0
self.stratThick[ids[tmpIDs], : self.step + 1] = 0.0
# Erode remaining stratal layers
if len(tmpIDs) < len(ids):
ero[tmpIDs] = 0.0
# Clear all stratigraphy points which are eroded
cumThick[cumThick < ero.reshape((len(ids), 1))] = 0
mask = (cumThick > 0).astype(int) == 0
tmpH = self.stratThick[ids, : self.step + 1]
tmpH[mask] = 0
self.stratThick[ids, : self.step + 1] = tmpH
# Update thickness of top stratigraphic layer
eroIDs = numpy.bincount(numpy.nonzero(cumThick)[0]) - 1
eroVals = cumThick[numpy.arange(len(ids)), eroIDs] - ero
eroVals[tmpIDs] = 0.0
self.stratThick[ids, eroIDs] = eroVals
return
[docs] def layerMesh(self, topsurf):
"""
Define stratigraphic layers mesh.
Args:
topsurf: elevation of the regular surface
"""
# Clear points with no stratigraphic layer
tmpIDs = numpy.where(self.stratIn == 0)[0]
# surf = numpy.tile(topsurf[tmpIDs].transpose(), (1, self.step+1)).reshape(self.step+1,len(tmpIDs)).transpose()
surf = numpy.array([topsurf[tmpIDs],] * int(self.step + 2)).transpose()
self.stratDepth[tmpIDs, : self.step + 2] = surf
self.stratThick[tmpIDs, : self.step + 2] = 0.0
# Find points with stratigraphic layers
tmpIDs = numpy.where(self.stratIn == 1)[0]
if len(tmpIDs) == 0:
return
# Compute cumulative stratal thicknesses
cumThick = numpy.cumsum(self.stratThick[tmpIDs, self.step + 1 :: -1], axis=1)[
:, ::-1
]
# Updata stratal depth
surf = numpy.array([topsurf[tmpIDs],] * int(self.step + 2)).transpose()
self.stratDepth[tmpIDs, : self.step + 2] = surf - cumThick
return
[docs] def write_hdf5_stratal(self, outstep):
"""
This function writes for each processor the HDF5 file containing sub-surface information.
Args:
outstep: output time step.
"""
sh5file = self.folder + "/" + self.h5file + str(outstep) + ".hdf5"
with h5py.File(sh5file, "w") as f:
# Write node coordinates
f.create_dataset(
"coords", shape=(self.ptsNb, 2), dtype="float64", compression="gzip"
)
f["coords"][:, :2] = self.xyi[self.ids]
# Write stratal layers depth per cells
f.create_dataset(
"layDepth",
shape=(self.ptsNb, self.step + 2),
dtype="float64",
compression="gzip",
)
f["layDepth"][:, : self.step + 2] = self.stratDepth[
self.ids, : self.step + 2
]
# Write stratal layers elevations per cells
f.create_dataset(
"layElev",
shape=(self.ptsNb, self.step + 2),
dtype="float64",
compression="gzip",
)
f["layElev"][:, : self.step + 2] = self.stratElev[self.ids, : self.step + 2]
# Write stratal layers thicknesses per cells
f.create_dataset(
"layThick",
shape=(self.ptsNb, self.step + 2),
dtype="float64",
compression="gzip",
)
f["layThick"][:, : self.step + 2] = self.stratThick[
self.ids, : self.step + 2
]
# Write stratal layers porosity per cells
f.create_dataset(
"layPoro",
shape=(self.ptsNb, self.step + 2),
dtype="float64",
compression="gzip",
)
f["layPoro"][:, : self.step + 2] = self.stratPoro[self.ids, : self.step + 2]
return