Skip to content

Commit

Permalink
Split a variable into two, to indicate different purposes.
Browse files Browse the repository at this point in the history
  • Loading branch information
@bangerth
bangerth committed Jul 14, 2023
1 parent ca748ee commit e96a62b
Show file tree
Hide file tree
Showing 2 changed files with 66 additions and 56 deletions.
52 changes: 29 additions & 23 deletions include/aspect/simulator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -750,49 +750,55 @@ namespace aspect

/**
* Assemble and solve the temperature equation.
* This function returns the residual after solving
* and can optionally compute and store an initial
* residual before solving the equation in the argument,
* if the latter is not a `nullptr`.
* This function returns the residual after solving.
*
* If the `redidual` argument is not a `nullptr`, the function computes
* the residual and puts it into this variable. The function returns
* the current residual divided by the initial residual given as the
* first argument. The two arguments may point to the same variable,
* in which case the function first computes the residual and at
* the end scales that residual by itself, thus returning 1.0.
*
* This function is implemented in
* <code>source/simulator/solver_schemes.cc</code>.
*/
double assemble_and_solve_temperature (double *initial_residual = nullptr);
double assemble_and_solve_temperature (const double &initial_residual = 0,
double *residual = nullptr);

/**
* Solve the composition equations with whatever method is selected
* (fields or particles). This function returns the residuals for
* all fields after solving
* and can optionally compute and store the initial
* residuals before solving the equation in the argument,
* if the latter is not a `nullptr`. For lack of a definition
* the residuals of all compositional fields that are advected
* using particles are considered zero.
* all fields after solving.
*
* If the `redidual` argument is not a `nullptr`, the function computes
* the residual and puts it into this variable. The function returns
* the current residual divided by the initial residual given as the
* first argument. The two arguments may point to the same variable,
* in which case the function first computes the residual and at
* the end scales that residual by itself, thus returning 1.0.
*
* This function is implemented in
* <code>source/simulator/solver_schemes.cc</code>.
*/
std::vector<double> assemble_and_solve_composition (std::vector<double> *initial_residual = nullptr);
std::vector<double> assemble_and_solve_composition (const std::vector<double> &initial_residual = {},
std::vector<double> *residual = nullptr);

/**
* Assemble and solve the Stokes equation.
* This function returns the nonlinear residual after solving
* and can optionally compute and store an initial
* residual before solving the equation in the argument,
* if the latter is not a `nullptr`.
*
* The returned nonlinear residual is normalized by the initial
* residual, i.e., it is the nonlinear residual computed by
* solve_stokes() divided by the initial residual as either
* already stored in the second argument, or as computed
* at the top of the function.
* This function returns the nonlinear residual after solving.
*
* If the `redidual` argument is not a `nullptr`, the function computes
* the residual and puts it into this variable. The function returns
* the current residual divided by the initial residual given as the
* first argument. The two arguments may point to the same variable,
* in which case the function first computes the residual and at
* the end scales that residual by itself, thus returning 1.0.
*
* This function is implemented in
* <code>source/simulator/solver_schemes.cc</code>.
*/
double assemble_and_solve_stokes (double *initial_nonlinear_residual = nullptr);
double assemble_and_solve_stokes (const double &initial_nonlinear_residual = 0,
double *nonlinear_residual = nullptr);

