Skip to content

Commit 39f923e

Browse files
bobmyhillbobmyhill
committed
merge residual and jacobian calculation
1 parent 26cfcc0 commit 39f923e

File tree

2 files changed

+44
-73
lines changed

2 files changed

+44
-73
lines changed

‎include/aspect/material_model/rheology/diffusion_dislocation.h

+8-21
Original file line numberDiff line numberDiff line change
@@ -90,32 +90,19 @@ namespace aspect
9090
const double temperature,
9191
const SymmetricTensor<2,dim> &strain_rate) const;
9292

93-
/**
94-
* Compute the derivative of the natural logarithm of the strain rate
95-
* with respect to the natural logarithm of the stress
96-
* at the current evaluation point (log stress) estimated by KINSOL.
97-
*/
98-
void
99-
compute_log_strain_rate_deriv(const Vector<double> &evaluation_point,
100-
double pressure,
101-
double temperature,
102-
const Rheology::DiffusionCreepParameters diffusion_creep_parameters,
103-
const Rheology::DislocationCreepParameters dislocation_creep_parameters,
104-
double &log_strain_rate_deriv) const;
105-
10693
/**
10794
* Compute the residual of the natural logarithm of the strain rate
108-
* at the current evaluation point (log stress) estimated by KINSOL.
95+
* and the derivative with respect to the logarithm of the stress
96+
* at the current value of the logarithm of the stress.
10997
*/
110-
void
111-
compute_log_strain_rate_residual(
112-
const Vector<double> &evaluation_point,
113-
Vector<double> &residual,
114-
double pressure,
115-
double temperature,
98+
std::pair<double, double>
99+
compute_log_strain_rate_residual_and_derivative(
100+
const double current_log_stress_ii,
101+
const double pressure,
102+
const double temperature,
116103
const Rheology::DiffusionCreepParameters diffusion_creep_parameters,
117104
const Rheology::DislocationCreepParameters dislocation_creep_parameters,
118-
double log_edot_ii) const;
105+
const double log_edot_ii) const;
119106

120107
/**
121108
* Compute the viscosity based on the composite viscous creep law,

‎source/material_model/rheology/diffusion_dislocation.cc

+36-52
Original file line numberDiff line numberDiff line change
@@ -80,20 +80,12 @@ namespace aspect
8080
(constants::gas_constant*temperature));
8181

8282
// Because the ratios of the diffusion and dislocation strain rates are not known, stress is also unknown
83-
// We use KINSOL to find the second invariant of the stress tensor.
84-
// Start with the assumption that all strain is accommodated by diffusion creep:
85-
// If the diffusion creep prefactor is very small, that means that the diffusion viscosity is very large.
86-
// In this case, use the maximum viscosity instead to compute the starting guess.
87-
double stress_ii = (prefactor_stress_diffusion > (0.5 / maximum_viscosity)
88-
?
89-
edot_ii/prefactor_stress_diffusion
90-
:
91-
0.5 / maximum_viscosity);
92-
Vector<double> log_stress_ii(1);
93-
log_stress_ii[0] = std::log(stress_ii);
94-
95-
typename SUNDIALS::KINSOL<Vector<double>>::AdditionalData additional_data;
83+
// We use KINSOL to find the second invariant of the stress tensor that satisfies the total strain rate.
84+
SUNDIALS::KINSOL<Vector<double>>::AdditionalData additional_data;
85+
additional_data.strategy = dealii::SUNDIALS::KINSOL<>::AdditionalData::newton;
9686
additional_data.function_tolerance = log_strain_rate_residual_threshold;
87+
additional_data.maximum_non_linear_iterations = stress_max_iteration_number;
88+
additional_data.maximum_setup_calls = 1;
9789

9890
SUNDIALS::KINSOL<Vector<double>> nonlinear_solver(additional_data);
9991

@@ -103,74 +95,66 @@ namespace aspect
10395
};
10496

10597
nonlinear_solver.residual =
106-
[&](const Vector<double> &evaluation_point,
98+
[&](const Vector<double> &current_log_stress_ii,
10799
Vector<double> &residual)
108100
{
109-
compute_log_strain_rate_residual(evaluation_point, residual, pressure, temperature, diffusion_creep_parameters, dislocation_creep_parameters, log_edot_ii);
101+
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);
110102
};
111103

112104
nonlinear_solver.setup_jacobian =
113-
[&](const Vector<double> &current_u,
105+
[&](const Vector<double> & /*current_log_stress_ii*/,
114106
const Vector<double> & /*current_f*/)
115107
{
116-
compute_log_strain_rate_deriv(current_u, pressure, temperature, diffusion_creep_parameters, dislocation_creep_parameters, log_strain_rate_deriv);
108+
// Do nothing here, because we calculate the Jacobian in the residual function
117109
};
118110

