From 6bb07d7b3932c2905ed4482bb76b8c9df23154e0 Mon Sep 17 00:00:00 2001
From: Philipp Grete <gretephi@msu.edu>
Date: Mon, 11 Mar 2024 21:07:52 +0100
Subject: [PATCH] Use 2nd order time integrator for particles

---
 src/tracers/tracers.cpp | 33 ++++++++++++++-------
 src/utils/robust.hpp    | 66 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 88 insertions(+), 11 deletions(-)
 create mode 100644 src/utils/robust.hpp

diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp
index b9f518de..cca35afd 100644
--- a/src/tracers/tracers.cpp
+++ b/src/tracers/tracers.cpp
@@ -118,6 +118,9 @@ TaskStatus AdvectTracers(MeshBlockData<Real> *mbd, const Real dt) {
   auto &x = swarm->Get<Real>("x").Get();
   auto &y = swarm->Get<Real>("y").Get();
   auto &z = swarm->Get<Real>("z").Get();
+  auto &vel_x = swarm->Get<Real>("vel_x").Get();
+  auto &vel_y = swarm->Get<Real>("vel_y").Get();
+  auto &vel_z = swarm->Get<Real>("vel_z").Get();
 
   auto swarm_d = swarm->GetDeviceContext();
 
@@ -131,17 +134,25 @@ TaskStatus AdvectTracers(MeshBlockData<Real> *mbd, const Real dt) {
           int k, j, i;
           swarm_d.Xtoijk(x(n), y(n), z(n), i, j, k);
 
-          // Predictor/corrector will first make sense with non constant interpolation
-          // TODO(pgrete) add non-cnonst interpolation
-          // predictor
-          // const Real kx = x(n) + 0.5 * dt * rhs1;
-          // const Real ky = y(n) + 0.5 * dt * rhs2;
-          // const Real kz = z(n) + 0.5 * dt * rhs3;
-
-          // corrector
-          x(n) += prim_pack(IV1, k, j, i) * dt;
-          y(n) += prim_pack(IV2, k, j, i) * dt;
-          z(n) += prim_pack(IV3, k, j, i) * dt;
+          // RK2/Heun's method (as the default in Flash)
+          // https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p62/node130.html#SECTION06813000000000000000
+          // Intermediate position and velocities
+          // x^{*,n+1} = x^n + dt * v^n
+          const auto x_star = x(n) + dt * vel_x(n);
+          const auto y_star = y(n) + dt * vel_y(n);
+          const auto z_star = z(n) + dt * vel_z(n);
+
+          // v^{*,n+1} = v(x^{*,n+1}, t^{n+1})
+          // First parameter b=0 assume to operate on a pack of a single block and needs
+          // to be updated if this becomes a MeshData function
+          const auto vel_x_star = LCInterp::Do(0, x_star, y_star, z_star, prim_pack, IV1);
+          const auto vel_y_star = LCInterp::Do(0, x_star, y_star, z_star, prim_pack, IV2);
+          const auto vel_z_star = LCInterp::Do(0, x_star, y_star, z_star, prim_pack, IV3);
+
+          // Full update using mean velocity
+          x(n) += dt * 0.5 * (vel_x(n) + vel_x_star);
+          y(n) += dt * 0.5 * (vel_y(n) + vel_y_star);
+          z(n) += dt * 0.5 * (vel_z(n) + vel_z_star);
 
           // The following call is required as it updates the internal block id following
           // the advection. The internal id is used in the subsequent task to communicate
diff --git a/src/utils/robust.hpp b/src/utils/robust.hpp
new file mode 100644
index 00000000..7488abea
--- /dev/null
+++ b/src/utils/robust.hpp
@@ -0,0 +1,66 @@
+
+//========================================================================================
+// AthenaPK - a performance portable block structured AMR astrophysical MHD code.
+// Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved.
+// Licensed under the BSD 3-Clause License (the "LICENSE").
+//========================================================================================
+// Copied from https://github.com/lanl/phoebus
+//========================================================================================
+// © 2022. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
+// is operated by Triad National Security, LLC for the U.S.
+// Department of Energy/National Nuclear Security Administration. All
+// rights in the program are reserved by Triad National Security, LLC,
+// and the U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others
+// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works,
+// distribute copies to the public, perform publicly and display
+// publicly, and to permit others to do so.
+
+#ifndef UTILS_ROBUST_HPP_
+#define UTILS_ROBUST_HPP_
+
+#include "basic_types.hpp"
+#include <kokkos_abstraction.hpp>
+#include <limits>
+
+namespace robust {
+using parthenon::Real;
+
+template <typename T = Real>
+KOKKOS_FORCEINLINE_FUNCTION constexpr auto LARGE() {
+  return 0.1 * std::numeric_limits<T>::max();
+}
+template <typename T = Real>
+KOKKOS_FORCEINLINE_FUNCTION constexpr auto SMALL() {
+  return 10 * std::numeric_limits<T>::min();
+}
+template <typename T = Real>
+KOKKOS_FORCEINLINE_FUNCTION constexpr auto EPS() {
+  return 10 * std::numeric_limits<T>::epsilon();
+}
+
+template <typename T>
+KOKKOS_FORCEINLINE_FUNCTION auto make_positive(const T val) {
+  return std::max(val, EPS<T>());
+}
+
+KOKKOS_FORCEINLINE_FUNCTION
+Real make_bounded(const Real val, const Real vmin, const Real vmax) {
+  return std::min(std::max(val, vmin + EPS()), vmax * (1.0 - EPS()));
+}
+
+template <typename T>
+KOKKOS_INLINE_FUNCTION int sgn(const T &val) {
+  return (T(0) <= val) - (val < T(0));
+}
+template <typename A, typename B>
+KOKKOS_INLINE_FUNCTION auto ratio(const A &a, const B &b) {
+  const B sgn = b >= 0 ? 1 : -1;
+  return a / (b + sgn * SMALL<B>());
+}
+} // namespace robust
+
+#endif // UTILS_ROBUST_HPP_