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 9, 2018
1 parent a893829 commit 3763b96
Show file tree
Hide file tree
Showing 4 changed files with 107 additions and 104 deletions.
69 changes: 40 additions & 29 deletions include/aspect/newton.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,19 +62,36 @@ 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 NewtonStabilization
{
none = 0,
symmetric = 1,
PD = 2,
SPD = symmetric | PD
};

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

/**
* 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 +135,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 +189,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 +211,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
20 changes: 12 additions & 8 deletions source/simulator/assemblers/newton_stokes.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,14 +119,16 @@ 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();
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 ?
// use the spd factor when the stabilization is PD or SPD
const double alpha = (preconditioner_stabilization | NewtonHandler<dim>::NewtonStabilization::PD) != NewtonHandler<dim>::NewtonStabilization::none ?
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)
// symmetrize when the stabilization is symmetric or SPD
if ((preconditioner_stabilization | NewtonHandler<dim>::NewtonStabilization::symmetric) != NewtonHandler<dim>::NewtonStabilization::none)
{
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 +167,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 +306,16 @@ 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();
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 ?
// use the spd factor when the stabilization is PD or SPD
const double alpha = (velocity_block_stabilization | NewtonHandler<dim>::NewtonStabilization::PD) != NewtonHandler<dim>::NewtonStabilization::none ?
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)
// symmetrize when the stabilization is symmetric or SPD
if ((velocity_block_stabilization | NewtonHandler<dim>::NewtonStabilization::symmetric) != NewtonHandler<dim>::NewtonStabilization::none)
{
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 +353,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

0 comments on commit 3763b96

Please sign in to comment.