# [aspect-devel] Newton solver

Wolfgang Bangerth bangerth at tamu.edu
Wed Sep 23 05:54:08 PDT 2015

```Hi Menno,
let me put this onto the mailing list so it's archived.

> I am currently looking into the possibilities of implementing a newton type
> solver in ASPECT. I have watched your deal.ii video lectures on this topic
> with great interest and also looked at how other people implement newton
> solvers. I found that there are basically two ways  to implement this:
>
> 1. Calculate the true jacobian with a Gateaux differential, as you showed in
> your video lectures/deal.ii tutorials,
>
> 2. Use a jacobian free method where we approximate the jacobian times a vector
> as used by May et al, 2015 (where they refer to Knoll and Keyes, 2004).
>
> The first method has the advantage that it can use the standard GMRES
> implementation in ASPECT (and even the preconditioning can probably stay the
> same, as for example assumed in May et al, 2015). But for this method to work
> the jacobian has to be calculated without knowing the exact rheology, and I am
> not sure whether this is possible.
>
> The second method has the advantage of less computation and storage, but for
> this method to work I would need to add the interface to the nonlinear GMRES
> versions of trilinos and petsc (not sure if a nonlinear GMRES is present in
> trilinos, but I have seen it in petsc) to deal.ii, and then to ASPECT.
>
> Since I heard you had already given some thought towards the newton method, I
> was wondering which method you where thinking of using, and what your advice
> in general would be in this matter.

I'm definitely in favor of option 1.

Here are in essence the steps that I believe one will have to do to make this
happen (preferably in small increments that can be individually reviewed):

1/ Figure out the sparsity pattern of the Newton matrix. For simplicity, let's
assume that you're currently only solving for (u,p,T), then currently the
system matrix has no entries in the (u,T) block because we assume that the
temperature is known when we solve the Stokes equations. In fact, the sparsity
pattern of the matrix can currently be represented like this:

[* *  ]
[* *  ]
[    *]

In other words, we solve the (u,p) system independently of the T equation.

If you want to be fully nonlinear in a Newton step, then you have to take into
account the other blocks as well: if, for example, the viscosity is
temperature dependent, then when you solve for the increments (du, dp, dT),
you will get a temperature term in the velocity equation. You may in that case
also not simply want to use a frozen velocity in the temperature equation, but
treat this nonlinearly as well, in which case you get a velocity term in the
temperature equation:

[* * *]
[* *  ]
[*   *]

To figure out what terms are nonzero, you need to ask the material model for
which coefficient depends on what. You may recall that we already have a bunch
of functions in the material models that describe which parameter depends on
which solution variable (e.g., "viscosity depends on pressure"). These flags
were all combined into the MaterialModel::Interface:: get_model_dependence()
function at the hackathon.

With this, you can define the sparsity pattern because you now know which of
the coefficients that appear in any given equation depend on which variables
(and therefore which updates of that variable appear in the corresponding
Newton equation).

2/ You need to assemble the system matrix for the Newton system. For this, you
need to know the *derivatives* of each coefficient with regard to each
variable. You already know from step 1 above which derivatives are nonzero to
begin with.

There are currently no functions in the material model that can compute these
derivatives. A good first step for this project may in fact be to figure out
how such a function should look like. I could imagine that maybe a function

virtual
void
evaluate_derivatives (const MaterialModelInputs &in,
const Variable solution_variable,
MaterialModelOutputs &)

could work -- the second argument denotes with regard to *which variable* to
compute derivatives. This should definitely work for scalar variables. One
would have to think about what to do with regard to derivatives with respect
to the strain tensor. Maybe it would be simpler to declare a structure
MaterialModelDerivatives and have the function look like this:

virtual
void
evaluate_derivatives (const MaterialModelInputs &in,
MaterialModelDerivatives &)

(You need a function like this because Newton's method only converges fast if
you have good information about the derivatives of your coefficients with
regard to the solution variables. Only the material model can provide this. If
you don't have this kind of information, Newton's method is not substantially
faster than the fixed point iteration we currently already have.)

3/ One would need to compute the Newton rhs vector. This is simply the
residual -- should not be overly complicated.

4/ One would need to solve the equations. In particular, this requires a
preconditioner. I think I know how to do this, but it will require a good deal
of experimentation and trial/error, so let's postpone this until the other
pieces are in place.

If you're interested in tackling this, my suggestion would be to do it in
small increments that can be individually designed, implemented and tested.
It's going to be a significant change to many parts of the core of ASPECT, so
don't go off for half a year implementing things, but keep discussing your
thoughts on the mailing list before you start writing code.

My suggestion would be to start with point 2 above. It's relatively
self-contained and could be implemented by itself without touching too many
other things.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth               email:            bangerth at math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

```