import os
import gc
import sys
import h5py
import glob
import shutil
import petsc4py
import numpy as np
from mpi4py import MPI
from time import process_time
from gospl.tools.constants import DISCHARGE_FLOOR
MPIrank = petsc4py.PETSc.COMM_WORLD.Get_rank()
MPIsize = petsc4py.PETSc.COMM_WORLD.Get_size()
MPIcomm = MPI.COMM_WORLD
[docs]
class WriteMesh(object):
"""
Class for writing model outputs. The outputs are written as hdf5 files for each mesh partition.
.. note::
The model outputs are all located in an user-defined output folder (`dir` key in the inputfile documentation) and consist of a time series file named `gospl.xdmf` and 2 other folders (`h5` and `xmf`).
The `XDMF` file is the main entry point for visualising the output and should be sufficient for most users. This file can easely be opened within `Paraview <https://www.paraview.org/download/>`_.
"""
def __init__(self):
"""
Initialise model outputs parameters.
"""
self.step = 0
self.stratStep = 1
self.file = "gospl"
# HDF5 dataset compression options, applied to every create_dataset via
# `**self._h5opts`. Driven by the YAML `output: compression:` setting
# (parsed in inputparser._readOut -> self.outCompress):
# "gzip" (default) -> {"compression": "gzip"} (historical behaviour)
# int 0..9 -> gzip at that level (trade CPU for file size)
# "none"/False/None -> {} (uncompressed: fastest I/O, largest files)
# gzip is CPU-serial per dataset and can dominate wall-clock at scale
# and high output cadence; disabling it trades disk for speed.
comp = getattr(self, "outCompress", "gzip")
if comp in (None, False, "none", "None", ""):
self._h5opts = {}
elif isinstance(comp, bool): # True -> default gzip
self._h5opts = {"compression": "gzip"}
elif isinstance(comp, int):
self._h5opts = {"compression": "gzip", "compression_opts": comp}
else:
self._h5opts = {"compression": comp}
if MPIrank == 0:
self._createOutputDir()
# Sync the chosen output dir to all processors
self.outputDir = MPIcomm.bcast(self.outputDir, root=0)
# In case of a restarting simulation, get the corresponding time
# values according to the restarting step
if self.rStep > 0:
self.step = self.rStep + 1
if self.strat > 0:
n = int(self.tout / self.strat)
self.stratStep = self.rStep * n + 1
self.tNow = self.tStart + self.rStep * self.tout
self.rStart = self.tStart + self.rStep * self.tout
self.saveTime = self.tNow + self.tout
self.saveStrat = self.tEnd + self.tout
if self.strat > 0:
self.saveStrat = self.tNow + self.strat
else:
self.rStart = self.tStart
if self.strataFile is not None:
self.stratStep = self.initLay
# Petsc vectors
self.upsG = self.hGlobal.duplicate()
self.upsL = self.hLocal.duplicate()
return
[docs]
def visModel(self):
"""
Main function to write model outputs on disk.
"""
# Output time step for first step
if self.saveTime == self.tStart:
self._outputMesh()
self.saveTime += self.tout
# Output time step after start time
elif self.tNow >= self.saveTime:
self._outputMesh()
self.saveTime += self.tout
return
[docs]
def _createOutputDir(self):
"""
Create a directory to store outputs. By default the folder will be called `output`. If a folder name is specified in the YAML input file, this name will be used.
.. note::
The input option `makedir` gives the ability to delete any existing output folder with the same name (if set to `False`) or to create a new folder with the given dir name plus a number at the end (*e.g.* `outputDir_XX` if set to `True` with `XX` the run number). It prevents overwriting on top of previous runs.
"""
# Get output directory
if self.outputDir is not None:
self.outputDir = os.getcwd() + "/" + self.outputDir
else:
self.outputDir = os.getcwd() + "/output"
if self.rStep == 0:
if self.makedir:
if os.path.exists(self.outputDir):
# Find the first free numeric suffix (handles gaps in the
# numbering caused by manual deletion of intermediate dirs).
base = self.outputDir
n = 0
while os.path.exists(f"{base}_{n}"):
n += 1
self.outputDir = f"{base}_{n}"
else:
if os.path.exists(self.outputDir):
shutil.rmtree(self.outputDir, ignore_errors=True)
os.makedirs(self.outputDir)
os.makedirs(self.outputDir + "/h5")
os.makedirs(self.outputDir + "/xmf")
shutil.copy(self.finput, self.outputDir)
return
[docs]
def _outputStrat(self):
"""
Saves mesh local stratigraphic information stored in the DMPlex as HDF5 file. The following variables will be recorded:
- elevation at time of deposition, considered to be to the current elevation for the top stratigraphic layer `stratZ`.
- thickness of each stratigrapic layer `stratH` accounting for both erosion & deposition events.
- porosity of sediment `phiS` in each stratigraphic layer computed at center of each layer.
.. important::
It is worth mentioning that the stratigraphic architecture is only outputed as HDF5 files and does not record the XMF and XDMF files. A set of post-processing scripts are then required to extract the information and visualise the stratigraphic records at any specific time step.
"""
t = process_time()
h5file = (
self.outputDir
+ "/h5/stratal."
+ str(self.step)
+ ".p"
+ str(MPIrank)
+ ".h5"
)
with h5py.File(h5file, "w") as f:
# Write stratal layers elevations per layers
f.create_dataset(
"stratZ",
shape=(self.lpoints, self.stratStep + 1),
dtype="float32",
**self._h5opts,
)
f["stratZ"][:, : self.stratStep + 1] = self.stratZ[:, : self.stratStep + 1]
# Write stratal layers thicknesses per layers
f.create_dataset(
"stratH",
shape=(self.lpoints, self.stratStep + 1),
dtype="float64",
**self._h5opts,
)
f["stratH"][:, : self.stratStep + 1] = self.stratH[:, : self.stratStep + 1]
# Write porosity values for coarse sediments
f.create_dataset(
"phiS",
shape=(self.lpoints, self.stratStep + 1),
dtype="float64",
**self._h5opts,
)
f["phiS"][:, : self.stratStep + 1] = self.phiS[:, : self.stratStep + 1]
# Per-layer erodibility multiplier (applied on top of self.K
# in the SPL evaluation when the layer is exposed). Always 1.0
# for layers deposited by goSPL; can be non-uniform when the
# initial-strata file supplies a `stratK` field.
if self.stratK is not None:
f.create_dataset(
"stratK",
shape=(self.lpoints, self.stratStep + 1),
dtype="float64",
**self._h5opts,
)
f["stratK"][:, : self.stratStep + 1] = self.stratK[
:, : self.stratStep + 1
]
# Dual-lithology fine-fraction fields. Only written when the
# dual coarse/fine option is enabled (self.stratHf allocated in
# readStratLayers); single-fraction outputs are unchanged.
if self.stratHf is not None:
f.create_dataset(
"stratHf",
shape=(self.lpoints, self.stratStep + 1),
dtype="float64",
**self._h5opts,
)
f["stratHf"][:, : self.stratStep + 1] = self.stratHf[
:, : self.stratStep + 1
]
f.create_dataset(
"phiF",
shape=(self.lpoints, self.stratStep + 1),
dtype="float64",
**self._h5opts,
)
f["phiF"][:, : self.stratStep + 1] = self.phiF[:, : self.stratStep + 1]
# In-model provenance: per-layer per-class thickness (lpoints,
# layers, classes). Only written when provenance tracers are on.
if getattr(self, "provOn", False):
f.create_dataset(
"stratP",
shape=(self.lpoints, self.stratStep + 1, self.provNb),
dtype="float64",
**self._h5opts,
)
f["stratP"][:, : self.stratStep + 1, :] = self.stratP[
:, : self.stratStep + 1, :
]
MPIcomm.Barrier()
if MPIrank == 0 and self.verbose:
print(
"Creating stratal outputfile \
(%0.02f seconds)"
% (process_time() - t),
flush=True,
)
return
[docs]
def _outputMesh(self):
"""
Saves mesh local information stored in the DMPlex to HDF5 file. If the file already exists, it will be overwritten. Mesh characteristics are recorded for each partition. The following variables will be available:
- surface elevation `elev`.
- cumulative erosion & deposition values `erodep`.
- erosion & deposition rate values `EDrate` for the considered time step.
- flow accumulation `FA`.
- flow accumulation `fillFA` considering pit filling.
- river sediment load `sedLoad`.
- soil thickness `soilH` if soil production is considered.
- uplift subsidence values if vertical tectonic forcing is considered `uplift`.
- flexural isostasy rebound `flexIso` if flexure is considered.
- precipitation maps based on forcing conditions `rain` (could also correspond to the orographic rain if the functionality is turned on).
"""
t = process_time()
if self.step == 0:
topology = self.outputDir + "/h5/topology.p" + str(MPIrank) + ".h5"
with h5py.File(topology, "w") as f:
f.create_dataset(
"coords",
shape=(self.lpoints, 3),
dtype="float32",
**self._h5opts,
)
f["coords"][:, :] = self.lcoords
f.create_dataset(
"cells",
shape=(len(self.lcells[:, 0]), 3),
dtype="int32",
**self._h5opts,
)
f["cells"][:, :] = self.lcells + 1
self.elems = MPIcomm.gather(len(self.lcells[:, 0]), root=0)
self.nodes = MPIcomm.gather(self.lpoints, root=0)
h5file = (
self.outputDir
+ "/h5/"
+ self.file
+ "."
+ str(self.step)
+ ".p"
+ str(MPIrank)
+ ".h5"
)
with h5py.File(h5file, "w") as f:
f.create_dataset(
"elev",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
f["elev"][:, 0] = self.hLocal.getArray()
f.create_dataset(
"erodep",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
f["erodep"][:, 0] = self.cumEDLocal.getArray()
f.create_dataset(
"EDrate",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
data = self.EbLocal.getArray().copy()
f["EDrate"][:, 0] = data
if not self.fast:
f.create_dataset(
"waterFill",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
f["waterFill"][:, 0] = self.waterFilled
f.create_dataset(
"fillFA",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
data = self.fillFAL.getArray().copy()
data[data <= DISCHARGE_FLOOR] = DISCHARGE_FLOOR
if not self.fast:
data[self.seaID] = 1.0
f["fillFA"][:, 0] = data
f.create_dataset(
"FA",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
data = self.FAL.getArray().copy()
data[data <= DISCHARGE_FLOOR] = DISCHARGE_FLOOR
if not self.fast:
data[self.seaID] = 1.0
f["FA"][:, 0] = data
if self.iceOn:
f.create_dataset(
"iceH",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
data = self.iceHL.getArray().copy()
data[data <= DISCHARGE_FLOOR] = DISCHARGE_FLOOR
if not self.fast:
data[self.seaID] = 1.0
f["iceH"][:, 0] = data
# Basal sliding speed (m/yr): the abrasion driver and a
# diagnostic of ice dynamics.
f.create_dataset(
"iceUb",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
f["iceUb"][:, 0] = self.iceUbL.getArray().copy()
# Ablation meltwater (m^3/yr) re-injected into the rivers: the
# glacial contribution to downstream discharge.
f.create_dataset(
"iceMelt",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
f["iceMelt"][:, 0] = self.iceMeltRiverL.getArray().copy()
# Glacial abrasion rate E_g = Kg|u_b|^l (m/yr); zero where
# abrasion is off (Kg = 0).
f.create_dataset(
"iceAbr",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
f["iceAbr"][:, 0] = self.iceAbrL.getArray().copy()
# Ice discharge (m^3/yr) from the diagnostic 'mfd' flow model
# (the ELA accumulation routed downhill).
f.create_dataset(
"iceFA",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
f["iceFA"][:, 0] = self.iceFAL.getArray().copy()
if self.flexOn:
f.create_dataset(
"flexIso",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
f["flexIso"][:, 0] = self.localFlex
if self.cptSoil:
f.create_dataset(
"soilH",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
f["soilH"][:, 0] = self.Lsoil.getArray().copy()
f.create_dataset(
"sedLoad",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
data = self.vSedLocal.getArray().copy()
data[data <= DISCHARGE_FLOOR] = DISCHARGE_FLOOR
f["sedLoad"][:, 0] = data
# Dual-lithology fine sediment load (m3/yr). `sedLoad` is the TOTAL
# (coarse+fine) flux; `sedLoadF` is the fine sub-flux, so coarse =
# sedLoad - sedLoadF. Only written when dual lithology is enabled.
if self.stratLith:
f.create_dataset(
"sedLoadF",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
dataF = self.vSedFLocal.getArray().copy()
dataF[dataF <= DISCHARGE_FLOOR] = DISCHARGE_FLOOR
f["sedLoadF"][:, 0] = dataF
if self.upsub is not None:
f.create_dataset(
"uplift",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
f["uplift"][:, 0] = self.upsub
if self.rainVal is not None:
f.create_dataset(
"rain",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
f["rain"][:, 0] = self.rainVal
if self.evapdata is not None and self.evapVal is not None:
f.create_dataset(
"evap",
shape=(self.lpoints, 1),
dtype="float32",
**self._h5opts,
)
f["evap"][:, 0] = self.evapVal
if self.memclear:
del data
gc.collect()
if self.stratNb > 0 and self.stratStep > 0:
self._outputStrat()
if MPIrank == 0:
self._save_DMPlex_XMF()
self._save_XDMF()
MPIcomm.Barrier()
if MPIrank == 0 and self.verbose:
print(
"Creating outputfile (%0.02f seconds)" % (process_time() - t),
flush=True,
)
if MPIrank == 0:
print("+++ Output Simulation Time: %0.02f years" % (self.tNow), flush=True)
self.step += 1
return
[docs]
def readData(self):
"""
When a simulation restarts, variables from previous HDF5 output files are read and assigned to the restarting run.
The following variables are used:
- surface elevation `elev`.
- cumulative erosion & deposition values `erodep`.
- erosion & deposition values `EDrate` for the considered time step.
- flow accumulation `fillFA` considering pit filling.
- river sediment load `sedLoad`.
- flexural isostasy induced tectonics `flexIso`.
- soil thickness `soilH`.
.. note::
If stratigraphy is turned on, the function also reads underlying stratigraphic information.
"""
h5file = (
self.outputDir
+ "/h5/"
+ self.file
+ "."
+ str(self.step - 1)
+ ".p"
+ str(MPIrank)
+ ".h5"
)
if not os.path.exists(h5file):
raise ValueError("Restart file is missing...")
with h5py.File(h5file, "r") as hf:
self.hLocal.setArray(np.array(hf["/elev"])[:, 0])
self.dm.localToGlobal(self.hLocal, self.hGlobal)
self.cumEDLocal.setArray(np.array(hf["/erodep"])[:, 0])
self.dm.localToGlobal(self.cumEDLocal, self.cumED)
self.vSedLocal.setArray(np.array(hf["/sedLoad"])[:, 0])
self.dm.localToGlobal(self.vSedLocal, self.vSed)
# Restart-safe restore of the fine sediment load (dual lithology).
# Older single-fraction outputs have no /sedLoadF; the field is
# recomputed each step in _getSedFlux anyway, so a missing dataset
# just leaves the init zeros.
if self.stratLith and "/sedLoadF" in hf:
self.vSedFLocal.setArray(np.array(hf["/sedLoadF"])[:, 0])
self.dm.localToGlobal(self.vSedFLocal, self.vSedF)
self.EbLocal.setArray(np.array(hf["/EDrate"])[:, 0])
self.dm.localToGlobal(self.EbLocal, self.Eb)
self.FAL.setArray(np.array(hf["/FA"])[:, 0])
self.fillFAL.setArray(np.array(hf["/fillFA"])[:, 0])
self.dm.localToGlobal(self.FAL, self.FAG)
self.elems = MPIcomm.gather(len(self.lcells[:, 0]), root=0)
self.nodes = MPIcomm.gather(len(self.lcoords[:, 0]), root=0)
if self.iceOn:
if "/iceH" in hf:
self.iceHL.setArray(np.array(hf["/iceH"])[:, 0])
else:
self.iceHL.set(0.)
if self.flexOn:
if "/flexIso" in hf:
self.localFlex = np.array(hf["/flexIso"])[:, 0]
else:
self.localFlex = np.zeros(self.lpoints)
if self.cptSoil:
if "/soilH" in hf:
self.Lsoil.setArray(np.array(hf["/soilH"])[:, 0])
else:
self.Lsoil.set(0.)
self.dm.localToGlobal(self.Lsoil, self.Gsoil)
if self.stratNb > 0 and self.stratStep > 0:
h5file = (
self.outputDir
+ "/h5/stratal."
+ str(self.step - 1)
+ ".p"
+ str(MPIrank)
+ ".h5"
)
if not os.path.exists(h5file):
raise ValueError("Restart file is missing...")
with h5py.File(h5file, "r") as hf:
self.stratZ.fill(0.0)
self.stratZ[:, : self.stratStep] = np.array(hf["/stratZ"])
self.stratH.fill(0.0)
self.stratH[:, : self.stratStep] = np.array(hf["/stratH"])
self.phiS.fill(0.0)
self.phiS[:, : self.stratStep] = np.array(hf["/phiS"])
if "/stratK" in hf and self.stratK is not None:
# Restart-safe restore of per-layer K multipliers.
# Older outputs without this dataset fall back to the
# initialised 1.0 array (no scaling), preserving
# behaviour for any pre-existing restart files.
self.stratK.fill(1.0)
self.stratK[:, : self.stratStep] = np.array(hf["/stratK"])
# Dual-lithology fine-fraction fields. Restarting a dual run
# from a single-fraction output (no /stratHf dataset) falls
# back to all-coarse zeros, so the restore stays robust.
if self.stratHf is not None:
self.stratHf.fill(0.0)
if "/stratHf" in hf:
self.stratHf[:, : self.stratStep] = np.array(hf["/stratHf"])
self.phiF.fill(0.0)
if "/phiF" in hf:
self.phiF[:, : self.stratStep] = np.array(hf["/phiF"])
# In-model provenance: restore per-layer per-class thickness.
# Restarting a provenance run from an output without /stratP
# falls back to the bedrock-seeded stratP (set in
# readStratLayers), so the restore stays robust.
if getattr(self, "provOn", False) and "/stratP" in hf:
self.stratP.fill(0.0)
self.stratP[:, : self.stratStep, :] = np.array(hf["/stratP"])
return
[docs]
def _save_DMPlex_XMF(self):
"""
Saves mesh local information stored in the HDF5 to XmF file. The XMF files are XML schema explaining how to read `gospl` data files.
The XmF file is written by a single processor (rank 0) and contains each partition HDF5 files in blocks. The variables described for the HDF5 file (function `_outputMesh` above) are all accessible from this file.
"""
xmf_file = self.outputDir + "/xmf/" + self.file + str(self.step) + ".xmf"
f = open(xmf_file, "w")
# Header for xml file
f.write('<?xml version="1.0" encoding="UTF-8"?>\n')
f.write('<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd">\n')
f.write('<Xdmf Version="2.0" xmlns:xi="http://www.w3.org/2001/XInclude">\n')
f.write(" <Domain>\n")
f.write(' <Grid GridType="Collection" CollectionType="Spatial">\n')
f.write(' <Time Type="Single" Value="%0.02f"/>\n' % self.saveTime)
for p in range(MPIsize):
pfile = (
"h5/" + str(self.file) + "." + str(self.step) + ".p" + str(p) + ".h5"
)
tfile = "h5/topology.p" + str(p) + ".h5"
f.write(' <Grid Name="Block.%s">\n' % (str(p)))
f.write(
' <Topology Type="Triangle" NumberOfElements="%d" BaseOffset="1">\n'
% self.elems[p]
)
f.write(' <DataItem Format="HDF" DataType="Int" ')
f.write('Dimensions="%d 3">%s:/cells</DataItem>\n' % (self.elems[p], tfile))
f.write(" </Topology>\n")
f.write(' <Geometry Type="XYZ">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 3">%s:/coords</DataItem>\n' % (self.nodes[p], tfile)
)
f.write(" </Geometry>\n")
f.write(' <Attribute Type="Scalar" Center="Node" Name="Z">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write('Dimensions="%d 1">%s:/elev</DataItem>\n' % (self.nodes[p], pfile))
f.write(" </Attribute>\n")
f.write(' <Attribute Type="Scalar" Center="Node" Name="ED">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/erodep</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
if not self.fast:
f.write(
' <Attribute Type="Scalar" Center="Node" Name="waterFill">\n'
)
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/waterFill</DataItem>\n'
% (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
f.write(' <Attribute Type="Scalar" Center="Node" Name="EDrate">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/EDrate</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
f.write(' <Attribute Type="Scalar" Center="Node" Name="FA">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/FA</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
f.write(' <Attribute Type="Scalar" Center="Node" Name="fillFA">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/fillFA</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
if self.iceOn:
f.write(' <Attribute Type="Scalar" Center="Node" Name="iceUb">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/iceUb</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
f.write(' <Attribute Type="Scalar" Center="Node" Name="iceH">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/iceH</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
f.write(' <Attribute Type="Scalar" Center="Node" Name="iceMelt">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/iceMelt</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
f.write(' <Attribute Type="Scalar" Center="Node" Name="iceAbr">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/iceAbr</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
f.write(' <Attribute Type="Scalar" Center="Node" Name="iceFA">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/iceFA</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
if self.flexOn:
f.write(' <Attribute Type="Scalar" Center="Node" Name="flexIso">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/flexIso</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
if self.cptSoil:
f.write(' <Attribute Type="Scalar" Center="Node" Name="soilH">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/soilH</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
f.write(' <Attribute Type="Scalar" Center="Node" Name="SL">\n')
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/sedLoad</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
if self.stratLith:
f.write(
' <Attribute Type="Scalar" Center="Node" Name="SLf">\n'
)
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/sedLoadF</DataItem>\n'
% (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
if self.upsub is not None:
f.write(
' <Attribute Type="Scalar" Center="Node" Name="vTec">\n'
)
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/uplift</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
if self.rainVal is not None:
f.write(
' <Attribute Type="Scalar" Center="Node" Name="Rain">\n'
)
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/rain</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
if self.evapdata is not None and self.evapVal is not None:
f.write(
' <Attribute Type="Scalar" Center="Node" Name="Evap">\n'
)
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/evap</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </Attribute>\n")
f.write(' <Attribute Type="Scalar" Center="Node" Name="sea">\n')
f.write(
' <DataItem ItemType="Function" Function="$0 * 0.00000000000000001 + %f" Dimensions="%d 1">\n'
% (self.sealevel, self.nodes[p])
)
f.write(
' <DataItem Format="HDF" NumberType="Float" Precision="4" '
)
f.write(
'Dimensions="%d 1">%s:/erodep</DataItem>\n' % (self.nodes[p], pfile)
)
f.write(" </DataItem>\n")
f.write(" </Attribute>\n")
f.write(" </Grid>\n")
f.write(" </Grid>\n")
f.write(" </Domain>\n")
f.write("</Xdmf>\n")
f.close()
return
[docs]
def _save_XDMF(self):
"""
This function writes the XDmF file which is calling the XmF files above. The XDmF file represents the *time series* of the model outputs and can be directly loaded and visualised with `Paraview <https://www.paraview.org/download/>`_.
.. note::
For a brief overview of the approach used to record `gospl` outputs, user can read this `visit documentation <https://www.visitusers.org/index.php?title=Using_XDMF_to_read_HDF5>`_
"""
xdmf_file = self.outputDir + "/" + self.file + ".xdmf"
with open(xdmf_file, "w") as f:
f.write('<?xml version="1.0" encoding="UTF-8"?>\n')
f.write('<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd">\n')
f.write('<Xdmf Version="2.0" xmlns:xi="http://www.w3.org/2001/XInclude">\n')
f.write(" <Domain>\n")
f.write(' <Grid GridType="Collection" CollectionType="Temporal">\n')
for s in range(self.step + 1):
xmf_file = "xmf/" + self.file + str(s) + ".xmf"
f.write(
' <xi:include href="%s" xpointer="xpointer(//Xdmf/Domain/Grid)"/>\n'
% xmf_file
)
f.write(" </Grid>\n")
f.write(" </Domain>\n")
f.write("</Xdmf>\n")
return
