From: Wolfgang Bangerth Date: Thu, 6 Jun 2024 19:57:23 +0000 (-0600) Subject: Add another KINSOL test. X-Git-Tag: v9.6.0-rc1~198^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c3223d7dd131f1d0edf219d61b078ebb95c47448;p=dealii.git Add another KINSOL test. --- diff --git a/tests/sundials/kinsol_09.cc b/tests/sundials/kinsol_09.cc new file mode 100644 index 0000000000..21b508f01c --- /dev/null +++ b/tests/sundials/kinsol_09.cc @@ -0,0 +1,209 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2021 - 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + +#include + +#include +#include + +#include + +#include "../tests.h" + +// A test for a 1d nonlinear problem taken from work of Bob Myhill in +// ASPECT. It shows that we are calling the residual() function more +// than perhaps necessary. + + +const double pressure = 59448242.437; +const double temperature = 327.2685405; +const double log_edot_ii = -40.053535387; +const double grain_size = 0.001; + +const double gas_constant = 8.31446; + +using namespace dealii; + +namespace DiffusionCreepParameters +{ + const double prefactor = 1.5e-16; + const double grain_size_exponent = 4; + const double stress_exponent = 1; + const double activation_energy = 375000; + const double activation_volume = 6e-06; +}; // namespace DiffusionCreepParameters + +namespace DislocationCreepParameters +{ + const double prefactor = 2e-15; + const double stress_exponent = 3.5; + const double activation_energy = 480000; + const double activation_volume = 8e-06; +}; // namespace DislocationCreepParameters + + +std::pair +compute_diffusion_log_strain_rate_and_derivative(const double log_stress, + const double pressure, + const double temperature) +{ + const double log_strain_rate_diffusion = + std::log(DiffusionCreepParameters::prefactor) + log_stress - + DiffusionCreepParameters::grain_size_exponent * std::log(grain_size) - + (DiffusionCreepParameters::activation_energy + + pressure * DiffusionCreepParameters::activation_volume) / + (gas_constant * temperature); + + const double dlog_strain_rate_dlog_stress_diffusion = 1.0; + + return std::make_pair(log_strain_rate_diffusion, + dlog_strain_rate_dlog_stress_diffusion); +} + + + +std::pair +compute_dislocation_log_strain_rate_and_derivative(const double log_stress, + const double pressure, + const double temperature) +{ + const double log_strain_rate_dislocation = + std::log(DislocationCreepParameters::prefactor) + + DislocationCreepParameters::stress_exponent * log_stress - + (DislocationCreepParameters::activation_energy + + pressure * DislocationCreepParameters::activation_volume) / + (gas_constant * temperature); + + const double dlog_strain_rate_dlog_stress_dislocation = + DislocationCreepParameters::stress_exponent; + + return std::make_pair(log_strain_rate_dislocation, + dlog_strain_rate_dlog_stress_dislocation); +} + + + +std::pair +compute_log_strain_rate_residual_and_derivative( + const double current_log_stress_ii, + const double pressure, + const double temperature, + const double log_edot_ii) +{ + const std::pair log_diff_edot_and_deriv = + compute_diffusion_log_strain_rate_and_derivative(current_log_stress_ii, + pressure, + temperature); + const std::pair log_disl_edot_and_deriv = + compute_dislocation_log_strain_rate_and_derivative(current_log_stress_ii, + pressure, + temperature); + + 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); +} + + +int +main() +{ + initlog(); + + const double maximum_viscosity = 1.e30; + double log_strain_rate_deriv; + + // For diffusion creep, viscosity is grain size dependent + const double prefactor_stress_diffusion = + DiffusionCreepParameters::prefactor * + std::pow(grain_size, -DiffusionCreepParameters::grain_size_exponent) * + std::exp( + -(std::max(DiffusionCreepParameters::activation_energy + + pressure * DiffusionCreepParameters::activation_volume, + 0.0)) / + (gas_constant * temperature)); + + SUNDIALS::KINSOL>::AdditionalData additional_data; + additional_data.strategy = dealii::SUNDIALS::KINSOL<>::AdditionalData::newton; + additional_data.function_tolerance = 1e-10; + additional_data.maximum_non_linear_iterations = 200; + additional_data.maximum_setup_calls = 10; + + int n_residual_evaluations = 0; + SUNDIALS::KINSOL> nonlinear_solver(additional_data); + + + nonlinear_solver.reinit_vector = [&](Vector &x) { x.reinit(1); }; + + + nonlinear_solver.residual = [&](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, + log_edot_ii); + + deallog << std::setprecision(11) + << " Computing residual at x=" << current_log_stress_ii[0] + << ", f(x)=" << residual(0) << std::endl; + n_residual_evaluations += 1; + }; + + + nonlinear_solver.setup_jacobian = + [&](const Vector ¤t_log_stress_ii, + const Vector & /*current_f*/) { + // Do nothing here, because we calculate the Jacobian in the residual + // function + deallog << " Recomputing J at x=" << current_log_stress_ii[0] + << std::endl; + }; + + + nonlinear_solver.solve_with_jacobian = [&](const Vector &residual, + Vector &solution, + const double /*tolerance*/) { + deallog << " Solving for dx with residual=" << residual[0] << std::endl; + + 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. + const double stress_ii = + (prefactor_stress_diffusion > (0.5 / maximum_viscosity) ? + std::exp(log_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); + deallog << n_residual_evaluations << " residual evaluations" << std::endl; +} diff --git a/tests/sundials/kinsol_09.output b/tests/sundials/kinsol_09.output new file mode 100644 index 0000000000..9d8955dfdf --- /dev/null +++ b/tests/sundials/kinsol_09.output @@ -0,0 +1,11 @@ + +DEAL:: Computing residual at x=-69.770699970, f(x)=-176.46700446 +DEAL:: Recomputing J at x=-69.770699970 +DEAL:: Solving for dx with residual=176.46700446 +DEAL:: Computing residual at x=-69.770701010, f(x)=-176.46700550 +DEAL:: Computing residual at x=106.69630449, f(x)=203.06842599 +DEAL:: Recomputing J at x=106.69630449 +DEAL:: Solving for dx with residual=-203.06842599 +DEAL:: Computing residual at x=106.69630608, f(x)=203.06843156 +DEAL:: Computing residual at x=48.676754204, f(x)=1.4210854715e-14 +DEAL::5 residual evaluations