How to solve advective equation with source term depending on variable

Question

I have the following equation $$ \dfrac{\partial s}{\partial t} + \nabla \cdot \left( \vec{v} s\right) = f(s) $$

Where $f(s)$ is an explicit source term that depends on $s$, e.g., ($\sin(s)\;cos(s)$).

When solving with finite volume method (fvm), should I solve the above equation once until a final tolerance? Or should I keep on solving the equation until the initial residual is below a tolence?

Basically, should it be:

1:

Begin time step

Assemble linear system of Eqs and solve for $s$

End time step

or 2:

Begin time step

while(EqInitialResidual > tolerance)
    Assemble linear system of Eqs and solve for $s$
    (Here, solve once for $s$. Take the new value of $s$, compute $f(s)$, solve again for $s$ until the result is bellow a tolerance.
End time step

What is the best solution approach?

Answer 1

It depends on your temporal discretization. If you have an explicit time discretization, then $$ s_{n+1} = \mathcal{L}(s_n) $$ Since $s_n$ is known, you can explicitly compute $s_{n+1}$ using the operator $\mathcal{L}(s)$, and this can be evaluated in a bounded finite number of steps.

For example, suppose you only have a 1D problem and are using the forward Euler method and an upwinding flux (assuming that $v > 0$). Then $$ s_{n+1}(x_i) = s_n(x_i) + \Delta t \left(f(s_n(x_i)) + \frac{v}{\Delta x} (s_{n}(x_{i-1}) - s_{n}(x_{i}))\right) $$ Simply plug in $s_n$ to quickly compute $s_{n+1}$. Trying to iterate this multiple times will not improve your solution since you will get the same result every time.

If you have an implicit temporal discretization, this is no longer true as we now have $$ \mathcal{L}(s_n, s_{n+1}) = 0 $$ where generally there is no way to express $s_{n+1}$ purely in terms of $s_n$. The most common solution to this is to use the Newton-Raphson method to break your problem up into a bunch of quasi-linear problems, and iterate until you find a solution $u_{n+1}$ which sufficiently minimizes the residual (what you described as method 2).

Repeating our previous 1D example, but this time using Backwards Euler: $$ \mathcal{L}(s_n,s_{n+1}) = 0 = s_{n+1}(x_i) - s_n(x_i) - \Delta t \left(f(s_{n+1}(x_i)) + \frac{v}{\Delta x} (s_{n+1}(x_{i-1}) - s_{n+1}(x_{i}))\right) $$ This time it's not possible to explicitly solve for $s_{n+1}$ in terms of $s_n$ and constants (or if it is possible it will be relatively difficult), so we must rely on some non-linear root solving algorithm, which often iterate until a solution sufficiently close to the true solution is found.

For the Newton-Raphson method, this looks like: $$ \partial_{s'_{n+1}(x_i)}\left(\mathcal{L}(s_n,s'_{n+1})\right) \Delta s'_{n+1}(x_i) = s'_{n+1}(x_i) - s_n(x_i) - \Delta t \left(f(s'_{n+1}(x_i)) + \frac{v}{\Delta x} (s'_{n+1}(x_{i-1}) - s'_{n+1}(x_{i}))\right) $$ where $s'$ is the current guess for $s_{n+1}$, and the next guess for $s_{n+1}$ is $$ s''_{n+1} = s'_{n+1} + \Delta s'_{n+1} $$ Repeat this until $\|\mathcal{L}(s_n, s'_{n+1})\| < \epsilon$. This could take one iteration, or it may take hundreds.

Side-note: This process could be non-convergent, though explaining methods for improving the Newton-Raphson method to achieve global convergence is beyond the scope of this answer. Search for "line search" or "trust regions" in the context of nonlinear optimization for more information on this.