From: Wolfgang Bangerth Date: Wed, 5 Jun 2024 20:56:07 +0000 (-0600) Subject: Add another IDA test for DAEs. X-Git-Tag: v9.6.0-rc1~193^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6477b1f0d385371cfdb8d81e6d92afabfb975847;p=dealii.git Add another IDA test for DAEs. --- diff --git a/tests/sundials/ida_07.cc b/tests/sundials/ida_07.cc new file mode 100644 index 0000000000..654a67d561 --- /dev/null +++ b/tests/sundials/ida_07.cc @@ -0,0 +1,131 @@ +// ------------------------------------------------------------------------ +// +// 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^{1/p} + * 0 = x^p-y + * with initial conditions + * x(0) = 1 + * y(0) = 1 + * This is the same problem as in the _06 test except that we define + * the otherwise entirely unnecessary variable y(t) as + * y(t) = x(t)^p + * with p>=1 instead of + * y(t) = x(t). + * + * That is, with Y=[x, y]: + * F(Y', Y, t) = [x' -a y^{1/p} ; -x^p + y] + * Y(0) = [1 1] + * + * The exact solution is still + * x(t) = exp(a t) + * but now + * y(t) = x(t)^p = [exp(a t)]^p + * + * The Jacobian to assemble is the following: + * + * J = dF/dY + alpha dF/dY' + * = [0 -ay^{1/p-1}/p ; -px^{p-1} 1] + alpha [1 0 ; 0 0] + * = [alpha -ay^{1/p-1}/p ; -px^{p-1} 1] + */ + +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); + + const double a = 1.0; + const double p = 1.5; + deallog << "Exponential growth factor = " << a << std::endl; + + using VectorType = Vector; + + VectorType y(2); + VectorType y_dot(2); + FullMatrix J(2, 2); + FullMatrix A(2, 2); + FullMatrix Jinv(2, 2); + + SUNDIALS::IDA> time_stepper(data); + + 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) = [x' -a y^{1/p} ; -x^p + y] + res = 0; + res[0] = y_dot[0] - a * std::pow(y[1], 1. / p); + res[1] = -std::pow(y[0], p) + y[1]; + }; + + time_stepper.setup_jacobian = [&](const double, + const VectorType &y, + const VectorType &, + const double alpha) { + // J = [alpha -ay^{1/p-1}/p ; -px^{p-1} 1] + J(0, 0) = alpha; + J(0, 1) = -a * std::pow(y[1], 1. / p - 1) / p; + J(1, 0) = -p * std::pow(y[0], p - 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; + }; + + + y[0] = y[1] = 1; + y_dot[0] = a; + y_dot[1] = + p * + std::pow(a, p - 1); // y'(0) = d/dt [x(0)^p] = p [x'(0)]^[p-1] = p a^{p-1} + time_stepper.solve_dae(y, y_dot); +} diff --git a/tests/sundials/ida_07.output b/tests/sundials/ida_07.output new file mode 100644 index 0000000000..151299a759 --- /dev/null +++ b/tests/sundials/ida_07.output @@ -0,0 +1,54 @@ + +DEAL::Exponential growth factor = 1.000000000 +DEAL::0.000000000 1.000000000 1.000000000 1.000000000 1.500000000 +DEAL::0.2000000000 1.221402758 1.349858808 1.221402759 2.024788213 +DEAL::0.4000000000 1.491824698 1.822118801 1.491824698 2.733178203 +DEAL::0.6000000000 1.822118801 2.459603112 1.822118801 3.689404669 +DEAL::0.8000000000 2.225540929 3.320116924 2.225540929 4.980175384 +DEAL::1.000000000 2.718281829 4.481689072 2.718281830 6.722533613 +DEAL::1.200000000 3.320116924 6.049647468 3.320116925 9.074471208 +DEAL::1.400000000 4.055199969 8.166169918 4.055199970 12.24925488 +DEAL::1.600000000 4.953032427 11.02317639 4.953032429 16.53476459 +DEAL::1.800000000 6.049647468 14.87973174 6.049647469 22.31959760 +DEAL::2.000000000 7.389056104 20.08553694 7.389056107 30.12830544 +DEAL::2.200000000 9.025013506 27.11263895 9.025013510 40.66895846 +DEAL::2.400000000 11.02317639 36.59823449 11.02317640 54.89735177 +DEAL::2.600000000 13.46373805 49.40244917 13.46373805 74.10367383 +DEAL::2.800000000 16.44464679 66.68633113 16.44464679 100.0294967 +DEAL::3.000000000 20.08553694 90.01713143 20.08553694 135.0256971 +DEAL::3.200000000 24.53253022 121.5104177 24.53253022 182.2656265 +DEAL::3.400000000 29.96410008 164.0219075 29.96410008 246.0328613 +DEAL::3.600000000 36.59823448 221.4064165 36.59823448 332.1096248 +DEAL::3.800000000 44.70118454 298.8674014 44.70118454 448.3011021 +DEAL::4.000000000 54.59815008 403.4287941 54.59815009 605.1431911 +DEAL::4.200000000 66.68633110 544.5719109 66.68633111 816.8578664 +DEAL::4.400000000 81.45086874 735.0951903 81.45086874 1102.642785 +DEAL::4.600000000 99.48431574 992.2747170 99.48431574 1488.412076 +DEAL::4.800000000 121.5104176 1339.430766 121.5104176 2009.146150 +DEAL::5.000000000 148.4131592 1808.042417 148.4131592 2712.063626 +DEAL::5.200000000 181.2722421 2440.601981 181.2722421 3660.902976 +DEAL::5.400000000 221.4064164 3294.468081 221.4064165 4941.702123 +DEAL::5.600000000 270.4264077 4447.066755 270.4264077 6670.600125 +DEAL::5.800000000 330.2995603 6002.912228 330.2995603 9004.368340 +DEAL::6.000000000 403.4287940 8103.083942 403.4287942 12154.62593 +DEAL::6.200000000 492.7490417 10938.01923 492.7490417 16407.02885 +DEAL::6.400000000 601.8450386 14764.78159 601.8450387 22147.17239 +DEAL::6.600000000 735.0951902 19930.37048 735.0951902 29895.55571 +DEAL::6.800000000 897.8472928 26903.18613 897.8472929 40354.77918 +DEAL::7.000000000 1096.633160 36315.50275 1096.633160 54473.25410 +DEAL::7.200000000 1339.430766 49020.80124 1339.430766 73531.20183 +DEAL::7.400000000 1635.984432 66171.16030 1635.984432 99256.74048 +DEAL::7.600000000 1998.195898 89321.72355 1998.195898 133982.5853 +DEAL::7.800000000 2440.601981 120571.7152 2440.601981 180857.5729 +DEAL::8.000000000 2980.957991 162754.7918 2980.957992 244132.1877 +DEAL::8.200000000 3640.950313 219695.9892 3640.950313 329543.9838 +DEAL::8.400000000 4447.066754 296558.5660 4447.066755 444837.8491 +DEAL::8.600000000 5431.659600 400312.1922 5431.659600 600468.2884 +DEAL::8.800000000 6634.244017 540364.9385 6634.244015 810547.4075 +DEAL::9.000000000 8103.083940 729416.3715 8103.083940 1094124.553 +DEAL::9.200000000 9897.129074 984609.1135 9897.129075 1476913.671 +DEAL::9.400000000 12088.38075 1329083.284 12088.38075 1993624.927 +DEAL::9.600000000 14764.78159 1794074.777 14764.78159 2691112.166 +DEAL::9.800000000 18033.74496 2421747.639 18033.74496 3632621.458 +DEAL::10.00000000 22026.46583 3269017.381 22026.46583 4903526.073 +DEAL::10.00000000 22026.46583 3269017.381 22026.46583 4903526.073