diff --git a/ode/impl/KokkosODE_BDF_impl.hpp b/ode/impl/KokkosODE_BDF_impl.hpp new file mode 100644 index 0000000000..deee0956f9 --- /dev/null +++ b/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 +struct BDF_table {}; + +template <> +struct BDF_table<1> { + static constexpr int order = 1; + Kokkos::Array coefficients{{-1.0, 1.0}}; +}; + +template <> +struct BDF_table<2> { + static constexpr int order = 2; + Kokkos::Array 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 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 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 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 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 +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 + 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 + 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 +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 diff --git a/ode/impl/KokkosODE_Newton_impl.hpp b/ode/impl/KokkosODE_Newton_impl.hpp index f0cb90810e..06d2ca473d 100644 --- a/ode/impl/KokkosODE_Newton_impl.hpp +++ b/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::zero(); + sys.residual(y0, rhs); + const norm_type norm0 = KokkosBlas::serial_nrm2(rhs); + norm_type norm = Kokkos::ArithTraits::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; } diff --git a/ode/src/KokkosODE_BDF.hpp b/ode/src/KokkosODE_BDF.hpp new file mode 100644 index 0000000000..3697d4d0df --- /dev/null +++ b/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 (lberge@sandia.gov) +/// \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 +struct BDF_coeff_helper { + using table_type = void; +}; + +template <> +struct BDF_coeff_helper { + using table_type = KokkosODE::Impl::BDF_table<1>; +}; + +template <> +struct BDF_coeff_helper { + using table_type = KokkosODE::Impl::BDF_table<2>; +}; + +template <> +struct BDF_coeff_helper { + using table_type = KokkosODE::Impl::BDF_table<3>; +}; + +template <> +struct BDF_coeff_helper { + using table_type = KokkosODE::Impl::BDF_table<4>; +}; + +template <> +struct BDF_coeff_helper { + using table_type = KokkosODE::Impl::BDF_table<5>; +}; + +template <> +struct BDF_coeff_helper { + using table_type = KokkosODE::Impl::BDF_table<6>; +}; + +template +struct BDF { + using table_type = typename BDF_coeff_helper::table_type; + + template + 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::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 diff --git a/ode/src/KokkosODE_Types.hpp b/ode/src/KokkosODE_Types.hpp index 068c4b17ed..b18278ca63 100644 --- a/ode/src/KokkosODE_Types.hpp +++ b/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_) diff --git a/ode/unit_test/Test_ODE.hpp b/ode/unit_test/Test_ODE.hpp index 5d4861879b..1b55171581 100644 --- a/ode/unit_test/Test_ODE.hpp +++ b/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 diff --git a/ode/unit_test/Test_ODE_BDF.hpp b/ode/unit_test/Test_ODE_BDF.hpp new file mode 100644 index 0000000000..99111a1fad --- /dev/null +++ b/ode/unit_test/Test_ODE_BDF.hpp @@ -0,0 +1,502 @@ +//@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 + +#include +#include "KokkosKernels_TestUtils.hpp" + +#include "KokkosODE_BDF.hpp" + +namespace Test { + +// Logistic equation +// Used to model population growth +// it is a simple nonlinear ODE with +// a lot of literature. +// +// Equation: y'(t) = r*y(t)*(1-y(t)/K) +// Jacobian: df/dy = r - 2*r*y/K +// Solution: y = K / (1 + ((K - y0) / y0)*exp(-rt)) +struct Logistic { + static constexpr int neqs = 1; + + const double r, K; + + Logistic(double r_, double K_) : r(r_), K(K_){}; + + template + KOKKOS_FUNCTION void evaluate_function(const double /*t*/, + const double /*dt*/, + const vec_type1& y, + const vec_type2& f) const { + f(0) = r * y(0) * (1.0 - y(0) / K); + } + + template + KOKKOS_FUNCTION void evaluate_jacobian(const double /*t*/, + const double /*dt*/, const vec_type& y, + const mat_type& jac) const { + jac(0, 0) = r - 2 * r * y(0) / K; + } + + template + KOKKOS_FUNCTION void solution(const double t, const vec_type& y0, + const vec_type& y) const { + y(0) = K / (1 + (K - y0) / y0 * Kokkos::exp(-r * t)); + } + +}; // Logistic + +// Lotka-Volterra equation +// A predator-prey model that describe +// population dynamics when two species +// interact. +// +// Equations: y0'(t) = alpha*y0(t) - beta*y0(t)*y1(t) +// y1'(t) = delta*y0(t)*y1(t) - gamma*y1(t) +// Jacobian: df0/dy = [alpha-beta*y1(t); beta*y0(t)] +// df1/dy = [delta*y1(t); delta*y0(t)-gamma] +// Solution: y = K / (1 + ((K - y0) / y0)*exp(-rt)) +struct LotkaVolterra { + static constexpr int neqs = 2; + + const double alpha, beta, delta, gamma; + + LotkaVolterra(double alpha_, double beta_, double delta_, double gamma_) + : alpha(alpha_), beta(beta_), delta(delta_), gamma(gamma_){}; + + template + KOKKOS_FUNCTION void evaluate_function(const double /*t*/, + const double /*dt*/, + const vec_type1& y, + const vec_type2& f) const { + f(0) = alpha * y(0) - beta * y(0) * y(1); + f(1) = delta * y(0) * y(1) - gamma * y(1); + } + + template + KOKKOS_FUNCTION void evaluate_jacobian(const double /*t*/, + const double /*dt*/, const vec_type& y, + const mat_type& jac) const { + jac(0, 0) = alpha - beta * y(1); + jac(0, 1) = -beta * y(0); + jac(1, 0) = delta * y(1); + jac(1, 1) = delta * y(0) - gamma; + } + +}; // LotkaVolterra + +// Robertson's autocatalytic chemical reaction: +// H. H. Robertson, The solution of a set of reaction rate equations, +// in J. Walsh (Ed.), Numerical Analysis: An Introduction, +// pp. 178–182, Academic Press, London (1966). +// +// Equations: y0' = -0.04*y0 + 10e4*y1*y2 +// y1' = 0.04*y0 - 10e4*y1*y2 - 3e-7 * y1**2 +// y2' = 3e-7 * y1**2 +struct StiffChemestry { + static constexpr int neqs = 3; + + StiffChemestry() {} + + template + KOKKOS_FUNCTION void evaluate_function(const double /*t*/, + const double /*dt*/, + const vec_type1& y, + const vec_type2& f) const { + f(0) = -0.04 * y(0) + 1.e4 * y(1) * y(2); + f(1) = 0.04 * y(0) - 1.e4 * y(1) * y(2) - 3.e7 * y(1) * y(1); + f(2) = 3.e7 * y(1) * y(1); + } + + template + KOKKOS_FUNCTION void evaluate_jacobian(const double /*t*/, + const double /*dt*/, const vec_type& y, + const mat_type& jac) const { + jac(0, 0) = -0.04; + jac(0, 1) = 1.e4 * y(2); + jac(0, 2) = 1.e4 * y(1); + jac(1, 0) = 0.04; + jac(1, 1) = -1.e4 * y(2) - 3.e7 * 2.0 * y(1); + jac(1, 2) = -1.e4 * y(1); + jac(2, 0) = 0.00; + jac(2, 1) = 3.e7 * 2.0 * y(1); + jac(2, 2) = 0.0; + } +}; + +template +struct BDFSolve_wrapper { + ode_type my_ode; + scalar_type tstart, tend; + int num_steps; + vec_type y_old, y_new, rhs, update; + mv_type y_vecs; + mat_type temp, jac; + + BDFSolve_wrapper(const ode_type& my_ode_, const scalar_type tstart_, + const scalar_type tend_, const int num_steps_, + 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_) + : my_ode(my_ode_), + tstart(tstart_), + tend(tend_), + num_steps(num_steps_), + y_old(y_old_), + y_new(y_new_), + rhs(rhs_), + update(update_), + y_vecs(y_vecs_), + temp(temp_), + jac(jac_) {} + + KOKKOS_FUNCTION + void operator()(const int /*idx*/) const { + KokkosODE::Experimental::BDF::Solve(my_ode, tstart, tend, + num_steps, y_old, y_new, rhs, + update, y_vecs, temp, jac); + } +}; + +template +void test_BDF_Logistic() { + using vec_type = Kokkos::View; + using mv_type = Kokkos::View; + using mat_type = Kokkos::View; + + Kokkos::RangePolicy myPolicy(0, 1); + Logistic mySys(1, 1); + + constexpr int num_tests = 7; + int num_steps[num_tests] = {512, 256, 128, 64, 32, 16, 8}; + double errors[num_tests] = {0}; + const scalar_type t_start = 0.0, t_end = 6.0; + vec_type y0("initial conditions", mySys.neqs), y_new("solution", mySys.neqs); + vec_type rhs("rhs", mySys.neqs), update("update", mySys.neqs); + mat_type jac("jacobian", mySys.neqs, mySys.neqs), + temp("temp storage", mySys.neqs, mySys.neqs + 4); + + scalar_type measured_order; + + // Test BDF1 +#if defined(HAVE_KOKKOSKERNELS_DEBUG) + std::cout << "\nBDF1 convergence test" << std::endl; +#endif + for (int idx = 0; idx < num_tests; idx++) { + mv_type y_vecs("history vectors", mySys.neqs, 1); + + Kokkos::deep_copy(y0, 0.5); + Kokkos::deep_copy(y_vecs, 0.5); + + BDFSolve_wrapper + solve_wrapper(mySys, t_start, t_end, num_steps[idx], y0, y_new, rhs, + update, y_vecs, temp, jac); + Kokkos::parallel_for(myPolicy, solve_wrapper); + + auto y_new_h = Kokkos::create_mirror_view(y_new); + Kokkos::deep_copy(y_new_h, y_new); + + errors[idx] = Kokkos::abs(y_new_h(0) - 1 / (1 + Kokkos::exp(-t_end))) / + Kokkos::abs(1 / (1 + Kokkos::exp(-t_end))); + } + measured_order = + Kokkos::pow(errors[num_tests - 1] / errors[0], 1.0 / (num_tests - 1)); + EXPECT_NEAR_KK_REL(measured_order, 2.0, 0.15); +#if defined(HAVE_KOKKOSKERNELS_DEBUG) + std::cout << "expected ratio: 2, actual ratio: " << measured_order + << ", order error=" << Kokkos::abs(measured_order - 2.0) / 2.0 + << std::endl; +#endif + + // Test BDF2 +#if defined(HAVE_KOKKOSKERNELS_DEBUG) + std::cout << "\nBDF2 convergence test" << std::endl; +#endif + for (int idx = 0; idx < num_tests; idx++) { + mv_type y_vecs("history vectors", mySys.neqs, 2); + Kokkos::deep_copy(y0, 0.5); + + BDFSolve_wrapper + solve_wrapper(mySys, t_start, t_end, num_steps[idx], y0, y_new, rhs, + update, y_vecs, temp, jac); + Kokkos::parallel_for(myPolicy, solve_wrapper); + + auto y_new_h = Kokkos::create_mirror_view(y_new); + Kokkos::deep_copy(y_new_h, y_new); + + errors[idx] = Kokkos::abs(y_new_h(0) - 1 / (1 + Kokkos::exp(-t_end))) / + Kokkos::abs(1 / (1 + Kokkos::exp(-t_end))); + } + measured_order = + Kokkos::pow(errors[num_tests - 1] / errors[0], 1.0 / (num_tests - 1)); + EXPECT_NEAR_KK_REL(measured_order, 4.0, 0.15); +#if defined(HAVE_KOKKOSKERNELS_DEBUG) + std::cout << "expected ratio: 4, actual ratio: " << measured_order + << ", order error=" << Kokkos::abs(measured_order - 4.0) / 4.0 + << std::endl; +#endif + + // Test BDF3 +#if defined(HAVE_KOKKOSKERNELS_DEBUG) + std::cout << "\nBDF3 convergence test" << std::endl; +#endif + for (int idx = 0; idx < num_tests; idx++) { + mv_type y_vecs("history vectors", mySys.neqs, 3); + Kokkos::deep_copy(y0, 0.5); + + BDFSolve_wrapper + solve_wrapper(mySys, t_start, t_end, num_steps[idx], y0, y_new, rhs, + update, y_vecs, temp, jac); + Kokkos::parallel_for(myPolicy, solve_wrapper); + + auto y_new_h = Kokkos::create_mirror_view(y_new); + Kokkos::deep_copy(y_new_h, y_new); + + errors[idx] = Kokkos::abs(y_new_h(0) - 1 / (1 + Kokkos::exp(-t_end))) / + Kokkos::abs(1 / (1 + Kokkos::exp(-t_end))); + } + measured_order = + Kokkos::pow(errors[num_tests - 1] / errors[0], 1.0 / (num_tests - 1)); + EXPECT_NEAR_KK_REL(measured_order, 8.0, 0.15); +#if defined(HAVE_KOKKOSKERNELS_DEBUG) + std::cout << "expected ratio: 8, actual ratio: " << measured_order + << ", order error=" << Kokkos::abs(measured_order - 8.0) / 8.0 + << std::endl; +#endif + + // Test BDF4 +#if defined(HAVE_KOKKOSKERNELS_DEBUG) + std::cout << "\nBDF4 convergence test" << std::endl; +#endif + for (int idx = 0; idx < num_tests; idx++) { + mv_type y_vecs("history vectors", mySys.neqs, 4); + Kokkos::deep_copy(y0, 0.5); + + BDFSolve_wrapper + solve_wrapper(mySys, t_start, t_end, num_steps[idx], y0, y_new, rhs, + update, y_vecs, temp, jac); + Kokkos::parallel_for(myPolicy, solve_wrapper); + + auto y_new_h = Kokkos::create_mirror_view(y_new); + Kokkos::deep_copy(y_new_h, y_new); + + errors[idx] = Kokkos::abs(y_new_h(0) - 1 / (1 + Kokkos::exp(-t_end))) / + Kokkos::abs(1 / (1 + Kokkos::exp(-t_end))); + } + measured_order = + Kokkos::pow(errors[num_tests - 1] / errors[0], 1.0 / (num_tests - 1)); +#if defined(HAVE_KOKKOSKERNELS_DEBUG) + std::cout << "expected ratio: 16, actual ratio: " << measured_order + << ", order error=" << Kokkos::abs(measured_order - 16.0) / 16.0 + << std::endl; +#endif + + // Test BDF5 +#if defined(HAVE_KOKKOSKERNELS_DEBUG) + std::cout << "\nBDF5 convergence test" << std::endl; +#endif + for (int idx = 0; idx < num_tests; idx++) { + mv_type y_vecs("history vectors", mySys.neqs, 5); + Kokkos::deep_copy(y0, 0.5); + + BDFSolve_wrapper + solve_wrapper(mySys, t_start, t_end, num_steps[idx], y0, y_new, rhs, + update, y_vecs, temp, jac); + Kokkos::parallel_for(myPolicy, solve_wrapper); + + auto y_new_h = Kokkos::create_mirror_view(y_new); + Kokkos::deep_copy(y_new_h, y_new); + + errors[idx] = Kokkos::abs(y_new_h(0) - 1 / (1 + Kokkos::exp(-t_end))) / + Kokkos::abs(1 / (1 + Kokkos::exp(-t_end))); + } + measured_order = + Kokkos::pow(errors[num_tests - 1] / errors[0], 1.0 / (num_tests - 1)); +#if defined(HAVE_KOKKOSKERNELS_DEBUG) + std::cout << "expected ratio: 32, actual ratio: " << measured_order + << ", order error=" << Kokkos::abs(measured_order - 32.0) / 32.0 + << std::endl; +#endif + +} // test_BDF_Logistic + +template +void test_BDF_LotkaVolterra() { + using vec_type = Kokkos::View; + using mv_type = Kokkos::View; + using mat_type = Kokkos::View; + + LotkaVolterra mySys(1.1, 0.4, 0.1, 0.4); + + const scalar_type t_start = 0.0, t_end = 100.0; + vec_type y0("initial conditions", mySys.neqs), y_new("solution", mySys.neqs); + vec_type rhs("rhs", mySys.neqs), update("update", mySys.neqs); + mat_type jac("jacobian", mySys.neqs, mySys.neqs), + temp("temp storage", mySys.neqs, mySys.neqs + 4); + + // Test BDF5 + mv_type y_vecs("history vectors", mySys.neqs, 5); + + Kokkos::deep_copy(y0, 10.0); + Kokkos::deep_copy(y_vecs, 10.0); + + Kokkos::RangePolicy myPolicy(0, 1); + BDFSolve_wrapper + solve_wrapper(mySys, t_start, t_end, 1000, y0, y_new, rhs, update, y_vecs, + temp, jac); + Kokkos::parallel_for(myPolicy, solve_wrapper); +} + +template +void test_BDF_StiffChemestry() { + using vec_type = Kokkos::View; + using mv_type = Kokkos::View; + using mat_type = Kokkos::View; + + StiffChemestry mySys{}; + + const scalar_type t_start = 0.0, t_end = 500.0; + vec_type y0("initial conditions", mySys.neqs), y_new("solution", mySys.neqs); + vec_type rhs("rhs", mySys.neqs), update("update", mySys.neqs); + mat_type jac("jacobian", mySys.neqs, mySys.neqs), + temp("temp storage", mySys.neqs, mySys.neqs + 4); + + // Test BDF5 + mv_type y_vecs("history vectors", mySys.neqs, 5); + + auto y0_h = Kokkos::create_mirror_view(y0); + y0_h(0) = 1.0; + y0_h(1) = 0.0; + y0_h(2) = 0.0; + Kokkos::deep_copy(y0, y0_h); + Kokkos::deep_copy(y_vecs, 0.0); + + Kokkos::RangePolicy myPolicy(0, 1); + BDFSolve_wrapper + solve_wrapper(mySys, t_start, t_end, 200000, y0, y_new, rhs, update, + y_vecs, temp, jac); + Kokkos::parallel_for(myPolicy, solve_wrapper); +} + +template +struct BDFSolve_parallel { + ode_type my_ode; + scalar_type tstart, tend; + int num_steps; + vec_type y_old, y_new, rhs, update; + mv_type y_vecs; + mat_type temp, jac; + + BDFSolve_parallel(const ode_type& my_ode_, const scalar_type tstart_, + const scalar_type tend_, const int num_steps_, + 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_) + : my_ode(my_ode_), + tstart(tstart_), + tend(tend_), + num_steps(num_steps_), + y_old(y_old_), + y_new(y_new_), + rhs(rhs_), + update(update_), + y_vecs(y_vecs_), + temp(temp_), + jac(jac_) {} + + KOKKOS_FUNCTION + void operator()(const int idx) const { + auto local_y_old = Kokkos::subview( + y_old, Kokkos::pair(my_ode.neqs * idx, my_ode.neqs * (idx + 1))); + auto local_y_new = Kokkos::subview( + y_new, Kokkos::pair(my_ode.neqs * idx, my_ode.neqs * (idx + 1))); + auto local_rhs = Kokkos::subview( + rhs, Kokkos::pair(my_ode.neqs * idx, my_ode.neqs * (idx + 1))); + auto local_update = Kokkos::subview( + update, Kokkos::pair(my_ode.neqs * idx, my_ode.neqs * (idx + 1))); + + auto local_y_vecs = Kokkos::subview( + y_vecs, Kokkos::pair(my_ode.neqs * idx, my_ode.neqs * (idx + 1)), + Kokkos::ALL()); + auto local_temp = Kokkos::subview( + temp, Kokkos::pair(my_ode.neqs * idx, my_ode.neqs * (idx + 1)), + Kokkos::ALL()); + auto local_jac = Kokkos::subview( + jac, Kokkos::pair(my_ode.neqs * idx, my_ode.neqs * (idx + 1)), + Kokkos::ALL()); + + KokkosODE::Experimental::BDF::Solve( + my_ode, tstart, tend, num_steps, local_y_old, local_y_new, local_rhs, + local_update, local_y_vecs, local_temp, local_jac); + } +}; + +template +void test_BDF_parallel() { + using vec_type = Kokkos::View; + using mv_type = Kokkos::View; + using mat_type = Kokkos::View; + + LotkaVolterra mySys(1.1, 0.4, 0.1, 0.4); + constexpr int num_solves = 1000; + + const scalar_type t_start = 0.0, t_end = 100.0; + vec_type y0("initial conditions", mySys.neqs * num_solves), + y_new("solution", mySys.neqs * num_solves); + vec_type rhs("rhs", mySys.neqs * num_solves), + update("update", mySys.neqs * num_solves); + mat_type jac("jacobian", mySys.neqs * num_solves, mySys.neqs), + temp("temp storage", mySys.neqs * num_solves, mySys.neqs + 4); + + // Test BDF5 + mv_type y_vecs("history vectors", mySys.neqs * num_solves, 5); + + Kokkos::deep_copy(y0, 10.0); + Kokkos::deep_copy(y_vecs, 10.0); + + Kokkos::RangePolicy myPolicy(0, num_solves); + BDFSolve_parallel + solve_wrapper(mySys, t_start, t_end, 1000, y0, y_new, rhs, update, y_vecs, + temp, jac); + Kokkos::parallel_for(myPolicy, solve_wrapper); +} + +} // namespace Test + +TEST_F(TestCategory, BDF_Logistic_serial) { + ::Test::test_BDF_Logistic(); +} +TEST_F(TestCategory, BDF_LotkaVolterra_serial) { + ::Test::test_BDF_LotkaVolterra(); +} +TEST_F(TestCategory, BDF_StiffChemestry_serial) { + ::Test::test_BDF_StiffChemestry(); +} +TEST_F(TestCategory, BDF_parallel_serial) { + ::Test::test_BDF_parallel(); +}