From 3ab637933af09eb59fc548cdcebc7a1fb921043a Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Fri, 31 May 2024 18:06:14 -0600 Subject: [PATCH 1/6] use sundials kinsol for diffusion dislocation creep --- .../rheology/diffusion_dislocation.h | 27 +++ .../rheology/diffusion_dislocation.cc | 177 +++++++++--------- 2 files changed, 115 insertions(+), 89 deletions(-) diff --git a/include/aspect/material_model/rheology/diffusion_dislocation.h b/include/aspect/material_model/rheology/diffusion_dislocation.h index 49315dda652..b2ed7a001a5 100644 --- a/include/aspect/material_model/rheology/diffusion_dislocation.h +++ b/include/aspect/material_model/rheology/diffusion_dislocation.h @@ -90,6 +90,33 @@ namespace aspect const double temperature, const SymmetricTensor<2,dim> &strain_rate) const; + /** + * Compute the derivative of the natural logarithm of the strain rate + * with respect to the natural logarithm of the stress + * at the current evaluation point (log stress) estimated by KINSOL. + */ + void + compute_log_strain_rate_deriv(const Vector &evaluation_point, + double pressure, + double temperature, + const Rheology::DiffusionCreepParameters diffusion_creep_parameters, + const Rheology::DislocationCreepParameters dislocation_creep_parameters, + double &log_strain_rate_deriv) const; + + /** + * Compute the residual of the natural logarithm of the strain rate + * at the current evaluation point (log stress) estimated by KINSOL. + */ + void + compute_log_strain_rate_residual( + const Vector &evaluation_point, + Vector &residual, + double pressure, + double temperature, + const Rheology::DiffusionCreepParameters diffusion_creep_parameters, + const Rheology::DislocationCreepParameters dislocation_creep_parameters, + double log_edot_ii) const; + /** * Compute the viscosity based on the composite viscous creep law, * averaging over all compositional fields according to their diff --git a/source/material_model/rheology/diffusion_dislocation.cc b/source/material_model/rheology/diffusion_dislocation.cc index d86b36f491c..be7cf72837b 100644 --- a/source/material_model/rheology/diffusion_dislocation.cc +++ b/source/material_model/rheology/diffusion_dislocation.cc @@ -23,6 +23,7 @@ #include #include +#include #include #include @@ -54,10 +55,10 @@ namespace aspect const double edot_ii = std::max(std::sqrt(std::max(-second_invariant(deviator(strain_rate)), 0.)), min_strain_rate); const double log_edot_ii = std::log(edot_ii); - + double log_strain_rate_deriv; // 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 composition_viscosities(n_chemical_composition_fields); for (unsigned int j=0; j < n_chemical_composition_fields; ++j) { @@ -79,7 +80,7 @@ namespace aspect (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. + // We use KINSOL 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. @@ -88,99 +89,97 @@ namespace aspect 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) - { - const std::pair log_diff_edot_and_deriv = diffusion_creep.compute_log_strain_rate_and_derivative(log_stress_ii, pressure, temperature, diffusion_creep_parameters); - const std::pair 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::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::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; - } - } + Vector log_stress_ii(1); + log_stress_ii[0] = std::log(stress_ii); + + typename SUNDIALS::KINSOL>::AdditionalData additional_data; + additional_data.function_tolerance = log_strain_rate_residual_threshold; + + SUNDIALS::KINSOL> nonlinear_solver(additional_data); + + // Declare lambda functions + // (i) resize a vector to the correct size + nonlinear_solver.reinit_vector = [&](Vector &x) + { + x.reinit(1); + }; + + // (ii) compute the residual vector + nonlinear_solver.residual = + [&](const Vector &evaluation_point, + Vector &residual) + { + compute_log_strain_rate_residual(evaluation_point, residual, pressure, temperature, diffusion_creep_parameters, dislocation_creep_parameters, log_edot_ii); + }; + + // (iii) compute the Jacobian matrix and factorization + nonlinear_solver.setup_jacobian = + [&](const Vector ¤t_u, + const Vector & /*current_f*/) + { + compute_log_strain_rate_deriv(current_u, pressure, temperature, diffusion_creep_parameters, dislocation_creep_parameters, log_strain_rate_deriv); + }; + + // (iv) solve a linear system with the Jacobian. + nonlinear_solver.solve_with_jacobian = [&](const Vector &rhs, + Vector &solution, + const double /*tolerance*/) + { + solution(0) -= rhs(0)/log_strain_rate_deriv; + }; + + // solve the problem + nonlinear_solver.solve(log_stress_ii); // The effective viscosity, with minimum and maximum bounds + 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 + void + DiffusionDislocation::compute_log_strain_rate_deriv( + const Vector &evaluation_point, + double pressure, + double temperature, + const Rheology::DiffusionCreepParameters diffusion_creep_parameters, + const Rheology::DislocationCreepParameters dislocation_creep_parameters, + double &log_strain_rate_deriv) const + { + const std::pair log_diff_edot_and_deriv = diffusion_creep.compute_log_strain_rate_and_derivative(evaluation_point[0], pressure, temperature, diffusion_creep_parameters); + const std::pair log_disl_edot_and_deriv = dislocation_creep.compute_log_strain_rate_and_derivative(evaluation_point[0], 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_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); + } + + + + template + void + DiffusionDislocation::compute_log_strain_rate_residual( + const Vector &evaluation_point, + Vector &residual, + double pressure, + double temperature, + const Rheology::DiffusionCreepParameters diffusion_creep_parameters, + const Rheology::DislocationCreepParameters dislocation_creep_parameters, + double log_edot_ii) const + { + const std::pair log_diff_edot_and_deriv = diffusion_creep.compute_log_strain_rate_and_derivative(evaluation_point[0], pressure, temperature, diffusion_creep_parameters); + const std::pair log_disl_edot_and_deriv = dislocation_creep.compute_log_strain_rate_and_derivative(evaluation_point[0], 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); + const double log_strain_rate_iterate = std::log(strain_rate_diffusion + strain_rate_dislocation); + residual(0) = log_edot_ii - log_strain_rate_iterate; + } + template double DiffusionDislocation::compute_viscosity (const double pressure, From fa84fc70b5118fae5389240edb0ce85d38518089 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Fri, 31 May 2024 18:11:24 -0600 Subject: [PATCH 2/6] added changelog --- doc/modules/changes/20240531_myhill | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 doc/modules/changes/20240531_myhill diff --git a/doc/modules/changes/20240531_myhill b/doc/modules/changes/20240531_myhill new file mode 100644 index 00000000000..1e3ad61dfb9 --- /dev/null +++ b/doc/modules/changes/20240531_myhill @@ -0,0 +1,4 @@ +Changed: The DiffusionDislocation rheology now uses SUNDIALS KINSOL +routines to calculate the equilibrium partitioning of strain rates. +
+(Bob Myhill, 2024/05/31) From fb9cb6dc36feef9e9037023c6c3f8125e5db5b99 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Fri, 31 May 2024 20:00:42 -0600 Subject: [PATCH 3/6] updated tests --- tests/diffusion_dislocation/screen-output | 12 ++++++------ tests/diffusion_dislocation/statistics | 2 +- .../screen-output | 12 ++++++------ .../statistics | 2 +- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/tests/diffusion_dislocation/screen-output b/tests/diffusion_dislocation/screen-output index e5b62ca7ac3..f0279f5d75d 100644 --- a/tests/diffusion_dislocation/screen-output +++ b/tests/diffusion_dislocation/screen-output @@ -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 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: diff --git a/tests/diffusion_dislocation/statistics b/tests/diffusion_dislocation/statistics index 97067103310..45ee5689647 100644 --- a/tests/diffusion_dislocation/statistics +++ b/tests/diffusion_dislocation/statistics @@ -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 diff --git a/tests/diffusion_dislocation_reference_density_profile/screen-output b/tests/diffusion_dislocation_reference_density_profile/screen-output index 4763cddc06d..076ffef2eb8 100644 --- a/tests/diffusion_dislocation_reference_density_profile/screen-output +++ b/tests/diffusion_dislocation_reference_density_profile/screen-output @@ -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: diff --git a/tests/diffusion_dislocation_reference_density_profile/statistics b/tests/diffusion_dislocation_reference_density_profile/statistics index 97067103310..45ee5689647 100644 --- a/tests/diffusion_dislocation_reference_density_profile/statistics +++ b/tests/diffusion_dislocation_reference_density_profile/statistics @@ -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 From 26cfcc0d933030653209bf2451138e561df45991 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Fri, 31 May 2024 21:34:36 -0600 Subject: [PATCH 4/6] removed superfluous docstrings --- source/material_model/rheology/diffusion_dislocation.cc | 7 ------- 1 file changed, 7 deletions(-) diff --git a/source/material_model/rheology/diffusion_dislocation.cc b/source/material_model/rheology/diffusion_dislocation.cc index be7cf72837b..fd0bd8c86d3 100644 --- a/source/material_model/rheology/diffusion_dislocation.cc +++ b/source/material_model/rheology/diffusion_dislocation.cc @@ -97,14 +97,11 @@ namespace aspect SUNDIALS::KINSOL> nonlinear_solver(additional_data); - // Declare lambda functions - // (i) resize a vector to the correct size nonlinear_solver.reinit_vector = [&](Vector &x) { x.reinit(1); }; - // (ii) compute the residual vector nonlinear_solver.residual = [&](const Vector &evaluation_point, Vector &residual) @@ -112,7 +109,6 @@ namespace aspect compute_log_strain_rate_residual(evaluation_point, residual, pressure, temperature, diffusion_creep_parameters, dislocation_creep_parameters, log_edot_ii); }; - // (iii) compute the Jacobian matrix and factorization nonlinear_solver.setup_jacobian = [&](const Vector ¤t_u, const Vector & /*current_f*/) @@ -120,7 +116,6 @@ namespace aspect compute_log_strain_rate_deriv(current_u, pressure, temperature, diffusion_creep_parameters, dislocation_creep_parameters, log_strain_rate_deriv); }; - // (iv) solve a linear system with the Jacobian. nonlinear_solver.solve_with_jacobian = [&](const Vector &rhs, Vector &solution, const double /*tolerance*/) @@ -128,10 +123,8 @@ namespace aspect solution(0) -= rhs(0)/log_strain_rate_deriv; }; - // solve the problem nonlinear_solver.solve(log_stress_ii); - // The effective viscosity, with minimum and maximum bounds stress_ii = std::exp(log_stress_ii[0]); composition_viscosities[j] = std::min(std::max(stress_ii/edot_ii/2, minimum_viscosity), maximum_viscosity); } From 39f923ea0819e41cdc92389903547036b94a6ddc Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Sat, 1 Jun 2024 12:14:32 -0600 Subject: [PATCH 5/6] merge residual and jacobian calculation --- .../rheology/diffusion_dislocation.h | 29 ++---- .../rheology/diffusion_dislocation.cc | 88 ++++++++----------- 2 files changed, 44 insertions(+), 73 deletions(-) diff --git a/include/aspect/material_model/rheology/diffusion_dislocation.h b/include/aspect/material_model/rheology/diffusion_dislocation.h index b2ed7a001a5..8f464eabda0 100644 --- a/include/aspect/material_model/rheology/diffusion_dislocation.h +++ b/include/aspect/material_model/rheology/diffusion_dislocation.h @@ -90,32 +90,19 @@ namespace aspect const double temperature, const SymmetricTensor<2,dim> &strain_rate) const; - /** - * Compute the derivative of the natural logarithm of the strain rate - * with respect to the natural logarithm of the stress - * at the current evaluation point (log stress) estimated by KINSOL. - */ - void - compute_log_strain_rate_deriv(const Vector &evaluation_point, - double pressure, - double temperature, - const Rheology::DiffusionCreepParameters diffusion_creep_parameters, - const Rheology::DislocationCreepParameters dislocation_creep_parameters, - double &log_strain_rate_deriv) const; - /** * Compute the residual of the natural logarithm of the strain rate - * at the current evaluation point (log stress) estimated by KINSOL. + * and the derivative with respect to the logarithm of the stress + * at the current value of the logarithm of the stress. */ - void - compute_log_strain_rate_residual( - const Vector &evaluation_point, - Vector &residual, - double pressure, - double temperature, + std::pair + 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, - double log_edot_ii) const; + const double log_edot_ii) const; /** * Compute the viscosity based on the composite viscous creep law, diff --git a/source/material_model/rheology/diffusion_dislocation.cc b/source/material_model/rheology/diffusion_dislocation.cc index fd0bd8c86d3..61d6cd03325 100644 --- a/source/material_model/rheology/diffusion_dislocation.cc +++ b/source/material_model/rheology/diffusion_dislocation.cc @@ -80,20 +80,12 @@ namespace aspect (constants::gas_constant*temperature)); // 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. - // 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); - Vector log_stress_ii(1); - log_stress_ii[0] = std::log(stress_ii); - - typename SUNDIALS::KINSOL>::AdditionalData additional_data; + // We use KINSOL to find the second invariant of the stress tensor that satisfies the total strain rate. + SUNDIALS::KINSOL>::AdditionalData additional_data; + additional_data.strategy = dealii::SUNDIALS::KINSOL<>::AdditionalData::newton; additional_data.function_tolerance = log_strain_rate_residual_threshold; + additional_data.maximum_non_linear_iterations = stress_max_iteration_number; + additional_data.maximum_setup_calls = 1; SUNDIALS::KINSOL> nonlinear_solver(additional_data); @@ -103,74 +95,66 @@ namespace aspect }; nonlinear_solver.residual = - [&](const Vector &evaluation_point, + [&](const Vector ¤t_log_stress_ii, Vector &residual) { - compute_log_strain_rate_residual(evaluation_point, residual, pressure, temperature, diffusion_creep_parameters, dislocation_creep_parameters, log_edot_ii); + 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 ¤t_u, + [&](const Vector & /*current_log_stress_ii*/, const Vector & /*current_f*/) { - compute_log_strain_rate_deriv(current_u, pressure, temperature, diffusion_creep_parameters, dislocation_creep_parameters, log_strain_rate_deriv); + // Do nothing here, because we calculate the Jacobian in the residual function }; - nonlinear_solver.solve_with_jacobian = [&](const Vector &rhs, + nonlinear_solver.solve_with_jacobian = [&](const Vector &residual, Vector &solution, const double /*tolerance*/) { - solution(0) -= rhs(0)/log_strain_rate_deriv; + solution(0) = residual(0) / log_strain_rate_deriv; }; - nonlinear_solver.solve(log_stress_ii); + // 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); + // KINSOL works on vectors and so the scalar (log) stress is inserted into + // a vector of length 1 + Vector log_stress_ii(1); + log_stress_ii[0] = std::log(stress_ii); + + 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 - void - DiffusionDislocation::compute_log_strain_rate_deriv( - const Vector &evaluation_point, - double pressure, - double temperature, - const Rheology::DiffusionCreepParameters diffusion_creep_parameters, - const Rheology::DislocationCreepParameters dislocation_creep_parameters, - double &log_strain_rate_deriv) const - { - const std::pair log_diff_edot_and_deriv = diffusion_creep.compute_log_strain_rate_and_derivative(evaluation_point[0], pressure, temperature, diffusion_creep_parameters); - const std::pair log_disl_edot_and_deriv = dislocation_creep.compute_log_strain_rate_and_derivative(evaluation_point[0], 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_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); - } - - - - template - void - DiffusionDislocation::compute_log_strain_rate_residual( - const Vector &evaluation_point, - Vector &residual, - double pressure, - double temperature, + std::pair + DiffusionDislocation::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, - double log_edot_ii) const + const double log_edot_ii) const { - const std::pair log_diff_edot_and_deriv = diffusion_creep.compute_log_strain_rate_and_derivative(evaluation_point[0], pressure, temperature, diffusion_creep_parameters); - const std::pair log_disl_edot_and_deriv = dislocation_creep.compute_log_strain_rate_and_derivative(evaluation_point[0], pressure, temperature, dislocation_creep_parameters); + const std::pair 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 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); - residual(0) = log_edot_ii - log_strain_rate_iterate; + return std::make_pair(log_strain_rate_iterate - log_edot_ii, log_strain_rate_deriv); } template From 3b01b3afcec41de7c2d289acf320f6ed56da6fe4 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Thu, 6 Jun 2024 11:56:31 -0600 Subject: [PATCH 6/6] addressed comments --- .../rheology/diffusion_dislocation.h | 7 +- .../rheology/diffusion_dislocation.cc | 64 ++++++++++--------- 2 files changed, 40 insertions(+), 31 deletions(-) diff --git a/include/aspect/material_model/rheology/diffusion_dislocation.h b/include/aspect/material_model/rheology/diffusion_dislocation.h index 8f464eabda0..a12d5bbc251 100644 --- a/include/aspect/material_model/rheology/diffusion_dislocation.h +++ b/include/aspect/material_model/rheology/diffusion_dislocation.h @@ -28,6 +28,8 @@ #include #include +#include + namespace aspect { namespace MaterialModel @@ -100,8 +102,8 @@ namespace aspect 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 Rheology::DiffusionCreepParameters &diffusion_creep_parameters, + const Rheology::DislocationCreepParameters &dislocation_creep_parameters, const double log_edot_ii) const; /** @@ -144,6 +146,7 @@ namespace aspect unsigned int n_chemical_composition_fields; MaterialUtilities::CompositionalAveragingOperation viscosity_averaging; + SUNDIALS::KINSOL>::AdditionalData kinsol_settings; }; diff --git a/source/material_model/rheology/diffusion_dislocation.cc b/source/material_model/rheology/diffusion_dislocation.cc index 61d6cd03325..0a6a2a46560 100644 --- a/source/material_model/rheology/diffusion_dislocation.cc +++ b/source/material_model/rheology/diffusion_dislocation.cc @@ -23,7 +23,6 @@ #include #include -#include #include #include @@ -55,14 +54,16 @@ namespace aspect const double edot_ii = std::max(std::sqrt(std::max(-second_invariant(deviator(strain_rate)), 0.)), min_strain_rate); const double log_edot_ii = std::log(edot_ii); - double log_strain_rate_deriv; // Find effective viscosities for each of the individual phases // Viscosities should have the same number of entries as compositional fields std::vector 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 @@ -73,22 +74,9 @@ 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 - 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 KINSOL to find the second invariant of the stress tensor that satisfies the total strain rate. - SUNDIALS::KINSOL>::AdditionalData additional_data; - additional_data.strategy = dealii::SUNDIALS::KINSOL<>::AdditionalData::newton; - additional_data.function_tolerance = log_strain_rate_residual_threshold; - additional_data.maximum_non_linear_iterations = stress_max_iteration_number; - additional_data.maximum_setup_calls = 1; - - SUNDIALS::KINSOL> nonlinear_solver(additional_data); + SUNDIALS::KINSOL> nonlinear_solver(kinsol_settings); + double log_strain_rate_deriv; nonlinear_solver.reinit_vector = [&](Vector &x) { x.reinit(1); @@ -98,7 +86,9 @@ namespace aspect [&](const Vector ¤t_log_stress_ii, Vector &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); + 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 = @@ -115,14 +105,25 @@ namespace aspect solution(0) = residual(0) / log_strain_rate_deriv; }; - // 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); + + // 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)); + + 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 @@ -142,8 +143,8 @@ namespace aspect 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 Rheology::DiffusionCreepParameters &diffusion_creep_parameters, + const Rheology::DislocationCreepParameters &dislocation_creep_parameters, const double log_edot_ii) const { const std::pair log_diff_edot_and_deriv = diffusion_creep.compute_log_strain_rate_and_derivative(current_log_stress_ii, pressure, temperature, diffusion_creep_parameters); @@ -339,6 +340,11 @@ namespace aspect dislocation_creep.initialize_simulator (this->get_simulator()); dislocation_creep.parse_parameters(prm, std::make_unique>(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; } } }