1
$\begingroup$

After successful implementation of characteristic-wise finite-difference WENO method to 1D Euler equations, I'm moving to 2D equations on cartesian grid: $$ \frac{\partial U}{\partial t} + \frac{\partial F(U)}{\partial x} + \frac{\partial G(U)}{\partial y} = 0 $$ where $$ U = \left( \begin{array}{c} \rho \\ \rho u \\ \rho v \\ e\end{array} \right), \quad F(U) = \left( \begin{array}{c} \rho u \\ \rho u^2 + p \\ \rho uv \\ (e+p)u\end{array} \right), \quad G(U) = \left( \begin{array}{c} \rho v \\ \rho uv \\ \rho v^2 + p \\ (e+p)v\end{array} \right), \quad $$

According to Shu's guidebook, numerical method for 2D case is quite simply derived from 1D case: $$ \frac{dU_{ij}}{dt} = - \frac{1}{\Delta x} \left( \hat{f}_{i+\frac{1}{2}, j} - \hat{f}_{i-\frac{1}{2}, j} \right) - \frac{1}{\Delta y} \left( \hat{g}_{i, j+\frac{1}{2}} - \hat{g}_{i, j-\frac{1}{2}} \right) $$ where $\hat{f},\hat{g}$ are found by 1D WENO method using corresponding eigenmatrices for $F(U)$ and $G(U)$.

The matrices I found in literature are (notation used: $q = \frac{u^2+v^2}{2},\; a=\sqrt{\frac{\gamma p}{\rho}}, \; h=\frac{a^2}{\gamma-1} + q$): $$ J_f = \frac{\partial F(U)}{\partial U} = \left( \begin{array}{cccc} 0 & 1 & 0 & 0\\ (\gamma-1)q-u^2 & (3-\gamma)u & (1-\gamma)v & \gamma-1\\ -uv & v & u & 0\\ \Big((\gamma-1)q-h\Big)u & h-(\gamma-1)u^2 & (1-\gamma)uv & \gamma u \end{array}\right) $$ $$ R_f = \left( \begin{array}{cccc} 1 & 1 & 1 & 0\\ u-a & u & u+a & 0\\ v & v & v & -1\\ h-au & q & h+au & -v \end{array}\right), \quad \lambda_f=(u-a,\;u,\;u+a,\;u) $$ $$ L_f = R^{-1}_f = \left( \begin{array}{cccc} \frac{(\gamma-1)q+au}{2a^2} & \frac{(1-\gamma)u-a}{2a^2} & \frac{(1-\gamma)v}{2a^2} & \frac{(\gamma-1)}{2a^2}\\ \frac{a^2-(\gamma-1)q}{a^2} & \frac{(\gamma-1)u}{a^2} & \frac{(\gamma-1)v}{a^2} & \frac{(1-\gamma)}{a^2}\\ \frac{(\gamma-1)q-au}{2a^2} & \frac{(1-\gamma)u+a}{2a^2} & \frac{(1-\gamma)v}{2a^2} & \frac{(\gamma-1)}{2a^2}\\ v & 0 & -1 & 0 \end{array}\right) $$ These are used to calculate fluxes over x-direction ($\hat{f}$).

Also: $$ J_g = \frac{\partial G(U)}{\partial U} = \left( \begin{array}{cccc} 0 & 0 & 1 & 0\\ -uv & v & u & 0\\ (\gamma-1)q-v^2 & (1-\gamma)u & (3-\gamma)v & \gamma-1\\ \Big((\gamma-1)q-h\Big)v & (1-\gamma)uv & h-(\gamma-1)v^2 & \gamma v \end{array}\right) $$ $$ R_g = \left( \begin{array}{cccc} 1 & 1 & 1 & 0\\ u & u & u & 1\\ v-a & v & v+a & 0\\ h-av & q & h+av & u \end{array}\right), \quad \lambda_g=(v-a,\;v,\;v+a,\;v) $$ $$ L_g = R^{-1}_g = \left( \begin{array}{cccc} \frac{(\gamma-1)q+av}{2a^2} & \frac{(1-\gamma)u}{2a^2} & \frac{(1-\gamma)v-a}{2a^2} & \frac{(\gamma-1)}{2a^2}\\ \frac{a^2-(\gamma-1)q}{a^2} & \frac{(\gamma-1)u}{a^2} & \frac{(\gamma-1)v}{a^2} & \frac{(1-\gamma)}{a^2}\\ \frac{(\gamma-1)q-av}{2a^2} & \frac{(1-\gamma)u}{2a^2} & \frac{(1-\gamma)v+a}{2a^2} & \frac{(\gamma-1)}{2a^2}\\ -u & 1 & 0 & 0 \end{array}\right) $$ These are used to calculate fluxes over y-direction ($\hat{g}$).

The correctness of the matrices have been tested numerically: they satisfy $ LR=E$, $LJR = \Lambda = \text{diag}(\lambda_i)$ for both $F$ and $G$.

Obviously, the code for $x$-wise and $y$-wise WENO reconstruction is very symmetrical: just substitute correct matrices and use array differences along [j] index instead of [i].

Here is my problem: when I test the program on simple stationary shock normal to $x$-axis, everything is OK. But for similar test for $y$-axis, some instability grows in the uniform-state area ahead of the shock. I still can't find any asymmetry in my code; is it possible that I'm using matrices incorrectly? Maybe the order of eigenvalues and eigenvectors should be different $(\lambda_y=(v-a,\; v,\; v,\; v+a)$, for example)?

$\endgroup$
2
  • 2
    $\begingroup$ I'm voting to close this question as off-topic because I think nobody will be able to benefit from this as the answer was an "undescribed bug" $\endgroup$ – Anton Menshov Oct 29 '18 at 15:30
  • $\begingroup$ Will the question itself still be visible and thus "googlable"? If so, I can't have any objections to closing the question. $\endgroup$ – omican Oct 30 '18 at 16:11
0
$\begingroup$

After 10th or 20th code sweep I finally found the bug in it. Quick tests show that everything is OK now, method works nice. Let this question hang in here for as a help for future googlers of 2D Euler decomposition.

$\endgroup$
6
  • $\begingroup$ Appreciated, I'm just learning finite volume methods myself. $\endgroup$ – AtmosphericPrisonEscape Oct 25 '17 at 23:09
  • $\begingroup$ Does it work with bugs? $\endgroup$ – nicoguaro Dec 29 '17 at 16:29
  • $\begingroup$ @nicoguaro the bug was in code, not in formulas here - they are correct. $\endgroup$ – omican Dec 30 '17 at 17:49
  • $\begingroup$ So, why is this an useful answer? $\endgroup$ – nicoguaro Dec 30 '17 at 18:10
  • 2
    $\begingroup$ @nicoguaro At this moment, the thread looks like: "question: Am I using right formulas? answer: Yes, right ones". I'm not sure how should it be properly archived. Obviously, the information contained in the question itself is useful (see 1st commenter)... $\endgroup$ – omican Jan 14 '18 at 10:39

Not the answer you're looking for? Browse other questions tagged or ask your own question.