From d4d3a38fbef404afc7632d016ec28a986c40483e Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sat, 1 Jun 2024 17:07:17 -0600 Subject: [PATCH] Add more SUNDIALS::IDA tests. --- tests/sundials/ida_06.cc | 143 +++++++++++++++++++++++++++++++++++ tests/sundials/ida_06.output | 54 +++++++++++++ tests/sundials/ida_06_in.prm | 19 +++++ 3 files changed, 216 insertions(+) create mode 100644 tests/sundials/ida_06.cc create mode 100644 tests/sundials/ida_06.output create mode 100644 tests/sundials/ida_06_in.prm diff --git a/tests/sundials/ida_06.cc b/tests/sundials/ida_06.cc new file mode 100644 index 0000000000..c7e2ec589e --- /dev/null +++ b/tests/sundials/ida_06.cc @@ -0,0 +1,143 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2017 - 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" + + +/** + * Solve an ODE problem of exponential growth, written as a DAE in + * which one could easily eliminate one variable, using a direct + * solver for the jacobian system. + * + * The equation we want to solve here is + * x' = a y + * 0 = x-y + * with initial conditions + * x(0) = 1 + * y(0) = 1 + * + * That is, with Y=[x, y]: + * F(Y', Y, t) = [1 0 ; 0 0] Y' + [0 -a ; -1 1] Y + * Y(0) = [1 1] + * + * The exact solution is + * + * x(t) = y(t) = exp(a t) + * + * The Jacobian to assemble is the following: + * + * J = dF/dY + alpha dF/dY' = [0 -a ; -1 1] + alpha [1 0 ; 0 0] + * = [alpha -a ; -1 1] + */ +class ExponentialGrowth +{ +public: + ExponentialGrowth( + double a_, + const typename SUNDIALS::IDA>::AdditionalData &data) + : time_stepper(data) + , y(2) + , y_dot(2) + , J(2, 2) + , A(2, 2) + , Jinv(2, 2) + , a(a_) + { + using VectorType = Vector; + + deallog << "Exponential growth factor = " << a << std::endl; + + time_stepper.reinit_vector = [&](VectorType &v) { v.reinit(2); }; + + + time_stepper.residual = [&](const double t, + const VectorType &y, + const VectorType &y_dot, + VectorType &res) { + // F(Y', Y, t) = [1 0 ; 0 0] Y' + [0 -a ; -1 1] Y + res = 0; + res[0] = y_dot[0] - a * y[0]; + res[1] = -y[0] + y[1]; + }; + + time_stepper.setup_jacobian = [&](const double, + const VectorType &, + const VectorType &, + const double alpha) { + // J = [alpha -a ; -1 1] + J(0, 0) = alpha; + J(0, 1) = -a; + J(1, 0) = -1; + J(1, 1) = 1; + + Jinv.invert(J); + }; + + time_stepper.solve_with_jacobian = + [&](const VectorType &src, VectorType &dst, const double) { + Jinv.vmult(dst, src); + }; + + time_stepper.output_step = [&](const double t, + const VectorType &sol, + const VectorType &sol_dot, + const unsigned int step_number) { + deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol_dot[0] << ' ' + << sol_dot[1] << std::endl; + }; + } + + void + run() + { + y[0] = y[1] = 1; + y_dot[0] = y_dot[1] = a; + time_stepper.solve_dae(y, y_dot); + } + SUNDIALS::IDA> time_stepper; + +private: + Vector y; + Vector y_dot; + FullMatrix J; + FullMatrix A; + FullMatrix Jinv; + double a; +}; + + +int +main() +{ + initlog(); + deallog << std::setprecision(10); + + SUNDIALS::IDA>::AdditionalData data; + ParameterHandler prm; + data.add_parameters(prm); + + std::ifstream ifile(SOURCE_DIR "/ida_06_in.prm"); + prm.parse_input(ifile); + + + ExponentialGrowth ode(1.0, data); + ode.run(); +} diff --git a/tests/sundials/ida_06.output b/tests/sundials/ida_06.output new file mode 100644 index 0000000000..1e98254c32 --- /dev/null +++ b/tests/sundials/ida_06.output @@ -0,0 +1,54 @@ + +DEAL::Exponential growth factor = 1.000000000 +DEAL::0.000000000 1.000000000 1.000000000 1.000000000 1.000000000 +DEAL::0.2000000000 1.221402760 1.221402760 1.221402758 1.221402758 +DEAL::0.4000000000 1.491824699 1.491824699 1.491824694 1.491824694 +DEAL::0.6000000000 1.822118802 1.822118802 1.822118802 1.822118802 +DEAL::0.8000000000 2.225540931 2.225540931 2.225540923 2.225540923 +DEAL::1.000000000 2.718281831 2.718281831 2.718281834 2.718281834 +DEAL::1.200000000 3.320116927 3.320116927 3.320116928 3.320116928 +DEAL::1.400000000 4.055199972 4.055199972 4.055199968 4.055199968 +DEAL::1.600000000 4.953032430 4.953032430 4.953032433 4.953032433 +DEAL::1.800000000 6.049647471 6.049647471 6.049647488 6.049647488 +DEAL::2.000000000 7.389056107 7.389056107 7.389056106 7.389056106 +DEAL::2.200000000 9.025013510 9.025013510 9.025013488 9.025013488 +DEAL::2.400000000 11.02317639 11.02317639 11.02317640 11.02317640 +DEAL::2.600000000 13.46373805 13.46373805 13.46373808 13.46373808 +DEAL::2.800000000 16.44464679 16.44464679 16.44464677 16.44464677 +DEAL::3.000000000 20.08553694 20.08553694 20.08553698 20.08553698 +DEAL::3.200000000 24.53253022 24.53253022 24.53253023 24.53253023 +DEAL::3.400000000 29.96410008 29.96410008 29.96410001 29.96410001 +DEAL::3.600000000 36.59823448 36.59823448 36.59823450 36.59823450 +DEAL::3.800000000 44.70118453 44.70118453 44.70118452 44.70118452 +DEAL::4.000000000 54.59815008 54.59815008 54.59815004 54.59815004 +DEAL::4.200000000 66.68633110 66.68633110 66.68633116 66.68633116 +DEAL::4.400000000 81.45086873 81.45086873 81.45086894 81.45086894 +DEAL::4.600000000 99.48431571 99.48431571 99.48431557 99.48431557 +DEAL::4.800000000 121.5104176 121.5104176 121.5104174 121.5104174 +DEAL::5.000000000 148.4131592 148.4131592 148.4131592 148.4131592 +DEAL::5.200000000 181.2722420 181.2722420 181.2722420 181.2722420 +DEAL::5.400000000 221.4064163 221.4064163 221.4064158 221.4064158 +DEAL::5.600000000 270.4264076 270.4264076 270.4264076 270.4264076 +DEAL::5.800000000 330.2995601 330.2995601 330.2995600 330.2995600 +DEAL::6.000000000 403.4287937 403.4287937 403.4287934 403.4287934 +DEAL::6.200000000 492.7490413 492.7490413 492.7490401 492.7490401 +DEAL::6.400000000 601.8450381 601.8450381 601.8450396 601.8450396 +DEAL::6.600000000 735.0951895 735.0951895 735.0951884 735.0951884 +DEAL::6.800000000 897.8472919 897.8472919 897.8472909 897.8472909 +DEAL::7.000000000 1096.633159 1096.633159 1096.633156 1096.633156 +DEAL::7.200000000 1339.430765 1339.430765 1339.430765 1339.430765 +DEAL::7.400000000 1635.984430 1635.984430 1635.984427 1635.984427 +DEAL::7.600000000 1998.195895 1998.195895 1998.195896 1998.195896 +DEAL::7.800000000 2440.601978 2440.601978 2440.601973 2440.601973 +DEAL::8.000000000 2980.957987 2980.957987 2980.957987 2980.957987 +DEAL::8.200000000 3640.950307 3640.950307 3640.950299 3640.950299 +DEAL::8.400000000 4447.066747 4447.066747 4447.066754 4447.066754 +DEAL::8.600000000 5431.659591 5431.659591 5431.659579 5431.659579 +DEAL::8.800000000 6634.244005 6634.244005 6634.244004 6634.244004 +DEAL::9.000000000 8103.083926 8103.083926 8103.083908 8103.083908 +DEAL::9.200000000 9897.129056 9897.129056 9897.129058 9897.129058 +DEAL::9.400000000 12088.38073 12088.38073 12088.38068 12088.38068 +DEAL::9.600000000 14764.78156 14764.78156 14764.78157 14764.78157 +DEAL::9.800000000 18033.74492 18033.74492 18033.74489 18033.74489 +DEAL::10.00000000 22026.46579 22026.46579 22026.46577 22026.46577 +DEAL::10.00000000 22026.46579 22026.46579 22026.46577 22026.46577 diff --git a/tests/sundials/ida_06_in.prm b/tests/sundials/ida_06_in.prm new file mode 100644 index 0000000000..fc46a20f46 --- /dev/null +++ b/tests/sundials/ida_06_in.prm @@ -0,0 +1,19 @@ +set Final time = 10 +set Initial time = 0 +set Time interval between each output = 0.2 +subsection Error control + set Absolute error tolerance = 1e-10 + set Ignore algebraic terms for error computations = true + set Relative error tolerance = 1e-10 +end +subsection Initial condition correction parameters + set Correction type at initial time = none + set Correction type after restart = none + set Maximum number of nonlinear iterations = 10 +end +subsection Running parameters + set Initial step size = 1e-6 + set Maximum number of nonlinear iterations = 10 + set Maximum order of BDF = 5 + set Minimum step size = 1e-7 +end -- 2.39.5