[aspect-devel] Newton solver

Fri Sep 25 09:47:22 PDT 2015

Thanks for your detailed response. I really like that you structured it in a way in which it can be implemented in small pieces. For the moment I still have two more questions. 

My first question is if you could elaborate a bit on why the first option has your preference? 

My second question is if it would be possible the make evaluating the derivative optional, because otherwise if I would implement this, I would need to implement the derivative for all current material models of which some work with quite complex rheologies and cut of the viscosity at a upper and lower threshold, which can make calculating the derivative very hard.

Thanks,

Menno

> From: bangerth at tamu.edu
> To: menno.fraters at outlook.com; aspect-devel at geodynamics.org
> Subject: Re: [aspect-devel] Newton solver
> Date: Wed, 23 Sep 2015 12:54:08 +0000
> 
> 
> 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/
> 
> _______________________________________________
> Aspect-devel mailing list
> Aspect-devel at geodynamics.org
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20150925/b814a508/attachment.html>