1. Riemann Solver

All project authors contributed to this assignment in equal parts.

1.2 - Getting Started

1.2.2 - Visualizations

Roe solver with input 500 (cells in x direction)

1.3 - F-wave Solver

File: Fwave.cpp (and Fwave.h)

For explanations on function inputs and outputs, see Solvers.

computeEigenvalues()

First we compute the height and particle velocity

\[\begin{split}\begin{split}\begin{aligned} h^{\text{Roe}}(q_l, q_r) &= \frac{1}{2} (h_l + h_r), \\ u^{\text{Roe}}(q_l, q_r) &= \frac{u_l \sqrt{h_l} + u_r \sqrt{h_r}}{\sqrt{h_l}+\sqrt{h_r}}. \end{aligned}\end{split}\end{split}\]

t_real l_hRoe = t_real(0.5) * (i_hL + i_hR);
t_real l_uRoe = l_hSqrtL * i_uL + l_hSqrtR * i_uR;
l_uRoe /= l_hSqrtL + l_hSqrtR;

Using those values, we can then compute the eigenvalues

\[\begin{split}\begin{split}\begin{aligned} \lambda^{\text{Roe}}_{1}(q_l, q_r) &= u^{\text{Roe}}(q_l, q_r) - \sqrt{gh^{\text{Roe}}(q_l, q_r)}, \\ \lambda^{\text{Roe}}_{2}(q_l, q_r) &= u^{\text{Roe}}(q_l, q_r) + \sqrt{gh^{\text{Roe}}(q_l, q_r)}, \end{aligned}\end{split}\end{split}\]

t_real l_ghSqrtRoe = m_gSqrt * std::sqrt(l_hRoe);
eigenvalueRoe_1 = l_uRoe - l_ghSqrtRoe;
eigenvalueRoe_2 = l_uRoe + l_ghSqrtRoe;

Note

m_gSqrt is the given square root of gravity.

computeEigencoefficients()

The task of this function is to compute

\[\begin{split}\begin{split}\begin{bmatrix} \alpha_1 \\ \alpha_2 \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ \lambda^{\text{Roe}}_1 & \lambda^{\text{Roe}}_2 \end{bmatrix}^{-1} \Delta f.\end{split}\end{split}\]

In order to do so, we start off by inverting the right eigenvector-matrix

t_real l_detInv = 1 / (eigenvalueRoe_2 - eigenvalueRoe_1);
t_real l_rInv[2][2] = {{0}};
l_rInv[0][0] = l_detInv * eigenvalueRoe_2;
l_rInv[0][1] = -l_detInv;
l_rInv[1][0] = -l_detInv * eigenvalueRoe_1;
l_rInv[1][1] = l_detInv;

Next, we need to calculate the jump in the flux function f, \(\Delta f := f(q_r) - f(q_l)\).

Note

Remember that \(f := [hu, hu^2 + \frac{1}{2}gh^2]^T\).

t_real f_delta[2] = {0};
f_delta[0] = i_huR - i_huL;
f_delta[1] = (i_huR * l_uR + t_real(0.5) * m_g * i_hR * i_hR) - (i_huL * l_uL + t_real(0.5) * m_g * i_hL * i_hL);

Finally, we can derive the desired output vector \(\alpha\):

alpha_1 = l_rInv[0][0] * f_delta[0] + l_rInv[0][1] * f_delta[1];
alpha_2 = l_rInv[1][0] * f_delta[0] + l_rInv[1][1] * f_delta[1];

netUpdates()

With the help of the eigenvalues, we can derive the eigenvectors:

\[\begin{split}\begin{split}\begin{aligned} r_1^{\text{Roe}} &= \begin{bmatrix} 1 \\ \lambda^{\text{Roe}}_1 \end{bmatrix}, \\ r_2^{\text{Roe}} &= \begin{bmatrix} 1 \\ \lambda^{\text{Roe}}_2 \end{bmatrix}. \end{aligned}\end{split}\end{split}\]

t_real eigenvectorRoe_1[2] = {1, eigenvalueRoe_1};
t_real eigenvectorRoe_2[2] = {1, eigenvalueRoe_2};

Now that we have the eigencoefficients \(\alpha_{1/2}\) and eigenvectors \(r_{1/2}\), we can compute the waves \(Z_{1/2}\):

\[Z_1 = \alpha_1 r_1, Z_2 = \alpha_2 r_2\]

t_real z1[2] = {0};
z1[0] = eigencoefficientRoe_1 * eigenvectorRoe_1[0];
z1[1] = eigencoefficientRoe_1 * eigenvectorRoe_1[1];

t_real z2[2] = {0};
z2[0] = eigencoefficientRoe_2 * eigenvectorRoe_2[0];
z2[1] = eigencoefficientRoe_2 * eigenvectorRoe_2[1];

All that is left to do is to set the net-updates depending on the wave speeds

\[\begin{split}\begin{split}\begin{split} A^- \Delta Q := \sum_{p:\{ \lambda_p^\text{Roe} < 0 \}} Z_p \\ A^+ \Delta Q := \sum_{p:\{ \lambda_p^\text{Roe} > 0 \}} Z_p \end{split}\end{split}\end{split}\]

for (unsigned short l_qt = 0; l_qt < 2; l_qt++)
{
  //init
  o_netUpdateL[l_qt] = 0;
  o_netUpdateR[l_qt] = 0;

  //wave 1
  if (eigenvalueRoe_1 < 0) o_netUpdateL[l_qt] += z1[l_qt];
  else o_netUpdateR[l_qt] += z1[l_qt];

  //wave 2
  if (eigenvalueRoe_2 < 0) o_netUpdateL[l_qt] += z2[l_qt];
  else o_netUpdateR[l_qt] += z2[l_qt];
}