How to compute WENO-reconstructed flux in local characteristic field?

Question

I'm trying to apply finite-difference characteristic-wise WENO method to 1D Euler equations: $$ U_t + F(U)_x = 0, \; \;\text{where} \;\; U = (\rho, \rho u, e) \;\; \text{and} \;\;F(U) = (\rho u, \rho u^2+p, (e+p)u) $$ Physical values are set at grid points $x_j$, so we need to approximate fluxes at mid-points $x_{j\pm \frac{1}{2}}$: $$ U_t \approx \frac{1}{h} \left[ \widehat{f}_{j+\frac{1}{2}} - \widehat{f}_{j-\frac{1}{2}} \right] $$ where $\widehat{f}_{j\pm\frac{1}{2}}$ are sought-for approximation of physical fluxes $f(U(x_{j\pm\frac{1}{2}}))$.

My algorithm is the following:

1) For each $x_{j+\frac{1}{2}}$, we calculate simple average state $U_{j+\frac{1}{2}} = \frac{1}{2} (U_j + U_{j+1})$ and, using it, local eigenvalues $(u-a, u, u+a)_{j+\frac{1}{2}}$, right eigenvector matrix $R_{j+\frac{1}{2}}$ and it's left counterpart $R^{-1}_{j+\frac{1}{2}}$. (Subscript ${j+\frac{1}{2}}$ will be omitted below)

2) Now we transform $U$, its differences $\Delta U$ and flux differences $\Delta F(U)$ to local characteristic field: $$ W=R^{-1}U, \;\; \Delta W=R^{-1}\Delta U, \;\; \Delta F(W)=R^{-1}\Delta F(U) $$ It is done only for relevant grid points which in my case ($r=2$) are $x_{j-1}, x_j, x_{j+1}, x_{j+2}$.

3) Then we reconstruct characteristic variable values by WENO method using $W$, $\Delta W$ and $\Delta F(W)$. We get two "candidates": $$ \widehat{W}^- = R_j(x_{j+\frac{1}{2}}) \; \text{and} \; \widehat{W}^+ = R_{j+1}(x_{j+\frac{1}{2}}) $$ where $R_j(x), R_{j+1}(x)$ are a convex combinations of corresponding polynomial functions obtained by WENO method.

The next step would be to compute 2-point flux function $h(\widehat{W}^-,\widehat{W}^+)$ using these reconstructed values. I'm using Lax-Friedrich's function: $$ h(a,b) = \frac{1}{2} \left[ F(a) + F(b) + \alpha(b-a) \right] $$ where $\alpha$ is maximum value of eigenvalue (different for each of equations, since they are decoupled now).

But how to compute these $F(a)$ and $F(b)$ using $\widehat{W}^\pm$? In component-wise approach it is simple: we reconstruct base values $-$ $(\widehat{\rho}^\pm, \widehat{u}^\pm, \widehat{p}^\pm)$, for example $-$ and the use explicit formulas to compute $F(a) = (\widehat{\rho}^- \widehat{u}^-,\; \widehat{\rho}^-{\widehat{u}^-}^2 + \widehat{p}^-,\; ...)$ and $F(b) = (\widehat{\rho}^+ \widehat{u}^+,\; \widehat{\rho}^+{\widehat{u}^+}^2 + \widehat{p}^+,\; ...)$, but we don't have such formulas for $F(\widehat{W})$ because $\widehat{W}$ is expressed in local characteristic field...

Answer 1

Although @omican has found the answer, if one HAS to proceed with the characteristic decoupling of the conservative variable vector only (and not the flux vector) as asked in the original question, here is how it can be done:

1. Instead of the algebraic averages as shown in point $1$ of the question, find the 'Roe linearized state' $U_{j+\frac{1}{2}} = Roe(U_j, U_{j+1})$, where $Roe()$ is the function returning Roe averaged quantities. Any standard book on CFD can be referred for this. Then proceed to finding eigensystem.

2. Get $\widehat{W}^- = \left[ \widehat{W_1}^- \widehat{W_2}^- \widehat{W_3}^- \right] ^ T$, and $\widehat{W}^+ = \left[ \widehat{W_1}^+ \widehat{W_2}^+ \widehat{W_3}^+ \right] ^ T$. This remains the same as described in points $2$ and $3$ in the question.

3. Now, based on the eigenvalues of the linearized (averaged) Jacobian (remember to perform Roe linearization and not take the algebraic average), construct a vector of characteristic variables at the face $\left(j+\frac{1}{2}\right)$. i.e. Consider the system to be subsonic, so that one of the eigenvalues is negative (let it be the first eigenvalue) and the remaining two are positive. So, the assembled vector of characteristic variables reads as, $\widehat{W}^* = \left[ \widehat{W_1}^+ \widehat{W_2}^- \widehat{W_3}^- \right] ^ T$. Here, owing to the wave propagation direction, the first characteristic field is advected from right to left (-ve eigenvalue) and remaining to waves from left to right.

4. Now find the vector of conserved variables from $\widehat{W}^*$, i.e. $U^* = R (\widehat{W}^*)$ where $R$ is the matrix of right-eigenvectors.

5. The numerical flux function $\widehat{f}_{j+\frac{1}{2}}$ is simply $f(U^*)$

This is another way of finding the numerical flux function . But this requires entropy fix since the numerical dissipation at the contact points reduces to zero.

However, the answer posted by @omican is correct if one wants to go with the Local Lax Friedrich's flux function.

Answer 2

It seems that I was misunderstanding how the flux is reconstructed:

My idea was to reconstruct base values $-$ $U$ (or $W$ if we are performing characteristic decomposition) $-$ and then substitute them to formulas of $F(U)$ (or $F(W)$) to obtain reconstructed fluxes.

It looks like the correct way is to apply reconstruction procedure to $F(U)$ (or $F(W)$ which is obtained by transformation $F(W) = R^{-1}F(U)$) itself. Some quick tests I performed with this approach show very plausible results.