Numerical solution of an advection equation, $\frac{\partial P}{\partial t}+\frac{\partial}{\partial x}\left(P^{5/3}\right)=0$, with finite volume

I was trying to solve the following equation numerically, $$\frac{\partial P}{\partial t}+\frac{\partial}{\partial x}\left(P^{5/3}\right)=0$$ I adopted the Godunov approach for discretising the equation numerically, $$\frac{P_{i}^{n+1}-P_{i}^{n}}{\Delta t}=-\frac{\mathcal{F}_{i+\frac{1}{2}}^{n}-\mathcal{F}_{i-\frac{1}{2}}^{n}}{x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}}$$ where $$\mathcal{F}$$ represents the numerical flux.

Following first order upwind scheme, I approximated the flux as, $$\mathcal{F}_{i+\frac{1}{2}}^{n}=\begin{cases} P_{i}^{5/3} & \text{if P_{i}^{2/3}>0}\\ P_{i+1}^{5/3} & \text{if P_{i}^{2/3}<0} \end{cases}.$$ However, with such a discretization, I am getting oscillations in the numerical solution. I am not sure if the discretization is correct. I will be grateful if someone can help me in this regard.

• The equation clearly only makes sense if $P\ge 0$. But then the question remains under what circumstances you consider $P^{2/3}<0$? If $P\ge 0$, isn't $P^{2/3}$ always greater or equal to zero? Jan 19 at 20:39
• Yes it is. I wrote it for completeness. The flux will only going to take the left value. Jan 20 at 10:47
• Then have you tried to debug your code? Maybe by solving a linear problem at first? Jan 20 at 18:04
• For linear problem, it is giving the right output. This is leading me to think if the calculation for "dt" that I am performing is correct or not. I was calculating dt in the following manner,$$dt=\mathcal{C}dx/P_{i}^{2/3}$$. For the right calculation instead of $P_{i}^{2/3}$ one needs to take $P_{i+1/2}^{2/3}$ but I do not know what is the value of $P_{i+1/2}$. Could this be the problem ? Jan 21 at 16:56
• you should have a single dt for the entire domain, try using the largest value in the entire domain of $a \approx \partial (P^{5/3})/\partial P = 5/3 P^{2/3}$ for that timestep. Jan 21 at 21:37

Your discretization with the upwind scheme looks correct.

One reason why you get oscillations might be that you are choosing an incorrect time step. Another one is if you have some complicated and not well-behaved boundary/initial conditions, some instabilities might arise if you don't handle them properly.

To compute a fair time step, that satisfies the Courant-Friedrichs-Lewy necessary stability condition, you can just take the maximum of the flow speed in your domain. Don't worry too much about how to compute the speed on the cell boundaries (fractional $$i$$), just take a fair value... The maximum of $$P^{2/3}$$ on the cell centers will be just fine. It is also a good practice to apply a safety factor $$c \in ]0,1[$$ to the $$\Delta t$$ you obtain in such a way, like

$$\Delta t = c 3/5 \Delta x / P^{2/3}$$

I implemented your scheme on your equation for a trivial set of boundary conditions ($$P=2$$) and initial conditions ($$P=2$$ on the left boundary, and drops linearly reaching value $$P=1$$ on the right boundary), and the solution I get looks well behaved:

import matplotlib.pyplot as plt
import numpy as np

#%% Settings
# Physical mesh points (you have to sum 1 left ghost cell)
Nx = 10
xStart = 0.
xEnd = 1.
tEnd = 1.

# Grid
x = np.linspace(xStart, xEnd, Nx)

# Initial condition
ic = np.linspace(2., 1., Nx)
# Boundary condition
bc = 2.

# Safety factor
Ccfl = 0.8

#%% Functions
def godunov(P, F, x, dt):
return P[1:] - dt * (F[1:] - F[:-1]) / (x[1:] - x[:-1])

def flux(P):
return P**(5/3)

#%% Solution
P = ic.copy()

Ps = []; ts = []
t = 0.
while t < tEnd:
dt = Ccfl*3/5 * (xEnd - xStart)/Nx / max(P)**(2/3)
F = flux(P)
P = np.append([bc], godunov(P, F, x, dt))
t += dt
# Save results
Ps.append(P)
ts.append(t)

#%% Print results
fig, ax = plt.subplots()
for tt, t in enumerate(ts):
if tt%3 == 0:
ax.plot(x, Ps[tt], 'o-', label = f"t={t}")
ax.legend()


New contributor
Rigel is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.