]> https://gitweb.dealii.org/ - dealii.git/commitdiff
ARKode: single step solve option 11810/head
authorSebastian Proell <proell@lnm.mw.tum.de>
Thu, 25 Feb 2021 08:28:58 +0000 (09:28 +0100)
committerSebastian Proell <proell@lnm.mw.tum.de>
Wed, 3 Mar 2021 17:33:43 +0000 (18:33 +0100)
include/deal.II/sundials/arkode.h
source/sundials/arkode.cc
tests/sundials/arkode_repeated_solve.cc [new file with mode: 0644]
tests/sundials/arkode_repeated_solve.prm [new file with mode: 0644]
tests/sundials/arkode_repeated_solve.with_sundials.geq.5.4.0.output [new file with mode: 0644]

index a6b59183123ed9983a1629392391c0b8c9885cc0..27a91311cae52882d0937d9c8d1a601c1205b5de 100644 (file)
@@ -41,6 +41,8 @@
 #  ifdef DEAL_II_WITH_MPI
 #    include <nvector/nvector_parallel.h>
 #  endif
+#  include <deal.II/base/discrete_time.h>
+
 #  include <deal.II/sundials/n_vector.h>
 #  include <deal.II/sundials/sunlinsol_wrapper.h>
 
@@ -535,10 +537,44 @@ namespace SUNDIALS
     /**
      * Integrate the initial value problem. This function returns the final
      * number of computed steps.
+     *
+     * @param solution On input, this vector contains the initial condition. On
+     *   output, it contains the solution at the final time.
      */
     unsigned int
     solve_ode(VectorType &solution);
 
+    /**
+     * Integrate the initial value problem. Compared to the function above, this
+     * function allows to specify an @p intermediate_time for the next solution.
+     * Repeated calls of this function must use monotonously increasing values
+     * for @p intermediate_time. The last solution state is saved internally
+     * along with the @p intermediate_time and will be reused as initial
+     * condition for the next call.
+     *
+     * Users may find this function useful when integrating ARKode into an outer
+     * time loop of their own, especially when output_step() is too restrictive.
+     *
+     * @note @p intermediate_time may be larger than AdditionalData::final_time,
+     *   which is ignored by this function.
+     *
+     * @param solution The final solution. If the solver restarts, either
+     *   because it is the first ever solve or the flag @p reset_solver is
+     *   set, the vector is also used as initial condition.
+     * @param intermediate_time The time for the incremental solution step. Must
+     *   be greater than the last time that was used in a previous call to this
+     *   function.
+     * @param reset_solver Optional flag to recreate all internal objects which
+     *   may be desirable for spatial adaptivity methods. If set to `true`,
+     *   reset() is called before solving the ODE, which sets @p solution as
+     *   initial condition. This will *not* reset the stored time from previous
+     *   calls to this function.
+     */
+    unsigned int
+    solve_ode_incrementally(VectorType & solution,
+                            const double intermediate_time,
+                            const bool   reset_solver = false);
+
     /**
      * Clear internal memory and start with clean objects. This function is
      * called when the simulation starts and when the user returns true to a
@@ -549,12 +585,12 @@ namespace SUNDIALS
      * a different function to solver_should_restart() that performs all mesh
      * changes, transfers the solution to the new mesh, and returns true.
      *
-     * @param[in] t  The new starting time
-     * @param[in] h  The new starting time step
-     * @param[in,out] y   The new initial solution
+     * @param t  The new starting time
+     * @param h  The new starting time step
+     * @param y  The new initial solution
      */
     void