/**
* Assemble and solve the defect correction form of the Stokes equation.
Expand Down
70 changes: 37 additions & 33 deletions source/simulator/solver_schemes.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,8 @@ namespace aspect


template <int dim>
double Simulator<dim>::assemble_and_solve_temperature (double *initial_residual)
double Simulator<dim>::assemble_and_solve_temperature (const double &initial_residual,
double *residual)
{
double current_residual = 0.0;

Expand All @@ -166,11 +167,8 @@ namespace aspect

assemble_advection_system (adv_field);

if (initial_residual)
{
Assert(initial_residual != nullptr, ExcInternalError());
*initial_residual = system_rhs.block(introspection.block_indices.temperature).l2_norm();
}
if (residual)
*residual = system_rhs.block(introspection.block_indices.temperature).l2_norm();

current_residual = solve_advection(adv_field);
break;
Expand Down Expand Up @@ -213,8 +211,8 @@ namespace aspect
current_linearization_point.block(introspection.block_indices.temperature)
= solution.block(introspection.block_indices.temperature);

if ((initial_residual != nullptr) && (*initial_residual > 0))
return current_residual / *initial_residual;
if (initial_residual > 0)
return current_residual / initial_residual;

return 0.0;
}
Expand All @@ -223,7 +221,8 @@ namespace aspect

template <int dim>
std::vector<double>
Simulator<dim>::assemble_and_solve_composition (std::vector<double> *initial_residual)
Simulator<dim>::assemble_and_solve_composition (const std::vector<double> &initial_residual,
std::vector<double> *residual)
{
// Advect the particles before they are potentially used to
// set up the compositional fields.
Expand All @@ -242,8 +241,8 @@ namespace aspect

std::vector<double> current_residual(introspection.n_compositional_fields,0.0);

if (initial_residual)
Assert(initial_residual->size() == introspection.n_compositional_fields, ExcInternalError());
if (residual)
Assert(residual->size() == introspection.n_compositional_fields, ExcInternalError());

std::vector<AdvectionField> fields_advected_by_particles;

Expand Down Expand Up @@ -274,8 +273,8 @@ namespace aspect

assemble_advection_system (adv_field);

if (initial_residual)
(*initial_residual)[c] = system_rhs.block(introspection.block_indices.compositional_fields[c]).l2_norm();
if (residual)
(*residual)[c] = system_rhs.block(introspection.block_indices.compositional_fields[c]).l2_norm();

current_residual[c] = solve_advection(adv_field);

Expand Down Expand Up @@ -346,8 +345,8 @@ namespace aspect
current_linearization_point.block(introspection.block_indices.compositional_fields[c])
= solution.block(introspection.block_indices.compositional_fields[c]);

if ((initial_residual != nullptr) && (*initial_residual)[c] > 0)
current_residual[c] /= (*initial_residual)[c];
if ((initial_residual.size() > 0) && (initial_residual[c] > 0))
current_residual[c] /= initial_residual[c];
else
current_residual[c] = 0.0;
}
Expand All @@ -359,7 +358,8 @@ namespace aspect

template <int dim>
double
Simulator<dim>::assemble_and_solve_stokes (double *initial_nonlinear_residual)
Simulator<dim>::assemble_and_solve_stokes (const double &initial_nonlinear_residual,
double *nonlinear_residual)
{
// If the Stokes matrix depends on the solution, or we have active
// velocity boundary conditions, we need to re-assemble the system matrix
Expand Down Expand Up @@ -396,11 +396,8 @@ namespace aspect
else
build_stokes_preconditioner();

if (initial_nonlinear_residual)
{
Assert(initial_nonlinear_residual != nullptr, ExcInternalError());
*initial_nonlinear_residual = compute_initial_stokes_residual();
}
if (nonlinear_residual)
*nonlinear_residual = compute_initial_stokes_residual();

const double current_nonlinear_residual = solve_stokes().first;

Expand All @@ -421,8 +418,8 @@ namespace aspect
current_linearization_point.block(fluid_pressure_block) = solution.block(fluid_pressure_block);
}

if ((initial_nonlinear_residual != nullptr) && (*initial_nonlinear_residual > 0))
return current_nonlinear_residual / *initial_nonlinear_residual;
if (initial_nonlinear_residual > 0)
return current_nonlinear_residual / initial_nonlinear_residual;
else
return 0.0;
}
Expand Down Expand Up @@ -792,7 +789,8 @@ namespace aspect
do
{
relative_residual =
assemble_and_solve_stokes(nonlinear_iteration == 0 ? &initial_stokes_residual : nullptr);
assemble_and_solve_stokes(initial_stokes_residual,
nonlinear_iteration == 0 ? &initial_stokes_residual : nullptr);

pcout << " Relative nonlinear residual (Stokes system) after nonlinear iteration " << nonlinear_iteration+1
<< ": " << relative_residual
Expand Down Expand Up @@ -989,10 +987,12 @@ namespace aspect
particle_world->restore_particles();

const double relative_temperature_residual =
assemble_and_solve_temperature(nonlinear_iteration == 0 ? &initial_temperature_residual : 0);
assemble_and_solve_temperature(initial_temperature_residual,
nonlinear_iteration == 0 ? &initial_temperature_residual : 0);

const std::vector<double> relative_composition_residual =
assemble_and_solve_composition(nonlinear_iteration == 0 ? &initial_composition_residual : 0);
assemble_and_solve_composition(initial_composition_residual,
nonlinear_iteration == 0 ? &initial_composition_residual : 0);

// write the residual output in the same order as the solutions
pcout << " Relative nonlinear residuals (temperature, compositional fields): " << relative_temperature_residual;
Expand Down Expand Up @@ -1109,13 +1109,16 @@ namespace aspect
particle_world->restore_particles();

const double relative_temperature_residual =
assemble_and_solve_temperature(nonlinear_iteration == 0 ? &initial_temperature_residual : nullptr);
assemble_and_solve_temperature(initial_temperature_residual,
nonlinear_iteration == 0 ? &initial_temperature_residual : nullptr);

const std::vector<double> relative_composition_residual =
assemble_and_solve_composition(nonlinear_iteration == 0 ? &initial_composition_residual : nullptr);
assemble_and_solve_composition(initial_composition_residual,
nonlinear_iteration == 0 ? &initial_composition_residual : nullptr);

const double relative_nonlinear_stokes_residual =
assemble_and_solve_stokes(nonlinear_iteration == 0 ? &initial_stokes_residual : nullptr);
assemble_and_solve_stokes(initial_stokes_residual,
nonlinear_iteration == 0 ? &initial_stokes_residual : nullptr);

// write the residual output in the same order as the solutions
pcout << " Relative nonlinear residuals (temperature, compositional fields, Stokes system): " << relative_temperature_residual;
Expand Down Expand Up @@ -1186,7 +1189,8 @@ namespace aspect
do
{
relative_residual =
assemble_and_solve_stokes(nonlinear_iteration == 0 ? &initial_stokes_residual : nullptr);
assemble_and_solve_stokes(initial_stokes_residual,
nonlinear_iteration == 0 ? &initial_stokes_residual : nullptr);

pcout << " Relative nonlinear residual (Stokes system) after nonlinear iteration " << nonlinear_iteration+1
<< ": " << relative_residual
Expand Down Expand Up @@ -1491,9 +1495,9 @@ namespace aspect
namespace aspect
{
#define INSTANTIATE(dim) \
template double Simulator<dim>::assemble_and_solve_temperature(double*); \
template std::vector<double> Simulator<dim>::assemble_and_solve_composition(std::vector<double> *); \
template double Simulator<dim>::assemble_and_solve_stokes(double*); \
template double Simulator<dim>::assemble_and_solve_temperature(const double &, double*); \
template std::vector<double> Simulator<dim>::assemble_and_solve_composition(const std::vector<double> &, std::vector<double> *); \
template double Simulator<dim>::assemble_and_solve_stokes(const double &, double*); \
template void Simulator<dim>::solve_single_advection_single_stokes(); \
template void Simulator<dim>::solve_no_advection_iterated_stokes(); \
template void Simulator<dim>::solve_no_advection_single_stokes(); \
Expand Down

0 comments on commit e96a62b

Please sign in to comment.