# How to implement Lax-Friedrich flux splitting with WENO scheme

I'm working on modeling a shock wave using the Euler equation with an advanced Equation of state and the fifth order WENO scheme. The equation are set up on the form: $$\begin{equation} \frac{\partial U}{\partial t} + \frac{\partial F(U)}{\partial x} = 0 \end{equation}$$ where $$\begin{equation} U = \begin{bmatrix}\rho\\ \rho v\\E\end{bmatrix}, F(U) = \begin{bmatrix}\rho v\\ \rho v^2 + p\\v(E+p)\end{bmatrix} \end{equation}$$ I've formulated the WENO scheme, but have yet to get it to work. What I find difficult to implement is the Lax-Friedrich Flux splitting. Specifically determining the $$\alpha$$ in the equation:

$$\begin{equation} f^{\pm} = f(U) \pm \alpha U \end{equation}$$

I do know that $$\alpha = \max \left|\frac{\partial f(U)}{\partial U}\right|$$, but for what index range of $$U$$ does the equation strech when i want to determine:

$$\begin{equation} \hat{f}_{i+1/2} = \hat{f}_{i+1/2}^+ + \hat{f}_{i+1/2}^- \end{equation}$$

• The $\hat{f}$ denotes the numerical formulation. In regards to the Equation of state, it is a somewhat complicated EOS adapted to the LJS fluid. I utilize a function which computes the temperature and pressure based on the internal energy and density. This does makes it difficult to determine analytical eigenvalues, compared to the ideal case. Feb 15, 2022 at 10:50
• I solve the temporal discretization using 5th explicit Runge Kutta with the programmingfunction solve_ivp() in python. What I specifically want, is to compute the $f^{\pm}$ in order to utilize the WENO scheme, but descriptions in litterature of $\alpha$ doesn't indicate the range of $U$. The spatial formulation is correct. Feb 15, 2022 at 12:28
• Ah got it, so your semi-discrete formulation reads $$\frac{\mathrm d \boldsymbol{u}} {\mathrm d t } = \frac{1}{\Delta x} \Big( \hat{ \boldsymbol{f}}_{i-1/2} - \hat{ \boldsymbol{f}}_{i+1/2} \Big)$$? Feb 15, 2022 at 12:35
• That is correct, where the $\hat{f}_{i\pm1/2}$ values are computed using the WENO scheme, with the splitted fluxes ($\hat{f}^{\pm}$). Feb 15, 2022 at 12:41

I think you mean with Lax-Friedrichs the local Lax-Friedrichs or Rusanov scheme ,  where the wave speeds (left & right-going) are given by the maximum eigenvalue $$\lambda_\max := \max_i | \lambda_i(\boldsymbol{U}) | \in \sigma\Big( \boldsymbol{F}'(\boldsymbol{U}) \Big)$$ where $$\sigma$$ denotes the spectrum of the Jacobian $$\boldsymbol{F}'(\boldsymbol{U})$$.
To finally answer your question, you need to take only the trace values $$\boldsymbol{U}^-, \boldsymbol{U}^+$$ into account at a certain face. Then,
$$\alpha := \max \bigg\{ \max_i \big \vert \lambda_i (\boldsymbol{U}^-) \big \vert \in \sigma\Big( \boldsymbol{F}'(\boldsymbol{U}) \Big) \Big \vert_{\boldsymbol{U}^-} , \max_i \big \vert \lambda_i (\boldsymbol{U}^+) \big \vert \in \sigma\Big( \boldsymbol{F}'(\boldsymbol{U}) \Big) \Big \vert_{\boldsymbol{U}^+} \bigg\}.$$ In other words at each face $$i$$ you jest need the trace values $$\boldsymbol{U}_i^+, \boldsymbol{U}_i^-$$. In that case, your numerical flux (I think you lack the factor if $$1/2$$ in your formulation and confused the signs) is given by \begin{align} \hat {\boldsymbol{f}}_{i+1/2}& = \frac{1}{2} \Big( \hat {\boldsymbol{f}}_{i+1/2}^+ + \hat {\boldsymbol{f}}_{i+1/2}^- \Big) =\frac{1}{2} \Big( F(\boldsymbol{U}_i^+) - \alpha \boldsymbol{U}_i^+ + F(\boldsymbol{U}_i^-) + \alpha \boldsymbol{U}_i^- \Big) \\ &= \frac{1}{2} \Big( F(\boldsymbol{U}_i^+) + F(\boldsymbol{U}_i^-) + \alpha (\boldsymbol{U}_i^- - \boldsymbol{U}_i^+ \Big) \end{align}
• Thank you! A couple of follow-up questions: So by interface of an index $i$ I only take the maximum approximate eigenvalue between the $i \pm 1/2$ values? Do you have any advice for computing the jacobian i.e. approximate values for values connected to the EOS? Feb 15, 2022 at 14:30
• More precise, you take the maximum of the absolute values of the eigenvalues of the Jacobian $\frac{\mathrm d \boldsymbol{F}}{\mathrm d \boldsymbol{U}} = F'(\boldsymbol{U})$ which is evaluated at both sides of the interface $i+ 1/2$. Based on your reconstruction / interpolation method you have the value $\boldsymbol{U}_{i+1/2}^-$ on the left of the $i+1/2$ cell boundary and $\boldsymbol{U}_{i+1/2}^+$ on the right. The Jacobian / eigenvalues (if analytically known) are evaluated at this states. Feb 16, 2022 at 7:57