Skip to content

Commit

Permalink
Change two enums for stabilisation into one enum which is part of the…
Browse files Browse the repository at this point in the history 
… NewtonHandler.
  • Loading branch information
@MFraters
MFraters committed Feb 8, 2018
1 parent ccdcf64 commit 1b83961
Show file tree
Hide file tree
Showing 4 changed files with 94 additions and 104 deletions.
60 changes: 31 additions & 29 deletions include/aspect/newton.h
View file
Expand Up @@ -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 <int dim>
class NewtonHandler: public SimulatorAccess<dim>
{
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,
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down
16 changes: 8 additions & 8 deletions source/simulator/assemblers/newton_stokes.cc
View file
Expand Up @@ -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<dim>::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<dim>::NewtonStabilization::SPD || preconditioner_stabilization == NewtonHandler<dim>::NewtonStabilization::PD ?
Utilities::compute_spd_factor<dim>(eta, strain_rate, viscosity_derivative_wrt_strain_rate, 0.9)
:
1;

if (preconditioner_stabilization == PreconditionerStabilization::SPD || preconditioner_stabilization == PreconditionerStabilization::symmetric)
if (preconditioner_stabilization == NewtonHandler<dim>::NewtonStabilization::SPD || preconditioner_stabilization == NewtonHandler<dim>::NewtonStabilization::symmetric)
{
for (unsigned int i = 0; i < stokes_dofs_per_cell; ++i)
for (unsigned int j = 0; j < stokes_dofs_per_cell; ++j)
Expand Down Expand Up @@ -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<dim>::NewtonStabilization::SPD)
{
// regardless of whether we do or do not add the Newton
// linearization terms, we ought to test whether the top-left
Expand Down Expand Up @@ -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<dim>::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<dim>::NewtonStabilization::SPD || velocity_block_stabilization ==NewtonHandler<dim>::NewtonStabilization::PD ?
Utilities::compute_spd_factor<dim>(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<dim>::NewtonStabilization::SPD || velocity_block_stabilization == NewtonHandler<dim>::NewtonStabilization::symmetric)
{
for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
Expand Down Expand Up @@ -349,7 +349,7 @@ namespace aspect
#if DEBUG
if (scratch.rebuild_newton_stokes_matrix)
{
if (this->get_newton_handler().get_velocity_block_stabilization() == VelocityBlockStabilization::SPD)
if (this->get_newton_handler().get_velocity_block_stabilization() == NewtonHandler<dim>::NewtonStabilization::SPD)
{
// regardless of whether we do or do not add the Newton
// linearization terms, we ought to test whether the top-left
Expand Down
111 changes: 49 additions & 62 deletions source/simulator/newton.cc
View file
Expand Up @@ -163,7 +163,7 @@ namespace aspect
}

template <int dim>
PreconditionerStabilization
typename NewtonHandler<dim>::NewtonStabilization
NewtonHandler<dim>::
get_preconditioner_stabilization() const
{
Expand All @@ -174,13 +174,13 @@ namespace aspect
template <int dim>
void
NewtonHandler<dim>::
set_preconditioner_stabilization(const PreconditionerStabilization preconditioner_stabilization_)
set_preconditioner_stabilization(const NewtonStabilization preconditioner_stabilization_)
{
preconditioner_stabilization = preconditioner_stabilization_;
}

template <int dim>
VelocityBlockStabilization
typename NewtonHandler<dim>::NewtonStabilization
NewtonHandler<dim>::
get_velocity_block_stabilization() const
{
Expand All @@ -191,7 +191,7 @@ namespace aspect
template <int dim>
void
NewtonHandler<dim>::
set_velocity_block_stabilization(const VelocityBlockStabilization velocity_block_stabilization_)
set_velocity_block_stabilization(const NewtonStabilization velocity_block_stabilization_)
{
velocity_block_stabilization = velocity_block_stabilization_;
}
Expand Down Expand Up @@ -247,44 +247,21 @@ namespace aspect
template <int dim>
std::string
NewtonHandler<dim>::
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 <int dim>
std::string
NewtonHandler<dim>::
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 "";
}
}

Expand All @@ -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 ();
}
Expand All @@ -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");

Expand Down

0 comments on commit 1b83961

Please sign in to comment.