Skip to content

Refactored line search functionality for do_one_defect_correction_Stokes_step in solver_schemes.cc #6772

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
arnoldkk13 wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Refactored line search functionality for do_one_defect_correction_Stokes_step in solver_schemes.cc #6772

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 20 additions & 0 deletions include/aspect/simulator.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1337,6 +1337,26 @@ namespace aspect
void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
LinearAlgebra::Vector &vec) const;

/**
* Perform a Newton line search to determine the optimal step length
* along the search direction. After the update, the current_linearization_point
* is set to the old current_linearlization_point plus a suitable update.
*
* @param dcr The defect correction residuals associated with the current nonlinear
* iteration.
* @param use_picard Whether a Picard iteration was used to update the nonlinear
* iteration (true) or a Newton update (false).
* @param search_direction The proposed update direction for the solution vector.
*
* @return This function returns the updated residual after the line
* search is performed.
*
* This function is implemented in
* <code>source/simulator/helper_functions.cc</code>
*/
double perform_line_search(const DefectCorrectionResiduals &dcr,
const bool use_picard,
LinearAlgebra::BlockVector &search_direction);

/**
* Add constraints to the given @p constraints object that are required
Expand Down
90 changes: 90 additions & 0 deletions source/simulator/helper_functions.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2631,6 +2631,93 @@ namespace aspect



template <int dim>
double Simulator<dim>::perform_line_search(const DefectCorrectionResiduals &dcr,
const bool use_picard,
LinearAlgebra::BlockVector &search_direction)
{
LinearAlgebra::BlockVector backup_linearization_point(introspection.index_sets.stokes_partitioning, mpi_communicator);
double step_length_factor = 1.0;
unsigned int line_search_iteration = 0;

// Many parts of the solver depend on the block layout (velocity = 0,
// pressure = 1). For example the linearized_stokes_initial_guess vector or the StokesBlock matrix
// wrapper. Let us make sure that this holds (and shorten their names):
const unsigned int pressure_block_index = (parameters.include_melt_transport) ?
introspection.variable("fluid pressure").block_index
: introspection.block_indices.pressure;
const unsigned int velocity_block_index = introspection.block_indices.velocities;
Assert(velocity_block_index == 0, ExcNotImplemented());
Assert(pressure_block_index == 1, ExcNotImplemented());
(void) pressure_block_index;

// Do the loop for the line search at least once with the full step length.
// If line search is disabled we will exit the loop in the first iteration.
do
{
if (line_search_iteration == 0)
{
// backup our starting point for the line search
backup_linearization_point.block(pressure_block_index) = current_linearization_point.block(pressure_block_index);
backup_linearization_point.block(velocity_block_index) = current_linearization_point.block(velocity_block_index);
}
else
{
// undo the last iteration and try again with smaller step length
current_linearization_point.block(pressure_block_index) = backup_linearization_point.block(pressure_block_index);
current_linearization_point.block(velocity_block_index) = backup_linearization_point.block(velocity_block_index);
search_direction.block(pressure_block_index) *= step_length_factor;
search_direction.block(velocity_block_index) *= step_length_factor;
}

// Update the current linearization point with the search direction
current_linearization_point.block(pressure_block_index) += search_direction.block(pressure_block_index);
current_linearization_point.block(velocity_block_index) += search_direction.block(velocity_block_index);

// Rebuild the rhs to determine the new residual.
assemble_newton_stokes_matrix = rebuild_stokes_preconditioner = false;
rebuild_stokes_matrix = !boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().empty();

assemble_stokes_system();

const double test_velocity_residual = system_rhs.block(velocity_block_index).l2_norm();
const double test_pressure_residual = system_rhs.block(pressure_block_index).l2_norm();
const double test_residual = std::sqrt(test_velocity_residual * test_velocity_residual
+ test_pressure_residual * test_pressure_residual);

// Determine if the residual has decreased sufficiently.
const double alpha = 1e-4;
if (test_residual < (1.0 - alpha * step_length_factor) * dcr.residual
||
line_search_iteration >= newton_handler->parameters.max_newton_line_search_iterations
||
use_picard)
{
pcout << " Newton system information: Norm of the rhs: " << test_residual
<< ", Derivative scaling factor: " << newton_handler->parameters.newton_derivative_scaling_factor << std::endl;
return test_residual;
}
else
{
pcout << " Line search iteration " << line_search_iteration << ", with norm of the rhs "
<< test_residual << " and going to " << (1.0 - alpha * step_length_factor) * dcr.residual
<< ", relative residual: " << test_residual/dcr.initial_residual << std::endl;

// The current search direction has not decreased the residual
// enough, so we take a smaller step and try again.
// TODO: make a parameter out of this.
step_length_factor = 2.0/3.0;
}

++line_search_iteration;
Assert(line_search_iteration <= newton_handler->parameters.max_newton_line_search_iterations,
ExcInternalError());
}
while (true);
}



template <int dim>
double
Simulator<dim>::compute_Eisenstat_Walker_linear_tolerance(const bool EisenstatWalkerChoiceOne,
Expand Down Expand Up @@ -2779,6 +2866,9 @@ namespace aspect
template void Simulator<dim>::restore_outflow_boundary_ids(const unsigned int boundary_id_offset); \
template void Simulator<dim>::check_consistency_of_boundary_conditions() const; \
template double Simulator<dim>::compute_initial_newton_residual(); \
template double Simulator<dim>::perform_line_search(const DefectCorrectionResiduals &dcr, \
const bool use_picard, \
LinearAlgebra::BlockVector &search_direction); \
template double Simulator<dim>::compute_Eisenstat_Walker_linear_tolerance(const bool EisenstatWalkerChoiceOne, \
const double maximum_linear_stokes_solver_tolerance, \
const double linear_stokes_solver_tolerance, \
Expand Down
167 changes: 52 additions & 115 deletions source/simulator/solver_schemes.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -530,6 +530,7 @@ namespace aspect
}



template <int dim>
void Simulator<dim>::do_one_defect_correction_Stokes_step(DefectCorrectionResiduals &dcr,
const bool use_picard)
Expand All @@ -547,6 +548,49 @@ namespace aspect
Assert(!parameters.include_melt_transport
|| introspection.variable("compaction pressure").block_index == 1, ExcNotImplemented());

// Define some lambda functions here for readability
auto update_residuals = [&]()
{
dcr.velocity_residual = system_rhs.block(velocity_block_index).l2_norm();
dcr.pressure_residual = system_rhs.block(pressure_block_index).l2_norm();
dcr.residual = std::sqrt(dcr.velocity_residual * dcr.velocity_residual +
dcr.pressure_residual * dcr.pressure_residual);
};

// Eisenstat Walker method for determining the linear solver tolerance
auto update_Eisenstat_Walker_tolerance = [&]()
{
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
<< ". ";
if (!use_picard)
{
pcout << "Stabilization Preconditioner is "
<< Newton::to_string(newton_handler->parameters.preconditioner_stabilization)
<< " and A block is "
<< Newton::to_string(newton_handler->parameters.velocity_block_stabilization)
<< '.';
}
pcout << std::endl;
};

// build the preconditioner
auto build_preconditioner = [&]()
{
if (stokes_matrix_free)
stokes_matrix_free->build_preconditioner();
else
build_stokes_preconditioner();
};

// Re-compute the pressure scaling factor for the Stokes assembly
pressure_scaling = compute_pressure_scaling_factor();

Expand Down Expand Up @@ -588,10 +632,7 @@ namespace aspect

assemble_stokes_system();

// recompute rhs
dcr.velocity_residual = system_rhs.block(velocity_block_index).l2_norm();
dcr.pressure_residual = system_rhs.block(pressure_block_index).l2_norm();
dcr.residual = std::sqrt(dcr.velocity_residual * dcr.velocity_residual + dcr.pressure_residual * dcr.pressure_residual);
update_residuals();

// Test whether the rhs has dropped so much that we can assume that the iteration is done.
if (dcr.residual < dcr.residual_old * 1e-8)
Expand All @@ -607,38 +648,15 @@ namespace aspect
}


// Eisenstat Walker method for determining the linear solver tolerance
if (nonlinear_iteration > 1)
{
if (!use_picard || newton_handler->parameters.use_Eisenstat_Walker_method_for_Picard_iterations)
{
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
<< ". ";
if (!use_picard)
{
pcout << "Stabilization Preconditioner is "
<< Newton::to_string(newton_handler->parameters.preconditioner_stabilization)
<< " and A block is "
<< Newton::to_string(newton_handler->parameters.velocity_block_stabilization)
<< '.';
}
pcout << std::endl;
update_Eisenstat_Walker_tolerance();
}
}

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

// The solution of the defect correction solver is an update direction
LinearAlgebra::BlockVector search_direction;
Expand Down Expand Up @@ -694,29 +712,16 @@ namespace aspect
*/
if (nonlinear_iteration > 1)
{
dcr.velocity_residual = system_rhs.block(velocity_block_index).l2_norm();
dcr.pressure_residual = system_rhs.block(pressure_block_index).l2_norm();
dcr.residual = std::sqrt(dcr.velocity_residual * dcr.velocity_residual + dcr.pressure_residual * dcr.pressure_residual);
update_residuals();

// Eisenstat Walker method for determining the linear solver tolerance
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;
update_Eisenstat_Walker_tolerance();
}
}

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

