Skip to content

Commit

Permalink
Simplify logic in newton solver
Browse files Browse the repository at this point in the history
  • Loading branch information
@gassmoeller
gassmoeller committed Nov 6, 2024
1 parent 4d576b2 commit 2a805d7
Showing 1 changed file with 59 additions and 64 deletions.
123 changes: 59 additions & 64 deletions source/simulator/solver_schemes.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -529,80 +529,75 @@ namespace aspect
else
build_stokes_preconditioner();

if (newton_handler->parameters.use_Newton_failsafe == false)
try
{
dcr.stokes_residuals = solve_stokes();
}
else
catch (const std::exception &exc)
{
try
{
dcr.stokes_residuals = solve_stokes();
}
catch (const std::exception &exc)
{
// Test that the exception we got is one of the two documented by
// throw_linear_solver_failure_exception(). If not, we have a genuine
// problem here, and will need to get outta here right away:
if ((dynamic_cast<const ExcMessage *>(&exc)==nullptr) &&
(dynamic_cast<const QuietException *>(&exc)==nullptr))
throw;

// start the solve over again and try with a stabilized version
pcout << "failed, trying again with stabilization" << std::endl;
newton_handler->parameters.preconditioner_stabilization = Newton::Parameters::Stabilization::SPD;
newton_handler->parameters.velocity_block_stabilization = Newton::Parameters::Stabilization::SPD;

// If the Stokes matrix depends on the solution, or we have active
// velocity boundary conditions, we need to re-assemble the system matrix
// (and preconditioner) every time. If we have active boundary conditions,
// they could a) depend on the solution, or b) be inhomogeneous. In both
// cases, just assembling the RHS will be incorrect. If no active
// boundaries exist, we only have no-slip or free slip conditions, so we
// don't need to force assembly of the matrix.
if (stokes_matrix_depends_on_solution()
||
(nonlinear_iteration == 0 && boundary_velocity_manager.get_active_boundary_velocity_conditions().size() > 0))
rebuild_stokes_matrix = rebuild_stokes_preconditioner = assemble_newton_stokes_matrix = true;
else if (parameters.enable_prescribed_dilation)
// The dilation requires the Stokes matrix (which is on the rhs
// in the Newton solver) to be updated.
rebuild_stokes_matrix = true;
// Test that we are trying to handle exceptions and that the
// exception we got is one of the two documented by
// throw_linear_solver_failure_exception(). If not, we have a genuine
// problem here, and will need to get outta here right away:
if (newton_handler->parameters.use_Newton_failsafe == false ||
((dynamic_cast<const ExcMessage *>(&exc)==nullptr) &&
(dynamic_cast<const QuietException *>(&exc)==nullptr)))
throw;

// start the solve over again and try with a stabilized version
pcout << "failed, trying again with stabilization" << std::endl;
newton_handler->parameters.preconditioner_stabilization = Newton::Parameters::Stabilization::SPD;
newton_handler->parameters.velocity_block_stabilization = Newton::Parameters::Stabilization::SPD;

// If the Stokes matrix depends on the solution, or we have active
// velocity boundary conditions, we need to re-assemble the system matrix
// (and preconditioner) every time. If we have active boundary conditions,
// they could a) depend on the solution, or b) be inhomogeneous. In both
// cases, just assembling the RHS will be incorrect. If no active
// boundaries exist, we only have no-slip or free slip conditions, so we
// don't need to force assembly of the matrix.
if (stokes_matrix_depends_on_solution()
||
(nonlinear_iteration == 0 && boundary_velocity_manager.get_active_boundary_velocity_conditions().size() > 0))
rebuild_stokes_matrix = rebuild_stokes_preconditioner = assemble_newton_stokes_matrix = true;
else if (parameters.enable_prescribed_dilation)
// The dilation requires the Stokes matrix (which is on the rhs
// in the Newton solver) to be updated.
rebuild_stokes_matrix = true;

assemble_stokes_system();

assemble_stokes_system();
/**
* Eisenstat Walker method for determining the tolerance
*/
if (nonlinear_iteration > 1)
{
dcr.residual_old = dcr.residual;
dcr.velocity_residual = system_rhs.block(introspection.block_indices.velocities).l2_norm();
dcr.pressure_residual = system_rhs.block(introspection.block_indices.pressure).l2_norm();
dcr.residual = std::sqrt(dcr.velocity_residual * dcr.velocity_residual + dcr.pressure_residual * dcr.pressure_residual);

/**
* Eisenstat Walker method for determining the tolerance
*/
if (nonlinear_iteration > 1)
if (!use_picard)
{
dcr.residual_old = dcr.residual;
dcr.velocity_residual = system_rhs.block(introspection.block_indices.velocities).l2_norm();
dcr.pressure_residual = system_rhs.block(introspection.block_indices.pressure).l2_norm();
dcr.residual = std::sqrt(dcr.velocity_residual * dcr.velocity_residual + dcr.pressure_residual * dcr.pressure_residual);

if (!use_picard)
{
const bool EisenstatWalkerChoiceOne = true;
parameters.linear_stokes_solver_tolerance = compute_Eisenstat_Walker_linear_tolerance(EisenstatWalkerChoiceOne,
newton_handler->parameters.maximum_linear_stokes_solver_tolerance,
parameters.linear_stokes_solver_tolerance,
dcr.stokes_residuals.second,
dcr.residual,
dcr.residual_old);

pcout << " The linear solver tolerance is set to " << parameters.linear_stokes_solver_tolerance << std::endl;
}
const bool EisenstatWalkerChoiceOne = true;
parameters.linear_stokes_solver_tolerance = compute_Eisenstat_Walker_linear_tolerance(EisenstatWalkerChoiceOne,
newton_handler->parameters.maximum_linear_stokes_solver_tolerance,
parameters.linear_stokes_solver_tolerance,
dcr.stokes_residuals.second,
dcr.residual,
dcr.residual_old);

pcout << " The linear solver tolerance is set to " << parameters.linear_stokes_solver_tolerance << std::endl;
}
}

if (stokes_matrix_free)
stokes_matrix_free->build_preconditioner();
else
build_stokes_preconditioner();
if (stokes_matrix_free)
stokes_matrix_free->build_preconditioner();
else
build_stokes_preconditioner();

// Give this another try:
dcr.stokes_residuals = solve_stokes();
}
// Give this another try:
dcr.stokes_residuals = solve_stokes();
}

dcr.velocity_residual = system_rhs.block(introspection.block_indices.velocities).l2_norm();
Expand Down

0 comments on commit 2a805d7

Please sign in to comment.