119-
nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
111+
nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &residual,
120112
Vector<double> &solution,
121113
const double /*tolerance*/)
122114
{
123-
solution(0) -= rhs(0)/log_strain_rate_deriv;
115+
solution(0) = residual(0) / log_strain_rate_deriv;
124116
};
125117

126-
nonlinear_solver.solve(log_stress_ii);
118+
// Start with the assumption that all strain is accommodated by diffusion creep:
119+
// If the diffusion creep prefactor is very small, that means that the diffusion viscosity is very large.
120+
// In this case, use the maximum viscosity instead to compute the starting guess.
121+
double stress_ii = (prefactor_stress_diffusion > (0.5 / maximum_viscosity)
122+
?
123+
edot_ii/prefactor_stress_diffusion
124+
:
125+
0.5 / maximum_viscosity);
127126

127+
// KINSOL works on vectors and so the scalar (log) stress is inserted into
128+
// a vector of length 1
129+
Vector<double> log_stress_ii(1);
130+
log_stress_ii[0] = std::log(stress_ii);
131+
132+
nonlinear_solver.solve(log_stress_ii);
128133
stress_ii = std::exp(log_stress_ii[0]);
129134
composition_viscosities[j] = std::min(std::max(stress_ii/edot_ii/2, minimum_viscosity), maximum_viscosity);
130135
}
131136
return composition_viscosities;
132137
}
133138

134-
135139
template <int dim>
136-
void
137-
DiffusionDislocation<dim>::compute_log_strain_rate_deriv(
138-
const Vector<double> &evaluation_point,
139-
double pressure,
140-
double temperature,
141-
const Rheology::DiffusionCreepParameters diffusion_creep_parameters,
142-
const Rheology::DislocationCreepParameters dislocation_creep_parameters,
143-
double &log_strain_rate_deriv) const
144-
{
145-
const std::pair<double, double> log_diff_edot_and_deriv = diffusion_creep.compute_log_strain_rate_and_derivative(evaluation_point[0], pressure, temperature, diffusion_creep_parameters);
146-
const std::pair<double, double> log_disl_edot_and_deriv = dislocation_creep.compute_log_strain_rate_and_derivative(evaluation_point[0], pressure, temperature, dislocation_creep_parameters);
147-
148-
const double strain_rate_diffusion = std::exp(log_diff_edot_and_deriv.first);
149-
const double strain_rate_dislocation = std::exp(log_disl_edot_and_deriv.first);
150-
log_strain_rate_deriv = (strain_rate_diffusion * log_diff_edot_and_deriv.second + strain_rate_dislocation * log_disl_edot_and_deriv.second)/
151-
(strain_rate_diffusion + strain_rate_dislocation);
152-
}
153-
154-
155-
156-
template <int dim>
157-
void
158-
DiffusionDislocation<dim>::compute_log_strain_rate_residual(
159-
const Vector<double> &evaluation_point,
160-
Vector<double> &residual,
161-
double pressure,
162-
double temperature,
140+
std::pair<double, double>
141+
DiffusionDislocation<dim>::compute_log_strain_rate_residual_and_derivative(
142+
const double current_log_stress_ii,
143+
const double pressure,
144+
const double temperature,
163145
const Rheology::DiffusionCreepParameters diffusion_creep_parameters,
164146
const Rheology::DislocationCreepParameters dislocation_creep_parameters,
165-
double log_edot_ii) const
147+
const double log_edot_ii) const
166148
{
167-
const std::pair<double, double> log_diff_edot_and_deriv = diffusion_creep.compute_log_strain_rate_and_derivative(evaluation_point[0], pressure, temperature, diffusion_creep_parameters);
168-
const std::pair<double, double> log_disl_edot_and_deriv = dislocation_creep.compute_log_strain_rate_and_derivative(evaluation_point[0], pressure, temperature, dislocation_creep_parameters);
149+
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);
150+
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);
169151

170152
const double strain_rate_diffusion = std::exp(log_diff_edot_and_deriv.first);
171153
const double strain_rate_dislocation = std::exp(log_disl_edot_and_deriv.first);
154+
double log_strain_rate_deriv = (strain_rate_diffusion * log_diff_edot_and_deriv.second + strain_rate_dislocation * log_disl_edot_and_deriv.second)/
155+
(strain_rate_diffusion + strain_rate_dislocation);
172156
const double log_strain_rate_iterate = std::log(strain_rate_diffusion + strain_rate_dislocation);
173-
residual(0) = log_edot_ii - log_strain_rate_iterate;
157+
return std::make_pair(log_strain_rate_iterate - log_edot_ii, log_strain_rate_deriv);
174158
}
175159

176160
template <int dim>

0 commit comments

Comments
 (0)