3 Bathemetry & Boundary Conditions
Links
Individual Contributions
Justus Dreßler: all members contributed equally
Thorsten Kröhl: all members contributed equally
Julius Halank: all members contributed equally
3.1 Non-zero Source Term
3.1.1 Extend FWave Solver with Bathemetry
Now calculating the \(\Delta x \Psi _{i-1/2}\) term
void tsunami_lab::solvers::FWave::deltaXPsi(t_real i_bL,
t_real i_bR,
t_real i_hL,
t_real i_hR,
t_real &o_deltaXPsi)
{
// compute deltaXPsi
o_deltaXPsi = -m_g * (i_bR - i_bL) * (i_hL + i_hR) / 2;
}
and subtracting it from the flux jump.
deltaXPsi(i_bL, i_bR, i_hL, i_hR, l_deltaXPsi);
// compute jump in fluxes
t_real l_flux0Jump = l_flux0R - l_flux0L;
t_real l_flux1Jump = l_flux1R - l_flux1L - l_deltaXPsi;
3.1.2 Implement an example which illustrates the effect of bathymetry
Dam break with h_l=10 and h_r=2 and a bathymetry that increases from -2 to -1 and goes back to -2 at x=0, x=5 and x=10 respectively.
Dam break with h_l=10 and h_r=2 as reference.
3.2 Reflecting boundary conditions
3.2.1 Implement the reflecting boundary conditions
Implemented a reflecting boundary on a wet dry interface in the FWave Solver.
bool updateL = true;
bool updateR = true;
// if both dry do nothing
if (i_hL <= 0 && i_hR <= 0)
{
o_netUpdateL[0] = 0;
o_netUpdateL[1] = 0;
o_netUpdateR[0] = 0;
o_netUpdateR[1] = 0;
return;
} // if only left side is dry, apply reflecting boundary condition
else if (i_hL <= 0)
{
i_hL = i_hR;
i_huL = -i_huR;
i_bL = i_bR;
updateL = false;
} // if only right side is dry, apply reflecting boundary condition
else if (i_hR <= 0)
{
i_hR = i_hL;
i_huR = -i_huL;
i_bR = i_bL;
updateR = false;
}
updateL and updateR are used to determine if the cells should be updated or not (dry cells don’t change).
Added boundary conditions to the command line parameters as -b 'WALL OPEN' in the main function.
// boundary
case 'b':
{
std::string l_arg(optarg);
// convert to upper case
std::transform(l_arg.begin(), l_arg.end(), l_arg.begin(), ::toupper);
// split string by space
std::stringstream l_stream(l_arg);
std::string l_boundaryLName, l_boundaryRName;
l_stream >> l_boundaryLName >> l_boundaryRName;
std::cout << "using boundary conditions " << l_boundaryLName << " " << l_boundaryRName << std::endl;
// convert to t_boundary
getBoundary(l_boundaryLName, &l_boundaryL);
getBoundary(l_boundaryRName, &l_boundaryR);
break;
}
with a helper function that translates strings to t_boundary enum members
// converts a string to a boundary condition (tsunami_lab::t_boundary)
void getBoundary(std::string i_name, tsunami_lab::t_boundary *o_boundary)
{
if (i_name == "WALL")
{
*o_boundary = tsunami_lab::t_boundary::WALL;
}
else if (i_name == "OPEN")
{
*o_boundary = tsunami_lab::t_boundary::OPEN;
}
else
{
std::cerr << "unknown boundary condition " << i_name << std::endl;
exit(EXIT_FAILURE);
}
}
and switches the ghost cells depending on the boundary conditions in WavePropagation1d.
// set left boundary
switch (m_boundaryLeft)
{
case t_boundary::OPEN:
{
l_h[0] = l_h[1];
l_hu[0] = l_hu[1];
l_b[0] = l_b[1];
break;
}
case t_boundary::WALL:
{
l_h[0] = 0;
l_hu[0] = 0;
l_b[m_nCells + 1] = 20;
break;
}
}
3.2.2 Show the implementation with the shock shock setup
Added new setup to easier simulate tasks (with user controlled h_l h_r hu_l hu_r and middle position)
else if (l_setupName == "CUSTOM1D")
{
double l_arg3 = std::stof(l_arg3Str);
double l_arg4 = std::stof(l_arg4Str);
double l_arg5 = std::stof(l_arg5Str);
std::cout << "using Custom1d(" << l_arg1 << "," << l_arg2 << "," << l_arg3 << "," << l_arg4 << "," << l_arg5 << ") setup" << std::endl;
l_setup = new tsunami_lab::setups::Custom1d(l_arg1,
l_arg2,
l_arg3,
l_arg4,
l_arg5);
}
reflecting right boundary condition with open left boundary condition, h=10 and u=10
Shock-Shock problem with h=10 and u=10
3.3 Hydraulic Jumps
3.3.1 Compute the location and value of the maximum Froude number
In \(x \in [0,25]\) the maximum Froude number is given by
3.3.2 Implement both cases through the base class setup
Implemented both setups and changed their endtime in the main function.
tsunami_lab::t_real tsunami_lab::setups::Supercritical1d::getBathymetry(t_real i_x,
t_real) const
{
if (8 < i_x && i_x < 12)
{
return -0.13 - 0.05 * (i_x - 10) * (i_x - 10);
}
else
{
return -0.33;
}
}
else if (l_setupName == "SUPERCRIT1D")
{
l_width = 25.0; // 25 m domain
l_endTime = 200; // 200 s simulation time
std::cout << " using Supercritical1d() setup" << std::endl;
l_setup = new tsunami_lab::setups::Supercritical1d();
}
3.3.3 Determine the position of the hydraulic jump
The hydraulic jump occurs between \(x_{id}=45\) and \(x_{id}=47\), which would represent \(x=0.45 \cdot 25 = 11.25\) and \(x=0.47 \cdot 25 = 11.75\) respectively. You can see a distinct spike in momentum around \(x_{id}=46\) which is the failure of our f-wave solver to converge to the constant momentum.
3.4 Tsunami simulation
We will use a csvReader library rapidcsv in our reader. Is a header only library that you can include by just adding the header file to your project.
3.4.1 Extract bathymetry data with 250m sampling
after applying the following commands to the cut bathymetry grid we get the following csv (only excerpt shown)
gmt grdtrack -GGeco.nc -E141.024949/37.316569/146.0/37.316569+i250e+d -Ar > data.csv
cat data.csv | tr -s '[:blank:]'' ',' > data.csv
3.4.2 Extend CSV class with reader for bathymetry data
Use rapidcsv to read the csv file and access the fourth column of it.
The param rapidcsv::LabelParams(-1, -1) is used to tell the reader that the csv file has no header and no index column.
void tsunami_lab::io::Csv::openCSV(const std::string &i_filePath, rapidcsv::Document &o_doc, size_t &o_rowCount)
{
// assume headless csv
o_doc = rapidcsv::Document(i_filePath, rapidcsv::LabelParams(-1, -1));
o_rowCount = o_doc.GetRowCount();
}
tsunami_lab::t_real tsunami_lab::io::Csv::readLine(const rapidcsv::Document &i_doc, size_t i_row)
{
float o_row = i_doc.GetRow<float>(i_row)[3];
return o_row;
}
3.4.3 Implement a setup that initializes 1d tsunamis
Added a TsunamiEvent1d setup that uses a rapidcsv::document csv file to read the bathymetry data and a function to calculate the displacement.
tsunami_lab::setups::TsunamiEvent1d::TsunamiEvent1d(rapidcsv::Document i_doc, size_t i_rowCount)
{
m_doc = i_doc;
m_rowCount = i_rowCount;
}
The displacement function is just a simple sine function between 175km and 225km.
tsunami_lab::t_real tsunami_lab::setups::TsunamiEvent1d::getDisplacement(t_real i_x) const
{
if (175000 < i_x && i_x < 225000)
{
return 10 * std::sin(M_PI * (i_x - 175000) / 37500 + M_PI);
}
else
{
return 0;
}
}
And the bathymetry gets read of the csv file and the displacement is added to it. A \(\delta\) (minimum offset from 0m height) of 20m is used so we don’t run into numeric problems with cells wetting and drying.
tsunami_lab::t_real tsunami_lab::setups::TsunamiEvent1d::getBathymetry(t_real i_x,
t_real) const
{
t_real l_bin = getBathymetryBin(i_x);
if (l_bin < 0)
{
// min(bin, -delta) + d
if (l_bin < -m_delta)
return l_bin + getDisplacement(i_x);
else
return -m_delta + getDisplacement(i_x);
}
// max(bin, delta) + d
if (l_bin > m_delta)
return l_bin + getDisplacement(i_x);
else
return m_delta + getDisplacement(i_x);
}
tsunami_lab::t_real tsunami_lab::setups::TsunamiEvent1d::getBathymetryBin(t_real i_x) const
{
// convert i_x to cell index (assuming 250m cells)
int l_row = i_x / 250;
return io::Csv::readLine(m_doc, l_row);
}
3.4.4 Visualize the tsunami setup
The Tsunami Setup simulated over an hour of time. We observe a runup of roughly 1.8 to 2 meters.
3.4.5 Impact of different initial displacements
A version with twice the initial displacement. (instead of \(10 \sin(\frac{(x - 175000)}{37500}\pi+ \pi)\) we used \(20 \sin(\frac{(x - 175000)}{37500}\pi + \pi)\)). The momentum traveling to both sides of the simulations are roughly twice as high. Maybe a linear relationship between the initial displacement and the momentum is present?
A version with a static 20 meter displacement in between 175km and 225km and a left reflective boundary. It seems to travel as a single big wave towards japan mainland hitting it with roughly 15m height and getting reflected to roughly half the height. This is basically our dambreak problem in 2 directions without an infinite water source. The lower video is just the height of the water without the momentum and bathymetry.