ODE: adding BDF algorithms

Implementing BDF formula for stiff ODEs.
Orders 1 to 5 are available and tested.
The integrators can be called on GPU to
solve multiple systems in parallel.

ode/impl/KokkosODE_BDF_impl.hpp

@@ -0,0 +1,149 @@
+//@HEADER
+// ************************************************************************
+//
+//                        Kokkos v. 4.0
+//       Copyright (2022) National Technology & Engineering
+//               Solutions of Sandia, LLC (NTESS).
+//
+// Under the terms of Contract DE-NA0003525 with NTESS,
+// the U.S. Government retains certain rights in this software.
+//
+// Part of Kokkos, under the Apache License v2.0 with LLVM Exceptions.
+// See https://kokkos.org/LICENSE for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//@HEADER
+
+#ifndef KOKKOSBLAS_BDF_IMPL_HPP
+#define KOKKOSBLAS_BDF_IMPL_HPP
+
+#include "Kokkos_Core.hpp"
+
+#include "KokkosODE_Newton.hpp"
+
+namespace KokkosODE {
+namespace Impl {
+
+template <int order>
+struct BDF_table {};
+
+template <>
+struct BDF_table<1> {
+  static constexpr int order = 1;
+  Kokkos::Array<double, 2> coefficients{{-1.0, 1.0}};
+};
+
+template <>
+struct BDF_table<2> {
+  static constexpr int order = 2;
+  Kokkos::Array<double, 3> coefficients{{-4.0 / 3.0, 1.0 / 3.0, 2.0 / 3.0}};
+};
+
+template <>
+struct BDF_table<3> {
+  static constexpr int order = 3;
+  Kokkos::Array<double, 4> coefficients{
+      {-18.0 / 11.0, 9.0 / 11.0, -2.0 / 11.0, 6.0 / 11.0}};
+};
+
+template <>
+struct BDF_table<4> {
+  static constexpr int order = 4;
+  Kokkos::Array<double, 5> coefficients{
+      {-48.0 / 25.0, 36.0 / 25.0, -16.0 / 25.0, 3.0 / 25.0, 12.0 / 25.0}};
+};
+
+template <>
+struct BDF_table<5> {
+  static constexpr int order = 5;
+  Kokkos::Array<double, 6> coefficients{{-300.0 / 137.0, 300.0 / 137.0,
+                                         -200.0 / 137.0, 75.0 / 137.0,
+                                         -12.0 / 137.0, 60.0 / 137.0}};
+};
+
+template <>
+struct BDF_table<6> {
+  static constexpr int order = 6;
+  Kokkos::Array<double, 7> coefficients{
+      {-360.0 / 147.0, 450.0 / 147.0, -400.0 / 147.0, 225.0 / 147.0,
+       -72.0 / 147.0, 10.0 / 147.0, 60.0 / 147.0}};
+};
+
+template <class system_type, class table_type, class mv_type>
+struct BDF_system_wrapper {
+  const system_type mySys;
+  const int neqs;
+  const table_type table;
+  const int order = table.order;
+
+  double t, dt;
+  mv_type yn;
+
+  KOKKOS_FUNCTION
+  BDF_system_wrapper(const system_type& mySys_, const table_type& table_,
+                     const double t_, const double dt_, const mv_type& yn_)
+      : mySys(mySys_),
+        neqs(mySys_.neqs),
+        table(table_),
+        t(t_),
+        dt(dt_),
+        yn(yn_) {}
+
+  template <class vec_type>
+  KOKKOS_FUNCTION void residual(const vec_type& y, const vec_type& f) const {
+    // f = f(t+dt, y)
+    mySys.evaluate_function(t, dt, y, f);
+
+    for (int eqIdx = 0; eqIdx < neqs; ++eqIdx) {
+      f(eqIdx) = y(eqIdx) - table.coefficients[order] * dt * f(eqIdx);
+      for (int orderIdx = 0; orderIdx < order; ++orderIdx) {
+        f(eqIdx) +=
+            table.coefficients[order - 1 - orderIdx] * yn(eqIdx, orderIdx);
+      }
+    }
+  }
+
+  template <class vec_type, class mat_type>
+  KOKKOS_FUNCTION void jacobian(const vec_type& y, const mat_type& jac) const {
+    mySys.evaluate_jacobian(t, dt, y, jac);
+
+    for (int rowIdx = 0; rowIdx < neqs; ++rowIdx) {
+      for (int colIdx = 0; colIdx < neqs; ++colIdx) {
+        jac(rowIdx, colIdx) =
+            -table.coefficients[order] * dt * jac(rowIdx, colIdx);
+      }
+      jac(rowIdx, rowIdx) += 1.0;
+    }
+  }
+};
+
+template <class ode_type, class table_type, class vec_type, class mv_type,
+          class mat_type, class scalar_type>
+KOKKOS_FUNCTION void BDFStep(ode_type& ode, const table_type& table,
+                             scalar_type t, scalar_type dt,
+                             const vec_type& y_old, const vec_type& y_new,
+                             const vec_type& rhs, const vec_type& update,
+                             const mv_type& y_vecs, const mat_type& temp,
+                             const mat_type& jac) {
+  using newton_params = KokkosODE::Experimental::Newton_params;
+
+  BDF_system_wrapper sys(ode, table, t, dt, y_vecs);
+  const newton_params param(50, 1e-14, 1e-12);
+
+  // first set y_new = y_old
+  for (int eqIdx = 0; eqIdx < sys.neqs; ++eqIdx) {
+    y_new(eqIdx) = y_old(eqIdx);
+  }
+
+  // solver the nonlinear problem
+  {
+    KokkosODE::Experimental::Newton::Solve(sys, param, jac, temp, y_new, rhs,
+                                           update);
+  }
+
+}  // BDFStep
+
+}  // namespace Impl
+}  // namespace KokkosODE
+
+#endif  // KOKKOSBLAS_BDF_IMPL_HPP