-    reset(const double t, const double h, VectorType &y);
+    reset(const double t, const double h, const VectorType &y);
 
     /**
      * Provides user access to the internally used ARKODE memory.
@@ -1224,6 +1260,14 @@ namespace SUNDIALS
                    << "Please provide an implementation for the function \""
                    << arg1 << "\"");
 
+    /**
+     * Internal routine to call ARKode repeatedly.
+     */
+    int
+    do_evolve_time(VectorType &          solution,
+                   dealii::DiscreteTime &time,
+                   const bool            do_reset);
+
 #  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
 
     /**
@@ -1238,9 +1282,11 @@ namespace SUNDIALS
     /**
      * Set up the solver and preconditioner for a non-identity mass matrix in
      * the ARKODE memory object based on the user-specified functions.
+     * @param solution The solution vector which is used as a template to create
+     *   new vectors.
      */
     void
-    setup_mass_solver();
+    setup_mass_solver(const VectorType &solution);
 
 #  endif
 
@@ -1262,16 +1308,6 @@ namespace SUNDIALS
      */
     void *arkode_mem;
 
-    /**
-     * ARKode solution vector.
-     */
-    internal::NVectorView<VectorType> yy;
-
-    /**
-     * ARKode absolute tolerances vector.
-     */
-    internal::NVectorView<VectorType> abs_tolls;
-
     /**
      * MPI communicator. SUNDIALS solver runs happily in
      * parallel. Note that if the library is compiled without MPI
@@ -1280,9 +1316,9 @@ namespace SUNDIALS
     MPI_Comm communicator;
 
     /**
-     * Memory pool of vectors.
+     * The final time in the last call to solve_ode().
      */
-    GrowingVectorMemory<VectorType> mem;
+    double last_end_time;
 
 #  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
     std::unique_ptr<internal::LinearSolverWrapper<VectorType>> linear_solver;
index 31d970927c4249a55932ea3fb76b323d6030fb62..f143a2d43bb60bdd420fb3dd76eb51d78d214669 100644 (file)
@@ -226,7 +226,6 @@ namespace SUNDIALS
       Assert(user_data != nullptr, ExcInternalError());
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(user_data);
-      GrowingVectorMemory<VectorType> mem;
 
       auto *src_v  = internal::unwrap_nvector_const<VectorType>(v);
       auto *src_y  = internal::unwrap_nvector_const<VectorType>(y);
@@ -385,6 +384,7 @@ namespace SUNDIALS
     , communicator(is_serial_vector<VectorType>::value ?
                      MPI_COMM_SELF :
                      Utilities::MPI::duplicate_communicator(mpi_comm))
+    , last_end_time(data.initial_time)
   {
     set_functions_to_trigger_an_assert();
   }
@@ -418,21 +418,64 @@ namespace SUNDIALS
     DiscreteTime time(data.initial_time,
                       data.final_time,
                       data.initial_step_size);
-    reset(time.get_current_time(), time.get_next_step_size(), solution);
 
-    if (output_step)
-      output_step(time.get_current_time(), solution, time.get_step_number());
+    return do_evolve_time(solution, time, /* force_solver_restart = */ true);
+  }
+
 
-    while (time.is_at_end() == false)
+
+  template <typename VectorType>
+  unsigned int
+  ARKode<VectorType>::solve_ode_incrementally(VectorType & solution,
+                                              const double intermediate_time,
+                                              const bool   reset_solver)
+  {
+    AssertThrow(
+      intermediate_time > last_end_time,
+      dealii::ExcMessage(
+        "The requested intermediate time is smaller than the last requested "
+        "intermediate time."));
+
+    const bool   do_reset = reset_solver || arkode_mem == nullptr;
+    DiscreteTime time(last_end_time, intermediate_time, data.initial_step_size);
+    return do_evolve_time(solution, time, do_reset);
+  }
+
+
+
+  template <typename VectorType>
+  int
+  ARKode<VectorType>::do_evolve_time(VectorType &  solution,
+                                     DiscreteTime &time,
+                                     const bool    do_reset)
+  {
+    if (do_reset)
+      {
+        reset(time.get_current_time(), time.get_next_step_size(), solution);
+        if (output_step)
+          output_step(time.get_current_time(),
+                      solution,
+                      time.get_step_number());
+      }
+
+    auto solution_nvector = internal::make_nvector_view(solution);
+
+    while (!time.is_at_end())
       {
         time.set_desired_next_step_size(data.output_period);
         double actual_next_time;
 #  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-        const auto status = SundialsARKode(
-          arkode_mem, time.get_next_time(), yy, &actual_next_time, ARK_NORMAL);
+        const auto status = SundialsARKode(arkode_mem,
+                                           time.get_next_time(),
+                                           solution_nvector,
+                                           &actual_next_time,
+                                           ARK_NORMAL);
 #  else
-        const auto status = ARKStepEvolve(
-          arkode_mem, time.get_next_time(), yy, &actual_next_time, ARK_NORMAL);
+        const auto status = ARKStepEvolve(arkode_mem,
+                                          time.get_next_time(),
+                                          solution_nvector,
+                                          &actual_next_time,
+                                          ARK_NORMAL);
 #  endif
         (void)status;
         AssertARKode(status);
@@ -450,17 +493,20 @@ namespace SUNDIALS
                       solution,
                       time.get_step_number());
       }
