From 3aaa6a734518c98d1b2da4eafef69ddbdbaafba0 Mon Sep 17 00:00:00 2001
From: Lucas Timotheo Sanches <lucas.t.s.carneiro@gmail.com>
Date: Wed, 8 Jan 2025 23:16:11 -0600
Subject: [PATCH] ODESolvers: Added mechanism for detecting if previous RHSs
 need to be refilled after initialization and regridding

---
 ODESolvers/interface.ccl |  2 ++
 ODESolvers/schedule.ccl  | 16 ++++++++++++++++
 ODESolvers/src/solve.cxx | 33 +++++++++++++++++++++++++++++----
 3 files changed, 47 insertions(+), 4 deletions(-)

diff --git a/ODESolvers/interface.ccl b/ODESolvers/interface.ccl
index a2576e42d..c27450a2e 100644
--- a/ODESolvers/interface.ccl
+++ b/ODESolvers/interface.ccl
@@ -10,3 +10,5 @@ void FUNCTION CallScheduleGroup(
   CCTK_POINTER IN cctkGH,
   CCTK_STRING IN groupname)
 REQUIRES FUNCTION CallScheduleGroup
+
+CCTK_REAL refill_prev_rhss TYPE=scalar "If true and using hybrid methods, this variable causes the stored prev. RHS to be refilled via regular RK step"
\ No newline at end of file
diff --git a/ODESolvers/schedule.ccl b/ODESolvers/schedule.ccl
index ad0fdf648..24aba4d41 100644
--- a/ODESolvers/schedule.ccl
+++ b/ODESolvers/schedule.ccl
@@ -4,8 +4,24 @@ SCHEDULE ODESolvers_Solve AT evol
 {
   LANG: C
   OPTIONS: level
+  READS: refill_prev_rhss
+  WRITES: refill_prev_rhss
 } "Solve ODEs"
 
+SCHEDULE ODESolvers_SetRefillPrevRHSS AT initial
+{
+  LANG: C
+  OPTIONS: global
+  WRITES: refill_prev_rhss
+} "Sets refill_prev_rhss to true, thus having previous RHS information recomputed via standard RK methods"
+
+SCHEDULE ODESolvers_SetRefillPrevRHSS AT postregrid
+{
+  LANG: C
+  OPTIONS: global
+  WRITES: refill_prev_rhss
+} "Sets refill_prev_rhss to true, thus having previous RHS information recomputed via standard RK methods"
+
 
 
 # CarpetX scheduled groups:
diff --git a/ODESolvers/src/solve.cxx b/ODESolvers/src/solve.cxx
index c9286b7d5..ce505e892 100644
--- a/ODESolvers/src/solve.cxx
+++ b/ODESolvers/src/solve.cxx
@@ -727,7 +727,7 @@ void mark_invalid(const std::vector<int> &groups) {
 ///////////////////////////////////////////////////////////////////////////////
 
 extern "C" void ODESolvers_Solve(CCTK_ARGUMENTS) {
-  DECLARE_CCTK_ARGUMENTS_ODESolvers_Solve;
+  DECLARE_CCTK_ARGUMENTSX_ODESolvers_Solve;
   DECLARE_CCTK_PARAMETERS;
 
   static bool did_output = false;
@@ -1018,7 +1018,13 @@ extern "C" void ODESolvers_Solve(CCTK_ARGUMENTS) {
   } else if (CCTK_EQUALS(method, "BMS422")) {
 
     // RK4 Bootstrapping
-    if (cctkGH->cctk_iteration <= 1) {
+    if (static_cast<CCTK_INT>(*refill_prev_rhss) != 0) {
+      if (verbose) {
+        CCTK_VINFO("  Taking RK4 step to fill prev. RHS");
+      }
+
+      *refill_prev_rhss -= -1;
+
       // k1 = f(y0)
       // k2 = f(y0 + h/2 k1)
       // k3 = f(y0 + h/2 k2)
@@ -1052,8 +1058,8 @@ extern "C" void ODESolvers_Solve(CCTK_ARGUMENTS) {
       calcrhs(4);
       calcupdate(4, dt, 0.0, reals<3>{1.0, dt / 6, dt / 6},
                  states<3>{&old, &kaccum, &rhs});
-
     } else {
+
       // k0 = f(y(t - h))
       // k1 = f(y(t))
       // k2 = f(y(t) + h * (a20 * k0 + a21 * k1))
@@ -1104,7 +1110,13 @@ extern "C" void ODESolvers_Solve(CCTK_ARGUMENTS) {
   } else if (CCTK_EQUALS(method, "BMS431")) {
 
     // RK4 Bootstrapping
-    if (cctkGH->cctk_iteration <= 2) {
+    if (static_cast<CCTK_INT>(*refill_prev_rhss) != 0) {
+      if (verbose) {
+        CCTK_VINFO("  Taking RK4 step to fill prev. RHS");
+      }
+
+      *refill_prev_rhss -= -1;
+
       // k1 = f(y0)
       // k2 = f(y0 + h/2 k1)
       // k3 = f(y0 + h/2 k2)
@@ -1484,4 +1496,17 @@ extern "C" void ODESolvers_Solve(CCTK_ARGUMENTS) {
   // TODO: Update time here, and not during time level cycling in the driver
 }
 
+extern "C" void ODESolvers_SetRefillPrevRHSS(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_ODESolvers_SetRefillPrevRHSS;
+  DECLARE_CCTK_PARAMETERS;
+
+  if (CCTK_EQUALS(method, "BMS422")) {
+    *refill_prev_rhss = 1;
+  } else if (CCTK_EQUALS(method, "BMS431")) {
+    *refill_prev_rhss = 2;
+  } else {
+    *refill_prev_rhss = 0;
+  }
+}
+
 } // namespace ODESolvers