Skip to content
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

use sundials kinsol for diffusion dislocation creep #5698

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open

use sundials kinsol for diffusion dislocation creep #5698

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
3ab6379
use sundials kinsol for diffusion dislocation creep
bobmyhill Jun 1, 2024
fa84fc7
added changelog
bobmyhill Jun 1, 2024
fb9cb6d
updated tests
bobmyhill Jun 1, 2024
26cfcc0
removed superfluous docstrings
bobmyhill Jun 1, 2024
39f923e
merge residual and jacobian calculation
bobmyhill Jun 1, 2024
3b01b3a
addressed comments
bobmyhill Jun 6, 2024
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
4 changes: 4 additions & 0 deletions doc/modules/changes/20240531_myhill
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
Changed: The DiffusionDislocation rheology now uses SUNDIALS KINSOL
routines to calculate the equilibrium partitioning of strain rates.
<br>
(Bob Myhill, 2024/05/31)
17 changes: 17 additions & 0 deletions include/aspect/material_model/rheology/diffusion_dislocation.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,8 @@
#include <aspect/material_model/rheology/diffusion_creep.h>
#include <aspect/material_model/rheology/dislocation_creep.h>

#include <deal.II/sundials/kinsol.h>

namespace aspect
{
namespace MaterialModel
Expand Down Expand Up @@ -90,6 +92,20 @@ namespace aspect
const double temperature,
const SymmetricTensor<2,dim> &strain_rate) const;

/**
* Compute the residual of the natural logarithm of the strain rate
* and the derivative with respect to the logarithm of the stress
* at the current value of the logarithm of the stress.
*/
std::pair<double, double>
compute_log_strain_rate_residual_and_derivative(
const double current_log_stress_ii,
const double pressure,
const double temperature,
const Rheology::DiffusionCreepParameters &diffusion_creep_parameters,
const Rheology::DislocationCreepParameters &dislocation_creep_parameters,
const double log_edot_ii) const;

/**
* Compute the viscosity based on the composite viscous creep law,
* averaging over all compositional fields according to their
Expand Down Expand Up @@ -130,6 +146,7 @@ namespace aspect
unsigned int n_chemical_composition_fields;

MaterialUtilities::CompositionalAveragingOperation viscosity_averaging;
SUNDIALS::KINSOL<Vector<double>>::AdditionalData kinsol_settings;

};

Expand Down
186 changes: 84 additions & 102 deletions source/material_model/rheology/diffusion_dislocation.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,13 +55,15 @@ namespace aspect
min_strain_rate);
const double log_edot_ii = std::log(edot_ii);


// Find effective viscosities for each of the individual phases
// Viscosities should have same number of entries as compositional fields
// Viscosities should have the same number of entries as compositional fields
std::vector<double> composition_viscosities(n_chemical_composition_fields);
for (unsigned int j=0; j < n_chemical_composition_fields; ++j)
{
// Power law creep equation
// Because the ratios of the diffusion and dislocation strain rates are not known, stress is also unknown
// We use KINSOL to find the second invariant of the stress tensor that satisfies the total strain rate.

// The power law creep equation is
// edot_ii_i = A_i * stress_ii_i^{n_i} * d^{-m} \exp\left(-\frac{E_i^\ast + PV_i^\ast}{n_iRT}\right)
// where ii indicates the square root of the second invariant and
// i corresponds to diffusion or dislocation creep
Expand All @@ -72,115 +74,90 @@ namespace aspect
// For dislocation creep, viscosity is grain size independent (m=0)
const Rheology::DislocationCreepParameters dislocation_creep_parameters = dislocation_creep.compute_creep_parameters(j);

// For diffusion creep, viscosity is grain size dependent

SUNDIALS::KINSOL<Vector<double>> nonlinear_solver(kinsol_settings);
double log_strain_rate_deriv;
nonlinear_solver.reinit_vector = [&](Vector<double> &x)
{
x.reinit(1);
};

nonlinear_solver.residual =
[&](const Vector<double> &current_log_stress_ii,
Vector<double> &residual)
{
std::tie(residual(0), log_strain_rate_deriv) = compute_log_strain_rate_residual_and_derivative(
current_log_stress_ii[0], pressure, temperature, diffusion_creep_parameters, dislocation_creep_parameters, log_edot_ii
);
};

nonlinear_solver.setup_jacobian =
[&](const Vector<double> & /*current_log_stress_ii*/,
const Vector<double> & /*current_f*/)
{
// Do nothing here, because we calculate the Jacobian in the residual function
};

nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &residual,
Vector<double> &solution,
const double /*tolerance*/)
{
solution(0) = residual(0) / log_strain_rate_deriv;
};


// Our starting guess for the stress assumes that all strain is
// accommodated by either diffusion creep, dislocation creep,
// or by the maximum viscosity limiter.
const double prefactor_stress_diffusion = diffusion_creep_parameters.prefactor *
std::pow(grain_size, -diffusion_creep_parameters.grain_size_exponent) *
std::exp(-(std::max(diffusion_creep_parameters.activation_energy + pressure*diffusion_creep_parameters.activation_volume,0.0))/
(constants::gas_constant*temperature));

// Because the ratios of the diffusion and dislocation strain rates are not known, stress is also unknown
// We use Newton's method to find the second invariant of the stress tensor.
// Start with the assumption that all strain is accommodated by diffusion creep:
// If the diffusion creep prefactor is very small, that means that the diffusion viscosity is very large.
// In this case, use the maximum viscosity instead to compute the starting guess.
double stress_ii = (prefactor_stress_diffusion > (0.5 / maximum_viscosity)
?
edot_ii/prefactor_stress_diffusion
:
0.5 / maximum_viscosity);
double log_stress_ii = std::log(stress_ii);
double log_strain_rate_residual = 2 * log_strain_rate_residual_threshold;
double log_strain_rate_deriv = 0;
unsigned int stress_iteration = 0;
while (std::abs(log_strain_rate_residual) > log_strain_rate_residual_threshold
&& stress_iteration < stress_max_iteration_number)
Comment on lines -95 to -96
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is your current stopping criterion. Could you check how many iterations you are using in this loop (in total or on average) in one of your benchmarks, and compare that against the number of iterations you're spending in KINSOLS? I think that might clarify why the code is so much slower.

You can also directly time individual sections by querying simulator_access.get_computing_timer(). You can take what's in particles/world.cc as an example of how you can put a timing section around a piece of code.

