# Exact Riemann solver for perfect gas mixture: problem with Newton's method convergence

I'm trying to solve multicomponent Euler equations for perfect gas mixture with Godunov-like scheme using exact Riemann solver. Of course, some approximate solver would probably be more cost-effective, but for now I'd like to make it work with exact solver.

Here is the problematic part of algorithm:

We have left and right conditions, for which we calculated relevant values, namely, density, velocity, pressure, heat capacity ratio and sound speed: $$\rho_l, u_l, p_l, \gamma_l, a_l$$ and $$\rho_r, u_r, p_r, \gamma_r, a_r$$.

we try to find the velocity and pressure $$v, p$$ in the middle state that may be reached through single shock wave or rarefaction wave from left and right states at the same time.

We construct functions $$\phi_l(p), \phi_r(p)$$: $$\phi_l(p) = \begin{cases} u_l + \frac{2a_l}{\sqrt{2\gamma_l(\gamma_l-1)}}\frac{\left(1-\frac{p}{p_l}\right)}{\sqrt{1+\beta\frac{p}{p_l}}} & \text{if} \;\; p \ge p_l \;\; \text{(shock wave)} \\ u_l + \frac{2a_l}{\gamma_l-1}\left( 1-\left(\frac{p}{p_l}\right)^{\alpha} \right) & \text{if} \;\; p < p_l \;\; \text{(rarefaction)} \end{cases}$$ $$\phi_r(p) = \begin{cases} u_r - \frac{2a_r}{\sqrt{2\gamma_r(\gamma_r-1)}}\frac{\left(1-\frac{p}{p_r}\right)}{\sqrt{1+\frac{\gamma_r+1}{\gamma_r-1}\frac{p}{p_r}}} & \text{if} \;\; p \ge p_r \;\; \text{(shock wave)} \\ u_r - \frac{2a_r}{\gamma_r-1}\left( 1-\left(\frac{p}{p_r}\right)^{\alpha} \right) & \text{if} \;\; p < p_r \;\; \text{(rarefaction)} \end{cases}$$ where $$\alpha = \frac{\gamma-1}{2\gamma}$$, $$\beta = \frac{\gamma+1}{\gamma-1}$$;

and their derivatives:

$$\phi'_l(p) = \begin{cases} \frac{2a_l}{\sqrt{2\gamma_l(\gamma_l-1)}} \frac{ \frac{-\sqrt{1+\beta\frac{p}{p_l}}}{p_l} - \frac{ \beta \left(1-\frac{p}{p_l}\right)} {2p_l \sqrt{1+\beta\frac{p}{p_l}}}} {1+\beta\frac{p}{p_l}} & \text{if} \;\; p \ge p_l \;\; \text{(shock wave)} \\ - \alpha \frac{1}{p_l^{\alpha}} \frac{2a_l}{\gamma_l-1}\left(\frac{p}{p_l}\right)^{\alpha - 1} & \text{if} \;\; p < p_l \;\; \text{(rarefaction)} \end{cases}$$ $$\phi'_r(p) = \begin{cases} -\frac{2a_r}{\sqrt{2\gamma_r(\gamma_r-1)}} \frac{ \frac{-\sqrt{1+\beta\frac{p}{p_r}}}{p_r} - \frac{ \beta \left(1-\frac{p}{p_r}\right)} {2p_r \sqrt{1+\beta\frac{p}{p_r}}}} {1+\beta\frac{p}{p_r}} & \text{if} \;\; p \ge p_r \;\; \text{(shock wave)} \\ \alpha \frac{1}{p_r^{\alpha}} \frac{2a_r}{\gamma_r-1}\left(\frac{p}{p_r}\right)^{\alpha - 1} & \text{if} \;\; p < p_r \;\; \text{(rarefaction)} \end{cases}$$

and try to find $$p$$ so that $$\phi(p) = \phi_l(p) - \phi_r(p) = 0$$

For that we use Newton's method: $$p_{n+1} = p_{n} - \frac{\phi(p_n)}{\phi'(p_n)}$$ starting with $$p_0 = (p_l+p_r)/2$$ until it converges.

For single-component Euler equations it almost always converges extremely fast: typical Sod's problem $$(\rho_l,u_l,p_l = 1,0,1$$, $$\rho_r,u_r,p_r = 0.125,0,0.1)$$ converges in 5-6 steps. But for my multi-component problems typical values are something like $$p \approx 10^5-10^7 Pa$$, while $$u,a \approx 10^2-10^3 m/s$$ and also $$\phi(p), \phi'(p) \approx 10^2-10^3$$. Thus the pressure changes very little at every step of the Newton's algorithm (for example only by $$\Delta p \approx 0.3 Pa$$ while $$p_0 = 550000 Pa$$) and the convergence is extremely slow.

In some cases, I was able to achieve convergence within hundred of steps by artificially multiplying $$\Delta p$$ by something like $$10^4$$ and it gave the correct solution. But I wasn't able to construct stable algorithm for dynamic multiplier for $$\Delta p$$ that would give at least reasonable convergence speed in general case.

Is there some more general way to speed up the convergence here? Or may some different way to determine middle state pressure?

• Line search can speed up a bit. en.wikipedia.org/wiki/Line_search But I am just starting to learn numerical optimization on my own. en.wikipedia.org/wiki/Halley%27s_method maybe faster than Newton if you can evaluate the 2nd derivative. By the way, better factor out the common parts of the equations. That makes everything much easier to read and check. – R zu Oct 21 '18 at 16:59
• I mean the Armijo-Goldstein condition and backtracking line search at en.wikipedia.org/wiki/Backtracking_line_search – R zu Oct 21 '18 at 17:06
• Use a library if possible. – R zu Oct 21 '18 at 17:07
• Also watch out for Maratos effect if you are using Newton's method with constraint: mcs.anl.gov/~anitescu/CLASSES/2011/LECTURES/… – R zu Oct 22 '18 at 0:13
• How are you computing the quantities for the mixture? Also, as indicated in the citation I posted on your other question, multicomponent Euler equations often suffers from pressure oscillations. I am not entirely sure if this is what your are seeing but I would be suspect this may be part of the problem. – Kyle Mandli Oct 22 '18 at 17:20

I have found very simple solution! At least, it seem to work well with my test case. Instead of explicit calculation of $$\phi'(p_n)$$, I use it's simple approximation: $$\phi'(p_n) \approx \frac{\phi(p_n)-\phi(p_{n-1})}{p_n-p_{n-1}}$$ To use this algorithm, we need to have two starting values ($$p_0, p_1$$) instead of one ($$p_0$$). Currently I calculate $$p_1$$ via single step of original Newton's method with explicit $$\phi'(p_0)$$. As stated above, it gives very small difference between $$p_1$$ and $$p_0$$, but following steps with approximate $$\phi'(p_n)$$ produce very good increments and the method converges in usual 4-5 steps.
Update: As Praveen Chandrashekar mentioned below, I should have checked my derivatives again. There is a pretty simple error in derivative for rarefaction wave functions, it should contain $$\frac{1}{p_l}$$ instead of $$\frac{1}{p_l^\alpha}$$ (and similar thing for $$p_r$$). With correct derivatives, original Newton's method converges well.
• What you are doing is essentially secant method which has slower convergence than Newton. Of course Newton does better only if sufficiently close to the root, otherwise could be slow. Have you checked that your derivatives are correct which you use in Newton method ? Plotting your function $\phi(p)$ might also tell you something useful. – cfdlab Oct 25 '18 at 3:44