// Give this another try:
dcr.stokes_residuals = solve_stokes(search_direction);
Expand Down Expand Up @@ -744,79 +749,11 @@ namespace aspect
dcr.residual_old = dcr.residual;

// We have computed an update. Prepare for a line search.
LinearAlgebra::BlockVector backup_linearization_point(introspection.index_sets.stokes_partitioning, mpi_communicator);

double step_length_factor = 1.0;
unsigned int line_search_iteration = 0;

// Do the loop for the line search at least once with the full step length.
// If line search is disabled we will exit the loop in the first iteration.
do
{
if (line_search_iteration == 0)
{
// backup our starting point for the line search
backup_linearization_point.block(pressure_block_index) = current_linearization_point.block(pressure_block_index);
backup_linearization_point.block(velocity_block_index) = current_linearization_point.block(velocity_block_index);
}
else
{
// undo the last iteration and try again with smaller step length
current_linearization_point.block(pressure_block_index) = backup_linearization_point.block(pressure_block_index);
current_linearization_point.block(velocity_block_index) = backup_linearization_point.block(velocity_block_index);
search_direction.block(pressure_block_index) *= step_length_factor;
search_direction.block(velocity_block_index) *= step_length_factor;
}

// Update the current linearization point with the search direction
current_linearization_point.block(pressure_block_index) += search_direction.block(pressure_block_index);
current_linearization_point.block(velocity_block_index) += search_direction.block(velocity_block_index);

// Rebuild the rhs to determine the new residual.
assemble_newton_stokes_matrix = rebuild_stokes_preconditioner = false;
rebuild_stokes_matrix = !boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().empty();

assemble_stokes_system();

const double test_velocity_residual = system_rhs.block(velocity_block_index).l2_norm();
const double test_pressure_residual = system_rhs.block(pressure_block_index).l2_norm();
const double test_residual = std::sqrt(test_velocity_residual * test_velocity_residual
+ test_pressure_residual * test_pressure_residual);

// Determine if the residual has decreased sufficiently.
const double alpha = 1e-4;
if (test_residual < (1.0 - alpha * step_length_factor) * dcr.residual
||
line_search_iteration >= newton_handler->parameters.max_newton_line_search_iterations
||
use_picard)
{
pcout << " Newton system information: Norm of the rhs: " << test_residual
<< ", Derivative scaling factor: " << newton_handler->parameters.newton_derivative_scaling_factor << std::endl;
dcr.residual = test_residual;
break;
}
else
{
pcout << " Line search iteration " << line_search_iteration << ", with norm of the rhs "
<< test_residual << " and going to " << (1.0 - alpha * step_length_factor) * dcr.residual
<< ", relative residual: " << test_residual/dcr.initial_residual << std::endl;

// The current search direction has not decreased the residual
// enough, so we take a smaller step and try again.
// TODO: make a parameter out of this.
step_length_factor = 2.0/3.0;
}

++line_search_iteration;
Assert(line_search_iteration <= newton_handler->parameters.max_newton_line_search_iterations,
ExcInternalError());
}
while (true);
dcr.residual = perform_line_search(dcr, use_picard, search_direction);
}


if (use_picard == true)
if (use_picard)
{
// When we are using (defect corrected) Picard, keep the
// newton_derivative_scaling_factor at zero. The newton_derivative_scaling_factor
Expand Down