Add another IDA test for DAEs.

diff --git a/tests/sundials/ida_07.cc b/tests/sundials/ida_07.cc
new file mode 100644 (file)
index 0000000..654a67d
--- /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 <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/sundials/ida.h>
+
+#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<Vector<double>>::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<double>;
+
+  VectorType         y(2);
+  VectorType         y_dot(2);
+  FullMatrix<double> J(2, 2);
+  FullMatrix<double> A(2, 2);
+  FullMatrix<double> Jinv(2, 2);
+
+  SUNDIALS::IDA<Vector<double>> 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 (file)
index 0000000..151299a
--- /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