Split newton stabilisation as options and add a failsafe #2006

MFraters · 2017-11-17T11:51:44Z

This is something which I have been working on and with for a while now, and since #1996 is (hopefully) nearly merged, it seems me like a good time to discuss this now. This pull request actually add two features at the same time, because they used to be intertwined, but I untangled them just before this pull request. If you want I can make two pull requests from this one.

The first feature is that it will now be possible to chose whether to use the stabilisation or not (and which of the two) for both the preconditioner and the A block.

My tests with the Spiegelman benchmark (I only tested with preconditioner being SPD) seem to indicate that if the solver doesn't fail, the convergence rate is by far the best if you only symmetrise the the A block (symmetric), then much closer together SPD, none, PD. But both the symmetric and the none fail when problems get tougher. This is why I introduce the second feature.

The second feature is that when the solver doesn't converge, we try again, but set everything to SPD. This way I can get the best of both world (being symmetric and SPD).

If you all agree with this approach, then will actually clean up the code, add more documentation and for example change the strings of "SPD" and "none" into enums.

gassmoeller

In principle I think this is a nice approach, and we should pursue it. I am a bit hesitant about all those Newton parameters ending up in the main Parameters class though. I would be very happy if those could be moved into the NewtonHandler class like for the melt settings. Let me know when I should take another look (I realize this is something that could take a bit of time).

gassmoeller · 2018-01-30T07:06:10Z

include/aspect/parameters.h

@@ -308,6 +308,9 @@ namespace aspect

    double                         nonlinear_tolerance;
    double                         nonlinear_switch_tolerance;
+    std::string                    use_Newton_stabilisation_preconditioner;


I feel like it is time to give the Newton handler its own declare_parameters and parse_parameters functions. This way we do not clutter the main parameters with things that are specific to the Newton handler. The melt handler already does this, and probably you can copy the same approach (see parameters.cc:1728, core.cc:544, melt.cc:1367). Also as you mentioned these can be enum

gassmoeller · 2018-01-30T07:07:30Z

include/aspect/newton.h

@@ -115,6 +126,7 @@ namespace aspect
       * explanation of the purpose of this factor.
       */
      double newton_derivative_scaling_factor;
+      std::pair<std::string,std::string> Newton_stabilisation;


Keep this as two separate variables and document what they do. You can return them as a pair in get_Newton_stabilization, but set them as two separate input arguments in set_Newton_stabilization

gassmoeller · 2018-01-30T07:08:15Z

include/aspect/newton.h

+       * Gets the Newton derivative scaling factor used for scaling the
+       * derivative part of the Newton Stokes solver in the assembly.
+       */
+      std::pair<std::string,std::string> get_Newton_stabilisation() const;


stabilization please. Otherwise we end up in an American vs. British English war 😄.

