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)

.. raw:: html

    <video width="100%" height="auto" controls>
      <source src="../../_static/assets/roe_solver_visualization.mp4" type="video/mp4">
    </video> 

1.3 - F-wave Solver
============================

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

For explanations on function inputs and outputs, see :ref:`ns-solvers`.  

computeEigenvalues()
^^^^^^^^^^^^^^^^^^^^^^^^^
First we compute the height and particle velocity

.. math::
    \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}

.. code-block:: cpp

  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

.. math::
    \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}

.. code-block:: cpp

  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

.. math::
    \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}

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

.. code-block:: cpp
      
  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, :math:`\Delta f := f(q_r) - f(q_l)`. 

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

.. code-block:: cpp

  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 :math:`\alpha`:

.. code-block:: cpp

  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:

.. math::
    \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}

.. code-block:: cpp

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

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

.. math:: 
   Z_1 = \alpha_1 r_1, Z_2 = \alpha_2 r_2

.. code-block:: cpp

  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

.. math::

  \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}

.. code-block:: cpp

  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];
  }