Source code for simulation.checkPoints

##~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~##
##                                                                                   ##
##  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 file saves simulation information to allow for checkpointing. It not only stores
some data for restarting a new simulation but also is the main entry to write the
simulation output.

**Badlands** uses `hdf5`_ library for outputs generation.
To read the **hdf5** files, the code creates a *XML schema* (**.xmf** files) describing
how the data is stored in the **hdf5** files.

Then `xdmf`_ files provides support to visualise temporal evolution of the output with
applications such as **Paraview** or Visit.

.. image:: img/output.png
   :scale: 50 %
   :alt: Badlands output folder
   :align: center

Above is a typical structure that you will find in **badlands** output folder:

- **h5** folder contains the **hdf5** data, all the information computed by the model are stored in these files. You will have at least the *tin* (surface) and *flow* (stream network) dataset and also the *sed* (stratigraphy) data if the stratal structure is computed in your simulation.
- **xmf** folder contains the XML files used to read the **hdf5** files contained in the **h5** folder.
- the **.xml** input file used to build this specific model.
- two **.xdmf** files for the surface (**tin_series.xdmf**) and the flow network (**flow_series.xdmf**) that read the **xmf** files through time.

"""

import time
import numpy as np
import os

if "READTHEDOCS" not in os.environ:
    from badlands import visualiseFlow, visualiseTIN, eroMesh


[docs]def write_checkpoints( input, recGrid, lGIDs, inIDs, tNow, FVmesh, force, flow, rain, elevation, fillH, cumdiff, cumhill, cumfail, wavediff, step, prop, mapero=None, cumflex=None, ): """ Create the checkpoint files (used for HDF5 output). Args: input: class containing XML input file parameters. recGrid: class describing the regular grid characteristics. lGIDs: global nodes indices. inIDs: local nodes indices. tNow: current time step FVmesh: finite volume mesh force: forcing parameters flow: flow parameters rain: rain value elevation: elevation mesh fillH: filled elevation mesh cumdiff: cumulative total erosion/deposition changes cumhill: cumulative hillslope erosion/deposition changes cumfail: cumulative failure induced erosion/deposition changes wavediff: cumulative wave induced erosion/deposition step: iterative step for output prop: proportion of each sediment type mapero: imposed erodibility maps cumflex: cumlative changes induced by flexural isostasy """ deepb = input.deepbasin if input.erolays and input.erolays >= 0: eroOn = True else: eroOn = False out_time = time.process_time() visXlim = np.zeros(2) visYlim = np.zeros(2) visXlim[0] = recGrid.rectX.min() visXlim[1] = recGrid.rectX.max() visYlim[0] = recGrid.rectY.min() visYlim[1] = recGrid.rectY.max() # Done when TIN has been built/rebuilt if FVmesh.outPts is None and FVmesh.outCells is None: FVmesh.outPts, FVmesh.outCells = visualiseTIN.output_cellsIDs( lGIDs, inIDs, visXlim, visYlim, FVmesh.node_coords[:, :2], FVmesh.cells ) tcells = np.zeros(1) tcells[0] = len(FVmesh.outCells) tnodes = np.zeros(1) tnodes[0] = len(lGIDs) # Done for every visualisation step flowIDs, polylines = visualiseFlow.output_Polylines( FVmesh.outPts, flow.receivers[FVmesh.outPts], visXlim, visYlim, FVmesh.node_coords[:, :2], ) fnodes = np.zeros(1) fnodes[0] = len(flowIDs) fline = np.zeros(1) fline[0] = len(polylines[:, 0]) # Compute flow parameters if deepb >= 5000.0: deepb = force.sealevel flow.view_receivers( fillH, elevation, FVmesh.neighbours, FVmesh.vor_edges, FVmesh.edge_length, lGIDs, deepb, ) # force.sealevel) flow.compute_parameters(elevation, force.sealevel) visdis = np.copy(flow.discharge) seaIDs = np.where(elevation < deepb)[0] # force.sealevel)[0] if len(seaIDs) > 0: visdis[seaIDs] = 1.0 flow.basinID[seaIDs] = -1 visdis[visdis < 1.0] = 1.0 rockOn = False if input.carbonate or input.pelagic: rockOn = True # Write HDF5 files if input.waveSed and tNow > input.tStart: waveOn = True meanH = force.meanH[lGIDs] meanS = force.meanS[lGIDs] wdiff = wavediff[lGIDs] else: waveOn = False meanH = None meanS = None wdiff = None if input.flexure: visualiseTIN.write_hdf5_flexure( input.outDir, input.th5file, step, FVmesh.node_coords[:, :2], elevation[lGIDs], fillH[lGIDs], rain[lGIDs], visdis[lGIDs], cumdiff[lGIDs], cumhill[lGIDs], cumfail[lGIDs], cumflex[lGIDs], FVmesh.outCells, input.oroRain, eroOn, flow.erodibility[lGIDs], FVmesh.control_volumes[lGIDs], waveOn, meanH, meanS, wdiff, rockOn, prop[lGIDs, :], force.sealevel, ) else: visualiseTIN.write_hdf5( input.outDir, input.th5file, step, FVmesh.node_coords[:, :2], elevation[lGIDs], fillH[lGIDs], rain[lGIDs], visdis[lGIDs], cumdiff[lGIDs], cumhill[lGIDs], cumfail[lGIDs], FVmesh.outCells, input.oroRain, eroOn, flow.erodibility[lGIDs], FVmesh.control_volumes[lGIDs], waveOn, meanH, meanS, wdiff, rockOn, prop[lGIDs, :], force.sealevel, ) if flow.sedload is not None: if flow.flowdensity is not None: visualiseFlow.write_hdf5( input.outDir, input.fh5file, step, FVmesh.node_coords[flowIDs, :2], elevation[flowIDs], visdis[flowIDs], flow.chi[flowIDs], flow.sedload[flowIDs], flow.flowdensity[flowIDs], flow.basinID[flowIDs], polylines, ) else: zeros = np.zeros(len(flowIDs), dtype=float) visualiseFlow.write_hdf5( input.outDir, input.fh5file, step, FVmesh.node_coords[flowIDs, :2], elevation[flowIDs], visdis[flowIDs], flow.chi[flowIDs], flow.sedload[flowIDs], zeros, flow.basinID[flowIDs], polylines, ) else: zeros = np.zeros(len(flowIDs), dtype=float) if flow.flowdensity is not None: visualiseFlow.write_hdf5( input.outDir, input.fh5file, step, FVmesh.node_coords[flowIDs, :2], elevation[flowIDs], visdis[flowIDs], flow.chi[flowIDs], zeros, flow.flowdensity[flowIDs], flow.basinID[flowIDs], polylines, ) else: visualiseFlow.write_hdf5( input.outDir, input.fh5file, step, FVmesh.node_coords[flowIDs, :2], elevation[flowIDs], visdis[flowIDs], flow.chi[flowIDs], zeros, zeros, flow.basinID[flowIDs], polylines, ) # Combine HDF5 files and write time series visualiseTIN.write_xmf( input.outDir, input.txmffile, input.txdmffile, step, tNow, tcells, tnodes, input.th5file, force.sealevel, input.flexure, input.oroRain, eroOn, waveOn, rockOn, prop.shape[1], ) visualiseFlow.write_xmf( input.outDir, input.fxmffile, input.fxdmffile, step, tNow, fline, fnodes, input.fh5file, ) print( " - Writing outputs (%0.02f seconds; tNow = %s)" % (time.process_time() - out_time, tNow) ) # Record erodibility maps if input.erolays and input.erolays >= 0: mapero.write_hdf5_erolay(step) return