6. Large Data Input and Output
6.1. netCDF Output
6.1.1. Install netCDF
Please install netCDF from here:
or
If you’re working with Debian-based distros, execute the provided commands to install netCDF.
sudo add-apt-repository universe
sudo apt-get update
sudo apt-get install zlib1g-dev libhdf5-dev libnetcdf-dev
now lets change our SConstruct:
netcdf_lib = env.Append(LIBS=['netcdf'])
zlib_lib = env.Append(LIBS=['z'])
hdf5_lib = env.Append(LIBS=['hdf5_serial'])
if 'LD_LIBRARY_PATH' not in env['ENV']:
env['ENV']['LD_LIBRARY_PATH'] = ''
env['ENV']['LD_LIBRARY_PATH'] = '/usr/lib/x86_64-linux-gnu' + env['ENV']['LD_LIBRARY_PATH']
6.1.2. Implementation of a class NetCdf
Now, let’s implement a NetCdf class by creating the following files: : NetCdf.cpp , NetCdf.h and NetCdf.test.cpp
Let’s proceed with the implementation of the header file named``NetCdf.h``
#ifndef TSUNAMI_LAB_IO_NETCDF
#define TSUNAMI_LAB_IO_NETCDF
#include "../../constants.h"
#include <cstring>
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <netcdf.h>
namespace tsunami_lab {
namespace io {
class NetCdf;
}
}
class tsunami_lab::io::NetCdf {
private:
int m_varIdX, m_varIdY, m_varIdHeight, m_varIdTime, m_varIdBathymetry, m_varIdImpulseX, m_varIdImpulseY;
t_idx m_time_step = 0;
static void checkNcErr( int i_err);
public:
void updateFile(t_idx i_nx,
t_idx i_ny,
t_idx i_stride,
t_real i_time,
t_real const * i_h,
t_real const * i_hu,
t_real const * i_hv,
const char* filename);
NetCdf(t_real l_nx,t_real l_ny, const char* filename);
void fillConstants(t_idx i_nx,
t_idx i_ny,
t_real i_dxy,
t_real i_domainstart_x,
t_real i_domainstart_y,
t_real i_stride,
t_real const * i_b,
const char* filename);
static int read( const char* filename,
const char* varname,
std::vector<t_real> &data);
};
#endif
In the document, you’ll find two functions: fillConstants, updateFile, and a constructor named netCDF.
netCDF:This constructer is responsible for creating the NetCDF output file.
fillConstants:The purpose of this function is to write the constants into the NetCDF output file.
updateFile:This function is designed to update the output file after each time step.
After incorporating the header file, we can proceed to develop the write function in the
NetCdf.cppfile.
void tsunami_lab::io::NetCdf::fillConstants(t_idx i_nx,
t_idx i_ny,
t_real i_dxy,
t_real i_domainstart_x,
t_real i_domainstart_y,
t_real i_stride,
t_real const * i_b,
const char* filename){
int l_ncId,l_err;
l_err = nc_open(filename,NC_WRITE, &l_ncId);
t_real *l_data0 = new t_real[i_ny];
for( t_idx l_iy = 0; l_iy < i_ny; l_iy++ )
{
l_data0[l_iy] = ((l_iy + 0.5) * i_dxy )+ i_domainstart_y;
}
// put y coordinates
l_err = nc_put_var_float(l_ncId, m_varIdY, l_data0);
checkNcErr(l_err);
delete[] l_data0;
t_real *l_data1 = new t_real[i_nx];
for(t_idx l_ix = 0; l_ix < i_nx; l_ix++)
{
l_data1[l_ix] = ((l_ix + 0.5) * i_dxy )+ i_domainstart_x;
}
// put x coordinates
l_err = nc_put_var_float(l_ncId, m_varIdX, l_data1);
checkNcErr(l_err);
delete[] l_data1;
t_real *l_data2 = new t_real[i_nx*i_ny];
for( t_idx l_iy = 0; l_iy < i_ny; l_iy++)
{
for( t_idx l_ix = 0; l_ix < i_nx; l_ix++)
{
t_idx l_id = (l_iy+1) * i_stride + (l_ix+1);
l_data2[(l_iy) * i_nx + (l_ix) ] = i_b[l_id];
}
}
// put bathymetry values
l_err = nc_put_var_float(l_ncId, m_varIdBathymetry, l_data2);
checkNcErr(l_err);
delete[] l_data2;
//close file for now
l_err = nc_close(l_ncId);
checkNcErr(l_err);
}
void tsunami_lab::io::NetCdf::updateFile(t_idx i_nx,
t_idx i_ny,
t_idx i_stride,
t_real i_time,
t_real const * i_h,
t_real const * i_hu,
t_real const * i_hv,
const char* filename){
int l_ncId, l_err;
//open file
l_err = nc_open(filename,NC_WRITE,&l_ncId);
checkNcErr(l_err);
//create memory
std::vector<t_real> l_temp_data_height(i_ny * i_nx);
std::vector<t_real> l_temp_data_momentum_x(i_ny * i_nx);
std::vector<t_real> l_temp_data_momentum_y(i_ny * i_nx);
//
for( t_idx l_iy = 1; l_iy < i_ny+1; l_iy++ ) {
for( t_idx l_ix = 1; l_ix < i_nx+1; l_ix++ ){
t_idx l_id = l_iy * i_stride + l_ix;
l_temp_data_height[(l_iy-1) * i_nx + (l_ix-1)] = i_h[l_id];
l_temp_data_momentum_x[(l_iy-1) * i_nx + (l_ix-1)] = i_hu[l_id];
l_temp_data_momentum_y[(l_iy-1) * i_nx + (l_ix-1)] = i_hv[l_id];
}
}
std::vector<size_t> l_startp = {m_time_step,0,0};
std::vector<size_t> l_endp = {1,i_ny,i_nx};
std::vector<ptrdiff_t> l_stridep = {1,1,1}; //no elements get skipped
l_err = nc_put_vars_float(l_ncId, m_varIdHeight, l_startp.data(), l_endp.data(), l_stridep.data(), l_temp_data_height.data());
checkNcErr(l_err);
l_err = nc_put_vars_float(l_ncId, m_varIdImpulseX, l_startp.data(), l_endp.data(), l_stridep.data(), l_temp_data_momentum_x.data());
checkNcErr(l_err);
l_err = nc_put_vars_float(l_ncId, m_varIdImpulseY, l_startp.data(), l_endp.data(), l_stridep.data(), l_temp_data_momentum_y.data());
checkNcErr(l_err);
//time step is how many timesteps there are and i_time what simtime it is
l_err = nc_put_var1_float(l_ncId, m_varIdTime, &m_time_step, &i_time);
checkNcErr(l_err);
l_err = nc_close(l_ncId);
checkNcErr(l_err);
m_time_step ++;
}
tsunami_lab::io::NetCdf::NetCdf(t_real l_nx,t_real l_ny,const char* filename) {
int l_ncId,l_err;
// Dimensions x, y, time
int l_dimXId, l_dimYId, l_dimTimeId;
int l_dim2Ids[2];
int l_dim3Ids[3];
std::cout << "generating netcdf-file ' " <<filename<< " '"<< std::endl;
l_err = nc_create(filename,NC_CLOBBER, &l_ncId);
checkNcErr(l_err);
//dimensions
l_err = nc_def_dim(l_ncId, "x", l_nx, &l_dimXId);
checkNcErr(l_err);
l_err = nc_def_dim(l_ncId, "y", l_ny, &l_dimYId);
checkNcErr(l_err);
l_err = nc_def_dim(l_ncId, "time", NC_UNLIMITED, &l_dimTimeId);
checkNcErr(l_err);
//variable x
l_err = nc_def_var(l_ncId, "x", NC_FLOAT, 1, &l_dimXId, &m_varIdX);
checkNcErr(l_err);
const char* units_attribute_x = "meter";
nc_put_att_text(l_ncId, m_varIdX, "units", strlen(units_attribute_x), units_attribute_x);
//variable y
l_err = nc_def_var(l_ncId, "y", NC_FLOAT, 1, &l_dimYId, &m_varIdY);
checkNcErr(l_err);
const char* units_attribute_y = "meter";
nc_put_att_text(l_ncId, m_varIdY, "units", strlen(units_attribute_y), units_attribute_y);
//variable time
l_err = nc_def_var(l_ncId, "time", NC_FLOAT, 1, &l_dimTimeId, &m_varIdTime);
checkNcErr(l_err);
const char* units_attribute_time = "seconds";
nc_put_att_text(l_ncId, m_varIdTime, "units", strlen(units_attribute_time), units_attribute_time);
//variable bathymetry
l_dim2Ids[0] =l_dimYId;
l_dim2Ids[1] =l_dimXId;
l_err = nc_def_var(l_ncId, "b", NC_FLOAT, 2, l_dim2Ids, &m_varIdBathymetry);
checkNcErr(l_err);
const char* units_attribute_Bathymetry = "meters";
nc_put_att_text(l_ncId, m_varIdBathymetry, "units", strlen(units_attribute_Bathymetry), units_attribute_Bathymetry);
//variable height
l_dim3Ids[0] =l_dimTimeId;
l_dim3Ids[1] =l_dimYId;
l_dim3Ids[2] =l_dimXId;
l_err = nc_def_var(l_ncId, "h", NC_FLOAT, 3, l_dim3Ids, &m_varIdHeight);
checkNcErr(l_err);
const char* units_attribute_height = "meters";
nc_put_att_text(l_ncId, m_varIdHeight, "units", strlen(units_attribute_height), units_attribute_height);
//variable momentum in x direction
l_err = nc_def_var(l_ncId, "hu", NC_FLOAT, 3, l_dim3Ids, &m_varIdImpulseX);
checkNcErr(l_err);
const char* units_attribute_impulseX = "meters";
nc_put_att_text(l_ncId, m_varIdImpulseX, "units", strlen(units_attribute_impulseX), units_attribute_impulseX);
//variable momentum in y direction
l_err = nc_def_var(l_ncId, "hv", NC_FLOAT, 3, l_dim3Ids, &m_varIdImpulseY);
checkNcErr(l_err);
const char* units_attribute_impulseY = "meters";
nc_put_att_text(l_ncId, m_varIdImpulseY, "units", strlen(units_attribute_impulseY), units_attribute_impulseY);
//enddef
l_err = nc_enddef( l_ncId );
checkNcErr( l_err );
}
lastly lets implement a test unit for the
NetCdf.test.cpp.
#include <catch2/catch.hpp>
#include "../../constants.h"
#include <sstream>
#include <vector>
#define private public
#include "NetCdf.h"
#include <filesystem>
#undef public
TEST_CASE( "Test the NetCdf-writer", "[NetCdfWrite]" ) {
tsunami_lab::t_real l_h[10] = { 0, 1, 2, 3, 4, 5, 6 ,7, 8 ,5};
tsunami_lab::t_real l_hu[10] = { 6, 5, 4, 3, 2, 1, 0 ,7, 8 ,5};
tsunami_lab::t_real l_hv[10] = { 6, 5, 4, 3, 2, 1, 0 ,7, 8 ,5};
tsunami_lab::t_real l_b[10] = { 0, 0, 0, 0, 0, 0, 0 ,7, 8 ,5};
tsunami_lab::io::NetCdf* l_netCdf = new tsunami_lab::io::NetCdf(5,2,"testsFiles/test.nc");
l_netCdf->fillConstants(5,2,1,-50,-50,0,l_b,"testsFiles/test.nc");
l_netCdf->updateFile(5,2,1,1,l_h,l_hu,l_hv,"testsFiles/test.nc");
int l_err;
int l_ncidp;
int l_dimXId,l_dimYId,l_TimeId;
size_t l_nx,l_ny , l_nt;
int l_vHId, l_vXId,l_vHuId, l_vYId, l_vTimeId, l_vHvId ,l_vBId;
REQUIRE(std::filesystem::exists("testsFiles/test.nc"));
l_err = nc_open("testsFiles/test.nc", NC_NOWRITE, &l_ncidp);
l_netCdf->checkNcErr(l_err);
l_err = nc_inq_dimid(l_ncidp,"x",&l_dimXId);
l_netCdf->checkNcErr(l_err);
l_err = nc_inq_dimid(l_ncidp,"y",&l_dimYId);
l_netCdf->checkNcErr(l_err);
l_err= nc_inq_dimid(l_ncidp, "time", &l_TimeId);
l_netCdf->checkNcErr(l_err);
l_err = nc_inq_dimlen(l_ncidp, l_dimXId, &l_nx);
l_netCdf->checkNcErr(l_err);
l_err = nc_inq_dimlen(l_ncidp, l_dimYId, &l_ny);
l_netCdf->checkNcErr(l_err);
l_err = nc_inq_dimlen(l_ncidp, l_TimeId, &l_nt);
l_netCdf->checkNcErr(l_err);
REQUIRE(l_nx == 5);
REQUIRE(l_ny == 2);
REQUIRE(l_nt == 1);
l_err = nc_inq_varid(l_ncidp, "x", &l_vXId);
l_netCdf->checkNcErr(l_err);
l_err = nc_inq_varid(l_ncidp, "y", &l_vYId);
l_netCdf->checkNcErr(l_err);
l_err = nc_inq_varid(l_ncidp, "time", &l_vTimeId);
l_netCdf->checkNcErr(l_err);
l_err = nc_inq_varid(l_ncidp, "h", &l_vHId);
l_netCdf->checkNcErr(l_err);
l_err = nc_inq_varid(l_ncidp, "hu", &l_vHuId);
l_netCdf->checkNcErr(l_err);
l_err = nc_inq_varid(l_ncidp, "hv", &l_vHvId);
l_netCdf->checkNcErr(l_err);
l_err = nc_inq_varid(l_ncidp, "b", &l_vBId);
l_netCdf->checkNcErr(l_err);
std::vector<tsunami_lab::t_real> x(10);
std::vector<tsunami_lab::t_real> y(2);
std::vector<tsunami_lab::t_real> time(1);
std::vector<tsunami_lab::t_real> height(10);
std::vector<tsunami_lab::t_real> momentum_x(10);
std::vector<tsunami_lab::t_real> momentum_y(10);
std::vector<tsunami_lab::t_real> bathymetry(10);
l_err = nc_get_var_float(l_ncidp, l_vXId, &x[0]);
l_netCdf->checkNcErr(l_err);
l_err = nc_get_var_float(l_ncidp, l_vYId, &y[0]);
l_netCdf->checkNcErr(l_err);
l_err = nc_get_var_float(l_ncidp, l_vTimeId, &time[0]);
l_netCdf->checkNcErr(l_err);
l_err = nc_get_var_float(l_ncidp, l_vHId, &height[0]);
l_netCdf->checkNcErr(l_err);
l_err = nc_get_var_float(l_ncidp, l_vHuId, &momentum_x[0]);
l_netCdf->checkNcErr(l_err);
l_err = nc_get_var_float(l_ncidp, l_vHvId, &momentum_y[0]);
l_netCdf->checkNcErr(l_err);
l_err = nc_get_var_float(l_ncidp, l_vBId, &bathymetry[0]);
l_netCdf->checkNcErr(l_err);
REQUIRE(x[0] == Approx(-49.5));
REQUIRE(x[1] == Approx(-48.5));
REQUIRE(x[2] == Approx(-47.5));
REQUIRE(x[3] == Approx(-46.5));
REQUIRE(x[4] == Approx(-45.5));
REQUIRE(y[0] == Approx(-49.5));
REQUIRE(y[1] == Approx(-48.5));
REQUIRE(time[0] == Approx(1));
REQUIRE(momentum_x[0] == Approx(4));
REQUIRE(momentum_x[1] == Approx(3));
REQUIRE(momentum_x[2] == Approx(2));
REQUIRE(momentum_x[3] == Approx(1));
REQUIRE(momentum_x[4] == Approx(0));
REQUIRE(momentum_x[5] == Approx(3));
REQUIRE(momentum_x[6] == Approx(2));
REQUIRE(momentum_x[7] == Approx(1));
REQUIRE(momentum_x[8] == Approx(0));
REQUIRE(momentum_x[9] == Approx(7));
REQUIRE(momentum_y[0] == Approx(4));
REQUIRE(momentum_y[1] == Approx(3));
REQUIRE(momentum_y[2] == Approx(2));
REQUIRE(momentum_y[3] == Approx(1));
REQUIRE(momentum_y[4] == Approx(0));
REQUIRE(momentum_y[5] == Approx(3));
REQUIRE(momentum_y[6] == Approx(2));
REQUIRE(momentum_y[7] == Approx(1));
REQUIRE(momentum_y[8] == Approx(0));
REQUIRE(momentum_y[9] == Approx(7));
for (size_t i = 0; i < 10; i++)
{
REQUIRE(bathymetry[i] == Approx(0));
}
delete l_netCdf;
}
6.2. NetCDF Input
6.2.1. Implementation and test for the artificial tsunami setup
Now, let’s proceed to implement the following files ArtificialTsunami2d.cpp , ArtificialTsunami2d.h and ArtificialTsunami2d.test.cpp utilizing the provided data.
\[\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}\]
and
\[\text{d}: [-500, +500]\times [-500, +500] \rightarrow \mathbb{R}\]
(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}\]
Let’s proceed with the implementation of the initial document
ArtificialTsunami2d.h:
#ifndef TSUNAMI_LAB_SETUPS_ARTIFICIALTSUNAMI2D_H
#define TSUNAMI_LAB_SETUPS_ARTIFICIALTSUNAMI2D_H
#include "../Setup.h"
#include <fstream>
#include <sstream>
#include <vector>
namespace tsunami_lab {
namespace setups {
class ArtificialTsunami2d;
}
}
/**
* ArtificialTsunami2d setup.
**/
class tsunami_lab::setups::ArtificialTsunami2d: public Setup {
private:
/**
* @param m_delta avoids running into numerical issues due to missing support for wetting and drying in our solver.
*/
t_real m_delta = 0;
/**
* @param m_Bin The swimming pool has a constant bathymetry of -100m before the “earthquake”
*/
t_real m_Bin = -100;
/**
* @brief the method adds the vertical displacement, typically caused by a subduction-zone earthquake.
* @param i_x is the distance from the Fukushima Daini Nuclear Power Plant.
*/
t_real displacement( t_real i_x ,t_real i_y) const;
t_real gFunction(t_real i_y) const;
t_real fFunction(t_real i_x) const;
public:
/**
* @brief The constructor for ArtificialTsunami2d.
* @param i_delta to avoid running into numerical issues (small value)
**/
ArtificialTsunami2d( t_real i_delta);
/**
* @brief is an arificial simulation so no inputs needed
* @return height is hardcoded to 100;
**/
t_real getHeight( t_real ,
t_real ) const;
/**
* @brief Gets the momentum in x-direction.
* @return Momentum in x-direction.
**/
t_real getMomentumX( t_real,
t_real) const;
/**
* @brief Gets the momentum in y-direction.
* @return Momentum in y-direction.
**/
t_real getMomentumY( t_real,
t_real ) const;
/**
* @brief Computes the bathymetry explained in <a href="https://scalable.uni-jena.de/opt/tsunami/chapters/assignment_5.html#equation-eq-artificial">5.2.1</a>.
* @param i_x position in the swimming pool.
* @param i_x position in the swimming pool.
* @return artificial computed Bathymetry.
*/
t_real getBathymetry( t_real i_x,
t_real i_y) const ;
};
#endif
Now, we will proceed with the implementation of the
ArtificialTsunami2d.cppfile :
#include "ArtificialTsunami2d.h"
#include <cmath>
#include <cstddef>
tsunami_lab::setups::ArtificialTsunami2d::ArtificialTsunami2d(t_real i_delta){
m_delta = i_delta;
}
tsunami_lab::t_real tsunami_lab::setups::ArtificialTsunami2d::getBathymetry( t_real i_x,
t_real i_y) const {
//m_Bin is hardcoded to -100
if(m_Bin < 0 ){
if( m_Bin < -m_delta){
return m_Bin + displacement(i_x,i_y);
}else{
return -m_delta + displacement(i_x,i_y);
}
}else{
if( m_Bin > m_delta){
return m_Bin + displacement(i_x,i_y);
}else{
return m_delta + displacement(i_x,i_y);
}
}
}
tsunami_lab::t_real tsunami_lab::setups::ArtificialTsunami2d::getHeight( t_real,
t_real )const{
//-100 m bathymetry and water level is at 0 -> 100m height
return 100;
}
tsunami_lab::t_real tsunami_lab::setups::ArtificialTsunami2d::displacement( t_real i_x,t_real i_y) const{
if((i_x > -500) && (i_y > -500) && (i_x < 500) && (i_y < 500) ){
return 5 * fFunction(i_x) * gFunction(i_y);
}else{
return 0;
}
}
tsunami_lab::t_real tsunami_lab::setups::ArtificialTsunami2d::fFunction(t_real i_x) const{
return sin(((i_x/500)+1)* M_PI );
}
tsunami_lab::t_real tsunami_lab::setups::ArtificialTsunami2d::gFunction(t_real i_y) const{
return -((i_y/500)*(i_y/500))+1;
}
tsunami_lab::t_real tsunami_lab::setups::ArtificialTsunami2d::getMomentumY( t_real ,
t_real )const{
return 0;
}
tsunami_lab::t_real tsunami_lab::setups::ArtificialTsunami2d::getMomentumX( t_real,
t_real )const{
return 0;
}
Finally, let’s implement the test unit for the ArtificialTsunami2d in the
ArtificialTsunami2d.test.cppfile:
#include <catch2/catch.hpp>
#include "ArtificialTsunami2d.h"
#include "string"
TEST_CASE("Test the two-dimensional ArtificialTsunami2d setup.", "[ArtificialTsunami2d]")
{
tsunami_lab::setups::ArtificialTsunami2d* l_ArtificialTsunami2d = new tsunami_lab::setups::ArtificialTsunami2d(20);
tsunami_lab::t_real l_momentumX = l_ArtificialTsunami2d->getMomentumX(0,0);
REQUIRE(l_momentumX == Approx(0));
tsunami_lab::t_real l_momentumY = l_ArtificialTsunami2d->getMomentumY(0,0);
REQUIRE(l_momentumY == Approx(0));
tsunami_lab::t_real l_bathymetryValue = l_ArtificialTsunami2d->getBathymetry(0,0);
REQUIRE(l_bathymetryValue == Approx(-99.84296f));
l_bathymetryValue = l_ArtificialTsunami2d->getBathymetry(9,0);
REQUIRE(l_bathymetryValue == Approx(-99.84296f));
l_bathymetryValue = l_ArtificialTsunami2d->getBathymetry(0,4);
REQUIRE(l_bathymetryValue == Approx(-100.15704f));
l_bathymetryValue = l_ArtificialTsunami2d->getBathymetry(9,4);
REQUIRE(l_bathymetryValue == Approx( -100.15704f));
tsunami_lab::t_real l_heightValue = l_ArtificialTsunami2d->getHeight(2,1);
l_heightValue = l_ArtificialTsunami2d->getHeight(0,0);
REQUIRE(l_heightValue == 100.0f);
l_heightValue = l_ArtificialTsunami2d->getHeight(0,4);
REQUIRE(l_heightValue == 100.0f);
l_heightValue = l_ArtificialTsunami2d->getHeight(9,0);
REQUIRE(l_heightValue == 100.0f);
l_heightValue = l_ArtificialTsunami2d->getHeight(9,4);
REQUIRE(l_heightValue == 100.0f);
}
//within the boundary [-500,500]x[-500,500]->R
//x = 300 , y = 300
/*(300 > -500) && (300 > -500) && (300 < 500) && (300 < 500) -> true
= m_Bin + displacement(i_x,i_y)
= -100 + 5 * fFunction(i_x) * gFunction(i_y);
= -100 + 5 * sin(((300/500)+1)* M_PI ) * (-((300/500)*(300/500))+1)
ca ~ 100 + -3.04338
= -103.04338f
*/
REQUIRE( l_ArtificialTsunami2d.getHeight( 300, 300 ) == 100.0f);
REQUIRE( l_ArtificialTsunami2d.getBathymetry( 300, 300 ) == -103.04338f );
REQUIRE( l_ArtificialTsunami2d.getMomentumX( 300, 300 ) == 0);
REQUIRE( l_ArtificialTsunami2d.getMomentumY( 300, 300 ) == 0 );
/*
x = 10 , y = 100;
-> both values in the domain ->
/*(10 > -500) && (10 > -500) && (100 < 500) && (100 < 500) -> true
= m_Bin + displacement(i_x,i_y)
= -100 + 5 * fFunction(i_x) * gFunction(i_y);
= -100 + 5 * sin(((10/500)+1)* M_PI ) * (-((100/500)*(100/500))+1)
ca ~ -100 + -0.30139449
= -100.30139449f
*/
REQUIRE( l_ArtificialTsunami2d.getHeight( 10, 100 ) == 100.0f );
REQUIRE( l_ArtificialTsunami2d.getBathymetry( 10, 100 ) == -100.30139449f );
REQUIRE( l_ArtificialTsunami2d.getMomentumX( 10, 100 ) == 0);
REQUIRE( l_ArtificialTsunami2d.getMomentumY( 10, 100 ) == 0 );
}
6.2.2. adding support for reading netCDF files
To enable reading from our NetCDF files, we need to implement the read function:
To implement the function, we’ll utilize the filename, variablename (which may represent a dimension or a variable), and a vector to store the data. The function is expected to return an integer indicating the state. If the function returns -1, it should result in the program’s exit.
now lests implement the following fucntion in the NetCdf.cpp the fowllowing function:
int tsunami_lab::io::NetCdf::read( const char* filename,
const char* varname,
std::vector<t_real> & data){
int l_ncId;
//if an error occurs it return -1
if (nc_open(filename, NC_NOWRITE, &l_ncId) != NC_NOERR) {
std::cerr << "Error opening NetCDF file: " << filename << std::endl;
return -1;
}
int varid;
//get the variable identifier
if (nc_inq_varid(l_ncId, varname, &varid) != NC_NOERR) {
std::cerr << "Error getting variable ID for variable: " << varname << std::endl;
nc_close(l_ncId);
return -1;
}
int ndims;
int dimids[NC_MAX_VAR_DIMS];
// get amount of dimensions (ndims) and the identifiers of the dimensions (dimids)
if (nc_inq_var(l_ncId, varid, nullptr, nullptr, &ndims, dimids, nullptr) != NC_NOERR) {
std::cerr << "Error getting variable information for variable: " << varname << std::endl;
nc_close(l_ncId);
return -1;
}
//get the length of the dimension
t_idx dataSize = 1;
for (int i = 0; i < ndims; ++i) {
char dimname[NC_MAX_NAME + 1];
size_t dimlen;
if (nc_inq_dim(l_ncId, dimids[i], dimname, &dimlen) != NC_NOERR) {
std::cerr << "Error getting dimension information." << std::endl;
nc_close(l_ncId);
}
dataSize *= dimlen;
}
//write the data in the vector
data.resize(dataSize);
int l_err = nc_get_var_float(l_ncId, varid , &data[0]);
checkNcErr(l_err);
if (nc_close(l_ncId) != NC_NOERR) {
std::cerr << "Error closing NetCDF file: " << filename << std::endl;
}
return 1;
}
now we have to add the function to the header file NetCdf.h:
int read( const char* filename,
const char* varname,
std::vector<t_real> &data);
Let’s now create a test unit for the reader function within the``NetCdf.test.cpp`` file:
Now, locate the file we utilized either in our project folder or download it from the provided source testfile
TEST_CASE( "Test the NetCdf-reader ", "[NetCdfreader]" ) {
std::vector<tsunami_lab::t_real> data;
tsunami_lab::io::NetCdf* l_netCdf = new tsunami_lab::io::NetCdf(10,10,"testsFiles/reader.nc");
l_netCdf->read("testsFiles/testFileForTheReader.nc" , "data", data ) ;
for (tsunami_lab::t_real i = 0; i < 71; i++)
{
REQUIRE(i == data[i]);
}
delete l_netCdf;
}
6.2.3. Integration the new class TsunamiEvent2d
lets implement TsunamiEvent2d setup. for the fowllowing setup we will need to implement the fowllowing files :
TsunamiEvent2d.cpp , TsunamiEvent2d.h and TsunamiEvent2d.test.cpp.
Prior to incorporating TsunamiEvent2d, it is essential to develop the following function, which will determine the nearest available value in the grid file in case the exact coordinates are unavailable:
tsunami_lab::t_idx tsunami_lab::setups::TsunamiEvent2d::findClosestIndex(const std::vector<t_real>& vec, t_real value) const {
t_idx closestIndex = 0;
for (t_idx index = 0; index < vec.size(); ++index) {
if (vec[index] > value) {
if (value - vec[index - 1] > vec[index] - value){
closestIndex = index ;
}
else{
closestIndex = index- 1;
}
break;
}
}
return closestIndex;
}
Now, let’s commence the implementation of the following header file,``TsunamiEvent2d.h`` :
/**
* @author Ward Tammaa
*
* @section DESCRIPTION
**/
#ifndef TSUNAMI_LAB_SETUPS_TSUNAMIEVENT2D_H
#define TSUNAMI_LAB_SETUPS_TSUNAMIEVENT2D_H
#include "../Setup.h"
#include <fstream>
#include <sstream>
#include <vector>
namespace tsunami_lab {
namespace setups {
class TsunamiEvent2d;
}
}
/**
* TsunamiEvent2d setup.
**/
class tsunami_lab::setups::TsunamiEvent2d: public Setup {
private:
/**
* @param m_delta avoids running into numerical issues due to missing support for wetting and drying in our solver.
*/
t_real m_delta = 0;
t_real m_width_bathymetry = 0;
t_real m_length_bathymetry = 0;
t_real m_width_displacement = 0;
t_real m_length_displacement = 0;
/**
* @param m_bathymetry_values contains all bathymetry values
*/
std::vector<t_real> m_bathymetry_values;
std::vector<t_real> m_bathymetry_x_values;
std::vector<t_real> m_bathymetry_y_values;
std::vector<t_real> m_displacement_values;
std::vector<t_real> m_displacement_x_values;
std::vector<t_real> m_displacement_y_values;
t_real getBathymetryNetCdf(t_real i_x, t_real i_y) const;
t_idx findClosestIndex(const std::vector<t_real>& vec, t_real value) const;
public:
/**
* @brief the method adds the vertical displacement, typically caused by a subduction-zone earthquake.
* @param i_x is the distance from the Fukushima Daini Nuclear Power Plant.
*/
t_real displacement( t_real i_x,t_real i_y) const;
/**
* @brief The constructor for TsunamiEvent2d.
* @param i_delta to avoid running into numerical issues (small value)
**/
TsunamiEvent2d( t_real i_delta);
/**
* @brief Gets the water height at a given point.
*
* @param i_x is the distance from the Fukushima Daini Nuclear Power Plant
* @return Height at the given point.
**/
t_real getHeight( t_real i_x,
t_real i_y) const;
/**
* @brief Gets the momentum in x-direction.
*
* @return Momentum in x-direction.
**/
t_real getMomentumX( t_real ,
t_real ) const;
/**
* @brief Gets the momentum in y-direction.
* @return Momentum in y-direction.
**/
t_real getMomentumY( t_real,
t_real ) const;
/**
* @brief Computes the bathymetry as explained in <a href="https://scalable.uni-jena.de/opt/tsunami/chapters/assignment_3.html#equation-eq-tsunami-event-1d">3.4.1</a>.
* @param i_x is the distance from the Fukushima Daini Nuclear Power Plant.
* @return Bathymetry.
*/
t_real getBathymetry( t_real i_x,
t_real i_y) const ;
/**
* @brief divide the distance i_x by 250 to determine the index for the bathymetry in the csv file.
* @param i_x is the distance from the Fukushima Daini Nuclear Power Plant.
* @return Bathymetry in the csv file (not the value that we use).
*/
};
#endif
Let’s now proceed with the implementation of the TsunamiEvent2d.cpp file. Before we implement the class, we need to read the bathymetry and displacement data from the files given in the task.
#include "TsunamiEvent2d.h"
#include "../../io/NetCdf/NetCdf.h"
#include <cmath>
#include <cstddef>
tsunami_lab::setups::TsunamiEvent2d::TsunamiEvent2d(t_real i_delta)
{
tsunami_lab::io::NetCdf::read("data/artificialtsunami_bathymetry_1000.nc","z",m_bathymetry_values);
tsunami_lab::io::NetCdf::read("data/artificialtsunami_bathymetry_1000.nc","x",m_bathymetry_x_values);
tsunami_lab::io::NetCdf::read("data/artificialtsunami_bathymetry_1000.nc","y",m_bathymetry_y_values);
tsunami_lab::io::NetCdf::read("data/artificialtsunami_displ_1000.nc","z",m_displacement_values);
tsunami_lab::io::NetCdf::read("data/artificialtsunami_displ_1000.nc","x",m_displacement_x_values);
tsunami_lab::io::NetCdf::read("data/artificialtsunami_displ_1000.nc","y",m_displacement_y_values);
m_delta = i_delta;
}
tsunami_lab::t_real tsunami_lab::setups::TsunamiEvent2d::getBathymetry( t_real i_x,
t_real i_y) const {
t_real l_batNetCdf = getBathymetryNetCdf(i_x, i_y);
if(l_batNetCdf < 0 )
{
if(l_batNetCdf < -m_delta)
{
return l_batNetCdf + displacement(i_x, i_y);
}
else
{
return displacement(i_x, i_y) - m_delta;
}
}
else
{
if(l_batNetCdf < m_delta)
{
return displacement(i_x, i_y) + m_delta;
}
else
{
return l_batNetCdf + displacement(i_x, i_y);
}
}
}
tsunami_lab::t_real tsunami_lab::setups::TsunamiEvent2d::displacement( t_real i_x,t_real i_y) const {
if (i_x < m_displacement_x_values[0] || i_x > m_displacement_x_values[m_displacement_x_values.size() - 1] ||
i_y < m_displacement_y_values[0] || i_y > m_displacement_y_values[m_displacement_y_values.size() - 1])
{
return 0;
}
//nachdem man schon eigentlich den index raushat wollen wir ihn nochmal berechnen?????
t_idx l_x = findClosestIndex(m_displacement_x_values, i_x);
t_idx l_y = findClosestIndex(m_displacement_y_values, i_y);
return m_displacement_values[l_x * m_displacement_y_values.size() + l_y];
}
tsunami_lab::t_real tsunami_lab::setups::TsunamiEvent2d::getHeight( t_real i_x,
t_real i_y)const{
t_real l_bin = getBathymetryNetCdf(i_x, i_y);
if (l_bin < 0) {
return (-l_bin > m_delta) ? -l_bin : m_delta;
}
return 0;
}
tsunami_lab::t_real tsunami_lab::setups::TsunamiEvent2d::getBathymetryNetCdf(t_real i_x, t_real i_y) const {
//check whether the position is within our domain
if (i_x < m_bathymetry_x_values[0] || i_x > m_bathymetry_x_values[m_bathymetry_x_values.size() - 1] ||
i_y < m_bathymetry_y_values[0] || i_y > m_bathymetry_y_values[m_bathymetry_y_values.size() - 1])
{
return 0;
}
t_idx l_x = findClosestIndex(m_bathymetry_x_values, i_x);
t_idx l_y = findClosestIndex(m_bathymetry_y_values, i_y);
return m_bathymetry_values[l_y * m_bathymetry_x_values.size() + l_x];
}
tsunami_lab::t_idx tsunami_lab::setups::TsunamiEvent2d::findClosestIndex(const std::vector<t_real>& vec, t_real value) const {
t_idx closestIndex = 0;
for (t_idx index = 0; index < vec.size(); ++index) {
if (vec[index] > value) {
if (value - vec[index - 1] > vec[index] - value){
closestIndex = index ;
}
else{
closestIndex = index- 1;
}
break;
}
}
return closestIndex;
}
tsunami_lab::t_real tsunami_lab::setups::TsunamiEvent2d::getMomentumY( t_real,
t_real)const{
return 0;
}
tsunami_lab::t_real tsunami_lab::setups::TsunamiEvent2d::getMomentumX( t_real ,
t_real)const{
return 0;
}
Finally, let’s proceed with implementing the unit tests in the TsunamiEvent2d.test.cpp file:
#include <catch2/catch.hpp>
#include "TsunamiEvent2d.h"
#include "string"
TEST_CASE("Test the two-dimensional tsunamiEvent setup.", "[TsunamiEvent2d]")
{
tsunami_lab::setups::TsunamiEvent2d* l_tsunamiEvent2d = new tsunami_lab::setups::TsunamiEvent2d(20);
tsunami_lab::t_real l_momentumX = l_tsunamiEvent2d->getMomentumX(0,0);
REQUIRE(l_momentumX == Approx(0));
tsunami_lab::t_real l_momentumY = l_tsunamiEvent2d->getMomentumY(0,0);
REQUIRE(l_momentumY == Approx(0));
tsunami_lab::t_real l_bathymetryValue = l_tsunamiEvent2d->getBathymetry(0,0);
REQUIRE(l_bathymetryValue == Approx(-99.84296f));
l_bathymetryValue = l_tsunamiEvent2d->getBathymetry(9,0);
REQUIRE(l_bathymetryValue == Approx(-99.84296f));
l_bathymetryValue = l_tsunamiEvent2d->getBathymetry(0,4);
REQUIRE(l_bathymetryValue == Approx(-100.15704f));
l_bathymetryValue = l_tsunamiEvent2d->getBathymetry(9,4);
REQUIRE(l_bathymetryValue == Approx( -100.15704f));
tsunami_lab::t_real l_heightValue = l_tsunamiEvent2d->getHeight(2,1);
l_heightValue = l_tsunamiEvent2d->getHeight(0,0);
REQUIRE(l_heightValue == 100.0f);
l_heightValue = l_tsunamiEvent2d->getHeight(0,4);
REQUIRE(l_heightValue == 100.0f);
l_heightValue = l_tsunamiEvent2d->getHeight(9,0);
REQUIRE(l_heightValue == 100.0f);
l_heightValue = l_tsunamiEvent2d->getHeight(9,4);
REQUIRE(l_heightValue == 100.0f);
}
6.2.4. Check the correctness
The last question was to compare the ArtificialTsunami2d with the TsunamiEvent2d the first animation is from the
ArtificialTsunami2d:
And now the TsunamiEvent2d animation
The second video had these settings : “domain_start_x” : -5050,”domain_start_y” : -5050 ,”dimension_x” : 10100,”dimension_y” : 10100,
As you can see both simulation are very similar but the best way to test the correctness is by having the same tests for both setups. In case they haven’t already noticed :D
6.3. Personal Contribution
Ward Tammaa, Daniel Schicker Doxygen Documentation
Mohamad Khaled Minawe, Ward Tammaa, Daniel Schicker Sphnix Documentation
Daniel Schicker, Mohamad Khaled Minawe , Ward Tammaa functions implementation
Mohamad Khaled Minawe, Daniel Schicker, Ward Tammaa Unit Testing
Mohamad Khaled Minawe, Daniel Schicker Geogebra Datei(Calculations for the Unit Tests)
Ward Tammaa Hosting the code , Action runner