Methods of solving non-linear advection-diffusion systems beyond Newton-Raphson?

Question

I'm working on a project where I have two adv-diff coupled domains through their respective source terms (one domain adds mass, the other subtracts mass). For brevity, I'm modeling them in steady state. The equations are your standard advection-diffusion transport equation with a source term look like this:

$$ \frac{\partial c_1}{\partial t} = 0 = \mathcal{F}_1 + \mathcal{Q}_1(c_1,c_2) \\ \frac{\partial c_2}{\partial t} = 0 = \mathcal{F}_2 + \mathcal{Q}_2(c_1,c_2) $$

Where $\mathcal{F}_i$ is diffusive and advective flux for species $i$, and $\mathcal{Q}_i$ is the source term for species $i$.

I have been able to write a solver for my problem using the Newton-Raphson method, and have completely coupled the two domains using a block mass matrix, ie:

$$ F_{coupled} = \left[\begin{array}{c c} A_1 & 0 \\ 0 & A_2 \\ \end{array}\right]\underbrace{ \left[\begin{array}{c} c_{1,i} \\ c_{2,i} \\ \end{array}\right] }_{x_i} - \left[\begin{array}{c} b_1(c_{1,i}, c_{2,i}) \\ b_2(c_{1,i}, c_{2,i}) \\ \end{array}\right] $$

The term $F_{coupled}$ is used to determine the Jacobian matrix and update both $c_1$ and $c_2$:

$$ \mathcal{J}(x_i) \left[x_{i+1} - x_i \right] = -\mathcal{F}_{coupled} $$

or

$$ x_{i+1} = x_i - \left(\mathcal{J}(x_i)\right)^{-1} \mathcal{F}_{coupled} $$

To speed things up, I don't calculate the Jacobian every iteration - right now I'm playing with every five iterations, which has seemed to work well enough and keep the solution steady.

The problem is: I'm going to be moving to a larger system where both domains are in 2D/2.5D, and calculating the Jacobian matrix is going to quickly deplete my available computer resources. I'm building this model to be used in an optimization setting later on, so I also can't be behind the wheel at every iteration tuning the damping factor, etc..

Am I right to be looking elsewhere for a more robust and algorithm for my problem, or is this is as good as it gets? I've looked a bit into Quasi-linearization, but am not sure how applicable it is w.r.t. to my system.

Are there any other slick algorithms that I may have missed that can solve a system of nonlinear equations without resorting to re-calculating the Jacobian as offen?

Answer 1

I'm assuming the limitation in 2D and 3D is storing the Jacobian.

One option is to retain the time derivatives and use an explicit "pseudo" time-stepping to iterate to steady state. Normally the CFL number you need for diffusive and reactive systems might get prohibitively small. You could try nonlinear multigrid (it's also called Full Approximation Storage multigrid) and local time-stepping to speed up convergence.

The other option is to use a fully implicit scheme as you're doing now, but not store the global Jacobian. You could use a matrix-free implicit scheme. $$ DF(u^n)\, \delta u^n = -F(u^n) $$ (where $DF$ is the Jacobian) can be solved with a Krylov subspace solver like GMRES and BiCGStab using the fact that $$ DF(u^n)\, \delta u \approx \frac{F(u^n+\epsilon\frac{\delta u}{\Vert\delta u\Vert}) - F(u^n)}{\epsilon}. $$ This is because GMRES and BiCGStab don't require a LHS matrix $A$, they only need to be able to compute its product $Ax$ given a vector $x$.

Now with a proper value of $\epsilon$ (usually about $10^{-7}$ for double-precision floats) you can execute a Newton loop without ever computing or storing a Jacobian. I know for a fact that this technique is used to solve some non-trivial cases in computational fluid dynamics. Note, however, that the number of evaluations of the function $F$ will be more than in a matrix-storage technique, instead of requiring a matrix-vector product.

Another thing to note is that if your system is such that a powerful preconditioner is needed (ie. Jacobi or block-Jacobi will not suffice), you might want to try using the above-mentioned method as a smoother in a multigrid scheme. If you want to try a point- or block-Jacobi preconditioner, you could compute and store only the diagonal elements or diagonal blocks of the Jacobian, which is not much. I would also mention that a Gauss-Seidel or SSOR preconditioner may be possible to implement without explicitly storing a Jacobian. This paper describes the implementation of a matrix-free GMRES preconditioned with matrix-free symmetric Gauss-Seidel in the context of computational fluid dynamics.

Answer 2

From my experience with Navier-Stokes equations, one can do very well without fully implicit schemes.

If you just want a fast numerical scheme for the solution of the time evolution, take a look at IMEX (implicit-explicit) schemes; see e.g. this paper by Ascher, Ruuth, Spiteri Implicit-Explicit Runge-Kutta Methods for Time-Dependent Partial Differential Equations.

You may even try to just use an explicit high-order time integration scheme with step size control (like Matlab's ODE45). However, you may run into trouble because of the stiffness of your system, which comes from the diffusive part. Luckily, the diffusive part is linear so that it can be treated implicitly (which is the idea of the IMEX schemes).

Answer 3

At first I considered to add only a remark, but the space was not enough, so I add some brief description of my experiences with this topic.

First, looking to your notation of $F_{coupled}$ I do not see the coupled form, I suppose that $b_1$ and $b_2$ shall be both dependent on $c_{1,i}$ and $c_{2,i}$. Moreover if $A_1$ and $A_2$ are matrix representations of approximations of $\mathcal{F}_1$ and $\mathcal{F}_2$ then they should not depend only on $c_i$, but also on neighbor values, but this may be only a wrong understanding of your notation.

As a general comment I would like to add that using analytical Jacobian seems to be the only way to obtain a quadratic convergence of nonlinear iterative solver (i.e. the Newton-Raphson solver in your case). Did you observe it in your case? It is quite important, because otherwise there could be some misunderstanding in your approximations (linearization).

In all applications I was involved in (some of them included large scale computations) we had never issue with time consumption of assembling the Jacobian, the most time consuming issue was always applying a linear solver. The analytical Jacobian (if available) has been always in applications I was working on the preferred choice because of quadratic convergence. In few cases such nonlinear solver produces a matrix that causes problems to the convergence of iterative linear solver, so we tried to use a simpler linearization than analytical Jacobian to help the linear solver. Such a trade off between behavior of nonlinear and linear algebraic solver depending on the linearization of nonlinear algebraic system was always tricky and I could not give a general recommendation.

But you are right that the drawback (or property) of analytical Jacobian for system of PDEs is that it produces coupled algebraic system, so if you decoupled such system in any way, e.g. solving separately the approximation of each PDE by, say, iterative splitting method, then again you loose the quadratic convergence of global solver. But then at least if you solve separately each discretizied (decoupled) PDE, you can again speed up the solution for this particular problem by using Newton-Raphson method.