{
const std::pair<double, double> log_diff_edot_and_deriv = diffusion_creep.compute_log_strain_rate_and_derivative(log_stress_ii, pressure, temperature, diffusion_creep_parameters);
const std::pair<double, double> log_disl_edot_and_deriv = dislocation_creep.compute_log_strain_rate_and_derivative(log_stress_ii, pressure, temperature, dislocation_creep_parameters);

const double strain_rate_diffusion = std::exp(log_diff_edot_and_deriv.first);
const double strain_rate_dislocation = std::exp(log_disl_edot_and_deriv.first);
log_strain_rate_residual = std::log(strain_rate_diffusion + strain_rate_dislocation) - log_edot_ii;
log_strain_rate_deriv = (strain_rate_diffusion * log_diff_edot_and_deriv.second + strain_rate_dislocation * log_disl_edot_and_deriv.second)/
(strain_rate_diffusion + strain_rate_dislocation);

// If the log strain rate derivative is zero, we catch it below.
if (log_strain_rate_deriv>std::numeric_limits<double>::min())
log_stress_ii -= log_strain_rate_residual/log_strain_rate_deriv;
stress_ii = std::exp(log_stress_ii);
stress_iteration += 1;

// In case the Newton iteration does not succeed, we do a fixpoint iteration.
// This allows us to bound both the diffusion and dislocation viscosity
// between a minimum and maximum value, so that we can compute the correct
// viscosity values even if the parameters lead to one or both of the
// viscosities being essentially zero or infinity.
// If anything that would be used in the next iteration is not finite, the
// Newton iteration would trigger an exception and we want to do the fixpoint
// iteration instead.
const bool abort_newton_iteration = !numbers::is_finite(stress_ii)
|| !numbers::is_finite(log_strain_rate_residual)
|| !numbers::is_finite(log_strain_rate_deriv)
|| log_strain_rate_deriv < std::numeric_limits<double>::min()
|| !numbers::is_finite(std::pow(stress_ii, diffusion_creep_parameters.stress_exponent-1))
|| !numbers::is_finite(std::pow(stress_ii, dislocation_creep_parameters.stress_exponent-1))
|| stress_iteration == stress_max_iteration_number;
if (abort_newton_iteration)
{
double diffusion_strain_rate = edot_ii;
double dislocation_strain_rate = min_strain_rate;
stress_iteration = 0;

do
{
const double old_diffusion_strain_rate = diffusion_strain_rate;

const double diffusion_prefactor = 0.5 * std::pow(diffusion_creep_parameters.prefactor,-1.0/diffusion_creep_parameters.stress_exponent);
const double diffusion_grain_size_dependence = std::pow(grain_size, diffusion_creep_parameters.grain_size_exponent/diffusion_creep_parameters.stress_exponent);
const double diffusion_strain_rate_dependence = std::pow(diffusion_strain_rate, (1.-diffusion_creep_parameters.stress_exponent)/diffusion_creep_parameters.stress_exponent);
const double diffusion_T_and_P_dependence = std::exp(std::max(diffusion_creep_parameters.activation_energy + pressure*diffusion_creep_parameters.activation_volume,0.0)/
(constants::gas_constant*temperature));

const double diffusion_viscosity = std::min(std::max(diffusion_prefactor * diffusion_grain_size_dependence
* diffusion_strain_rate_dependence * diffusion_T_and_P_dependence,
minimum_viscosity), maximum_viscosity);

const double dislocation_prefactor = 0.5 * std::pow(dislocation_creep_parameters.prefactor,-1.0/dislocation_creep_parameters.stress_exponent);
const double dislocation_strain_rate_dependence = std::pow(dislocation_strain_rate, (1.-dislocation_creep_parameters.stress_exponent)/dislocation_creep_parameters.stress_exponent);
const double dislocation_T_and_P_dependence = std::exp(std::max(dislocation_creep_parameters.activation_energy + pressure*dislocation_creep_parameters.activation_volume,0.0)/
(dislocation_creep_parameters.stress_exponent*constants::gas_constant*temperature));

const double dislocation_viscosity = std::min(std::max(dislocation_prefactor * dislocation_strain_rate_dependence
* dislocation_T_and_P_dependence,
minimum_viscosity), maximum_viscosity);

diffusion_strain_rate = dislocation_viscosity / (diffusion_viscosity + dislocation_viscosity) * edot_ii;
dislocation_strain_rate = diffusion_viscosity / (diffusion_viscosity + dislocation_viscosity) * edot_ii;

++stress_iteration;
AssertThrow(stress_iteration < stress_max_iteration_number,
ExcMessage("No convergence has been reached in the loop that determines "
"the ratio of diffusion/dislocation viscosity. Aborting! "
"Residual is " + Utilities::to_string(log_strain_rate_residual) +
" after " + Utilities::to_string(stress_iteration) + " iterations. "
"You can increase the number of iterations by adapting the "
"parameter 'Maximum strain rate ratio iterations'."));

log_strain_rate_residual = std::abs(std::log(diffusion_strain_rate/old_diffusion_strain_rate));
stress_ii = 2.0 * edot_ii * 1./(1./diffusion_viscosity + 1./dislocation_viscosity);
}
while (log_strain_rate_residual > log_strain_rate_residual_threshold);

break;
}
}

// The effective viscosity, with minimum and maximum bounds
const double prefactor_stress_dislocation = dislocation_creep_parameters.prefactor *
std::exp(-(std::max(dislocation_creep_parameters.activation_energy + pressure*dislocation_creep_parameters.activation_volume,0.0))/
(constants::gas_constant*temperature));

double stress_ii = std::min(
{
std::pow(edot_ii/prefactor_stress_diffusion, 1./diffusion_creep_parameters.stress_exponent),
std::pow(edot_ii/prefactor_stress_dislocation, 1./dislocation_creep_parameters.stress_exponent),
0.5 / maximum_viscosity
});

// KINSOL works on vectors and so the scalar (log) stress is inserted into
// a vector of length 1
Vector<double> log_stress_ii(1);
log_stress_ii[0] = std::log(stress_ii);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is your starting guess. How is it chosen? Based on the current log stress? It might be worth trying the output of the previous loop iteration, under the assumption that what you found in the previous iteration is a good starting guess for the current iteration. Not sure whether that's true, since you're not looping over quadrature points, for example, but a thought.


nonlinear_solver.solve(log_stress_ii);
stress_ii = std::exp(log_stress_ii[0]);
composition_viscosities[j] = std::min(std::max(stress_ii/edot_ii/2, minimum_viscosity), maximum_viscosity);
}
return composition_viscosities;
}