ode/impl/KokkosODE_Newton_impl.hpp

@@ -41,7 +41,9 @@ KOKKOS_FUNCTION KokkosODE::Experimental::newton_solver_status NewtonSolve(
   // the norm of the residual.
   using norm_type = typename Kokkos::Details::InnerProductSpaceTraits<
       typename vec_type::non_const_value_type>::mag_type;
-  norm_type norm = Kokkos::ArithTraits<norm_type>::zero();
+  sys.residual(y0, rhs);
+  const norm_type norm0 = KokkosBlas::serial_nrm2(rhs);
+  norm_type norm        = Kokkos::ArithTraits<norm_type>::zero();
 
   // LBV - 07/24/2023: for now assume that we take
   // a full Newton step. Eventually this value can
@@ -59,7 +61,7 @@ KOKKOS_FUNCTION KokkosODE::Experimental::newton_solver_status NewtonSolve(
     // with J=du/dx, rhs=f(u_n+update)-f(u_n)
     norm = KokkosBlas::serial_nrm2(rhs);
 
-    if ((norm < params.rel_tol) ||
+    if ((norm < (params.rel_tol * norm0)) ||
         (it > 0 ? KokkosBlas::serial_nrm2(update) < params.abs_tol : false)) {
       return newton_solver_status::NLS_SUCCESS;
     }

ode/src/KokkosODE_BDF.hpp

@@ -0,0 +1,135 @@
+//@HEADER
+// ************************************************************************
+//
+//                        Kokkos v. 4.0
+//       Copyright (2022) National Technology & Engineering
+//               Solutions of Sandia, LLC (NTESS).
+//
+// Under the terms of Contract DE-NA0003525 with NTESS,
+// the U.S. Government retains certain rights in this software.
+//
+// Part of Kokkos, under the Apache License v2.0 with LLVM Exceptions.
+// See https://kokkos.org/LICENSE for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//@HEADER
+
+#ifndef KOKKOSODE_BDF_HPP
+#define KOKKOSODE_BDF_HPP
+
+/// \author Luc Berger-Vergiat ([email protected])
+/// \file KokkosODE_BDF.hpp
+
+#include "Kokkos_Core.hpp"
+#include "KokkosODE_Types.hpp"
+#include "KokkosODE_RungeKutta.hpp"
+
+#include "KokkosODE_BDF_impl.hpp"
+
+namespace KokkosODE {
+namespace Experimental {
+
+enum BDF_type : int {
+  BDF1 = 0,
+  BDF2 = 1,
+  BDF3 = 2,
+  BDF4 = 3,
+  BDF5 = 4,
+  BDF6 = 5
+};
+
+template <BDF_type T>
+struct BDF_coeff_helper {
+  using table_type = void;
+};
+
+template <>
+struct BDF_coeff_helper<BDF_type::BDF1> {
+  using table_type = KokkosODE::Impl::BDF_table<1>;
+};
+
+template <>
+struct BDF_coeff_helper<BDF_type::BDF2> {
+  using table_type = KokkosODE::Impl::BDF_table<2>;
+};
+
+template <>
+struct BDF_coeff_helper<BDF_type::BDF3> {
+  using table_type = KokkosODE::Impl::BDF_table<3>;
+};
+
+template <>
+struct BDF_coeff_helper<BDF_type::BDF4> {
+  using table_type = KokkosODE::Impl::BDF_table<4>;
+};
+
+template <>
+struct BDF_coeff_helper<BDF_type::BDF5> {
+  using table_type = KokkosODE::Impl::BDF_table<5>;
+};
+
+template <>
+struct BDF_coeff_helper<BDF_type::BDF6> {
+  using table_type = KokkosODE::Impl::BDF_table<6>;
+};
+
+template <BDF_type T>
+struct BDF {
+  using table_type = typename BDF_coeff_helper<T>::table_type;
+
+  template <class ode_type, class vec_type, class mv_type, class mat_type,
+            class scalar_type>
+  KOKKOS_FUNCTION static void Solve(
+      const ode_type& ode, const scalar_type t_start, const scalar_type t_end,
+      const int num_steps, const vec_type& y0, const vec_type& y,
+      const vec_type& rhs, const vec_type& update, const mv_type& y_vecs,
+      const mat_type& temp, const mat_type& jac) {
+    const table_type table;
+
+    const double dt = (t_end - t_start) / num_steps;
+    double t        = t_start;
+
+    // Load y0 into y_vecs(:, 0)
+    for (int eqIdx = 0; eqIdx < ode.neqs; ++eqIdx) {
+      y_vecs(eqIdx, 0) = y0(eqIdx);
+    }
+
+    // Compute initial start-up history vectors
+    // Using a non adaptive explicit method.
+    const int init_steps = table.order - 1;
+    if (num_steps < init_steps) {
+      return;
+    }
+    KokkosODE::Experimental::ODE_params params(table.order - 1);
+    for (int stepIdx = 0; stepIdx < init_steps; ++stepIdx) {
+      KokkosODE::Experimental::RungeKutta<RK_type::RKF45>::Solve(
+          ode, params, t, t + dt, y0, y, update, temp);
+
+      for (int eqIdx = 0; eqIdx < ode.neqs; ++eqIdx) {
+        y_vecs(eqIdx, stepIdx + 1) = y(eqIdx);
+        y0(eqIdx)                  = y(eqIdx);
+      }
+      t += dt;
+    }
+
+    for (int stepIdx = init_steps; stepIdx < num_steps; ++stepIdx) {
+      KokkosODE::Impl::BDFStep(ode, table, t, dt, y0, y, rhs, update, y_vecs,
+                               temp, jac);
+
+      // Update history
+      for (int eqIdx = 0; eqIdx < ode.neqs; ++eqIdx) {
+        y0(eqIdx) = y(eqIdx);
+        for (int orderIdx = 0; orderIdx < table.order - 1; ++orderIdx) {
+          y_vecs(eqIdx, orderIdx) = y_vecs(eqIdx, orderIdx + 1);
+        }
+        y_vecs(eqIdx, table.order - 1) = y(eqIdx);
+      }
+      t += dt;
+    }
+  }
+};
+
+}  // namespace Experimental
+}  // namespace KokkosODE
+
+#endif  // KOKKOSODE_BDF_HPP

ode/src/KokkosODE_Types.hpp

@@ -57,9 +57,11 @@ struct Newton_params {
   int max_iters;
   double abs_tol, rel_tol;
 
-  // Constructor that only specify the desired number of steps.
-  // In this case no adaptivity is provided, the time step will
-  // be constant such that dt = (tend - tstart) / num_steps;
+  // Constructor that sets basic solver parameters
+  // used while solving the nonlinear system
+  // int max_iters_  [in]: maximum number of iterations allowed
+  // double abs_tol_ [in]: absolute tolerance to reach for successful solve
+  // double rel_tol_ [in]: relative tolerance to reach for successful solve
   KOKKOS_FUNCTION
   Newton_params(const int max_iters_, const double abs_tol_,
                 const double rel_tol_)

ode/unit_test/Test_ODE.hpp

@@ -22,5 +22,6 @@
 
 // Implicit integrators
 #include "Test_ODE_Newton.hpp"
+#include "Test_ODE_BDF.hpp"
 
 #endif  // TEST_ODE_HPP