From c6fed550913cc477cd964e8d137502da5555090c Mon Sep 17 00:00:00 2001
From: Michael B Kuhn <31661049+mbkuhn@users.noreply.github.com>
Date: Tue, 29 Aug 2023 10:57:42 -0600
Subject: [PATCH] nodal projection: allow compgrad to be called on supplied phi
 (#119)

---
 Projections/hydro_NodalProjector.H   |  2 ++
 Projections/hydro_NodalProjector.cpp | 33 ++++++++++++++++++++++++++++
 2 files changed, 35 insertions(+)

diff --git a/Projections/hydro_NodalProjector.H b/Projections/hydro_NodalProjector.H
index f7356fe4..6f89a933 100644
--- a/Projections/hydro_NodalProjector.H
+++ b/Projections/hydro_NodalProjector.H
@@ -81,6 +81,8 @@ public:
     amrex::Vector<       amrex::MultiFab* > getPhi          ()       {return GetVecOfPtrs(m_phi);}
     amrex::Vector< const amrex::MultiFab* > getPhiConst     () const {return GetVecOfConstPtrs(m_phi);}
 
+    amrex::Vector<       amrex::MultiFab* > calcGradPhi     (const amrex::Vector<amrex::MultiFab*>& a_phi);
+
     void computeRHS ( const amrex::Vector<amrex::MultiFab*>&       a_rhs,
                       const amrex::Vector<amrex::MultiFab*>&       a_vel,
                       const amrex::Vector<amrex::MultiFab*>&       a_S_cc = {},
diff --git a/Projections/hydro_NodalProjector.cpp b/Projections/hydro_NodalProjector.cpp
index 86619992..577a446b 100644
--- a/Projections/hydro_NodalProjector.cpp
+++ b/Projections/hydro_NodalProjector.cpp
@@ -360,6 +360,39 @@ NodalProjector::project ( const Vector<MultiFab*>& a_phi, Real a_rtol, Real a_at
     }
 }
 
+amrex::Vector<amrex::MultiFab*>
+NodalProjector::calcGradPhi ( const Vector<MultiFab*>& a_phi )
+{
+    BL_PROFILE("NodalProjector::calcGradPhi");
+    AMREX_ALWAYS_ASSERT(!m_need_bcs);
+
+    // Copy input
+    AMREX_ALWAYS_ASSERT(a_phi.size()==m_phi.size());
+    for (int lev=0; lev < m_phi.size(); ++lev )
+    {
+        MultiFab::Copy(m_phi[lev],*a_phi[lev],0,0,1,m_phi[lev].nGrow());
+    }
+
+    // Set coeffs involved with fluxes
+    for (int lev = 0; lev < m_sigma.size(); ++lev)
+    {
+        m_linop -> setSigma(lev, *m_sigma[lev]);
+    }
+
+    // Perform projection
+    for (int lev(0); lev < m_phi.size(); ++lev)
+    {
+        // set m_fluxes = grad(phi)
+        m_linop->compGrad(lev,m_fluxes[lev],m_phi[lev]);
+    }
+
+    // Average down
+    averageDown(GetVecOfPtrs(m_fluxes));
+
+    // Output fluxes
+    return GetVecOfPtrs(m_fluxes);
+}
+
 
 //
 // Compute RHS: div(u) + S_nd + S_cc