[aspect-devel] Newton solver
Menno Fraters
menno.fraters at outlook.com
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>
More information about the Aspect-devel
mailing list