From 1b83961ed4ecda43dd07a3abc1143d92e113b7dc Mon Sep 17 00:00:00 2001 From: MFraters Date: Thu, 8 Feb 2018 01:50:06 +0100 Subject: [PATCH] Change two enums for stabilisation into one enum which is part of the NewtonHandler. --- include/aspect/newton.h | 60 +++++----- source/simulator/assemblers/newton_stokes.cc | 16 +-- source/simulator/newton.cc | 111 ++++++++----------- source/simulator/solver_schemes.cc | 11 +- 4 files changed, 94 insertions(+), 104 deletions(-) diff --git a/include/aspect/newton.h b/include/aspect/newton.h index 0cbbb5c6237..f1a15c76072 100644 --- a/include/aspect/newton.h +++ b/include/aspect/newton.h @@ -62,19 +62,27 @@ namespace aspect * material model outputs for the Newton solver. */ - enum class PreconditionerStabilization - { - SPD, PD, symmetric, none - }; - enum class VelocityBlockStabilization - { - SPD, PD, symmetric, none - }; - template class NewtonHandler: public SimulatorAccess { public: + + /** + * 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 + * Positive Definite, symmetric represents that the matrix is + * only symmetrized, and PD represents that we do the same as + * what we do for SPD, but without the symmetrization. + */ + enum class NewtonStabilization + { + none = 0, + symmetric = 1, + PD = 2, + SPD = symmetric | PD + }; + /** * Determine, based on the run-time parameters of the current simulation, * which functions need to be called in order to assemble linear systems, @@ -118,22 +126,22 @@ namespace aspect /** * Get the stabilization type used in the preconditioner. */ - PreconditionerStabilization get_preconditioner_stabilization() const; + NewtonStabilization get_preconditioner_stabilization() const; /** * Set the stabilization type used in the preconditioner. */ - void set_preconditioner_stabilization(const PreconditionerStabilization preconditioner_stabilization); + void set_preconditioner_stabilization(const NewtonStabilization preconditioner_stabilization); /** * Get the stabilization type used in the velocity block. */ - VelocityBlockStabilization get_velocity_block_stabilization() const; + NewtonStabilization get_velocity_block_stabilization() const; /** * Sets the stabilization type used in the velocity block. */ - void set_velocity_block_stabilization(const VelocityBlockStabilization velocity_block_stabilization); + void set_velocity_block_stabilization(const NewtonStabilization velocity_block_stabilization); /** * Get whether to use the Newton failsafe. If the failsafe is used, a failure @@ -172,13 +180,7 @@ namespace aspect /** * get a std::string describing the stabilization type used for the preconditioner. */ - std::string get_preconditioner_stabilization_string(PreconditionerStabilization preconditioner_stabilization); - - /** - * get a std::string describing the stabilization type used for the velocity block. - */ - std::string get_velocity_block_stabilization_string(VelocityBlockStabilization velocity_block_stabilization); - + std::string get_newton_stabilization_string(const NewtonStabilization preconditioner_stabilization) const; /** * Declare additional parameters that are needed for the Newton. @@ -200,15 +202,15 @@ namespace aspect * See the get_newton_derivative_scaling_factor() function for an * explanation of the purpose of this factor. */ - double newton_derivative_scaling_factor; - PreconditionerStabilization preconditioner_stabilization; - VelocityBlockStabilization velocity_block_stabilization; - bool use_Newton_failsafe; - double nonlinear_switch_tolerance; - unsigned int max_pre_newton_nonlinear_iterations; - unsigned int max_newton_line_search_iterations; - bool use_newton_residual_scaling_method; - double maximum_linear_stokes_solver_tolerance; + double newton_derivative_scaling_factor; + NewtonStabilization preconditioner_stabilization; + NewtonStabilization velocity_block_stabilization; + bool use_Newton_failsafe; + double nonlinear_switch_tolerance; + unsigned int max_pre_newton_nonlinear_iterations; + unsigned int max_newton_line_search_iterations; + bool use_newton_residual_scaling_method; + double maximum_linear_stokes_solver_tolerance; }; namespace Assemblers diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc index 14e450d7287..92c80180f98 100644 --- a/source/simulator/assemblers/newton_stokes.cc +++ b/source/simulator/assemblers/newton_stokes.cc @@ -119,14 +119,14 @@ namespace aspect const SymmetricTensor<2,dim> viscosity_derivative_wrt_strain_rate = derivatives->viscosity_derivative_wrt_strain_rate[q]; const SymmetricTensor<2,dim> strain_rate = scratch.material_model_inputs.strain_rate[q]; - const PreconditionerStabilization preconditioner_stabilization = this->get_newton_handler().get_preconditioner_stabilization(); + const typename NewtonHandler::NewtonStabilization preconditioner_stabilization = this->get_newton_handler().get_preconditioner_stabilization(); // todo: make this 0.9 into a global input parameter - const double alpha = preconditioner_stabilization == PreconditionerStabilization::SPD || preconditioner_stabilization == PreconditionerStabilization::PD ? + const double alpha = preconditioner_stabilization == NewtonHandler::NewtonStabilization::SPD || preconditioner_stabilization == NewtonHandler::NewtonStabilization::PD ? Utilities::compute_spd_factor(eta, strain_rate, viscosity_derivative_wrt_strain_rate, 0.9) : 1; - if (preconditioner_stabilization == PreconditionerStabilization::SPD || preconditioner_stabilization == PreconditionerStabilization::symmetric) + if (preconditioner_stabilization == NewtonHandler::NewtonStabilization::SPD || preconditioner_stabilization == NewtonHandler::NewtonStabilization::symmetric) { for (unsigned int i = 0; i < stokes_dofs_per_cell; ++i) for (unsigned int j = 0; j < stokes_dofs_per_cell; ++j) @@ -165,7 +165,7 @@ namespace aspect } #if DEBUG { - if (this->get_newton_handler().get_preconditioner_stabilization() == PreconditionerStabilization::SPD) + if (this->get_newton_handler().get_preconditioner_stabilization() == NewtonHandler::NewtonStabilization::SPD) { // regardless of whether we do or do not add the Newton // linearization terms, we ought to test whether the top-left @@ -304,14 +304,14 @@ namespace aspect const SymmetricTensor<2,dim> viscosity_derivative_wrt_strain_rate = derivatives->viscosity_derivative_wrt_strain_rate[q]; const double viscosity_derivative_wrt_pressure = derivatives->viscosity_derivative_wrt_pressure[q]; - const VelocityBlockStabilization velocity_block_stabilization = this->get_newton_handler().get_velocity_block_stabilization(); + const typename NewtonHandler::NewtonStabilization velocity_block_stabilization = this->get_newton_handler().get_velocity_block_stabilization(); // todo: make this 0.9 into a global input parameter - const double alpha = velocity_block_stabilization == VelocityBlockStabilization::SPD || velocity_block_stabilization ==VelocityBlockStabilization::PD ? + const double alpha = velocity_block_stabilization == NewtonHandler::NewtonStabilization::SPD || velocity_block_stabilization ==NewtonHandler::NewtonStabilization::PD ? Utilities::compute_spd_factor(eta, strain_rate, viscosity_derivative_wrt_strain_rate, 0.9) : 1; - if (velocity_block_stabilization == VelocityBlockStabilization::SPD || velocity_block_stabilization == VelocityBlockStabilization::symmetric) + if (velocity_block_stabilization == NewtonHandler::NewtonStabilization::SPD || velocity_block_stabilization == NewtonHandler::NewtonStabilization::symmetric) { for (unsigned int i=0; iget_newton_handler().get_velocity_block_stabilization() == VelocityBlockStabilization::SPD) + if (this->get_newton_handler().get_velocity_block_stabilization() == NewtonHandler::NewtonStabilization::SPD) { // regardless of whether we do or do not add the Newton // linearization terms, we ought to test whether the top-left diff --git a/source/simulator/newton.cc b/source/simulator/newton.cc index 2bf7a65d4f7..4b89bd14f4b 100644 --- a/source/simulator/newton.cc +++ b/source/simulator/newton.cc @@ -163,7 +163,7 @@ namespace aspect } template - PreconditionerStabilization + typename NewtonHandler::NewtonStabilization NewtonHandler:: get_preconditioner_stabilization() const { @@ -174,13 +174,13 @@ namespace aspect template void NewtonHandler:: - set_preconditioner_stabilization(const PreconditionerStabilization preconditioner_stabilization_) + set_preconditioner_stabilization(const NewtonStabilization preconditioner_stabilization_) { preconditioner_stabilization = preconditioner_stabilization_; } template - VelocityBlockStabilization + typename NewtonHandler::NewtonStabilization NewtonHandler:: get_velocity_block_stabilization() const { @@ -191,7 +191,7 @@ namespace aspect template void NewtonHandler:: - set_velocity_block_stabilization(const VelocityBlockStabilization velocity_block_stabilization_) + set_velocity_block_stabilization(const NewtonStabilization velocity_block_stabilization_) { velocity_block_stabilization = velocity_block_stabilization_; } @@ -247,44 +247,21 @@ namespace aspect template std::string NewtonHandler:: - get_preconditioner_stabilization_string(PreconditionerStabilization preconditioner_stabilization) + get_newton_stabilization_string(const NewtonStabilization preconditioner_stabilization) const { switch (preconditioner_stabilization) { - case PreconditionerStabilization::SPD: + case NewtonStabilization::SPD: return "SPD"; - break; - case PreconditionerStabilization::PD: + case NewtonStabilization::PD: return "PD"; - break; - case PreconditionerStabilization::symmetric: + case NewtonStabilization::symmetric: return "symmetric"; - break; - case PreconditionerStabilization::none: + case NewtonStabilization::none: return "none"; - break; - } - } - - template - std::string - NewtonHandler:: - get_velocity_block_stabilization_string(VelocityBlockStabilization velocity_block_stabilization) - { - switch (velocity_block_stabilization) - { - case VelocityBlockStabilization::SPD: - return "SPD"; - break; - case VelocityBlockStabilization::PD: - return "PD"; - break; - case VelocityBlockStabilization::symmetric: - return "symmetric"; - break; - case VelocityBlockStabilization::none: - return "none"; - break; + default: + Assert(false,ExcNotImplemented()); + return ""; } } @@ -301,49 +278,59 @@ namespace aspect Patterns::Double(0,1), "A relative tolerance with respect to the residual of the first " "iteration, up to which the nonlinear Picard solver will iterate, " - "before changing to the newton solver."); + "before changing to the Newton solver."); prm.declare_entry ("Max pre-Newton nonlinear iterations", "10", Patterns::Integer (0), - "The maximum number of Picard nonlinear iterations to be performed " - "before switching to Newton iterations."); + "If the 'Nonlinear Newton solver switch tolerance' is reached before the " + "maximal number of Picard iteration, then the solver switches to Newton " + "solves anyway."); prm.declare_entry ("Max Newton line search iterations", "5", Patterns::Integer (0), "The maximum number of line search iterations allowed. If the " - "criterion is not reached after this iteration, we apply the scaled " - "increment to the solution and continue."); + "criterion is not reached after this number of iteration, we apply " + "the scaled increment even though it does not satisfy the necessary " + "criteria and simply continue with the next Newton iteration."); prm.declare_entry ("Use Newton residual scaling method", "false", Patterns::Bool (), "This method allows to slowly introduce the derivatives based on the improvement " "of the residual. If set to false, the scaling factor for the Newton derivatives " - "is set to one immediately when switching on the Newton solver."); + "is set to one immediately when switching on the Newton solver. When this is set to " + "true, the derivatives are slowly introduced by the following equation: max(0.0, " + "(1.0-(residual/switch_initial_residual))), where switch_initial_residual is the " + "residual at the time when the Newton solver is switched on."); prm.declare_entry ("Maximum linear Stokes solver tolerance", "0.9", Patterns::Double (0,1), - "When the linear Stokes solver tolerance is dynamically chosen, this defines " - "the most loose tolerance allowed."); + "The linear Stokes solver tolerance is dynamically chosen for the Newton solver, based " + "on the Eisenstat walker 1994 paper (https://doi.org/10.1137/0917003), equation 2.2. " + "Because this value can become larger then one, we limit this value by this parameter."); prm.declare_entry ("Stabilization preconditioner", "SPD", Patterns::Selection ("SPD|PD|symmetric|none"), - "This parameters allows for the stabilisation of the preconditioner. By default, the " - "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 " + "This parameters allows for the stabilization of the preconditioner. If one derives the Newton " + "method without any modifications, the 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 preconditioner 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 guaranteed Symmetric " "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 " - "matrix created for the velocity block is not necessarily Symmetric Positive Definite. " - "This is problematic (see Fraters et al, in prep). When none is chosen, the velocity block " - "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 " + "This parameters allows for the stabilization of the velocity block. If one derives the Newton " + "method without any modifications, the matrix created for the velocity block is not necessarily " + "Symmetric Positive Definite. This is problematic (see Fraters et al., in prep). When `none' is chosen, " + "the velocity block 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 guaranteed Symmetric " "Positive Definite."); + prm.declare_entry ("Use Newton failsafe", "false", Patterns::Bool (), - "Switches on SPD stabilization when solver fails."); + "When this parameter is true and the linear solver fails, we try again, but now with SPD stabilization " + "for both the preconditioner and the velocity block. The SPD stabilization will remain active untill " + "the next timestep, when the default values are restored."); } prm.leave_subsection (); } @@ -366,23 +353,23 @@ namespace aspect maximum_linear_stokes_solver_tolerance = prm.get_double("Maximum linear Stokes solver tolerance"); std::string preconditioner_stabilization_string = prm.get("Stabilization preconditioner"); if (preconditioner_stabilization_string == "SPD") - preconditioner_stabilization = PreconditionerStabilization::SPD; + preconditioner_stabilization = NewtonStabilization::SPD; else if (preconditioner_stabilization_string == "PD") - preconditioner_stabilization = PreconditionerStabilization::PD; + preconditioner_stabilization = NewtonStabilization::PD; else if (preconditioner_stabilization_string == "symmetric") - preconditioner_stabilization = PreconditionerStabilization::symmetric; + preconditioner_stabilization = NewtonStabilization::symmetric; else if (preconditioner_stabilization_string == "none") - preconditioner_stabilization = PreconditionerStabilization::none; + preconditioner_stabilization = NewtonStabilization::none; std::string velocity_block_stabilization_string = prm.get("Stabilization velocity block"); if (velocity_block_stabilization_string == "SPD") - velocity_block_stabilization = VelocityBlockStabilization::SPD; + velocity_block_stabilization = NewtonStabilization::SPD; else if (velocity_block_stabilization_string == "PD") - velocity_block_stabilization = VelocityBlockStabilization::PD; + velocity_block_stabilization = NewtonStabilization::PD; else if (velocity_block_stabilization_string == "symmetric") - velocity_block_stabilization = VelocityBlockStabilization::symmetric; + velocity_block_stabilization = NewtonStabilization::symmetric; else if (velocity_block_stabilization_string == "none") - velocity_block_stabilization = VelocityBlockStabilization::none; + velocity_block_stabilization = NewtonStabilization::none; use_Newton_failsafe = prm.get_bool("Use Newton failsafe"); diff --git a/source/simulator/solver_schemes.cc b/source/simulator/solver_schemes.cc index ad8e8518155..1bd08b08790 100644 --- a/source/simulator/solver_schemes.cc +++ b/source/simulator/solver_schemes.cc @@ -579,9 +579,9 @@ namespace aspect residual, residual_old); - 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; + const std::string preconditioner_stabilization = newton_handler->get_newton_stabilization_string(newton_handler->get_preconditioner_stabilization()); + const std::string velocity_block_stabilization = newton_handler->get_newton_stabilization_string(newton_handler->get_velocity_block_stabilization()); + pcout << " The linear solver tolerance is set to " << parameters.linear_stokes_solver_tolerance << ". Stabilization Preconditioner is " << preconditioner_stabilization << " and A block is " << velocity_block_stabilization << "." << std::endl; } } @@ -599,10 +599,11 @@ namespace aspect } catch (...) { + // start the solve over again and try with a stabilized version computing_timer.exit_section(); pcout << "Solve failed and catched, try again with stabilisation" << std::endl; - newton_handler->set_preconditioner_stabilization(PreconditionerStabilization::SPD); - newton_handler->set_velocity_block_stabilization(VelocityBlockStabilization::SPD); + newton_handler->set_preconditioner_stabilization(NewtonHandler::NewtonStabilization::SPD); + newton_handler->set_velocity_block_stabilization(NewtonHandler::NewtonStabilization::SPD); rebuild_stokes_matrix = rebuild_stokes_preconditioner = assemble_newton_stokes_matrix = true; assemble_stokes_system();