From 19279b101bae1830e742f90489f2058903825269 Mon Sep 17 00:00:00 2001
From: dmontgomeryNREL <David.Montgomery@nrel.gov>
Date: Mon, 21 Oct 2024 16:19:12 -0600
Subject: [PATCH] Now passes in MultiFabs for state_old and state_new. This
 provides access to state_old.Factory() etc.

---
 .../EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp  |  7 ++++---
 Source/PeleLMeX_Forces.cpp                           |  8 ++++----
 Source/PeleLMeX_ProblemSpecificFunctions.H           |  6 +++---
 Source/PeleLMeX_ProblemSpecificFunctions.cpp         | 12 +++++++-----
 4 files changed, 18 insertions(+), 15 deletions(-)

diff --git a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp
index 5c06c6cd..d7669544 100644
--- a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp
+++ b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp
@@ -28,11 +28,11 @@ problem_modify_ext_sources(
   Real /*time*/,
   Real /*dt*/,
   int lev,
-  MultiArray4<const Real> const& state_old_arr,
-  MultiArray4<const Real> const& /*state_new_arr*/,
+  const MultiFab& state_old,
+  const MultiFab& /*state_new*/,
   Vector<std::unique_ptr<MultiFab>>& a_extSource,
   const GeometryData& /*geomdata*/,
-  ProbParm const& prob_parm)
+  const ProbParm& prob_parm)
 {
   /*
   Notes:
@@ -42,6 +42,7 @@ problem_modify_ext_sources(
   */
 
   auto ext_source_arr = a_extSource[lev]->arrays();
+  auto const& state_old_arr = state_old.const_arrays();
 
   ParallelFor(
     *a_extSource[lev],
diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp
index 0db93119..164fac2f 100644
--- a/Source/PeleLMeX_Forces.cpp
+++ b/Source/PeleLMeX_Forces.cpp
@@ -271,11 +271,11 @@ PeleLM::getExternalSources(
   // User defined external sources
   if (m_user_defined_ext_sources) {
     for (int lev = 0; lev <= finest_level; lev++) {
+      auto* ldata_p_old = getLevelDataPtr(lev, a_timestamp_old);
+      auto* ldata_p_new = getLevelDataPtr(lev, a_timestamp_new);
       problem_modify_ext_sources(
-        getTime(lev, a_timestamp_new), m_dt, lev,
-        getLevelDataPtr(lev, a_timestamp_old)->state.const_arrays(),
-        getLevelDataPtr(lev, a_timestamp_new)->state.const_arrays(),
-        m_extSource, geom[lev].data(), *prob_parm_d);
+        getTime(lev, a_timestamp_new), m_dt, lev, ldata_p_old->state,
+        ldata_p_new->state, m_extSource, geom[lev].data(), *prob_parm_d);
     }
   }
 }
\ No newline at end of file
diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.H b/Source/PeleLMeX_ProblemSpecificFunctions.H
index 697064b3..534c1934 100644
--- a/Source/PeleLMeX_ProblemSpecificFunctions.H
+++ b/Source/PeleLMeX_ProblemSpecificFunctions.H
@@ -11,9 +11,9 @@ void problem_modify_ext_sources(
   Real time,
   Real dt,
   int lev,
-  MultiArray4<const Real> const& state_old,
-  MultiArray4<const Real> const& state_new,
+  const MultiFab& state_old,
+  const MultiFab& state_new,
   Vector<std::unique_ptr<MultiFab>>& a_extSource,
   const GeometryData& geomdata,
-  ProbParm const& prob_parm);
+  const ProbParm& prob_parm);
 #endif
\ No newline at end of file
diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.cpp b/Source/PeleLMeX_ProblemSpecificFunctions.cpp
index ff2d4f78..e228a088 100644
--- a/Source/PeleLMeX_ProblemSpecificFunctions.cpp
+++ b/Source/PeleLMeX_ProblemSpecificFunctions.cpp
@@ -28,11 +28,11 @@ problem_modify_ext_sources(
   Real /*time*/,
   Real /*dt*/,
   int /*lev*/,
-  MultiArray4<const Real> const& /*state_old_arr*/,
-  MultiArray4<const Real> const& /*state_new_arr*/,
+  const MultiFab& /*state_old*/,
+  const MultiFab& /*state_new*/,
   Vector<std::unique_ptr<MultiFab>>& /*a_extSource*/,
   const GeometryData& /*geomdata*/,
-  ProbParm const& /*prob_parm*/)
+  const ProbParm& /*prob_parm*/)
 {
   /*
   Notes:
@@ -42,13 +42,15 @@ problem_modify_ext_sources(
 
   // Example: Exponential decay ode quantity
   auto ext_source_arr = a_extSource[lev]->arrays();
+  auto const& state_old_arr = state_old.const_arrays();
+
   ParallelFor(
     *a_extSource[lev],
     [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept {
-      for (int n = 0; n < NUM_ODE; n++){
+      for (int n = 0; n < NUM_ODE; n++) {
         Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n);
         Real src_strength = -1.0 * pow(10.0,n+1);
-        ext_source_arr[box_no](i, j, k, FIRSTODE + n) += src_strength * B_n;
+        ext_source_arr[box_no](i, j, k, FIRSTODE + n) += src;
       }
     });
   Gpu::streamSynchronize();