Add a test.

diff --git a/tests/sundials/ida_03.cc b/tests/sundials/ida_03.cc
new file mode 100644 (file)
index 0000000..1273517
--- /dev/null
+++ b/tests/sundials/ida_03.cc
@@ -0,0 +1,138 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2017 - 2022 by the deal.II authors
+//
+//    This file is part of the deal.II library.
+//
+//    The deal.II library is free software; you can use it, redistribute
+//    it, and/or modify it under the terms of the GNU Lesser General
+//    Public License as published by the Free Software Foundation; either
+//    version 2.1 of the License, or (at your option) any later version.
+//    The full text of the license can be found in the file LICENSE.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 the exponential decay problem:
+ *
+ * y' = -k y
+ * y (0) = 1
+ *
+ * with k=2. That is
+ *
+ * F(y', y, t) = y' + k y = 0
+ * y_0  = 1
+ *
+ * The exact solution is
+ *
+ * y(t) = exp(-k*t)
+ *
+ * The Jacobian we then need to assemble is the following:
+ *
+ * J = dF/dy + alpha dF/dy'
+ *   = k + alpha 1
+ *   = (k+alpha)
+ */
+class HarmonicOscillator
+{
+public:
+  HarmonicOscillator(
+    double                                                        _kappa,
+    const typename SUNDIALS::IDA<Vector<double>>::AdditionalData &data)
+    : time_stepper(data)
+    , y(1)
+    , y_dot(1)
+    , kappa(_kappa)
+  {
+    using VectorType = Vector<double>;
+
+    time_stepper.reinit_vector = [&](VectorType &v) { v.reinit(1); };
+
+
+    time_stepper.residual = [&](const double      t,
+                                const VectorType &y,
+                                const VectorType &y_dot,
+                                VectorType &      res) {
+      res[0] = y_dot[0] + kappa * y[0];
+    };
+
+    time_stepper.setup_jacobian = [&](const double,
+                                      const VectorType &,
+                                      const VectorType &,
+                                      const double alpha) {
+      J = kappa + alpha;
+    };
+
+    // Used only in ver < 4.0.0
+    time_stepper.solve_jacobian_system =
+      [&](const VectorType &src, VectorType &dst) { dst[0] = src[0] / J; };
+
+    // Used in ver >= 4.0.0
+    time_stepper.solve_with_jacobian =
+      [&](const VectorType &src, VectorType &dst, const double) {
+        dst[0] = src[0] / J;
+      };
+
+    time_stepper.output_step = [&](const double       t,
+                                   const VectorType & sol,
+                                   const VectorType & sol_dot,
+                                   const unsigned int step_number) {
+      deallog << "Intermediate output:" << std::endl;
+      deallog << "  t =" << t << std::endl;
+      deallog << "  y =" << sol[0] << "  (exact: " << std::exp(-kappa * t)
+              << ')' << std::endl;
+      deallog << "  y'=" << sol_dot[0]
+              << "  (exact: " << -kappa * std::exp(-kappa * t) << ')'
+              << std::endl;
+    };
+  }
+
+  void
+  run()
+  {
+    y[0]     = 1;
+    y_dot[0] = -kappa;
+    time_stepper.solve_dae(y, y_dot);
+  }
+  SUNDIALS::IDA<Vector<double>> time_stepper;
+
+private:
+  Vector<double> y;
+  Vector<double> y_dot;
+  double         J;
+  double         kappa;
+};
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  SUNDIALS::IDA<Vector<double>>::AdditionalData data;
+  ParameterHandler                              prm;
+  data.add_parameters(prm);
+
+  // std::ofstream ofile(SOURCE_DIR "/ida_03.prm");
+  // prm.print_parameters(ofile, ParameterHandler::ShortText);
+  // ofile.close();
+
+  std::ifstream ifile(SOURCE_DIR "/ida_03_in.prm");
+  prm.parse_input(ifile);
+
+
+  HarmonicOscillator ode(2.0, data);
+  ode.run();
+}

diff --git a/tests/sundials/ida_03.output b/tests/sundials/ida_03.output
new file mode 100644 (file)
index 0000000..2b3ff26
--- /dev/null
+++ b/tests/sundials/ida_03.output
@@ -0,0 +1,49 @@
+
+DEAL::Intermediate output:
+DEAL::  t =0.000000000
+DEAL::  y =1.000000000  (exact: 1.000000000)
+DEAL::  y'=-2.000000000  (exact: -2.000000000)
+DEAL::Intermediate output:
+DEAL::  t =0.2000000000
+DEAL::  y =0.6703200457  (exact: 0.6703200460)
+DEAL::  y'=-1.340640088  (exact: -1.340640092)
+DEAL::Intermediate output:
+DEAL::  t =0.4000000000
+DEAL::  y =0.4493289638  (exact: 0.4493289641)
+DEAL::  y'=-0.8986579256  (exact: -0.8986579282)
+DEAL::Intermediate output:
+DEAL::  t =0.6000000000
+DEAL::  y =0.3011942117  (exact: 0.3011942119)
+DEAL::  y'=-0.6023884233  (exact: -0.6023884238)
+DEAL::Intermediate output:
+DEAL::  t =0.8000000000
+DEAL::  y =0.2018965179  (exact: 0.2018965180)
+DEAL::  y'=-0.4037930380  (exact: -0.4037930360)
+DEAL::Intermediate output:
+DEAL::  t =1.000000000
+DEAL::  y =0.1353352831  (exact: 0.1353352832)
+DEAL::  y'=-0.2706705663  (exact: -0.2706705665)
+DEAL::Intermediate output:
+DEAL::  t =1.200000000
+DEAL::  y =0.09071795314  (exact: 0.09071795329)
+DEAL::  y'=-0.1814359026  (exact: -0.1814359066)
+DEAL::Intermediate output:
+DEAL::  t =1.400000000
+DEAL::  y =0.06081006262  (exact: 0.06081006263)
+DEAL::  y'=-0.1216201252  (exact: -0.1216201253)
+DEAL::Intermediate output:
+DEAL::  t =1.600000000
+DEAL::  y =0.04076220398  (exact: 0.04076220398)
+DEAL::  y'=-0.08152441186  (exact: -0.08152440796)
+DEAL::Intermediate output:
+DEAL::  t =1.800000000
+DEAL::  y =0.02732372240  (exact: 0.02732372245)
+DEAL::  y'=-0.05464744495  (exact: -0.05464744489)
+DEAL::Intermediate output:
+DEAL::  t =2.000000000
+DEAL::  y =0.01831563874  (exact: 0.01831563889)
+DEAL::  y'=-0.03663127757  (exact: -0.03663127778)
+DEAL::Intermediate output:
+DEAL::  t =2.000000000
+DEAL::  y =0.01831563874  (exact: 0.01831563889)
+DEAL::  y'=-0.03663127757  (exact: -0.03663127778)

diff --git a/tests/sundials/ida_03_in.prm b/tests/sundials/ida_03_in.prm
new file mode 100644 (file)
index 0000000..e6fa795
--- /dev/null
+++ b/tests/sundials/ida_03_in.prm
@@ -0,0 +1,19 @@
+set Final time                        = 2
+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