Add another IDA test, this time for lacking initial velocities and one initial condition.

diff --git a/tests/sundials/ida_09.cc b/tests/sundials/ida_09.cc
new file mode 100644 (file)
index 0000000..4bc43ca
--- /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 <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"
+
+
+// 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<Vector<double>>::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<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;
+  };
+
+  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 (file)
index 0000000..01fa9e1
--- /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