[aspect-devel] Newton solver

[Wolfgang Bangerth]

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/