5. Large Data Input and Output

This part of the project extends our current solver, already capable of simulating artificial two-dimensional setups, with the most important data input for tsunami simulations: bathymetry and vertical displacements of the sea floor generated by earthquakes.

5.1. NetCDF Output

Up to this point we used simple ASCII-CSV files for our file output. This chapter adds support for binary netCDF (network Common Data Form) to our software. An introduction to netCDF and a description of its data model are given in the NetCDF Users Guide.

In addition to using a binary format over an ASCII one which results in faster file accesses, netCDF is also a single-file and self-describing format: We can store all time steps in one file and provide additional information. This removes overhead introduced by numerous files.

Tasks

  1. Install netCDF on your system and provide a mechanism which shares netCDF’s files with your solver and allows you to link with netCDF.

    Hints:

    • You’ll have to consider the netCDF dependencies zlib and HDF5.

    • SCons provides the courtesy methods CheckLib and CheckLibWithHeader to aid your discovery of available libraries.

  2. Implement a class io::NetCdf which acts as your interface to netCDF files. Implement a function which allows you to write your current time step to a netCDF file. Make sure to follow the COARDS convention (see Conventions for the standardization of NetCDF files) up to the point that you can postprocess your data easily, e.g., through visualization in ParaView.

    Note:

    • Values in ghost cells should not be written by default.

    • The keyword “since” in the attribute description, e.g., “seconds since the earthquake event”, triggers the time domain in ParaView.

    • NetCDF supports unlimited dimensions which is helpful to realize the time dimension.

5.2. NetCDF Input

In the previous parts of the project, we mostly used hardcoded and artificial setups for our initial values. Our CSV-based bathymetry input for the 1D tsunami simulations in Ch. 3.4 is a notable exception. This chapter adds a new setup for 2D simulations to our software which handles data input for the bathymetry and vertical displacements. For this we use two netCDF files following the COARDS conventions:

  • The first file contains space-dependent bathymetry \(b_\text{in}(x,y)\).

  • The second file contains space-dependent vertical displacements \(d(x,y)\). Ignoring temporary ground shanking, these displacements summarize the vertical and permanent displacements of the sea floor caused by an earthquake.

../_images/gebco20_ucsb3_topo_displ_annotated.png

Fig. 5.2.1 Illustration of the input data for the March 11, 2011, M 9.1 Tohoku earthquake and tsunami event. The bathymetry data was extracted from the GEBCO_2020 grid. The displacements have been derived from the UCSB preliminary model III. Shown is a cut in West-East direction with \(x \in [-200\,\text{km}, 500\,\text{km}]\) cut through the epicenter (\(y=0\)). The given bathymetry data is shown in black. The vertical displacements, scaled by a factor of 100 \(\times\), in red.

As an example Fig. 5.2.1 shows a vertical cut through bathymetry and displacement data for the March 11, 2011, M 9.1 Tohoku event. Each file contains the data, either the bathymetry or the displacement, using a regular grid. Thus, we have:

  1. Grid coordinates given through the coordinate variables (see COARDS convention) \(x\) and \(y\),

  2. An associated value in a third variable \(z\).

The computational domain of the setups covers the extend of the given bathymetry grid. However, the given displacement data may only cover smaller regions. In other parts of the computational domain the displacement is zero.

In the following, we will implement a setup setups::TsunamiEvent2d which is able to handle data input. Similar to Ch. 3.4’s one-dimensional setups::TsunamiEvent1d, we use the following initial values for the water height \(h\), momenta \(hu\) and \(hv\), and bathymetry \(b\):

\[\begin{split}\begin{split} h &= \begin{cases} \max( -b_\text{in}, \delta), &\text{if } b_\text{in} < 0 \\ 0, &\text{else} \end{cases}\\ hu &= 0\\ hv &= 0\\ b &= \begin{cases} \min(b_\text{in}, -\delta) + d, & \text{ if } b_\text{in} < 0\\ \max(b_\text{in}, \delta) + d, & \text{ else}. \end{cases} \end{split}\end{split}\]

Depending on the resolution of our simulations we have to request bathymetry and displacement data at coordinates not necessarily matching the grid coordinates of our input files. Thus, we fall back to using the closest available value in the grid file if the exact coordinates are not available.

../_images/query.svg

Fig. 5.2.2 The solver queries the value at point \(\times (30,7.5)\). The scenario returns point \({\scriptsize \blacksquare} (25,0)\).

An example of such a query is illustrated in Fig. 5.2.2. Note, that this might lead to aliasing issues which we ignore for now.

To test our netCDF reader, we simulate a “tsunami in a swimming pool”. The swimming pool has a constant bathymetry of -100m before the “earthquake” and a size of 10km \(\times\) 10km. The “earthquake” generates a displacement in a 1km \(\times\) 1km square in the center of the swimming pool.

../_images/sin.png

Fig. 5.2.3 Visualization of the vertical displacement of the artificial “tsunami in a swimming pool”.

The displacement is illustrated in Fig. 5.2.3 and can be described via the function \(\text{d}(x,y)\) where \((0,0)\) is the center of the computational domain:

\[\text{d}: [-500, +500]\times [-500, +500] \rightarrow \mathbb{R}\]
(5.2.1)\[\begin{split}\begin{aligned} \text{d}(x, y) & = & 5 \cdot f(x)g(y) \\ \text{f}(x) & = & \sin\left(\left(\frac{x}{500}+1\right) \cdot \pi\right) \\ \text{g}(y) & = & -\left(\frac{y}{500}\right)^2 + 1 \end{aligned}\end{split}\]

Tasks

  1. Implement and test the artificial tsunami setup by hard-coding the functions in Eq. 5.2.1. (no file input) in a new class setups::ArtificialTsunami2d.

  2. Extend your class io:NetCdf of Ch. 5.1 with support for reading netCDF files as described above.

    Hints:

    • The tool ncgen can be used to generate small netCDF files for your unit tests.

    • You can use the command ncdump -h <file.nc> to display available variables and dimensions in a netCDF file

  3. Integrate the new class setups::TsunamiEvent2d into your code. Users should be able to set the total simulation time and resolution through runtime options, parsed from the command line or a configuration file.

  4. Check the correctness of your file input-based class setups::TsunamiEvent2d by comparing it to your artificial setup in setups::ArtificialTsunami2d. Use the linked input data of the artificial tsunami for bathymetry and displacements in your tests.