After reading through the rest of the PR, I think it would be better to split both functions into separate functions that could be named set_preconditioner_stabilization and set_velocity_block_stabilization respectively (and the same for the get_... functions. This is more intuitive than returning or setting a pair of parameters. Otherwise it would get real complicated if somebody added a third stabilization option.

gassmoeller · 2018-01-30T07:09:48Z

include/aspect/parameters.h

@@ -308,6 +308,9 @@ namespace aspect

    double                         nonlinear_tolerance;
    double                         nonlinear_switch_tolerance;
+    std::string                    use_Newton_stabilisation_preconditioner;
+    std::string                    use_Newton_stabilisation_A_block;
+    bool                           use_Newton_failsafe;


If we decide this is a good way, we might even omit this variable and always fall back to the failsafe.

gassmoeller · 2018-01-30T07:13:14Z

source/simulator/assemblers/newton_stokes.cc

@@ -96,6 +96,7 @@ namespace aspect

          const double JxW = scratch.finite_element_values.JxW(q);

+


Do you need this line?

bangerth · 2018-02-05T18:39:30Z

@MFraters -- can you address @gassmoeller 's comments first? Once you're done with those, let me know and I'll take another look.

I'd say make this a priority today, if you can!

MFraters · 2018-02-05T22:58:05Z

I addressed all comments except the enum one. If I make it an enum, is there an easy way to output the enum name to the screen without a bunch of if statements?

bangerth

A couple comments on the header file. The rest to come later today.

bangerth · 2018-02-07T15:00:31Z

include/aspect/newton.h

+      /**
+       * Set the stabilization type used in the preconditioner.
+       */
+      void set_preconditioner_stabilization(const std::string preconditioner_stabilization);


I see what you mean by using these strings. I think it would indeed be cleaner if you defined enums here and in similar places that these functions could take/return

bangerth · 2018-02-07T15:00:40Z

include/aspect/newton.h

+      /**
+       * Sets the stabilization type used in the velocity block.
+       */
+      void set_velocity_block_stabilization(const std::string velocity_block_stabilization);


bangerth · 2018-02-07T15:01:09Z

include/aspect/newton.h

+      /**
+       * Get whether to use the Newton failsafe
+       */
+      bool get_use_Newton_failsafe();


can you add another sentence to the documentation that explains what the "failsafe" is?

Also end the sentence with a period.

bangerth · 2018-02-07T15:01:30Z

include/aspect/newton.h

+
+      /**
+       * Get the nonlinear tolerance at which to switch the
+       * nonlinear solver from defect corrected Piccard to


-> defect correction
-> Picard (one 'c')

bangerth · 2018-02-07T15:02:59Z

include/aspect/newton.h

+      double get_nonlinear_switch_tolerance();
+
+      /**
+       * Get the maximum amount of pre-Newton nonlinear iterations


amount -> number
(maybe also in other places; "amount" is used for uncountable things, such as milk, flour, love; "number" is used for countable things, such as iterations, words per doc string, software bugs, etc.)

bangerth · 2018-02-07T15:05:44Z

include/aspect/newton.h

-      double newton_derivative_scaling_factor;
+      double       newton_derivative_scaling_factor;
+      std::string  preconditioner_stabilization;
+      std::string  velocity_block_stabilization;


These two variables would simply become enums then, possibly of different enum types. Of course you'd have to have a way to translate from a string to an enum and possibly back, but that should be easy to write in the parse_parameters() function.

The advantage is that enums have only finitely many possible values, and so you can easily cover them in a switch statement without much trouble. For std::string objects, there are infinitely many possibilities that, at least in principle, you'd have to check at one point or other. You also can't use them in a switch statement.

bangerth · 2018-02-07T22:02:29Z

include/aspect/newton.h

@@ -61,6 +61,16 @@ namespace aspect
   * A Class which can declare and parse parameters and creates
   * material model outputs for the Newton solver.
   */
+
+  enum class PreconditionerStabilization


The comment above this enum is still about the NewtonHandler class.

I'd say move these enums into the NewtonHandler class. Also add documentation for both.

bangerth · 2018-02-07T22:03:34Z

include/aspect/newton.h

+
+  enum class PreconditionerStabilization
+  {
+    SPD, PD, symmetric, none


I think there is an easier way: declare it like this:

none = 0, symmetric = 1; PD = 2; SPD = symmetric | PD;

You can then in the appropriate place ask whether a particular bit is set.

Same for the enum below.

bangerth · 2018-02-07T22:04:06Z

include/aspect/newton.h

+  };
+  enum class VelocityBlockStabilization
+  {
+    SPD, PD, symmetric, none


Thinking about it, since the two are the same, could you not use the same enum for both purposes, and simply have two variables of this enum type?

bangerth · 2018-02-07T22:06:18Z

source/simulator/core.cc

@@ -646,6 +646,7 @@ namespace aspect
      {
        assemble_newton_stokes_system = true;
        newton_handler->initialize_simulator(*this);
+        newton_handler->parse_parameters(prm);


Could you investigate whether this whole if block could be moved to the same location where we currently initialize the melt handler and the free surface handler (around line 530)?

bangerth · 2018-02-07T22:07:15Z

source/simulator/newton.cc

+  template <int dim>
+  std::string
+  NewtonHandler<dim>::
+  get_velocity_block_stabilization_string(VelocityBlockStabilization velocity_block_stabilization)


make this argument const. same in the declaration of the function if it's not yet there. same also for the previous function.

bangerth · 2018-02-07T22:18:57Z

source/simulator/newton.cc

+                           "matrix created for the preconditioning is not necessarily Symmetric Positive Definite. "
+                           "This is problematic (see Fraters et al, in prep). When none is chosen, the perconitioner "
+                           "is not stabilized. The Symmetric parameters symmetrizes the matrix, and PD makes the matrix "
+                           "Positive Definite. SPD is the full stabilization, where the matrix is garanteed Symmetric "


-> guaranteed

bangerth · 2018-02-07T22:20:25Z

source/simulator/newton.cc

+                           "Positive Definite.");
+        prm.declare_entry ("Stabilization velocity block", "SPD",
+                           Patterns::Selection ("SPD|PD|symmetric|none"),
+                           "This parameters allows for the stabilisation of the velocity block. By default, the "


all the same issues here. You may also just reference the preconditioner parameter (or the other way around) and just say something like "Stabilization for the velocity block of the system matrix used in the Newton method. Possible values for this parameter are the same as for ..."

bangerth · 2018-02-07T22:21:27Z

source/simulator/newton.cc

+                           "Positive Definite.");
+        prm.declare_entry ("Use Newton failsafe", "false",
+                           Patterns::Bool (),
+                           "Switches on SPD stabilization when solver fails.");


Can you elaborate on this in a couple more sentences? I.e., what solver needs to fail, and maybe reference what other parameter values are then changed? And for how long -- all future time steps, or just the current time step?

bangerth · 2018-02-07T22:22:47Z

source/simulator/solver_schemes.cc

-                pcout << "   The linear solver tolerance is set to " << parameters.linear_stokes_solver_tolerance << std::endl;
+                const std::string preconditioner_stabilization = newton_handler->get_preconditioner_stabilization_string(newton_handler->get_preconditioner_stabilization());
+                const std::string velocity_block_stabilization = newton_handler->get_velocity_block_stabilization_string(newton_handler->get_velocity_block_stabilization());
+                pcout << "   The linear solver tolerance is set to " << parameters.linear_stokes_solver_tolerance << ". Stabilisation Preconditioner is " << preconditioner_stabilization << " and A block is " << velocity_block_stabilization << "." << std::endl;


-> stabilization (with a z)

I'd also show the value for the matrix first, before the value for the preconditioner

Why do you think that that order is more logical? In my mind the preconditioner comes first.

I don't know -- I think as preconditioners as an ancillary tool to solving the linear system.

I don't feel strongly about the issue if you do :-)

bangerth · 2018-02-07T22:23:48Z

source/simulator/solver_schemes.cc

+              }
+            catch (...)
+              {
+                computing_timer.exit_section();


leave a comment here of the form

// start the solve over again and try with a stabilized version

About line 602: you are exiting a timer block here, but what timer block does the remainder of this block correspond to? Don't you have to start a timer block again somewhere?

So this exit solver.cc:525 computing_timer.enter_section (" Solve Stokes system"); (at the beginning of solve_stokes(). When the solver fails the computing_timer.exit_section(); on line 921 is never executed, so I do it here.

I don't think I understand this.

You have code like this:

timer.enter_section (...); try { ...do something... } catch (...) { timer.leave_section(); ...code now not associated with any timer section... } ...code that is now under the timer section again, but only if the try succeeded; if the try failed, this whole code block is not timed...

(I can't compare with the current code any more, because the line numbers you cite don't seem to match any more.)

Or are you saying that the solve_stokes() function enters a section at the top but the leaves through a throw statement with leaving the section? That would be a bug in that function -- it needs to leave the timer section before throwing the exception. If you don't right away see how to do this, show me the place tomorrow.

Or are you saying that the solve_stokes() function enters a section at the top but the leaves through a throw statement with leaving the section? That would be a bug in that function -- it needs to leave the timer section before throwing the exception. If you don't right away see how to do this, show me the place tomorrow.

Yes. I opened a separate pull request for this #2084 and removed it from this one.

bangerth

I think these are the only comments for the file I haven't looked through in any of the other reviews you've already gotten.

bangerth · 2018-02-08T04:39:07Z

source/simulator/assemblers/newton_stokes.cc

              // todo: make this 0.9 into a global input parameter
-              double alpha = Utilities::compute_spd_factor<dim>(eta, strain_rate, viscosity_derivative_wrt_strain_rate, 0.9);
+              const double alpha = preconditioner_stabilization == NewtonHandler<dim>::NewtonStabilization::SPD || preconditioner_stabilization == NewtonHandler<dim>::NewtonStabilization::PD ?


Break long lines.

If you do the trick with the bit masks for the enum, you can just write this as

const double alpha = (preconditioner_stabilization | NewtonHandler<dim>::NewtonStabilization::SPD != 0 ? ...);

bangerth · 2018-02-08T04:39:27Z

source/simulator/assemblers/newton_stokes.cc

+                                   :
+                                   1;
+
+              if (preconditioner_stabilization == NewtonHandler<dim>::NewtonStabilization::SPD || preconditioner_stabilization == NewtonHandler<dim>::NewtonStabilization::symmetric)


same in all other places where you have to make this sort of comparison

i.e., here you'll have to test preconditioner_stabilization | ...::symmetric

…re control over the way both the preconditioner and the A block are stabelized.

…e aplha to the unsymetrisized version of the equation.

…abilization velocity block

MFraters · 2018-02-08T19:43:40Z

include/aspect/newton.h

+       * only symmetrized, and PD represents that we do the same as
+       * what we do for SPD, but without the symmetrization.
+       */
+      enum NewtonStabilization


I had change this form enum class NewtonStabilization to enum NewtonStabilization` to make sure I could use the operator|. The downside of this is that is isn't strongly typed anymore. An other option I tried was to add a below this the function:

inline NewtonStabilization operator|(NewtonStabilization a, NewtonStabilization b) {return static_cast<NewtonStabilization>(static_cast<int>(a) | static_cast<int>(b));}

but then I got the error:

error: ‘aspect::NewtonHandler<dim>::NewtonStabilization aspect::NewtonHandler<dim>::operator|(aspect::NewtonHandler<dim>::NewtonStabilization, aspect::NewtonHandler<dim>::NewtonStabilization)’ must take exactly one argument inline NewtonStabilization operator|(NewtonStabilization a, NewtonStabilization b)

You can't make this function a member of the surrounding NewtonHandler class. It needs to be a function at namespace scope. Just declare it below the end of the NewtonHandler class as

inline typename NewtonHandler<dim>::NewtonStabilization operator|(typename NewtonHandler<dim>::NewtonStabilization a, typename NewtonHandler<dim>::NewtonStabilization b) { return static_cast<typename NewtonHandler<dim>::NewtonStabilization>( static_cast<int>(a) | static_cast<int>(b)); }

hmm, this somehow doesn't work. It doesn't seem to recognise that the inline opertor is pressent, but it compiles newton.h fine. So when I put:

template <int dim> inline typename NewtonHandler<dim>::NewtonStabilization operator|(typename NewtonHandler<dim>::NewtonStabilization a, typename NewtonHandler<dim>::NewtonStabilization b) { return static_cast<typename NewtonHandler<dim>::NewtonStabilization>( static_cast<int>(a) | static_cast<int>(b)); }

I get as error:

error: no match for ‘operator|’ (operand types are ‘aspect::NewtonHandler<2>::NewtonStabilization’ and ‘aspect::NewtonHandler<2>::NewtonStabilization’) const double alpha = preconditioner_stabilization | NewtonHandler<dim>::NewtonStabilization::PD ?

So I have tried to add NewtonHandler<dim>::NewtonStabilization operator|(typename NewtonHandler<dim>::NewtonStabilization a, typename NewtonHandler<dim>::NewtonStabilization b); \ to the bottom of newton.cc, but that doesn't help. I also tried to just explicity instanciate them:

typename NewtonHandler<2>::NewtonStabilization operator|(typename NewtonHandler<2>::NewtonStabilization a, typename NewtonHandler<2>::NewtonStabilization b) { return static_cast<typename NewtonHandler<2>::NewtonStabilization>( static_cast<int>(a) | static_cast<int>(b)); } typename NewtonHandler<3>::NewtonStabilization operator|(typename NewtonHandler<3>::NewtonStabilization a, typename NewtonHandler<3>::NewtonStabilization b) { return static_cast<typename NewtonHandler<3>::NewtonStabilization>( static_cast<int>(a) | static_cast<int>(b)); }

This works, but then I will have to specifically convert it into an int or bool again because it complains:

error: could not convert ‘aspect::operator|(velocity_block_stabilization, (aspect::NewtonHandler<2>::NewtonStabilization)1)’ from ‘aspect::NewtonHandler<2>::NewtonStabilization’ to ‘bool’

So if we really want this approach, we would also need to implement a way to covert this class to a bool.

I think I understand why it doesn't work.

Providing two specializations is ok with me. If you don't like this approach: What happens if you put the operator back into the NewtonHandler class as you had it initially, but prefix the whole declaration with the friend keyword? I.e.,

friend NewtonStabilization operator| (const NewtonStabilization a, const NewtonStabilization b) { return static_cast<NewtonStabilization>( static_cast<int>(a) | static_cast<int>(b)); }

(But if you like the explicit specialization you show, then no need to try this approach.)

As for the conversion to bool: Just do the comparison via

if ((value | ...::SPD) != ...::none)

...::none has numerical value zero, but has type equal to the enum, so you can compare with it.

ok, I got it working now with the class enum. I do not have a strong preference for either way, but I guess it is nicer to not have explicit specializations.

bangerth

Almost good. Please fix the two issues marked up and merge yourself!

Thanks for sticking with me through all of these changes!

bangerth · 2018-02-09T14:50:43Z

include/aspect/newton.h

+      /**
+       * This enum describes the type of stabilization is used
+       * for the Newton solver. None represents no stabilization,
+       * SPD represent that the resulting matrix is make Symmetric


make -> made

bangerth · 2018-02-09T14:51:27Z

include/aspect/newton.h

+       * only symmetrized, and PD represents that we do the same as
+       * what we do for SPD, but without the symmetrization.
+       */
+      enum class NewtonStabilization


I blanked -- enum class is only C++11, but we can't presume that yet for ASPECT. Can you just make it enum? I think everything should continue to work anyway without further modifications.

… NewtonHandler.

MFraters force-pushed the add_Newton_stabilisation_as_failsafe branch from 8ad1385 to aad87cf Compare November 23, 2017 12:27

MFraters force-pushed the add_Newton_stabilisation_as_failsafe branch from e87ec4a to 9b68a88 Compare January 14, 2018 16:00

gassmoeller reviewed Jan 30, 2018

View reviewed changes

MFraters changed the title ~~WIP: Split newton stabilisation as options and add a failsafe~~ Split newton stabilisation as options and add a failsafe Feb 6, 2018

bangerth reviewed Feb 7, 2018

View reviewed changes

MFraters force-pushed the add_Newton_stabilisation_as_failsafe branch 2 times, most recently from cfec7ed to 4e2d658 Compare February 7, 2018 18:44

bangerth mentioned this pull request Feb 7, 2018

Move Newton parameters into a separate struct #2083

Closed

bangerth reviewed Feb 7, 2018

View reviewed changes

MFraters force-pushed the add_Newton_stabilisation_as_failsafe branch 2 times, most recently from 1b83961 to b1b4d6a Compare February 8, 2018 03:01

bangerth reviewed Feb 8, 2018

View reviewed changes

MFraters mentioned this pull request Feb 8, 2018

add the exit timer in the linear solver catch #2084

Merged

MFraters added 11 commits February 8, 2018 18:58

Add interface for the failsafe option

cca3cf7

Make use of the failsafe option in newton_stokes.cc

adc655f

set and use the failsafe in solver_schemes.cc

cebba3d

update newton_stokes.cc for failsafe.

0ecae4a

implement the failsafe in a more flexible way togheter with giving mo…

644e335

…re control over the way both the preconditioner and the A block are stabelized.

fix the Newton failsafe for deal.ii versions before 9.0.0, and add th…

f462490

…e aplha to the unsymetrisized version of the equation.

update after rebase and fix for deal.ii 8.5.

c013d8b

bugfixes for none symmetric assembly.

ff76e40

adress comments. Except for the enum.

37e91b9

add desciption for the parameters Stabilization preconditioner and St…

e6c3a35

…abilization velocity block

Make stabilization an enum.

a893829

MFraters force-pushed the add_Newton_stabilisation_as_failsafe branch 3 times, most recently from 037c654 to 7037fad Compare February 8, 2018 19:38

MFraters commented Feb 8, 2018

View reviewed changes

MFraters force-pushed the add_Newton_stabilisation_as_failsafe branch from 7037fad to 0781284 Compare February 9, 2018 04:44

bangerth reviewed Feb 9, 2018

View reviewed changes

MFraters force-pushed the add_Newton_stabilisation_as_failsafe branch from 0781284 to 3d62c9e Compare February 9, 2018 15:16

MFraters added 2 commits February 9, 2018 16:22

Change two enums for stabilisation into one enum which is part of the…

2b30ab6

… NewtonHandler.

update tests.

8cc6994

MFraters force-pushed the add_Newton_stabilisation_as_failsafe branch from 3d62c9e to 8cc6994 Compare February 9, 2018 15:25

MFraters merged commit 4e0a02a into geodynamics:master Feb 9, 2018

MFraters mentioned this pull request Feb 13, 2018

fix | bug in newton solver stablization. #2092

Merged

		@@ -96,6 +96,7 @@ namespace aspect

		const double JxW = scratch.finite_element_values.JxW(q);

Split newton stabilisation as options and add a failsafe #2006

Split newton stabilisation as options and add a failsafe #2006

Conversation

MFraters commented Nov 17, 2017

gassmoeller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bangerth commented Feb 5, 2018

MFraters commented Feb 5, 2018

bangerth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bangerth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bangerth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment