3. Finite Volume Discretization
Finite Volume Discretization is a powerful and widely used numerical technique in the field of computational fluid dynamics (CFD) and other scientific and engineering disciplines. It plays a crucial role in the simulation and analysis of physical processes, particularly those governed by partial differential equations (PDEs)
3.1. Integrating our f-wave solver into the one-dimensional implementation of a wave propagation
There are several ways to integrate our F-wave solver into the one-dimensional wave propagation implementation, including the following options:
Using the command prompt.
Utilizing a configuration file, and so on.
Today, we will focus on using the command prompt to integrate our solver with the following files: WavePropagation1d.h, WavePropagation1d.cpp, WavePropagation1d.test.cpp. Let’s proceed step by step:
3.1.1. Integration of the F-wave solver
In the
WavePropagation1d.h, we will define a boolean member to choose the two solver. Update the code in “WavePropagation1d.h” as follows:
class tsunami_lab::patches::WavePropagation1d: public WavePropagation {
private:
bool m_choice = true;
}
We will modify the constructor in the
WavePropagation1d.cppfile to include the boolean variable. This change enables us to select between the Roe solver and our F-wave solver within the test files and the command prompt.
tsunami_lab::patches::WavePropagation1d::WavePropagation1d( t_idx i_nCells,bool i_choice ) {
m_choice = i_choice;
m_nCells = i_nCells;
// allocate memory including a single ghost cell on each side
for( unsigned short l_st = 0; l_st < 2; l_st++ ) {
m_h[l_st] = new t_real[ m_nCells + 2 ];
m_hu[l_st] = new t_real[ m_nCells + 2 ];
}
// init to zero
for( unsigned short l_st = 0; l_st < 2; l_st++ ) {
for( t_idx l_ce = 0; l_ce < m_nCells; l_ce++ ) {
m_h[l_st][l_ce] = 0;
m_hu[l_st][l_ce] = 0;
}
}
}
In the same file, we will navigate to the
timeStepfunction and add a call tonetupdatesbased on the user’s choice..
void tsunami_lab::patches::WavePropagation1d::timeStep( t_real i_scaling ) {
// compute net-updates
if(m_choice){
solvers::Roe::netUpdates( l_hOld[l_ceL],
l_hOld[l_ceR],
l_huOld[l_ceL],
l_huOld[l_ceR],
l_netUpdates[0],
l_netUpdates[1] );
}else{
solvers::fwave::netUpdates( l_hOld[l_ceL],
l_hOld[l_ceR],
l_huOld[l_ceL],
l_huOld[l_ceR],
l_netUpdates[0],
l_netUpdates[1]);
}
}
To conclude, we navigate to the “main.cpp” file, where we need to modify the solver constructor to permit the selection between our F-wave and Roe solver options.
// construct solver
tsunami_lab::patches::WavePropagation *l_waveProp;
if (i_argc > 2) {
for (char* ptr = i_argv[2]; *ptr; ++ptr) {
*ptr = std::tolower(static_cast<unsigned char>(*ptr));
}
if (std::strstr(i_argv[2], "roe") != nullptr) {
std::cout << "ROE Will be Used" << std::endl;
l_waveProp = new tsunami_lab::patches::WavePropagation1d( l_nx , true );
} else {
std::cout << "FWAVE Will be Used" << std::endl;
l_waveProp = new tsunami_lab::patches::WavePropagation1d( l_nx , false );
}
}else{
std::cout << "FWAVE Will be Used" << std::endl;
l_waveProp = new tsunami_lab::patches::WavePropagation1d( l_nx , false );
}
Important
To select between the F-wave and Roe solver, we execute the following command.
./build/tsunami_lab 200 Roe
or
./build/tsunami_lab 200 fwave
if the user does not choose any solver, the F-wave solver will be choosen automatically.
To create a unit test for our solvers, we need to include a
trueorfalseargument in the constructor of the unit test:TEST_CASE( "Test the 1d wave propagation solver.", "[WaveProp1d]" ) { // constructs the 1d wave propagation solver. tsunami_lab::patches::WavePropagation1d m_waveProp( 100 , true ); ... }
Important
truewill use theRoe solverandfalsewill use theFwave solver
3.1.2. Constant middle states
Now, we are going to check for some middle states. This file, middle_states.csv , contains a collection of constant middle states that arise immediately in the Riemann solution at the initial discontinuity. We will use it as a sanity test. Let’s implement a unit test for the middle state.”
We access a CSV file and search for suitable values to use in a unit test. In this example, we will select the initial values from the file.
Note
hLeft: 8899.326826472694,hRight: 8899.326826472694,huLeft: 122.0337839252433,huRight: -122.0337839252433,hStar: 8899.739847378269Important
hStarrepresents the height in the middle state.
we construct a test unit like that
TEST_CASE( "Test the 1d wave propagation wave for the middle states.", "[middle states]" ){ tsunami_lab::patches::WavePropagation1d middle_states1(100, true); for (std::size_t l_ce = 0; l_ce < 50; l_ce++) { middle_states1.setHeight(l_ce, 0, 8899.326826472694); middle_states1.setMomentumX(l_ce, 0, 122.0337839252433); } for (std::size_t l_ce = 50; l_ce < 100; l_ce++) { middle_states1.setHeight(l_ce, 0, 8899.326826472694); middle_states1.setMomentumX(l_ce, 0, -122.0337839252433); } middle_states1.setGhostOutflow(); for (int i = 0; i < 17; i++) { middle_states1.timeStep(0.001); } REQUIRE(middle_states1.getHeight()[49] == Approx(8899.74)); REQUIRE(middle_states1.getHeight()[50] == Approx(8899.74)); }
Important
We will observe that the height after 17 time steps is approaching the initial height of 8899.326826.
3.1.3. GitHub Actions
We’ve established a workflow within GitHub that currently operates as a self-hosted background process. We’ve implemented a “tests-runner.yml” file to automatically test everything following each push and pull request.
3.2. Shock and Rarefaction Waves
Next, we will delve into the setup implementation of both the shock shock problem and the rare_reare problem:
3.2.1. Shock-Shock problem
We’ll apply our solver to solve situations where two streams of water move in opposite directions and crash into each other \(x_\text{dis}\) . The situation is defined by the following configuration:
\[\begin{split}\begin{cases} Q_i = q_{l} \quad &\text{if } x_i \le x_\text{dis} \\ Q_i = q_{r} \quad &\text{if } x_i > x_\text{dis} \end{cases} \qquad q_l \in \mathbb{R}^+ \times \mathbb{R}^+, \; q_r \in \mathbb{R}^+ \times \mathbb{R}^-,\end{split}\]
Before we begin, let’s navigate to the setup folder and generate three files:
ShockShock.cpp,ShockShock.hand for testing purposes, createShockShock.test.cpp.Now, let’s start implementing each file, beginning with
ShockShock.h.
#ifndef TSUNAMI_LAB_SETUPS_SHOCK_SHOCK_H
#define TSUNAMI_LAB_SETUPS_SHOCK_SHOCK_H
#include "Setup.h"
namespace tsunami_lab {
namespace setups {
class ShockShock;
}
}
class tsunami_lab::setups::ShockShock: public Setup {
private:
//! height
t_real m_height = 0;
//! impulse
t_real m_hu = 0;
//! location of the dam
t_real m_locationDam = 0;
public:
/**
* Constructor.
*
* @param i_height water height
* @param i_hu water impulse
* @param i_locationDam location (x-coordinate) of the dam.
**/
ShockShock( t_real i_height,
t_real i_hu,
t_real i_locationDam);
/**
* Gets the water height at a given point.
*
* @return height at the given point.
**/
t_real getHeight( t_real,
t_real) const;
/**
* Gets the momentum in x-direction.
* @param i_x x-coordinate of the queried point.
* @return momentum in x-direction.
**/
t_real getMomentumX( t_real i_x,
t_real ) const;
/**
* Gets the momentum in y-direction.
* @return momentum in y-direction.
**/
t_real getMomentumY( t_real,
t_real ) const;
};
We will implement the configuration for the
Shock Shock Problemin theShockShock.hfile :
#include "ShockShock.h"
tsunami_lab::setups::ShockShock::ShockShock(t_real i_height,
t_real i_hu,
t_real i_locationDam){
m_height = i_height;
m_hu = i_hu;
m_locationDam = i_locationDam;
}
tsunami_lab::t_real tsunami_lab::setups::ShockShock::getHeight( t_real,
t_real)const{
return m_height;
}
tsunami_lab::t_real tsunami_lab::setups::ShockShock::getMomentumX(t_real i_x,
t_real)const{
if( i_x <= m_locationDam ) {
return m_hu;
}
else {
return -m_hu;
}
}
tsunami_lab::t_real tsunami_lab::setups::ShockShock::getMomentumY(t_real,
t_real)const{
return 0;
}
Lastly, to verify that the implemented setup functions correctly, we’ll implement the
ShockShock.test.cppfile :
#include <catch2/catch.hpp>
#include "ShockShock.h"
TEST_CASE( "Test the shock-shock setup.", "[ShockShock]" ) {
tsunami_lab::setups::ShockShock l_shockShock( 25, 55, 3 );
// left side
REQUIRE( l_shockShock.getHeight( 2, 0 ) == 25 );
REQUIRE( l_shockShock.getMomentumX( 2, 0 ) == 55 );
REQUIRE( l_shockShock.getMomentumY( 2, 0 ) == 0 );
REQUIRE( l_shockShock.getHeight( 3, 0 ) == 25 );
REQUIRE( l_shockShock.getMomentumX( 3, 0 ) == 55 );
REQUIRE( l_shockShock.getMomentumY( 3, 0 ) == 0 );
// right side
REQUIRE( l_shockShock.getHeight( 5, 0 ) == 25 );
REQUIRE( l_shockShock.getMomentumX( 5, 0 ) == -55 );
REQUIRE( l_shockShock.getMomentumY( 5, 0 ) == 0 );
REQUIRE( l_shockShock.getHeight( 60, 0 ) == 25 );
REQUIRE( l_shockShock.getMomentumX( 60, 0 ) == -55 );
REQUIRE( l_shockShock.getMomentumY( 60, 0 ) == 0 );
}
3.2.2. Rare Rare Porblems
We can configure rare-rare Riemann problems with two streams of water moving away from each other at a specific position, denoted as \(x_\text{dis}\) . The scenario is outlined as follows:
and
Before we begin, let’s navigate to the setup folder again and generate three files:
RareRare.cpp,RareRare.hand for testing purposes, createRareRare.test.cpp.Now, let’s start implementing each file, beginning with
RareRare.h:
#ifndef TSUNAMI_LAB_SETUPS_Rare_Rare_H
#define TSUNAMI_LAB_SETUPS_Rare_Rare_H
#include "Setup.h"
namespace tsunami_lab {
namespace setups {
class RareRare;
}
}
class tsunami_lab::setups::RareRare: public Setup {
private:
//! height
t_real m_height = 0;
//! impulse
t_real m_hu = 0;
//! location of the dam
t_real m_locationDam = 0;
public:
/**
* Constructor.
*
* @param i_height water height
* @param i_hu water impulse
* @param i_locationDam location (x-coordinate) of the dam.
**/
RareRare( t_real i_height,
t_real i_hu,
t_real i_locationDam);
/**
* Gets the water height at a given point.
* @return height at the given point.
**/
t_real getHeight( t_real ,
t_real ) const;
/**
* Gets the momentum in x-direction.
* @param i_x x-coordinate of the queried point.
* @return momentum in x-direction.
**/
t_real getMomentumX( t_real i_x,
t_real ) const;
/**
* Gets the momentum in y-direction.
* @return momentum in y-direction.
**/
t_real getMomentumY( t_real,
t_real ) const;
};
We will implement the configuration for the
Rare Rare Problemin the :
#include "RareRare.h"
tsunami_lab::setups::RareRare::RareRare(t_real i_height,
t_real i_hu,
t_real i_locationDam){
m_height = i_height;
m_hu = i_hu;
m_locationDam = i_locationDam;
}
tsunami_lab::t_real tsunami_lab::setups::RareRare::getHeight(t_real,
t_real)const{
return m_height;
}
tsunami_lab::t_real tsunami_lab::setups::RareRare::getMomentumX(t_real i_x,
t_real)const{
if( i_x <= m_locationDam ) {
return -m_hu;
}
else {
return m_hu;
}
}
tsunami_lab::t_real tsunami_lab::setups::RareRare::getMomentumY(t_real,
t_real)const{
return 0;
}
Lastly, to verify that the implemented setup functions correctly, we’ll implement the
RareRare.test.cppfile :
#include <catch2/catch.hpp>
#include "RareRare.h"
TEST_CASE( "Test the rare-rare setup.", "[RareRare]" ) {
tsunami_lab::setups::RareRare l_rareRare( 25,55,3);
// left side
REQUIRE( l_rareRare.getHeight( 2, 0 ) == 25 );
REQUIRE( l_rareRare.getMomentumX( 2, 0 ) == -55 );
REQUIRE( l_rareRare.getMomentumY( 2, 0 ) == 0 );
REQUIRE( l_rareRare.getHeight( 3, 0 ) == 25 );
REQUIRE( l_rareRare.getMomentumX( 3, 0 ) == -55 );
REQUIRE( l_rareRare.getMomentumY( 3, 0 ) == 0 );
// right side
REQUIRE( l_rareRare.getHeight( 5, 0 ) == 25 );
REQUIRE( l_rareRare.getMomentumX( 5, 0 ) == 55 );
REQUIRE( l_rareRare.getMomentumY( 5, 0 ) == 0 );
REQUIRE( l_rareRare.getHeight( 60, 0 ) == 25 );
REQUIRE( l_rareRare.getMomentumX( 60, 0 ) == 55 );
REQUIRE( l_rareRare.getMomentumY( 60, 0 ) == 0 );
}
3.2.3. The influence of \(u_l\) and \(h_l\) on the wave speed in both shock-shock and rare-rare setup
After experimenting with different sets of initial water heights \(h_l\) and particle velocities \(u_l\) . we observed that the primary variable affecting wave velocities is the initial height. Both \(u_l\) and \(h_l\) exert an influence on momentum. Specifically, as \(h_l`and :math:`u_l\) increase in magnitude, the momentum also increases, leading to higher middle state heights. However, the changes in particle velocities do not affect the wave speed; only the initial height has an impact on both \(\lambda_{1}\) and \(\lambda_{2}\). This relationship becomes evident through our Geogebra calculations, which demonstrate that only \(h_l\) affect the wave speed. This connection is further illuminated by the following equations for our approximate Riemann solver.
And this relationship is also evident through our generated CSV files.
1. Shock-Shock setup:
\(h_l\) = 10 , \(u_l\) = 15 :
\(h_l\) = 20 , \(u_l\) = 15 :
\(h_l\) = 40 , \(u_l\) = 15 :
\(h_l\) = 10 , \(u_l\) = 20 :
\(h_l\) = 20 , \(u_l\) = 20 :
2. Shock-Shock setup:
\(h_l\) = 10 , \(u_l\) = 5 :
\(h_l\) = 10 , \(u_l\) = 8 :
\(h_l\) = 20 , \(u_l\) = 5 :
It can be observed that simulations featuring greater initial height will consistently exhibit higher wave speeds for both setups.
3. Geogebra calculations:
both setup are outlined as follows:
so we do the calculations:
\(h_l\) = \(h_r\) , \(-(hu_l) =\) = \(hu_r\)
if \(h_l\) = \(h_r\) , \(hu_l =\) = \(-(hu_r)\) then \(u^{\text{Roe}}(q_l, q_r)\) will always equal zero.
This implies that \(\lambda_{1/2} = \mp \sqrt{gh}\)
For example, consider the situation illustrated in the picture.
3.3. Dam-Break
3.3.1. Evaluating how various factors influence the configuration of a dam-break scenario
To begin, we will explore a range of initial water height scenarios, aiming to formulate a conclusion based on our experimentation:
\(h_l\) = 14 , \(h_r\) = 4, \(hu\) = 0 :
\(h_l\) = 14 , \(h_r\) = 8, \(hu\) = 0 :
\(h_l\) = 50 , \(h_r\) = 12, \(hu\) = 0 :
\(h_l\) = 50 , \(h_r\) = 18, \(hu\) = 0 :
To understand the significance of particle velocity, it is necessary to make a temporary modification to the getMomentumX() function.
Note: \(hu_r = h_r * u_r\) so we need to modify the getMomentumX to return either hu_ur or hu_ul
Now, let’s conduct a few tests to further enhance our conclusion.
\(h_l\) = 18 , \(h_r\) = 5, \(hu_l\) = 0, \(hu_r\) = 0.5 :
\(h_l\) = 40 , \(h_r\) = 10, \(hu_l\) = 10, \(hu_r\) = 14 :
So we noticed these points:
The greater the initial height difference, the more significant the momentum observed within the middle state.
The speed of the shock wave seems largely unaffected by variations in the right height, whereas the rarefaction wave demonstrates reduced speed with increasing initial height differences.
The speed of the shock wave appears to be directly proportional to the square root of the left height.
Although a higher momentum on the right side does result in a marginal increase in the speed of the shock wave, the effect is generally minor.
3.3.2. Evacuation time
The provided information consists of the following data: \(q_l=[14, 0]^T\) , \(q_r=[3.5, 0.7]^T\) and \(s_{\text{village}}= 25 km\)
Now, let’s calculate the wave speed using this information.
so for \(\lambda_{1/2}\) we get \(\lambda_{1}= 8.71 \frac{m}{s}\) and \(\lambda_{2}= -7.97 \frac{m}{s}\).You can determine the time it takes for the wave to reach the village by using the formula: \(time= \frac{distance}{wave speed}\).
so we get :
\(time_{\text{evacuation}} = \frac{25000 m}{8.71 \frac{m}{s}} = 2870.26 s = 0.797 h = 47 min\)
This can also be calculated through a simulation:
The simulation has 118 CSV solutions, and the shock wave reaches 9.79m, which means it has moved around 5.14m. Therefore, we have \(t = \frac{41}{118} * 1,5s = 0.52 s\). the speed of the shock wave is given by: \(v_{\text{shockwave}} = \frac{5.14m}{0.52s} = 9,88 \frac{m}{s} = 35.56 \frac{km}{h}\). Hence, the evacuation time is: \(time_{\text{evacuation}} = \frac{25 km}{35.56 \frac{km}{h}} = 0.70h = 42 min\)
3.4. Personal Contribution
Ward Tammaa, Daniel Schicker Doxygen Documentation
Mohamad Khaled Minawe, Ward Tammaa 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(Calculation for the Unit Tests)
Ward Tammaa Hosting the code, Action runner