-
+    last_end_time = time.get_current_time();
     return time.get_step_number();
   }
 
+
+
 #  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
   template <typename VectorType>
   void
-  ARKode<VectorType>::reset(const double current_time,
-                            const double current_time_step,
-                            VectorType & solution)
+  ARKode<VectorType>::reset(const double      current_time,
+                            const double      current_time_step,
+                            const VectorType &solution)
   {
+    last_end_time = current_time;
     if (arkode_mem)
       ARKodeFree(&arkode_mem);
 
@@ -474,21 +520,21 @@ namespace SUNDIALS
 
     // just a view on the memory in solution, all write operations on yy by
     // ARKODE will automatically be mirrored to solution
-    yy = internal::make_nvector_view(solution);
+    auto initial_condition_nvector = internal::make_nvector_view(solution);
 
     status = ARKodeInit(
       arkode_mem,
       explicit_function ? &t_arkode_explicit_function<VectorType> : nullptr,
       implicit_function ? &t_arkode_implicit_function<VectorType> : nullptr,
       current_time,
-      yy);
+      initial_condition_nvector);
     AssertARKode(status);
 
     if (get_local_tolerances)
       {
-        abs_tolls = make_nvector_view(get_local_tolerances());
+        const auto abs_tols = make_nvector_view(get_local_tolerances());
         status =
-          ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls);
+          ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
         AssertARKode(status);
       }
     else
@@ -570,8 +616,9 @@ namespace SUNDIALS
   void
   ARKode<VectorType>::reset(const double current_time,
                             const double current_time_step,
-                            VectorType &solution)
+                            const VectorType &solution)
   {
+    last_end_time = current_time;
     if (arkode_mem)
       ARKStepFree(&arkode_mem);
 
@@ -580,7 +627,7 @@ namespace SUNDIALS
 
     // just a view on the memory in solution, all write operations on yy by
     // ARKODE will automatically be mirrored to solution
-    yy = internal::make_nvector_view(solution);
+    auto initial_condition_nvector = internal::make_nvector_view(solution);
 
     Assert(explicit_function || implicit_function,
            ExcFunctionNotProvided("explicit_function || implicit_function"));
@@ -589,15 +636,16 @@ namespace SUNDIALS
       explicit_function ? &t_arkode_explicit_function<VectorType> : nullptr,
       implicit_function ? &t_arkode_implicit_function<VectorType> : nullptr,
       current_time,
-      yy);
+      initial_condition_nvector);
 
     Assert(arkode_mem != nullptr, ExcInternalError());
 
     if (get_local_tolerances)
       {
-        abs_tolls = internal::make_nvector_view(get_local_tolerances());
+        const auto abs_tols =
+          internal::make_nvector_view(get_local_tolerances());
         status =
-          ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls);
+          ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
         AssertARKode(status);
       }
     else
@@ -614,7 +662,7 @@ namespace SUNDIALS
 
     setup_system_solver(solution);
 
