From 5454fe884e93e50f176edc2420f5228ce38e081f Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 5 Jun 2024 17:01:47 -0600 Subject: [PATCH] Add another IDA test, this time for lacking initial velocities and one initial condition. --- tests/sundials/ida_09.cc | 110 +++++++++++++++++++++++++++++++++++ tests/sundials/ida_09.output | 54 +++++++++++++++++ 2 files changed, 164 insertions(+) create mode 100644 tests/sundials/ida_09.cc create mode 100644 tests/sundials/ida_09.output diff --git a/tests/sundials/ida_09.cc b/tests/sundials/ida_09.cc new file mode 100644 index 0000000000..4bc43ca5a8 --- /dev/null +++ b/tests/sundials/ida_09.cc @@ -0,0 +1,110 @@ +// ------------------------------------------------------------------------ +// +// 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" + + +// Like the _07 test, but we only provide correct initial conditions +// for Y[0]', not for Y[1]' -- this is not uncommon in cases where we +// have a first order DAE where at best we know initial conditions for +// all variables (or could compute them), but definitely not initial +// velocities Y'. + +int +main() +{ + initlog(); + deallog << std::setprecision(10); + + SUNDIALS::IDA>::AdditionalData data; + ParameterHandler prm; + data.add_parameters(prm); + + std::ifstream ifile(SOURCE_DIR "/ida_08_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; + }; + + time_stepper.differential_components = []() { + IndexSet x(2); + x.add_index(0); + return x; + }; + + // Provide correct initial conditions x(0), but incorrect initial + // conditions for y(0) and derivatives x'(0), y'(0): + y[0] = 1; // correct + y[1] = 42; // wrong + y_dot[0] = 0; // wrong + y_dot[1] = 0; // wrong + time_stepper.solve_dae(y, y_dot); +} diff --git a/tests/sundials/ida_09.output b/tests/sundials/ida_09.output new file mode 100644 index 0000000000..01fa9e1943 --- /dev/null +++ b/tests/sundials/ida_09.output @@ -0,0 +1,54 @@ + +DEAL::Exponential growth factor = 1.000000000 +DEAL::0.000000000 1.000000000 1.000000000 0.9999999977 0.000000000 +DEAL::0.2000000000 1.221402759 1.349858809 1.221402759 2.024788213 +DEAL::0.4000000000 1.491824699 1.822118803 1.491824700 2.733178187 +DEAL::0.6000000000 1.822118803 2.459603116 1.822118807 3.689404694 +DEAL::0.8000000000 2.225540932 3.320116931 2.225540933 4.980175346 +DEAL::1.000000000 2.718281834 4.481689084 2.718281841 6.722533663 +DEAL::1.200000000 3.320116931 6.049647487 3.320116939 9.074471279 +DEAL::1.400000000 4.055199979 8.166169949 4.055199988 12.24925498 +DEAL::1.600000000 4.953032441 11.02317644 4.953032454 16.53476475 +DEAL::1.800000000 6.049647487 14.87973181 6.049647500 22.31959778 +DEAL::2.000000000 7.389056130 20.08553705 7.389056149 30.12830576 +DEAL::2.200000000 9.025013541 27.11263911 9.025013558 40.66895871 +DEAL::2.400000000 11.02317644 36.59823472 11.02317646 54.89735242 +DEAL::2.600000000 13.46373811 49.40244951 13.46373813 74.10367417 +DEAL::2.800000000 16.44464687 66.68633162 16.44464691 100.0294981 +DEAL::3.000000000 20.08553705 90.01713214 20.08553707 135.0256977 +DEAL::3.200000000 24.53253036 121.5104187 24.53253042 182.2656292 +DEAL::3.400000000 29.96410025 164.0219090 29.96410027 246.0328616 +DEAL::3.600000000 36.59823471 221.4064186 36.59823479 332.1096297 +DEAL::3.800000000 44.70118483 298.8674044 44.70118482 448.3011011 +DEAL::4.000000000 54.59815046 403.4287983 54.59815056 605.1432001 +DEAL::4.200000000 66.68633158 544.5719168 66.68633170 816.8578791 +DEAL::4.400000000 81.45086935 735.0951986 81.45086945 1102.642801 +DEAL::4.600000000 99.48431650 992.2747285 99.48431664 1488.412099 +DEAL::4.800000000 121.5104186 1339.430782 121.5104189 2009.146180 +DEAL::5.000000000 148.4131605 1808.042440 148.4131608 2712.063678 +DEAL::5.200000000 181.2722436 2440.602013 181.2722440 3660.903032 +DEAL::5.400000000 221.4064184 3294.468126 221.4064190 4941.702220 +DEAL::5.600000000 270.4264103 4447.066817 270.4264108 6670.600256 +DEAL::5.800000000 330.2995635 6002.912316 330.2995644 9004.368533 +DEAL::6.000000000 403.4287981 8103.084066 403.4287986 12154.62604 +DEAL::6.200000000 492.7490469 10938.01940 492.7490484 16407.02921 +DEAL::6.400000000 601.8450453 14764.78184 601.8450461 22147.17257 +DEAL::6.600000000 735.0951988 19930.37082 735.0952014 29895.55643 +DEAL::6.800000000 897.8473039 26903.18662 897.8473074 40354.78020 +DEAL::7.000000000 1096.633174 36315.50345 1096.633179 54473.25552 +DEAL::7.200000000 1339.430785 49020.80224 1339.430791 73531.20394 +DEAL::7.400000000 1635.984456 66171.16176 1635.984464 99256.74322 +DEAL::7.600000000 1998.195929 89321.72563 1998.195941 133982.5897 +DEAL::7.800000000 2440.602022 120571.7183 2440.602037 180857.5784 +DEAL::8.000000000 2980.958045 162754.7961 2980.958067 244132.1968 +DEAL::8.200000000 3640.950383 219695.9955 3640.950410 329543.9943 +DEAL::8.400000000 4447.066844 296558.5750 4447.066852 444837.8636 +DEAL::8.600000000 5431.659710 400312.2045 5431.659714 600468.3068 +DEAL::8.800000000 6634.244153 540364.9551 6634.244159 810547.4344 +DEAL::9.000000000 8103.084108 729416.3941 8103.084115 1094124.592 +DEAL::9.200000000 9897.129280 984609.1443 9897.129292 1476913.720 +DEAL::9.400000000 12088.38100 1329083.326 12088.38102 1993624.992 +DEAL::9.600000000 14764.78190 1794074.834 14764.78192 2691112.257 +DEAL::9.800000000 18033.74534 2421747.717 18033.74537 3632621.584 +DEAL::10.00000000 22026.46631 3269017.487 22026.46631 4903526.238 +DEAL::10.00000000 22026.46631 3269017.487 22026.46631 4903526.238 -- 2.39.5