template <int dim>
std::pair<double, double>
DiffusionDislocation<dim>::compute_log_strain_rate_residual_and_derivative(
const double current_log_stress_ii,
const double pressure,
const double temperature,
const Rheology::DiffusionCreepParameters &diffusion_creep_parameters,
const Rheology::DislocationCreepParameters &dislocation_creep_parameters,
const double log_edot_ii) const
{
const std::pair<double, double> log_diff_edot_and_deriv = diffusion_creep.compute_log_strain_rate_and_derivative(current_log_stress_ii, pressure, temperature, diffusion_creep_parameters);
const std::pair<double, double> log_disl_edot_and_deriv = dislocation_creep.compute_log_strain_rate_and_derivative(current_log_stress_ii, pressure, temperature, dislocation_creep_parameters);

const double strain_rate_diffusion = std::exp(log_diff_edot_and_deriv.first);
const double strain_rate_dislocation = std::exp(log_disl_edot_and_deriv.first);
double log_strain_rate_deriv = (strain_rate_diffusion * log_diff_edot_and_deriv.second + strain_rate_dislocation * log_disl_edot_and_deriv.second)/
(strain_rate_diffusion + strain_rate_dislocation);
const double log_strain_rate_iterate = std::log(strain_rate_diffusion + strain_rate_dislocation);
return std::make_pair(log_strain_rate_iterate - log_edot_ii, log_strain_rate_deriv);
}

template <int dim>
double
DiffusionDislocation<dim>::compute_viscosity (const double pressure,
Expand Down Expand Up @@ -363,6 +340,11 @@ namespace aspect
dislocation_creep.initialize_simulator (this->get_simulator());
dislocation_creep.parse_parameters(prm, std::make_unique<std::vector<unsigned int>>(n_chemical_composition_fields));

// KINSOL settings
kinsol_settings.strategy = dealii::SUNDIALS::KINSOL<>::AdditionalData::newton;
kinsol_settings.function_tolerance = log_strain_rate_residual_threshold;
kinsol_settings.maximum_non_linear_iterations = stress_max_iteration_number;
kinsol_settings.maximum_setup_calls = 1;
}
}
}
Expand Down
12 changes: 6 additions & 6 deletions tests/diffusion_dislocation/screen-output
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,16 +13,16 @@ Number of degrees of freedom: 1,237 (578+81+289+289)
Solving temperature system... 0 iterations.
Solving depleted_lithosphere system ... 0 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 16+0 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.68721e-16, 1.05609e-16, 0.000715698
Relative nonlinear residual (total system) after nonlinear iteration 2: 0.000715698
Solving Stokes system... 21+0 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.68721e-16, 1.05609e-16, 0.000716824
Relative nonlinear residual (total system) after nonlinear iteration 2: 0.000716824
Comment on lines -16 to +18
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's a bummer that (at the moment) you need more iterations, but a good sign that the results are essentially identical.


Solving temperature system... 0 iterations.
Solving depleted_lithosphere system ... 0 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 8+0 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.68721e-16, 1.05609e-16, 5.22177e-08
Relative nonlinear residual (total system) after nonlinear iteration 3: 5.22177e-08
Solving Stokes system... 7+0 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.68721e-16, 1.05609e-16, 8.28348e-07
Relative nonlinear residual (total system) after nonlinear iteration 3: 8.28348e-07


Postprocessing:
Expand Down
2 changes: 1 addition & 1 deletion tests/diffusion_dislocation/statistics
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,4 +11,4 @@
# 11: Iterations for Stokes solver
# 12: Velocity iterations in Stokes preconditioner
# 13: Schur complement iterations in Stokes preconditioner
0 0.000000000000e+00 0.000000000000e+00 64 659 289 289 3 0 0 30 36 35
0 0.000000000000e+00 0.000000000000e+00 64 659 289 289 3 0 0 34 40 39
12 changes: 6 additions & 6 deletions tests/diffusion_dislocation_reference_density_profile/screen-output
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,16 +13,16 @@ Number of degrees of freedom: 1,237 (578+81+289+289)
Solving temperature system... 0 iterations.
Solving depleted_lithosphere system ... 0 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 16+0 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.15502e-16, 1.05609e-16, 0.000715698
Relative nonlinear residual (total system) after nonlinear iteration 2: 0.000715698
Solving Stokes system... 21+0 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.15502e-16, 1.05609e-16, 0.000716824
Relative nonlinear residual (total system) after nonlinear iteration 2: 0.000716824

Solving temperature system... 0 iterations.
Solving depleted_lithosphere system ... 0 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 8+0 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.15502e-16, 1.05609e-16, 5.22177e-08
Relative nonlinear residual (total system) after nonlinear iteration 3: 5.22177e-08
Solving Stokes system... 7+0 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.15502e-16, 1.05609e-16, 8.28348e-07
Relative nonlinear residual (total system) after nonlinear iteration 3: 8.28348e-07


Postprocessing:
Expand Down
Loading
Loading