Add more SUNDIALS::IDA tests.

diff --git a/tests/sundials/ida_06.cc b/tests/sundials/ida_06.cc
new file mode 100644 (file)
index 0000000..c7e2ec5
--- /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 <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
+ *   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<Vector<double>>::AdditionalData &data)
+    : time_stepper(data)
+    , y(2)
+    , y_dot(2)
+    , J(2, 2)
+    , A(2, 2)
+    , Jinv(2, 2)
+    , a(a_)
+  {
+    using VectorType = Vector<double>;
+
+    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<Vector<double>> time_stepper;
+
+private:
+  Vector<double>     y;
+  Vector<double>     y_dot;
+  FullMatrix<double> J;
+  FullMatrix<double> A;
+  FullMatrix<double> Jinv;
+  double             a;
+};
+
+
+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);
+
+
+  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 (file)
index 0000000..1e98254
--- /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 (file)
index 0000000..fc46a20
--- /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