-    setup_mass_solver();
+    setup_mass_solver(solution);
 
     status = ARKStepSetInitStep(arkode_mem, current_time_step);
     AssertARKode(status);
@@ -657,8 +705,9 @@ namespace SUNDIALS
           {
             // use default solver from SUNDIALS
             // TODO give user options
+            auto y_template = internal::make_nvector_view(solution);
             sun_linear_solver =
-              SUNLinSol_SPGMR(yy,
+              SUNLinSol_SPGMR(y_template,
                               PREC_NONE,
                               0 /*krylov subvectors, 0 uses default*/);
           }
@@ -709,7 +758,7 @@ namespace SUNDIALS
 
   template <typename VectorType>
   void
-  ARKode<VectorType>::setup_mass_solver()
+  ARKode<VectorType>::setup_mass_solver(const VectorType &solution)
   {
     int status;
     (void)status;
@@ -726,8 +775,9 @@ namespace SUNDIALS
           }
         else
           {
+            auto y_template = internal::make_nvector_view(solution);
             sun_mass_linear_solver =
-              SUNLinSol_SPGMR(yy,
+              SUNLinSol_SPGMR(y_template,
                               PREC_NONE,
                               0 /*krylov subvectors, 0 uses default*/);
           }
diff --git a/tests/sundials/arkode_repeated_solve.cc b/tests/sundials/arkode_repeated_solve.cc
new file mode 100644 (file)
index 0000000..ad4a629
--- /dev/null
@@ -0,0 +1,129 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2021 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/arkode.h>
+
+#include "../tests.h"
+
+// Test repeated calls to solve_ode with harmonic oscillator problem
+
+using VectorType = Vector<double>;
+
+double kappa = 1.0;
+
+std::unique_ptr<SUNDIALS::ARKode<VectorType>>
+create_solver()
+{
+  ParameterHandler                             prm;
+  SUNDIALS::ARKode<VectorType>::AdditionalData data;
+  data.add_parameters(prm);
+
+  std::ifstream ifile(SOURCE_DIR "/arkode_repeated_solve.prm");
+  prm.parse_input(ifile);
+
+  auto ode = std::make_unique<SUNDIALS::ARKode<VectorType>>(data);
+
+  // will yield analytic solution y[0] = sin(kappa*t); y[1] = kappa*cos(kappa*t)
+  ode->explicit_function =
+    [&](double, const VectorType &y, VectorType &ydot) -> int {
+    ydot[0] = y[1];
+    ydot[1] = -kappa * kappa * y[0];
+    return 0;
+  };
+
+  ode->output_step = [&](const double       t,
+                         const VectorType & sol,
+                         const unsigned int step_number) -> int {
+    deallog << t << " " << sol[0] << " " << sol[1] << std::endl;
+    return 0;
+  };
+  return ode;
+}
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, numbers::invalid_unsigned_int);
+
+  auto ode = create_solver();
+  {
+    Vector<double> y(2);
+    y[0] = 0;
+    y[1] = kappa;
+
+    // solve the ode in one go
+    ode->solve_ode(y);
+
+    deallog << "First solve done" << std::endl;
+
+    // solve again -> will call reset automatically
+    y[0] = 0;
+    y[1] = kappa;
+    ode->solve_ode(y);
+    deallog << "Second solve done" << std::endl;
+
+    // Solve for a time greater than the originally requested final time
+    ode->solve_ode_incrementally(y, 3.0);
+    deallog << "Solve further than final time done" << std::endl;
+  }
+
+  // solve in two steps with new object
+  ode = create_solver();
+
+  {
+    Vector<double> y(2);
+    y[0] = 0;
+    y[1] = kappa;
+
+    ode->solve_ode_incrementally(y, 1.0);
+    deallog << "Split solve - part 1 done" << std::endl;
+  }
+
+  {
+    // N.B. initial conditions are not necessary as the internal state is reused
+    Vector<double> y(2);
+
+    ode->solve_ode_incrementally(y, 2.0);
+    deallog << "Split solve - part 2 done" << std::endl;
+  }
+
+  {
+    Vector<double> y(2);
+    y[0] = 0;
+    y[1] = kappa;
+    // solving with single-argument call reruns from start
+    ode->solve_ode(y);
+    deallog << "Rerun with single argument call done" << std::endl;
+  }
+
+  {
+    Vector<double> y0(2), y(2);
+    y0[0] = 0;
+    y0[1] = kappa;
+    // manually reset solver
+    ode->reset(0.0, 0.01, y0);
+    ode->solve_ode_incrementally(y, 2.0);
+  }
+
+  return 0;
+}
diff --git a/tests/sundials/arkode_repeated_solve.prm b/tests/sundials/arkode_repeated_solve.prm
new file mode 100644 (file)
index 0000000..2bdb712
--- /dev/null
@@ -0,0 +1,13 @@
+set Final time                        = 2.0
+set Initial time                      = 0
+set Time interval between each output = 0.1
+subsection Error control
+  set Absolute error tolerance = 0.000001
+  set Relative error tolerance = 0.000010
+end
+subsection Running parameters
+  set Initial step size                      = 0.010000
+  set Maximum number of nonlinear iterations = 10
+  set Maximum order of ARK                   = 5
+  set Minimum step size                      = 0.000001
+end
diff --git a/tests/sundials/arkode_repeated_solve.with_sundials.geq.5.4.0.output b/tests/sundials/arkode_repeated_solve.with_sundials.geq.5.4.0.output
new file mode 100644 (file)
index 0000000..3f6116f
--- /dev/null
@@ -0,0 +1,121 @@
+
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.100000 0.0998334 0.995004
+DEAL::0.200000 0.198669 0.980067
+DEAL::0.300000 0.295521 0.955336
+DEAL::0.400000 0.389420 0.921060
+DEAL::0.500000 0.479426 0.877583
+DEAL::0.600000 0.564643 0.825335
+DEAL::0.700000 0.644219 0.764841
+DEAL::0.800000 0.717356 0.696707
+DEAL::0.900000 0.783330 0.621606
+DEAL::1.00000 0.841478 0.540292
+DEAL::1.10000 0.891214 0.453587
+DEAL::1.20000 0.932042 0.362356
+DEAL::1.30000 0.963560 0.267498
+DEAL::1.40000 0.985452 0.169958
+DEAL::1.50000 0.997498 0.0707265
+DEAL::1.60000 0.999576 -0.0292043
+DEAL::1.70000 0.991667 -0.128845
+DEAL::1.80000 0.973849 -0.227205
+DEAL::1.90000 0.946302 -0.323292
+DEAL::2.00000 0.909300 -0.416148
+DEAL::First solve done
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.100000 0.0998334 0.995004
+DEAL::0.200000 0.198669 0.980067
+DEAL::0.300000 0.295521 0.955336
+DEAL::0.400000 0.389420 0.921060
+DEAL::0.500000 0.479426 0.877583
+DEAL::0.600000 0.564643 0.825335
+DEAL::0.700000 0.644219 0.764841
+DEAL::0.800000 0.717356 0.696707
+DEAL::0.900000 0.783330 0.621606
+DEAL::1.00000 0.841478 0.540292
+DEAL::1.10000 0.891214 0.453587
+DEAL::1.20000 0.932042 0.362356
+DEAL::1.30000 0.963560 0.267498
+DEAL::1.40000 0.985452 0.169958
+DEAL::1.50000 0.997498 0.0707265
+DEAL::1.60000 0.999576 -0.0292043
+DEAL::1.70000 0.991667 -0.128845
+DEAL::1.80000 0.973849 -0.227205
+DEAL::1.90000 0.946302 -0.323292
+DEAL::2.00000 0.909300 -0.416148
+DEAL::Second solve done
+DEAL::2.10000 0.863210 -0.504849
+DEAL::2.20000 0.808497 -0.588504
+DEAL::2.30000 0.745707 -0.666278
+DEAL::2.40000 0.675463 -0.737397
+DEAL::2.50000 0.598472 -0.801147
+DEAL::2.60000 0.515502 -0.856891
+DEAL::2.70000 0.427379 -0.904075
+DEAL::2.80000 0.334987 -0.942226
+DEAL::2.90000 0.239250 -0.970961
+DEAL::3.00000 0.141118 -0.989995
+DEAL::Solve further than final time done
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.100000 0.0998334 0.995004
+DEAL::0.200000 0.198669 0.980067
+DEAL::0.300000 0.295521 0.955336
+DEAL::0.400000 0.389420 0.921060
+DEAL::0.500000 0.479426 0.877583
+DEAL::0.600000 0.564643 0.825335
+DEAL::0.700000 0.644219 0.764841
+DEAL::0.800000 0.717356 0.696707
+DEAL::0.900000 0.783330 0.621606
+DEAL::1.00000 0.841478 0.540292
+DEAL::Split solve - part 1 done
+DEAL::1.10000 0.891214 0.453587
+DEAL::1.20000 0.932042 0.362356
+DEAL::1.30000 0.963560 0.267498
+DEAL::1.40000 0.985452 0.169958
+DEAL::1.50000 0.997498 0.0707265
+DEAL::1.60000 0.999576 -0.0292043
+DEAL::1.70000 0.991667 -0.128845
+DEAL::1.80000 0.973849 -0.227205
+DEAL::1.90000 0.946302 -0.323292
+DEAL::2.00000 0.909300 -0.416148
+DEAL::Split solve - part 2 done
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.100000 0.0998334 0.995004
+DEAL::0.200000 0.198669 0.980067
+DEAL::0.300000 0.295521 0.955336
+DEAL::0.400000 0.389420 0.921060
+DEAL::0.500000 0.479426 0.877583
+DEAL::0.600000 0.564643 0.825335
+DEAL::0.700000 0.644219 0.764841
+DEAL::0.800000 0.717356 0.696707
+DEAL::0.900000 0.783330 0.621606
+DEAL::1.00000 0.841478 0.540292
+DEAL::1.10000 0.891214 0.453587
+DEAL::1.20000 0.932042 0.362356
+DEAL::1.30000 0.963560 0.267498
+DEAL::1.40000 0.985452 0.169958
+DEAL::1.50000 0.997498 0.0707265
+DEAL::1.60000 0.999576 -0.0292043
+DEAL::1.70000 0.991667 -0.128845
+DEAL::1.80000 0.973849 -0.227205
+DEAL::1.90000 0.946302 -0.323292
+DEAL::2.00000 0.909300 -0.416148
+DEAL::Rerun with single argument call done
+DEAL::0.100000 0.0998334 0.995004
+DEAL::0.200000 0.198669 0.980067
+DEAL::0.300000 0.295521 0.955336
+DEAL::0.400000 0.389420 0.921060
+DEAL::0.500000 0.479426 0.877583
+DEAL::0.600000 0.564643 0.825335
+DEAL::0.700000 0.644219 0.764841
+DEAL::0.800000 0.717356 0.696707
+DEAL::0.900000 0.783330 0.621606
+DEAL::1.00000 0.841478 0.540292
+DEAL::1.10000 0.891214 0.453587
+DEAL::1.20000 0.932042 0.362356
+DEAL::1.30000 0.963560 0.267498
+DEAL::1.40000 0.985452 0.169958
+DEAL::1.50000 0.997498 0.0707265
+DEAL::1.60000 0.999576 -0.0292043
+DEAL::1.70000 0.991667 -0.128845
+DEAL::1.80000 0.973849 -0.227205
+DEAL::1.90000 0.946302 -0.323292
+DEAL::2.00000 0.909300 -0.416148

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.