Skip to content

Commit

Permalink
Merge pull request #6133 from gassmoeller/code_cleanup
Browse files Browse the repository at this point in the history 
Further newton code cleanup
  • Loading branch information
@gassmoeller
gassmoeller authored Nov 7, 2024
2 parents e7745ad + 1ebc54d commit c3b9737
Showing 1 changed file with 31 additions and 55 deletions.
86 changes: 31 additions & 55 deletions source/simulator/solver_schemes.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -435,7 +435,8 @@ namespace aspect
const unsigned int pressure_block_index = (parameters.include_melt_transport) ?
introspection.variable("fluid pressure").block_index
: introspection.block_indices.pressure;
Assert(introspection.block_indices.velocities == 0, ExcNotImplemented());
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;
Assert(!parameters.include_melt_transport
Expand All @@ -447,15 +448,14 @@ namespace aspect
dcr.switch_initial_residual = dcr.initial_residual;
dcr.residual_old = dcr.initial_residual;
dcr.residual = dcr.initial_residual;
}

assemble_newton_stokes_system = assemble_newton_stokes_matrix = true;

if (nonlinear_iteration == 0)
{
assemble_newton_stokes_system = assemble_newton_stokes_matrix = false;
}
else
assemble_newton_stokes_system = assemble_newton_stokes_matrix = true;

// Make sure we use the correct assemblers for first and all other iterations
// We solve the normal system in the first iteration, and the defect correction
// system in all other iterations.
if (nonlinear_iteration <= 1)
{
set_assemblers();
Expand All @@ -481,8 +481,8 @@ namespace aspect
assemble_stokes_system();

// recompute rhs
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.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);

// Test whether the rhs has dropped so much that we can assume that the iteration is done.
Expand Down Expand Up @@ -573,8 +573,8 @@ namespace aspect
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.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);

if (!use_picard)
Expand All @@ -600,22 +600,14 @@ namespace aspect
dcr.stokes_residuals = solve_stokes();
}

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);

double test_residual = dcr.residual;
if (nonlinear_iteration == 0)
{
/**
* The first nonlinear iteration we are computing the whole system in a non-defect corrected Picard way,
* to make sure that the boundary conditions are correct in combination with correct initial guesses.
*/
current_linearization_point.block(introspection.block_indices.velocities) = solution.block(introspection.block_indices.velocities);
current_linearization_point.block(introspection.block_indices.pressure) = solution.block(introspection.block_indices.pressure);
// The first nonlinear iteration we are computing the whole system in a non-defect corrected Picard way,
// to make sure that the boundary conditions are correct in combination with correct initial guesses.
current_linearization_point.block(velocity_block_index) = solution.block(velocity_block_index);
current_linearization_point.block(pressure_block_index) = solution.block(pressure_block_index);

dcr.residual = dcr.stokes_residuals.first;

pcout << " Newton system information: Norm of the rhs: " << dcr.stokes_residuals.first << std::endl;
}
else
Expand All @@ -629,32 +621,18 @@ namespace aspect
* is met, or the line search iteration limit is reached.
*/
LinearAlgebra::BlockVector backup_linearization_point(introspection.index_sets.stokes_partitioning, mpi_communicator);
backup_linearization_point.block(introspection.block_indices.pressure) = current_linearization_point.block(introspection.block_indices.pressure);
backup_linearization_point.block(introspection.block_indices.velocities) = current_linearization_point.block(introspection.block_indices.velocities);

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);

double test_velocity_residual = 0;
double test_pressure_residual = 0;
double step_length_factor = 1;
double alpha = 1e-4;
double step_length_factor = 1.0;
unsigned int line_search_iteration = 0;


/**
* Do the loop for the line search. Even when we
* don't do a line search we go into this loop
*/
// Do the loop for the line search at least once
do
{
// Reset the current linearization point and the search direction
current_linearization_point.block(introspection.block_indices.pressure) = backup_linearization_point.block(introspection.block_indices.pressure);
current_linearization_point.block(introspection.block_indices.velocities) = backup_linearization_point.block(introspection.block_indices.velocities);

LinearAlgebra::BlockVector search_direction = solution;
search_direction *= step_length_factor;

current_linearization_point.block(introspection.block_indices.pressure) += search_direction.block(introspection.block_indices.pressure);
current_linearization_point.block(introspection.block_indices.velocities) += search_direction.block(introspection.block_indices.velocities);
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);
current_linearization_point.block(pressure_block_index).add(step_length_factor, solution.block(pressure_block_index));
current_linearization_point.block(velocity_block_index).add(step_length_factor, solution.block(velocity_block_index));

// Rebuild the rhs to determine the new residual.
assemble_newton_stokes_matrix = rebuild_stokes_preconditioner = false;
Expand All @@ -663,12 +641,13 @@ namespace aspect

assemble_stokes_system();

test_velocity_residual = system_rhs.block(introspection.block_indices.velocities).l2_norm();
test_pressure_residual = system_rhs.block(introspection.block_indices.pressure).l2_norm();
test_residual = std::sqrt(test_velocity_residual * test_velocity_residual
+ test_pressure_residual * test_pressure_residual);
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 decrease is sufficient.
// 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
Expand All @@ -682,7 +661,6 @@ namespace aspect
}
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;
Expand All @@ -696,9 +674,7 @@ namespace aspect

++line_search_iteration;
Assert(line_search_iteration <= newton_handler->parameters.max_newton_line_search_iterations,
ExcMessage ("This tests the while condition. This condition should "
"actually never be false, because the break statement "
"above should have caught it."));
ExcInternalError());
}
while (true);
}
Expand All @@ -721,7 +697,7 @@ namespace aspect
* by Raid Hassani.
*/
if (newton_handler->parameters.use_newton_residual_scaling_method)
dcr.newton_residual_for_derivative_scaling_factor = test_residual;
dcr.newton_residual_for_derivative_scaling_factor = dcr.residual;
else
dcr.newton_residual_for_derivative_scaling_factor = 0;
}
Expand Down

0 comments on commit c3b9737

Please sign in to comment.