I'm trying to solve the following equation with FVM

$$u_t + f(u)_x = g(u)$$

where $g$ is some smooth function of $u$ and $f(u) = \frac{u^2}{2}$. This is really similar to Burger's equation, except for the source term. My main problem is that I don't know how to derive correctly the numerical scheme since I have a source term here.

More precisely, taking a space-time cell:

$\int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} u(t_{n+1},x) - u(t_{n},x) + \int_{t_n}^{t_{n+ 1}} f(u(x_{j+\frac{1}{2}}),t) - f(u(x_{j-\frac{1}{2}},t))dt = \int_{t_n}^{t_{n+1}} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} g(u)du$

Now I define as usual \begin{align} \bar{u}_j^n = \frac{1}{\Delta x} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} u(x,t_n)dx \end{align} and

\begin{align} F_{j+ \frac{1}{2}^n} = \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} f(u(x_{j+\frac{1}{2}},t))dt \end{align}

and hence the numerical scheme writes

$$\bar{u}_{j+1}^n = \bar{u}_j^n - \frac{\Delta t}{\Delta x} (F_{j+ \frac{1}{2}^n} - F_{j-\frac{1}{2}}^n) + g(\bar{u}_j^n)$$

where I replaced the average of $g(u)$ with its evaluation on the average, which is a second order approximation.

My question: in order to determine the numerical fluxes I usually solve Riemann problems at the interfaces $x=x_{j+\frac{1}{2}}$, like

$u_t + f(u)_x = 0$ with initial data $u(x,t_n)$ given by \begin{align} \bar{u}_{j-1}^n , \quad x< x_{j-\frac{1}{2}}\\ \bar{u}_j^n, \quad x> x< x_{j+\frac{1}{2}} \end{align}

My problem is that now I should also include the source term in the Riemann problem, so how can I handle it? I need to understand how to find an expression for the numerical fluxes in this case

  • 3
    $\begingroup$ You asked this question already: scicomp.stackexchange.com/questions/34691/…. People were trying to help you with this. Deleting your question and asking again really discourages anyone from trying to help. $\endgroup$ Mar 31, 2020 at 7:43
  • $\begingroup$ @DavidKetcheson yes but I realized that in my question was not asked in a good way, so I decided to consider a simpler problem, in order to understand. $\endgroup$
    – VoB
    Mar 31, 2020 at 8:32
  • $\begingroup$ @VoB In these situations, the question should be edited rather than deleted and asked again. This is a link to chat with a previous discussion and helpful comments. $\endgroup$
    – Anton Menshov
    Mar 31, 2020 at 13:34
  • $\begingroup$ @AntonMenshov I thought that the edited question was not so clear, I'm sorry, it wasn't my intention to spam on the site. Also, I think that this one is clearer because the equation should be a bit simpler $\endgroup$
    – VoB
    Mar 31, 2020 at 14:09
  • $\begingroup$ @VoB don't worry :) Your question seems to be pretty interesting, and well described now. $\endgroup$
    – Anton Menshov
    Mar 31, 2020 at 14:36

2 Answers 2


The usual techniques for solving Riemann problem rely on the self-similar structure of the solutions. Some general techniques can be developed which can be applied to any hyperbolic problem.

If you add a source term, the self-similarity is lost.

Whether you can solve the RP with source term depends on the PDE and the precise form of the source terms, so there may not be general methods. It may be possible to do it in some simple cases.

The simplest approach is to derive the fluxes ignoring the source term and then add the source term, leading to an unsplit method. This may work well or not, depends on the problem.

There may be other considering like stiff source terms, need for well-balancing, etc., which may necessitate special methods.

Operator splitting is another idea which has been mentioned by others. May be useful to deal with stiff problems, via an explicit-implicit approach.

I would suggest seeing Chapter 7 of [1] which is all about source terms, and references given there.

[1] LeVeque, Finite Volume Methods for Hyperbolic Problems.


The time evolution equation in hand is $\frac{\partial}{\partial{t}}{u} = L_1(u) + L_2(u)$ where the operators in the RHS are $L_1 = -\frac{\partial}{\partial{x}}{f(u)}$ and $L_2 = g(u)$. The operator splitting techniques consists of doing time steps for the PDE combining two substeps, one using only the first operator in the RHS to produce an intermediate update, the other using only the second operator, e.g., for the first order explicit Euler,

$u^{*} = u^j + \tau L_1 (u^j)$


$u^{j+1} = u^{*} + \tau L_2 (u^{*})$

The total time step obtained in such a manner has low order accuracy (even if the individual substeps were highly accurate) but there are variations of this method that improve the accuracy order, e.g., see Alternating Directions Implicit method and Strang splitting on Wikipedia.

What this means for the problem in hand, is that one can combine time steps with the flux divergence only (ignoring the source) and with the source only (ignoring flux divergence) to produce a numerical scheme for time evolution. So, if a Riemann solver is used for treating the flux divergence term, in this method it would be used in the standard way, without the source.

Here are some references to get acquainted with the operator splitting method (also called the fractional steps method)






https://www.cec.uchile.cl/cinetica/pcordero/MC_libros/NumericalRecipesinC.pdf (p. 856)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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