@@ -97,41 +97,34 @@ namespace aspect
97
97
98
98
SUNDIALS::KINSOL<Vector<double >> nonlinear_solver (additional_data);
99
99
100
- // Declare lambda functions
101
- // (i) resize a vector to the correct size
102
100
nonlinear_solver.reinit_vector = [&](Vector<double > &x)
103
101
{
104
102
x.reinit (1 );
105
103
};
106
104
107
- // (ii) compute the residual vector
108
105
nonlinear_solver.residual =
109
106
[&](const Vector<double > &evaluation_point,
110
107
Vector<double > &residual)
111
108
{
112
109
compute_log_strain_rate_residual (evaluation_point, residual, pressure, temperature, diffusion_creep_parameters, dislocation_creep_parameters, log_edot_ii);
113
110
};
114
111
115
- // (iii) compute the Jacobian matrix and factorization
116
112
nonlinear_solver.setup_jacobian =
117
113
[&](const Vector<double > ¤t_u,
118
114
const Vector<double > & /* current_f*/ )
119
115
{
120
116
compute_log_strain_rate_deriv (current_u, pressure, temperature, diffusion_creep_parameters, dislocation_creep_parameters, log_strain_rate_deriv);
121
117
};
122
118
123
- // (iv) solve a linear system with the Jacobian.
124
119
nonlinear_solver.solve_with_jacobian = [&](const Vector<double > &rhs,
125
120
Vector<double > &solution,
126
121
const double /* tolerance*/ )
127
122
{
128
123
solution (0 ) -= rhs (0 )/log_strain_rate_deriv;
129
124
};
130
125
131
- // solve the problem
132
126
nonlinear_solver.solve (log_stress_ii);
133
127
134
- // The effective viscosity, with minimum and maximum bounds
135
128
stress_ii = std::exp (log_stress_ii[0 ]);
136
129
composition_viscosities[j] = std::min (std::max (stress_ii/edot_ii/2 , minimum_viscosity), maximum_viscosity);
137
130
}
0 commit comments