From 48147044ca4dfd18f77472674a060ba13f5de962 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Tue, 12 Nov 2024 12:48:00 -0800
Subject: [PATCH 01/20] WIP GMRES example

---
 .../GMRES/Poisson/AMReX_GMRES_Poisson.H       | 175 ++++++++++++++
 .../LinearSolvers/GMRES/Poisson/GNUmakefile   |  17 ++
 .../LinearSolvers/GMRES/Poisson/Make.package  |   2 +
 .../LinearSolvers/GMRES/Poisson/inputs        |   7 +
 .../LinearSolvers/GMRES/Poisson/main.cpp      | 226 ++++++++++++++++++
 5 files changed, 427 insertions(+)
 create mode 100644 ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
 create mode 100644 ExampleCodes/LinearSolvers/GMRES/Poisson/GNUmakefile
 create mode 100644 ExampleCodes/LinearSolvers/GMRES/Poisson/Make.package
 create mode 100644 ExampleCodes/LinearSolvers/GMRES/Poisson/inputs
 create mode 100644 ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
new file mode 100644
index 00000000..f9ea8002
--- /dev/null
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -0,0 +1,175 @@
+#ifndef AMREX_GMRES_POISSON_H_
+#define AMREX_GMRES_POISSON_H_
+#include <AMReX_Config.H>
+
+#include <AMReX_GMRES.H>
+
+#include <utility>
+
+namespace amrex {
+
+/**
+ * \brief Solve Poisson equation using null preconditioner
+ */
+template <typename MF>
+class GMRESPoissonT
+{
+public:
+    using GM = GMRES<MF,GMRESPoissonT<MF>>;
+
+    explicit GMRESPoissonT ();
+
+    /**
+     * \brief Solve the linear system
+     *
+     * \param a_sol     unknowns, i.e., x in A x = b.
+     * \param a_rhs     RHS, i.e., b in A x = b.
+     * \param a_tol_rel relative tolerance.
+     * \param a_tol_abs absolute tolerance.
+     */
+    void solve (MF& a_sol, MF const& a_rhs, amrex::Real a_tol_rel, amrex::Real a_tol_abs);
+
+    //! Sets verbosity.
+    void setVerbose (int v) { m_gmres.setVerbose(v); }
+
+    //! Sets the max number of iterations
+    void setMaxIters (int niters) { m_gmres.setMaxIters(niters); }
+
+    //! Gets the number of iterations.
+    [[nodiscard]] int getNumIters () const { return m_gmres.getNumIters(); }
+
+    //! Gets the 2-norm of the residual.
+    [[nodiscard]] amrex::Real getResidualNorm () const { return m_gmres.getResidualNorm(); }
+
+    //! Get the GMRES object.
+    GM& getGMRES () { return m_gmres; }
+
+    //! Make MultiFab without ghost cells
+    MF makeVecRHS () const;
+
+    //! Make MultiFab with ghost cells and set ghost cells to zero
+    MF makeVecLHS () const;
+
+    amrex::Real norm2 (MF const& mf) const;
+
+    static void scale (MF& mf, amrex::Real scale_factor);
+
+    amrex::Real dotProduct (MF const& mf1, MF const& mf2) const;
+
+    //! lhs = 0
+    static void setToZero (MF& lhs);
+
+    //! lhs = rhs
+    static void assign (MF& lhs, MF const& rhs);
+
+    //! lhs += a*rhs
+    static void increment (MF& lhs, MF const& rhs, amrex::Real a);
+
+    //! lhs = a*rhs_a + b*rhs_b
+    static void linComb (MF& lhs, amrex::Real a, MF const& rhs_a, amrex::Real b, MF const& rhs_b);
+
+    //! lhs = L(rhs)
+    void apply (MF& lhs, MF const& rhs) const;
+
+    void precond (MF& lhs, MF const& rhs) const;
+
+private:
+    GM m_gmres;
+};
+
+template <typename MF>
+GMRESPoissonT<MF>::GMRESPoissonT ()
+{
+
+}
+
+template <typename MF>
+auto GMRESPoissonT<MF>::makeVecRHS () const -> MF
+{
+//    return m_linop->make(0, 0, IntVect(0));
+    MultiFab x;
+    return x;
+}
+
+template <typename MF>
+auto GMRESPoissonT<MF>::makeVecLHS () const -> MF
+{
+//    auto mf = m_linop->make(0, 0, IntVect(1));
+//    setBndry(mf, amrex::Real(0), 0, nComp(mf));
+//    return mf;
+    MultiFab x;
+    return x;
+}
+
+template <typename MF>
+auto GMRESPoissonT<MF>::norm2 (MF const& mf) const -> amrex::Real
+{
+//    auto r = m_linop->xdoty(0, 0, mf, mf, false);
+//    return std::sqrt(r);
+    return 0.;
+}
+
+template <typename MF>
+void GMRESPoissonT<MF>::scale (MF& mf, amrex::Real scale_factor)
+{
+    Scale(mf, scale_factor, 0, nComp(mf), 0);
+}
+
+template <typename MF>
+auto GMRESPoissonT<MF>::dotProduct (MF const& mf1, MF const& mf2) const -> amrex::Real
+{
+//    return m_linop->xdoty(0, 0, mf1, mf2, false);
+    return 0.;
+}
+
+template <typename MF>
+void GMRESPoissonT<MF>::setToZero (MF& lhs)
+{
+    setVal(lhs, amrex::Real(0.0));
+}
+
+template <typename MF>
+void GMRESPoissonT<MF>::assign (MF& lhs, MF const& rhs)
+{
+    LocalCopy(lhs, rhs, 0, 0, nComp(lhs), IntVect(0));
+}
+
+template <typename MF>
+void GMRESPoissonT<MF>::increment (MF& lhs, MF const& rhs, amrex::Real a)
+{
+    Saxpy(lhs, a, rhs, 0, 0, nComp(lhs), IntVect(0));
+}
+
+template <typename MF>
+void GMRESPoissonT<MF>::linComb (MF& lhs, amrex::Real a, MF const& rhs_a, amrex::Real b, MF const& rhs_b)
+{
+    LinComb(lhs, a, rhs_a, 0, b, rhs_b, 0, 0, nComp(lhs), IntVect(0));
+}
+
+template <typename MF>
+void GMRESPoissonT<MF>::apply (MF& lhs, MF const& rhs) const
+{
+    // fixme
+}
+
+template <typename MF>
+void GMRESPoissonT<MF>::precond (MF& lhs, MF const& rhs) const
+{
+    LocalCopy(lhs, rhs, 0, 0, nComp(lhs), IntVect(0));
+}
+
+template <typename MF>
+void GMRESPoissonT<MF>::solve (MF& a_sol, MF const& a_rhs, amrex::Real a_tol_rel, amrex::Real a_tol_abs)
+{
+    /*
+    auto res = makeVecRHS();
+    auto cor = makeVecLHS();
+    m_gmres.solve(cor, res, a_tol_rel, a_tol_abs); // L(cor) = res
+    */
+}
+
+using GMRESPoisson = GMRESPoissonT<MultiFab>;
+
+}
+
+#endif
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/GNUmakefile b/ExampleCodes/LinearSolvers/GMRES/Poisson/GNUmakefile
new file mode 100644
index 00000000..c216d09a
--- /dev/null
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/GNUmakefile
@@ -0,0 +1,17 @@
+# AMREX_HOME defines the directory in which we will find all the AMReX code.
+AMREX_HOME ?= ../../../../../amrex
+
+DEBUG        = FALSE
+USE_MPI      = FALSE
+USE_OMP      = FALSE
+COMP         = gnu
+DIM          = 3
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.defs
+
+include ./Make.package
+
+include $(AMREX_HOME)/Src/Base/Make.package
+include $(AMREX_HOME)/Src/LinearSolvers/Make.package
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.rules
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/Make.package b/ExampleCodes/LinearSolvers/GMRES/Poisson/Make.package
new file mode 100644
index 00000000..d533b57b
--- /dev/null
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/Make.package
@@ -0,0 +1,2 @@
+CEXE_sources += main.cpp
+CEXE_headers += AMReX_GMRES_Poisson.H
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs b/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs
new file mode 100644
index 00000000..5b098f82
--- /dev/null
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs
@@ -0,0 +1,7 @@
+n_cell = 32
+max_grid_size = 16
+
+nsteps = 1000
+plot_int = 100
+
+dt = 1.e-5
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
new file mode 100644
index 00000000..bf3173b5
--- /dev/null
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
@@ -0,0 +1,226 @@
+/*
+ * A simplified single file version of the HeatEquation_EX0_C exmaple.
+ * This code is designed to be used with Demo_Tutorial.rst.
+ *
+ */
+
+
+#include <AMReX.H>
+#include <AMReX_PlotFileUtil.H>
+#include <AMReX_ParmParse.H>
+#include <AMReX_GMRES_Poisson.H>
+
+int main (int argc, char* argv[])
+{
+    amrex::Initialize(argc,argv);
+    {
+
+    // **********************************
+    // DECLARE SIMULATION PARAMETERS
+    // **********************************
+
+    // number of cells on each side of the domain
+    int n_cell;
+
+    // size of each box (or grid)
+    int max_grid_size;
+
+    // total steps in simulation
+    int nsteps;
+
+    // how often to write a plotfile
+    int plot_int;
+
+    // time step
+    amrex::Real dt;
+
+    // **********************************
+    // READ PARAMETER VALUES FROM INPUT DATA
+    // **********************************
+    // inputs parameters
+    {
+        // ParmParse is way of reading inputs from the inputs file
+        // pp.get means we require the inputs file to have it
+        // pp.query means we optionally need the inputs file to have it - but we must supply a default here
+        amrex::ParmParse pp;
+
+        // We need to get n_cell from the inputs file - this is the number of cells on each side of
+        //   a square (or cubic) domain.
+        pp.get("n_cell",n_cell);
+
+        // The domain is broken into boxes of size max_grid_size
+        pp.get("max_grid_size",max_grid_size);
+
+        // Default nsteps to 10, allow us to set it to something else in the inputs file
+        nsteps = 10;
+        pp.query("nsteps",nsteps);
+
+        // Default plot_int to -1, allow us to set it to something else in the inputs file
+        //  If plot_int < 0 then no plot files will be written
+        plot_int = -1;
+        pp.query("plot_int",plot_int);
+
+        // time step
+        pp.get("dt",dt);
+    }
+
+    // **********************************
+    // DEFINE SIMULATION SETUP AND GEOMETRY
+    // **********************************
+
+    // make BoxArray and Geometry
+    // ba will contain a list of boxes that cover the domain
+    // geom contains information such as the physical domain size,
+    // number of points in the domain, and periodicity
+    amrex::BoxArray ba;
+    amrex::Geometry geom;
+
+    // define lower and upper indices
+    amrex::IntVect dom_lo(0,0,0);
+    amrex::IntVect dom_hi(n_cell-1, n_cell-1, n_cell-1);
+
+    // Make a single box that is the entire domain
+    amrex::Box domain(dom_lo, dom_hi);
+
+    // Initialize the boxarray "ba" from the single box "domain"
+    ba.define(domain);
+
+    // Break up boxarray "ba" into chunks no larger than "max_grid_size" along a direction
+    ba.maxSize(max_grid_size);
+
+    // This defines the physical box, [0,1] in each direction.
+    amrex::RealBox real_box({ 0., 0., 0.},
+                     { 1., 1., 1.});
+
+    // periodic in all direction
+    amrex::Array<int,3> is_periodic{1,1,1};
+
+    // This defines a Geometry object
+    geom.define(domain, real_box, amrex::CoordSys::cartesian, is_periodic);
+
+    // extract dx from the geometry object
+    amrex::GpuArray<amrex::Real,3> dx = geom.CellSizeArray();
+
+    // Nghost = number of ghost cells for each array
+    int Nghost = 1;
+
+    // Ncomp = number of components for each array
+    int Ncomp = 1;
+
+    // How Boxes are distrubuted among MPI processes
+    amrex::DistributionMapping dm(ba);
+
+    // we allocate two phi multifabs; one will store the old state, the other the new.
+    amrex::MultiFab phi_old(ba, dm, Ncomp, Nghost);
+    amrex::MultiFab phi_new(ba, dm, Ncomp, Nghost);
+
+    // time = starting time in the simulation
+    amrex::Real time = 0.0;
+
+    // **********************************
+    // INITIALIZE DATA LOOP
+    // **********************************
+
+    // loop over boxes
+    for (amrex::MFIter mfi(phi_old); mfi.isValid(); ++mfi)
+    {
+        const amrex::Box& bx = mfi.validbox();
+
+        const amrex::Array4<amrex::Real>& phiOld = phi_old.array(mfi);
+
+        // set phi = 1 + e^(-(r-0.5)^2)
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
+        {
+
+            // **********************************
+            // SET VALUES FOR EACH CELL
+            // **********************************
+
+            amrex::Real x = (i+0.5) * dx[0];
+            amrex::Real y = (j+0.5) * dx[1];
+            amrex::Real z = (k+0.5) * dx[2];
+            amrex::Real rsquared = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5))/0.01;
+            phiOld(i,j,k) = 1. + std::exp(-rsquared);
+        });
+    }
+
+    // **********************************
+    // WRITE INITIAL PLOT FILE
+    // **********************************
+
+    // Write a plotfile of the initial data if plot_int > 0
+    if (plot_int > 0)
+    {
+        int step = 0;
+        const std::string& pltfile = amrex::Concatenate("plt",step,5);
+        WriteSingleLevelPlotfile(pltfile, phi_old, {"phi"}, geom, time, 0);
+    }
+
+
+    // **********************************
+    // MAIN TIME EVOLUTION LOOP
+    // **********************************
+
+    for (int step = 1; step <= nsteps; ++step)
+    {
+        // fill periodic ghost cells
+        phi_old.FillBoundary(geom.periodicity());
+
+        // new_phi = old_phi + dt * Laplacian(old_phi)
+        // loop over boxes
+        for ( amrex::MFIter mfi(phi_old); mfi.isValid(); ++mfi )
+        {
+            const amrex::Box& bx = mfi.validbox();
+
+            const amrex::Array4<amrex::Real>& phiOld = phi_old.array(mfi);
+            const amrex::Array4<amrex::Real>& phiNew = phi_new.array(mfi);
+
+            // advance the data by dt
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+            {
+
+                // **********************************
+                // EVOLVE VALUES FOR EACH CELL
+                // **********************************
+
+                phiNew(i,j,k) = phiOld(i,j,k) + dt *
+                    ( (phiOld(i+1,j,k) - 2.*phiOld(i,j,k) + phiOld(i-1,j,k)) / (dx[0]*dx[0])
+                     +(phiOld(i,j+1,k) - 2.*phiOld(i,j,k) + phiOld(i,j-1,k)) / (dx[1]*dx[1])
+                     +(phiOld(i,j,k+1) - 2.*phiOld(i,j,k) + phiOld(i,j,k-1)) / (dx[2]*dx[2])
+                        );
+            });
+        }
+
+        // **********************************
+        // INCREMENT
+        // **********************************
+
+        // update time
+        time = time + dt;
+
+        // copy new solution into old solution
+        amrex::MultiFab::Copy(phi_old, phi_new, 0, 0, 1, 0);
+
+        // Tell the I/O Processor to write out which step we're doing
+        amrex::Print() << "Advanced step " << step << "\n";
+
+
+        // **********************************
+        // WRITE PLOTFILE AT GIVEN INTERVAL
+        // **********************************
+
+        // Write a plotfile of the current data (plot_int was defined in the inputs file)
+        if (plot_int > 0 && step%plot_int == 0)
+        {
+            const std::string& pltfile = amrex::Concatenate("plt",step,5);
+            WriteSingleLevelPlotfile(pltfile, phi_new, {"phi"}, geom, time, step);
+        }
+    }
+
+
+    }
+    amrex::Finalize();
+    return 0;
+}
+
+

From d44c77aee22e3b7aa6ee142cdf5cebb3ce9bd455 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Fri, 15 Nov 2024 12:57:57 -0800
Subject: [PATCH 02/20] progress

---
 .../GMRES/Poisson/AMReX_GMRES_Poisson.H       | 115 ++++++++++++------
 .../LinearSolvers/GMRES/Poisson/main.cpp      | 106 ++--------------
 2 files changed, 85 insertions(+), 136 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index f9ea8002..974ad14d 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -1,23 +1,27 @@
-#ifndef AMREX_GMRES_POISSON_H_
-#define AMREX_GMRES_POISSON_H_
+#ifndef AMREX_GMRES_MLMG_H_
+#define AMREX_GMRES_MLMG_H_
 #include <AMReX_Config.H>
 
 #include <AMReX_GMRES.H>
-
 #include <utility>
 
 namespace amrex {
 
 /**
- * \brief Solve Poisson equation using null preconditioner
+ * \brief Solve using GMRES with multigrid as preconditioner
+ *
+ * The linear system to solve is provided by MLMG, which is also being used
+ * as the preconditioner.
+ *
  */
 template <typename MF>
-class GMRESPoissonT
+class GMRESPOISSONT
 {
 public:
-    using GM = GMRES<MF,GMRESPoissonT<MF>>;
+    using RT = amrex::Real; // typename MF::RT; // double or float
+    using GM = GMRES<MF,GMRESPOISSONT<MF>>;
 
-    explicit GMRESPoissonT ();
+    explicit GMRESPOISSONT (const BoxArray& ba, const DistributionMapping& dm);
 
     /**
      * \brief Solve the linear system
@@ -27,7 +31,7 @@ public:
      * \param a_tol_rel relative tolerance.
      * \param a_tol_abs absolute tolerance.
      */
-    void solve (MF& a_sol, MF const& a_rhs, amrex::Real a_tol_rel, amrex::Real a_tol_abs);
+    void solve (MF& a_sol, MF const& a_rhs, RT a_tol_rel, RT a_tol_abs);
 
     //! Sets verbosity.
     void setVerbose (int v) { m_gmres.setVerbose(v); }
@@ -39,22 +43,34 @@ public:
     [[nodiscard]] int getNumIters () const { return m_gmres.getNumIters(); }
 
     //! Gets the 2-norm of the residual.
-    [[nodiscard]] amrex::Real getResidualNorm () const { return m_gmres.getResidualNorm(); }
+    [[nodiscard]] RT getResidualNorm () const { return m_gmres.getResidualNorm(); }
 
     //! Get the GMRES object.
     GM& getGMRES () { return m_gmres; }
 
+    /**
+     * \brief Set MLMG's multiplicative property of zero
+     *
+     * This should NOT be called unless MLMG has the multiplicative property
+     * of zero. MLMG is not a matrix, and usually does not have the
+     * properties of a matrix such as the multiplicative property of zero
+     * (i.e., M*0=0) because how domain boundary conditions are
+     * handled. However, if MLMG has the property of zero, calling this
+     * function with true can have a small performance benefit.
+     */
+    void setPropertyOfZero (bool b) { m_prop_zero = b; }
+
     //! Make MultiFab without ghost cells
     MF makeVecRHS () const;
 
     //! Make MultiFab with ghost cells and set ghost cells to zero
     MF makeVecLHS () const;
 
-    amrex::Real norm2 (MF const& mf) const;
+    RT norm2 (MF const& mf) const;
 
-    static void scale (MF& mf, amrex::Real scale_factor);
+    static void scale (MF& mf, RT scale_factor);
 
-    amrex::Real dotProduct (MF const& mf1, MF const& mf2) const;
+    RT dotProduct (MF const& mf1, MF const& mf2) const;
 
     //! lhs = 0
     static void setToZero (MF& lhs);
@@ -63,46 +79,54 @@ public:
     static void assign (MF& lhs, MF const& rhs);
 
     //! lhs += a*rhs
-    static void increment (MF& lhs, MF const& rhs, amrex::Real a);
+    static void increment (MF& lhs, MF const& rhs, RT a);
 
     //! lhs = a*rhs_a + b*rhs_b
-    static void linComb (MF& lhs, amrex::Real a, MF const& rhs_a, amrex::Real b, MF const& rhs_b);
+    static void linComb (MF& lhs, RT a, MF const& rhs_a, RT b, MF const& rhs_b);
 
     //! lhs = L(rhs)
     void apply (MF& lhs, MF const& rhs) const;
 
     void precond (MF& lhs, MF const& rhs) const;
 
+    //! Control whether or not to use MLMG as preconditioner.
+    bool usePrecond (bool new_flag) { return std::exchange(m_use_precond, new_flag); }
+
+    //! Set the number of MLMG preconditioner iterations per GMRES iteration.
+    void setPrecondNumIters (int precond_niters) { m_precond_niters = precond_niters; }
+
 private:
     GM m_gmres;
+    BoxArray m_ba;
+    DistributionMapping m_dm;
+    bool m_use_precond = true;
+    bool m_prop_zero = false;
+    int m_precond_niters = 1;
 };
 
 template <typename MF>
-GMRESPoissonT<MF>::GMRESPoissonT ()
+GMRESPOISSONT<MF>::GMRESPOISSONT (const BoxArray& ba, const DistributionMapping& dm)
+    : m_ba(ba), m_dm(dm)
 {
-
+    m_gmres.define(*this);
 }
 
 template <typename MF>
-auto GMRESPoissonT<MF>::makeVecRHS () const -> MF
+auto GMRESPOISSONT<MF>::makeVecRHS () const -> MF
 {
-//    return m_linop->make(0, 0, IntVect(0));
-    MultiFab x;
-    return x;
+    return MultiFab(m_ba, m_dm, 1, 0);
 }
 
 template <typename MF>
-auto GMRESPoissonT<MF>::makeVecLHS () const -> MF
+auto GMRESPOISSONT<MF>::makeVecLHS () const -> MF
 {
-//    auto mf = m_linop->make(0, 0, IntVect(1));
-//    setBndry(mf, amrex::Real(0), 0, nComp(mf));
-//    return mf;
-    MultiFab x;
-    return x;
+    MF mf(m_ba, m_dm, 1, 0);
+    mf.setVal(0.);
+    return mf;
 }
 
 template <typename MF>
-auto GMRESPoissonT<MF>::norm2 (MF const& mf) const -> amrex::Real
+auto GMRESPOISSONT<MF>::norm2 (MF const& mf) const -> RT
 {
 //    auto r = m_linop->xdoty(0, 0, mf, mf, false);
 //    return std::sqrt(r);
@@ -110,65 +134,76 @@ auto GMRESPoissonT<MF>::norm2 (MF const& mf) const -> amrex::Real
 }
 
 template <typename MF>
-void GMRESPoissonT<MF>::scale (MF& mf, amrex::Real scale_factor)
+void GMRESPOISSONT<MF>::scale (MF& mf, RT scale_factor)
 {
     Scale(mf, scale_factor, 0, nComp(mf), 0);
 }
 
 template <typename MF>
-auto GMRESPoissonT<MF>::dotProduct (MF const& mf1, MF const& mf2) const -> amrex::Real
+auto GMRESPOISSONT<MF>::dotProduct (MF const& mf1, MF const& mf2) const -> RT
 {
 //    return m_linop->xdoty(0, 0, mf1, mf2, false);
     return 0.;
 }
 
 template <typename MF>
-void GMRESPoissonT<MF>::setToZero (MF& lhs)
+void GMRESPOISSONT<MF>::setToZero (MF& lhs)
 {
-    setVal(lhs, amrex::Real(0.0));
+    setVal(lhs, RT(0.0));
 }
 
 template <typename MF>
-void GMRESPoissonT<MF>::assign (MF& lhs, MF const& rhs)
+void GMRESPOISSONT<MF>::assign (MF& lhs, MF const& rhs)
 {
     LocalCopy(lhs, rhs, 0, 0, nComp(lhs), IntVect(0));
 }
 
 template <typename MF>
-void GMRESPoissonT<MF>::increment (MF& lhs, MF const& rhs, amrex::Real a)
+void GMRESPOISSONT<MF>::increment (MF& lhs, MF const& rhs, RT a)
 {
     Saxpy(lhs, a, rhs, 0, 0, nComp(lhs), IntVect(0));
 }
 
 template <typename MF>
-void GMRESPoissonT<MF>::linComb (MF& lhs, amrex::Real a, MF const& rhs_a, amrex::Real b, MF const& rhs_b)
+void GMRESPOISSONT<MF>::linComb (MF& lhs, RT a, MF const& rhs_a, RT b, MF const& rhs_b)
 {
     LinComb(lhs, a, rhs_a, 0, b, rhs_b, 0, 0, nComp(lhs), IntVect(0));
 }
 
 template <typename MF>
-void GMRESPoissonT<MF>::apply (MF& lhs, MF const& rhs) const
+void GMRESPOISSONT<MF>::apply (MF& lhs, MF const& rhs) const
 {
-    // fixme
+//    m_linop->apply(0, 0, lhs, const_cast<MF&>(rhs),
+//                   MLLinOpT<MF>::BCMode::Homogeneous,
+//                   MLLinOpT<MF>::StateMode::Correction);
 }
 
 template <typename MF>
-void GMRESPoissonT<MF>::precond (MF& lhs, MF const& rhs) const
+void GMRESPOISSONT<MF>::precond (MF& lhs, MF const& rhs) const
 {
-    LocalCopy(lhs, rhs, 0, 0, nComp(lhs), IntVect(0));
+    if (m_use_precond) {
+
+
+    } else {
+        LocalCopy(lhs, rhs, 0, 0, nComp(lhs), IntVect(0));
+    }
 }
 
 template <typename MF>
-void GMRESPoissonT<MF>::solve (MF& a_sol, MF const& a_rhs, amrex::Real a_tol_rel, amrex::Real a_tol_abs)
+void GMRESPOISSONT<MF>::solve (MF& a_sol, MF const& a_rhs, RT a_tol_rel, RT a_tol_abs)
 {
     /*
     auto res = makeVecRHS();
+    m_mlmg->apply({&res}, {&a_sol}); // res = L(sol)
+    increment(res, a_rhs, RT(-1)); // res = L(sol) - rhs
     auto cor = makeVecLHS();
+    m_linop->setDirichletNodesToZero(0,0,res);
     m_gmres.solve(cor, res, a_tol_rel, a_tol_abs); // L(cor) = res
+    increment(a_sol, cor, RT(-1)); // sol = sol - cor
     */
 }
 
-using GMRESPoisson = GMRESPoissonT<MultiFab>;
+using GMRESPOISSON = GMRESPOISSONT<MultiFab>;
 
 }
 
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
index bf3173b5..8e90595a 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
@@ -50,18 +50,6 @@ int main (int argc, char* argv[])
 
         // The domain is broken into boxes of size max_grid_size
         pp.get("max_grid_size",max_grid_size);
-
-        // Default nsteps to 10, allow us to set it to something else in the inputs file
-        nsteps = 10;
-        pp.query("nsteps",nsteps);
-
-        // Default plot_int to -1, allow us to set it to something else in the inputs file
-        //  If plot_int < 0 then no plot files will be written
-        plot_int = -1;
-        pp.query("plot_int",plot_int);
-
-        // time step
-        pp.get("dt",dt);
     }
 
     // **********************************
@@ -101,46 +89,35 @@ int main (int argc, char* argv[])
     // extract dx from the geometry object
     amrex::GpuArray<amrex::Real,3> dx = geom.CellSizeArray();
 
-    // Nghost = number of ghost cells for each array
-    int Nghost = 1;
-
-    // Ncomp = number of components for each array
-    int Ncomp = 1;
-
     // How Boxes are distrubuted among MPI processes
     amrex::DistributionMapping dm(ba);
 
     // we allocate two phi multifabs; one will store the old state, the other the new.
-    amrex::MultiFab phi_old(ba, dm, Ncomp, Nghost);
-    amrex::MultiFab phi_new(ba, dm, Ncomp, Nghost);
-
-    // time = starting time in the simulation
-    amrex::Real time = 0.0;
+    amrex::MultiFab rhs(ba, dm, 1, 0);
+    amrex::MultiFab phi(ba, dm, 1, 1);
 
     // **********************************
     // INITIALIZE DATA LOOP
     // **********************************
 
     // loop over boxes
-    for (amrex::MFIter mfi(phi_old); mfi.isValid(); ++mfi)
+    for (amrex::MFIter mfi(rhs); mfi.isValid(); ++mfi)
     {
         const amrex::Box& bx = mfi.validbox();
 
-        const amrex::Array4<amrex::Real>& phiOld = phi_old.array(mfi);
+        const amrex::Array4<amrex::Real>& rhs_p = rhs.array(mfi);
 
-        // set phi = 1 + e^(-(r-0.5)^2)
+        // set rhs = 1 + e^(-(r-0.5)^2)
         amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
         {
-
             // **********************************
             // SET VALUES FOR EACH CELL
             // **********************************
-
             amrex::Real x = (i+0.5) * dx[0];
             amrex::Real y = (j+0.5) * dx[1];
             amrex::Real z = (k+0.5) * dx[2];
             amrex::Real rsquared = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5))/0.01;
-            phiOld(i,j,k) = 1. + std::exp(-rsquared);
+            rhs_p(i,j,k) = 1. + std::exp(-rsquared);
         });
     }
 
@@ -149,75 +126,12 @@ int main (int argc, char* argv[])
     // **********************************
 
     // Write a plotfile of the initial data if plot_int > 0
-    if (plot_int > 0)
-    {
-        int step = 0;
-        const std::string& pltfile = amrex::Concatenate("plt",step,5);
-        WriteSingleLevelPlotfile(pltfile, phi_old, {"phi"}, geom, time, 0);
-    }
-
-
-    // **********************************
-    // MAIN TIME EVOLUTION LOOP
-    // **********************************
-
-    for (int step = 1; step <= nsteps; ++step)
-    {
-        // fill periodic ghost cells
-        phi_old.FillBoundary(geom.periodicity());
-
-        // new_phi = old_phi + dt * Laplacian(old_phi)
-        // loop over boxes
-        for ( amrex::MFIter mfi(phi_old); mfi.isValid(); ++mfi )
-        {
-            const amrex::Box& bx = mfi.validbox();
-
-            const amrex::Array4<amrex::Real>& phiOld = phi_old.array(mfi);
-            const amrex::Array4<amrex::Real>& phiNew = phi_new.array(mfi);
-
-            // advance the data by dt
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-            {
-
-                // **********************************
-                // EVOLVE VALUES FOR EACH CELL
-                // **********************************
-
-                phiNew(i,j,k) = phiOld(i,j,k) + dt *
-                    ( (phiOld(i+1,j,k) - 2.*phiOld(i,j,k) + phiOld(i-1,j,k)) / (dx[0]*dx[0])
-                     +(phiOld(i,j+1,k) - 2.*phiOld(i,j,k) + phiOld(i,j-1,k)) / (dx[1]*dx[1])
-                     +(phiOld(i,j,k+1) - 2.*phiOld(i,j,k) + phiOld(i,j,k-1)) / (dx[2]*dx[2])
-                        );
-            });
-        }
-
-        // **********************************
-        // INCREMENT
-        // **********************************
-
-        // update time
-        time = time + dt;
-
-        // copy new solution into old solution
-        amrex::MultiFab::Copy(phi_old, phi_new, 0, 0, 1, 0);
-
-        // Tell the I/O Processor to write out which step we're doing
-        amrex::Print() << "Advanced step " << step << "\n";
-
-
-        // **********************************
-        // WRITE PLOTFILE AT GIVEN INTERVAL
-        // **********************************
-
-        // Write a plotfile of the current data (plot_int was defined in the inputs file)
-        if (plot_int > 0 && step%plot_int == 0)
-        {
-            const std::string& pltfile = amrex::Concatenate("plt",step,5);
-            WriteSingleLevelPlotfile(pltfile, phi_new, {"phi"}, geom, time, step);
-        }
-    }
+    WriteSingleLevelPlotfile("rhs", rhs, {"rhs"}, geom, 0., 0);
 
+    amrex::GMRESPOISSON gmres_poisson(ba,dm);
 
+    gmres_poisson.solve(phi, rhs, 1.e-12, 0.);
+    
     }
     amrex::Finalize();
     return 0;

From 7fd891b8a2e500707b3b7290247676cdf2ff5121 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Fri, 15 Nov 2024 12:58:20 -0800
Subject: [PATCH 03/20] unused variables

---
 ExampleCodes/LinearSolvers/GMRES/Poisson/inputs | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs b/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs
index 5b098f82..9b5aeb18 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs
@@ -1,7 +1,2 @@
 n_cell = 32
 max_grid_size = 16
-
-nsteps = 1000
-plot_int = 100
-
-dt = 1.e-5

From afee2f56e9fcb4271e06ff00680912f95ab29628 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Fri, 15 Nov 2024 14:00:21 -0800
Subject: [PATCH 04/20] compiles but generates NaNs

---
 .../GMRES/Poisson/AMReX_GMRES_Poisson.H       | 81 ++++++++++---------
 .../LinearSolvers/GMRES/Poisson/main.cpp      | 21 ++---
 2 files changed, 46 insertions(+), 56 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index 974ad14d..bd72c073 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -21,7 +21,7 @@ public:
     using RT = amrex::Real; // typename MF::RT; // double or float
     using GM = GMRES<MF,GMRESPOISSONT<MF>>;
 
-    explicit GMRESPOISSONT (const BoxArray& ba, const DistributionMapping& dm);
+    explicit GMRESPOISSONT (const BoxArray& ba, const DistributionMapping& dm, const Geometry& geom);
 
     /**
      * \brief Solve the linear system
@@ -36,30 +36,9 @@ public:
     //! Sets verbosity.
     void setVerbose (int v) { m_gmres.setVerbose(v); }
 
-    //! Sets the max number of iterations
-    void setMaxIters (int niters) { m_gmres.setMaxIters(niters); }
-
-    //! Gets the number of iterations.
-    [[nodiscard]] int getNumIters () const { return m_gmres.getNumIters(); }
-
-    //! Gets the 2-norm of the residual.
-    [[nodiscard]] RT getResidualNorm () const { return m_gmres.getResidualNorm(); }
-
     //! Get the GMRES object.
     GM& getGMRES () { return m_gmres; }
 
-    /**
-     * \brief Set MLMG's multiplicative property of zero
-     *
-     * This should NOT be called unless MLMG has the multiplicative property
-     * of zero. MLMG is not a matrix, and usually does not have the
-     * properties of a matrix such as the multiplicative property of zero
-     * (i.e., M*0=0) because how domain boundary conditions are
-     * handled. However, if MLMG has the property of zero, calling this
-     * function with true can have a small performance benefit.
-     */
-    void setPropertyOfZero (bool b) { m_prop_zero = b; }
-
     //! Make MultiFab without ghost cells
     MF makeVecRHS () const;
 
@@ -85,28 +64,24 @@ public:
     static void linComb (MF& lhs, RT a, MF const& rhs_a, RT b, MF const& rhs_b);
 
     //! lhs = L(rhs)
-    void apply (MF& lhs, MF const& rhs) const;
+    void apply (MF& lhs, MF& rhs) const;
 
     void precond (MF& lhs, MF const& rhs) const;
 
     //! Control whether or not to use MLMG as preconditioner.
     bool usePrecond (bool new_flag) { return std::exchange(m_use_precond, new_flag); }
 
-    //! Set the number of MLMG preconditioner iterations per GMRES iteration.
-    void setPrecondNumIters (int precond_niters) { m_precond_niters = precond_niters; }
-
 private:
     GM m_gmres;
     BoxArray m_ba;
     DistributionMapping m_dm;
-    bool m_use_precond = true;
-    bool m_prop_zero = false;
-    int m_precond_niters = 1;
+    Geometry m_geom;
+    bool m_use_precond = false;
 };
 
 template <typename MF>
-GMRESPOISSONT<MF>::GMRESPOISSONT (const BoxArray& ba, const DistributionMapping& dm)
-    : m_ba(ba), m_dm(dm)
+GMRESPOISSONT<MF>::GMRESPOISSONT (const BoxArray& ba, const DistributionMapping& dm, const Geometry& geom)
+    : m_ba(ba), m_dm(dm), m_geom(geom)
 {
     m_gmres.define(*this);
 }
@@ -128,8 +103,8 @@ auto GMRESPOISSONT<MF>::makeVecLHS () const -> MF
 template <typename MF>
 auto GMRESPOISSONT<MF>::norm2 (MF const& mf) const -> RT
 {
-//    auto r = m_linop->xdoty(0, 0, mf, mf, false);
-//    return std::sqrt(r);
+    auto r = MultiFab::Dot(mf,0,1,0);
+    return std::sqrt(r);
     return 0.;
 }
 
@@ -142,8 +117,7 @@ void GMRESPOISSONT<MF>::scale (MF& mf, RT scale_factor)
 template <typename MF>
 auto GMRESPOISSONT<MF>::dotProduct (MF const& mf1, MF const& mf2) const -> RT
 {
-//    return m_linop->xdoty(0, 0, mf1, mf2, false);
-    return 0.;
+    return MultiFab::Dot(mf1,0,mf2,0,1,0);
 }
 
 template <typename MF>
@@ -171,11 +145,33 @@ void GMRESPOISSONT<MF>::linComb (MF& lhs, RT a, MF const& rhs_a, RT b, MF const&
 }
 
 template <typename MF>
-void GMRESPOISSONT<MF>::apply (MF& lhs, MF const& rhs) const
+void GMRESPOISSONT<MF>::apply (MF& lhs, MF& rhs) const
 {
-//    m_linop->apply(0, 0, lhs, const_cast<MF&>(rhs),
-//                   MLLinOpT<MF>::BCMode::Homogeneous,
-//                   MLLinOpT<MF>::StateMode::Correction);
+    // apply matrix to rhs for output lhs
+    rhs.FillBoundary(m_geom.periodicity());
+    
+    const GpuArray<Real, AMREX_SPACEDIM> dx = m_geom.CellSizeArray();
+
+    for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+        
+        const Box& bx = mfi.tilebox();
+
+        const Array4<const Real> & rhs_p = rhs.array(mfi);
+        const Array4<      Real> & lhs_p = lhs.array(mfi);
+
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        {
+            lhs_p(i,j,k) = ( rhs_p(i+1,j,k) - 2.*rhs_p(i,j,k) + rhs_p(i-1,j,k) ) / (dx[0]*dx[0])
+                         + ( rhs_p(i,j+1,k) - 2.*rhs_p(i,j,k) + rhs_p(i,j-1,k) ) / (dx[1]*dx[1])
+#if (AMREX_SPACEDIM == 3)
+                         + ( rhs_p(i,j,k+1) - 2.*rhs_p(i,j,k) + rhs_p(i,j,k-1) ) / (dx[2]*dx[2])
+#endif
+                ;
+        });
+
+    }
+
+    
 }
 
 template <typename MF>
@@ -183,7 +179,6 @@ void GMRESPOISSONT<MF>::precond (MF& lhs, MF const& rhs) const
 {
     if (m_use_precond) {
 
-
     } else {
         LocalCopy(lhs, rhs, 0, 0, nComp(lhs), IntVect(0));
     }
@@ -201,6 +196,12 @@ void GMRESPOISSONT<MF>::solve (MF& a_sol, MF const& a_rhs, RT a_tol_rel, RT a_to
     m_gmres.solve(cor, res, a_tol_rel, a_tol_abs); // L(cor) = res
     increment(a_sol, cor, RT(-1)); // sol = sol - cor
     */
+    auto res = makeVecRHS();
+    apply(res,a_sol);
+    increment(res, a_rhs, RT(-1)); // res = L(sol) - rhs
+    auto cor = makeVecLHS();
+    m_gmres.solve(cor, res, a_tol_rel, a_tol_abs); // L(cor) = res
+    increment(a_sol, cor, RT(-1)); // sol = sol - cor
 }
 
 using GMRESPOISSON = GMRESPOISSONT<MultiFab>;
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
index 8e90595a..da66ba00 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
@@ -25,15 +25,6 @@ int main (int argc, char* argv[])
     // size of each box (or grid)
     int max_grid_size;
 
-    // total steps in simulation
-    int nsteps;
-
-    // how often to write a plotfile
-    int plot_int;
-
-    // time step
-    amrex::Real dt;
-
     // **********************************
     // READ PARAMETER VALUES FROM INPUT DATA
     // **********************************
@@ -78,7 +69,7 @@ int main (int argc, char* argv[])
 
     // This defines the physical box, [0,1] in each direction.
     amrex::RealBox real_box({ 0., 0., 0.},
-                     { 1., 1., 1.});
+                            { 1., 1., 1.});
 
     // periodic in all direction
     amrex::Array<int,3> is_periodic{1,1,1};
@@ -121,16 +112,14 @@ int main (int argc, char* argv[])
         });
     }
 
-    // **********************************
-    // WRITE INITIAL PLOT FILE
-    // **********************************
-
-    // Write a plotfile of the initial data if plot_int > 0
     WriteSingleLevelPlotfile("rhs", rhs, {"rhs"}, geom, 0., 0);
 
-    amrex::GMRESPOISSON gmres_poisson(ba,dm);
+    amrex::GMRESPOISSON gmres_poisson(ba,dm,geom);
 
+    gmres_poisson.setVerbose(2);
     gmres_poisson.solve(phi, rhs, 1.e-12, 0.);
+
+    WriteSingleLevelPlotfile("phi", phi, {"phi"}, geom, 0., 0);
     
     }
     amrex::Finalize();

From fff361685f84381e8d7106af40cc458233405c39 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Fri, 15 Nov 2024 14:13:56 -0800
Subject: [PATCH 05/20] fixes

---
 ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H | 3 +--
 ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp              | 2 +-
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index bd72c073..bd21ece2 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -95,7 +95,7 @@ auto GMRESPOISSONT<MF>::makeVecRHS () const -> MF
 template <typename MF>
 auto GMRESPOISSONT<MF>::makeVecLHS () const -> MF
 {
-    MF mf(m_ba, m_dm, 1, 0);
+    MF mf(m_ba, m_dm, 1, 1);
     mf.setVal(0.);
     return mf;
 }
@@ -105,7 +105,6 @@ auto GMRESPOISSONT<MF>::norm2 (MF const& mf) const -> RT
 {
     auto r = MultiFab::Dot(mf,0,1,0);
     return std::sqrt(r);
-    return 0.;
 }
 
 template <typename MF>
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
index da66ba00..96a2cba9 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
@@ -55,7 +55,7 @@ int main (int argc, char* argv[])
     amrex::Geometry geom;
 
     // define lower and upper indices
-    amrex::IntVect dom_lo(0,0,0);
+    amrex::IntVect dom_lo(       0,        0,        0);
     amrex::IntVect dom_hi(n_cell-1, n_cell-1, n_cell-1);
 
     // Make a single box that is the entire domain

From 5c678a60add62c52795fcd3b904d6288f6beeec9 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Fri, 15 Nov 2024 14:50:53 -0800
Subject: [PATCH 06/20] no more NaNs at least :shrug:

---
 ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H | 1 +
 ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp              | 3 +++
 2 files changed, 4 insertions(+)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index bd21ece2..5bcc21ef 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -199,6 +199,7 @@ void GMRESPOISSONT<MF>::solve (MF& a_sol, MF const& a_rhs, RT a_tol_rel, RT a_to
     apply(res,a_sol);
     increment(res, a_rhs, RT(-1)); // res = L(sol) - rhs
     auto cor = makeVecLHS();
+    res.FillBoundary(m_geom.periodicity()); // FIXME not sure this is needed
     m_gmres.solve(cor, res, a_tol_rel, a_tol_abs); // L(cor) = res
     increment(a_sol, cor, RT(-1)); // sol = sol - cor
 }
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
index 96a2cba9..315c7cac 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
@@ -116,6 +116,9 @@ int main (int argc, char* argv[])
 
     amrex::GMRESPOISSON gmres_poisson(ba,dm,geom);
 
+    // initial guess
+    phi.setVal(0.);
+
     gmres_poisson.setVerbose(2);
     gmres_poisson.solve(phi, rhs, 1.e-12, 0.);
 

From c1b8ae0e55fe920929bbd01b9abea08f6d8f8922 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Mon, 25 Nov 2024 10:38:11 -0800
Subject: [PATCH 07/20] multifab specific routines

---
 .../GMRES/Poisson/AMReX_GMRES_Poisson.H            | 14 ++++++--------
 1 file changed, 6 insertions(+), 8 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index 5bcc21ef..08bb4f7d 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -95,9 +95,7 @@ auto GMRESPOISSONT<MF>::makeVecRHS () const -> MF
 template <typename MF>
 auto GMRESPOISSONT<MF>::makeVecLHS () const -> MF
 {
-    MF mf(m_ba, m_dm, 1, 1);
-    mf.setVal(0.);
-    return mf;
+    return MultiFab(m_ba, m_dm, 1, 1);
 }
 
 template <typename MF>
@@ -110,7 +108,7 @@ auto GMRESPOISSONT<MF>::norm2 (MF const& mf) const -> RT
 template <typename MF>
 void GMRESPOISSONT<MF>::scale (MF& mf, RT scale_factor)
 {
-    Scale(mf, scale_factor, 0, nComp(mf), 0);
+    mf.mult(scale_factor);
 }
 
 template <typename MF>
@@ -122,25 +120,25 @@ auto GMRESPOISSONT<MF>::dotProduct (MF const& mf1, MF const& mf2) const -> RT
 template <typename MF>
 void GMRESPOISSONT<MF>::setToZero (MF& lhs)
 {
-    setVal(lhs, RT(0.0));
+    lhs.setVal(0.);
 }
 
 template <typename MF>
 void GMRESPOISSONT<MF>::assign (MF& lhs, MF const& rhs)
 {
-    LocalCopy(lhs, rhs, 0, 0, nComp(lhs), IntVect(0));
+    MultiFab::Copy(lhs,rhs,0,0,1,0);
 }
 
 template <typename MF>
 void GMRESPOISSONT<MF>::increment (MF& lhs, MF const& rhs, RT a)
 {
-    Saxpy(lhs, a, rhs, 0, 0, nComp(lhs), IntVect(0));
+    MultiFab::Saxpy(lhs,a,rhs,0,0,1,0);
 }
 
 template <typename MF>
 void GMRESPOISSONT<MF>::linComb (MF& lhs, RT a, MF const& rhs_a, RT b, MF const& rhs_b)
 {
-    LinComb(lhs, a, rhs_a, 0, b, rhs_b, 0, 0, nComp(lhs), IntVect(0));
+    MultiFab::LinComb(lhs,a,rhs_a,0,b,rhs_b,0,0,1,0);
 }
 
 template <typename MF>

From 62cac5e34accfa9f79f564b669e6b093276a5dd2 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Mon, 25 Nov 2024 11:02:08 -0800
Subject: [PATCH 08/20] seems to work

---
 .../GMRES/Poisson/AMReX_GMRES_Poisson.H       | 22 +++----------------
 .../LinearSolvers/GMRES/Poisson/main.cpp      |  3 +--
 2 files changed, 4 insertions(+), 21 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index 08bb4f7d..65f16501 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -101,8 +101,7 @@ auto GMRESPOISSONT<MF>::makeVecLHS () const -> MF
 template <typename MF>
 auto GMRESPOISSONT<MF>::norm2 (MF const& mf) const -> RT
 {
-    auto r = MultiFab::Dot(mf,0,1,0);
-    return std::sqrt(r);
+    return mf.norm2();
 }
 
 template <typename MF>
@@ -177,29 +176,14 @@ void GMRESPOISSONT<MF>::precond (MF& lhs, MF const& rhs) const
     if (m_use_precond) {
 
     } else {
-        LocalCopy(lhs, rhs, 0, 0, nComp(lhs), IntVect(0));
+        MultiFab::Copy(lhs,rhs,0,0,1,0);
     }
 }
 
 template <typename MF>
 void GMRESPOISSONT<MF>::solve (MF& a_sol, MF const& a_rhs, RT a_tol_rel, RT a_tol_abs)
 {
-    /*
-    auto res = makeVecRHS();
-    m_mlmg->apply({&res}, {&a_sol}); // res = L(sol)
-    increment(res, a_rhs, RT(-1)); // res = L(sol) - rhs
-    auto cor = makeVecLHS();
-    m_linop->setDirichletNodesToZero(0,0,res);
-    m_gmres.solve(cor, res, a_tol_rel, a_tol_abs); // L(cor) = res
-    increment(a_sol, cor, RT(-1)); // sol = sol - cor
-    */
-    auto res = makeVecRHS();
-    apply(res,a_sol);
-    increment(res, a_rhs, RT(-1)); // res = L(sol) - rhs
-    auto cor = makeVecLHS();
-    res.FillBoundary(m_geom.periodicity()); // FIXME not sure this is needed
-    m_gmres.solve(cor, res, a_tol_rel, a_tol_abs); // L(cor) = res
-    increment(a_sol, cor, RT(-1)); // sol = sol - cor
+    m_gmres.solve(a_sol, a_rhs, a_tol_rel, a_tol_abs);
 }
 
 using GMRESPOISSON = GMRESPOISSONT<MultiFab>;
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
index 315c7cac..684a042a 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
@@ -107,8 +107,7 @@ int main (int argc, char* argv[])
             amrex::Real x = (i+0.5) * dx[0];
             amrex::Real y = (j+0.5) * dx[1];
             amrex::Real z = (k+0.5) * dx[2];
-            amrex::Real rsquared = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5))/0.01;
-            rhs_p(i,j,k) = 1. + std::exp(-rsquared);
+            rhs_p(i,j,k) = sin(2.*M_PI*x) * sin(4.*M_PI*y) * sin(8.*M_PI*z);
         });
     }
 

From f0fd9bd246bee5a6594f10d13b7543549e33a939 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Mon, 25 Nov 2024 11:08:56 -0800
Subject: [PATCH 09/20] trailing whitespace

---
 .../LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H       | 6 +++---
 ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp           | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index 65f16501..73faaa86 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -145,11 +145,11 @@ void GMRESPOISSONT<MF>::apply (MF& lhs, MF& rhs) const
 {
     // apply matrix to rhs for output lhs
     rhs.FillBoundary(m_geom.periodicity());
-    
+
     const GpuArray<Real, AMREX_SPACEDIM> dx = m_geom.CellSizeArray();
 
     for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
-        
+
         const Box& bx = mfi.tilebox();
 
         const Array4<const Real> & rhs_p = rhs.array(mfi);
@@ -167,7 +167,7 @@ void GMRESPOISSONT<MF>::apply (MF& lhs, MF& rhs) const
 
     }
 
-    
+
 }
 
 template <typename MF>
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
index 684a042a..595006e3 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
@@ -122,7 +122,7 @@ int main (int argc, char* argv[])
     gmres_poisson.solve(phi, rhs, 1.e-12, 0.);
 
     WriteSingleLevelPlotfile("phi", phi, {"phi"}, geom, 0., 0);
-    
+
     }
     amrex::Finalize();
     return 0;

From 1676eb645c355fb1ac582de437b17e9d3fcf1690 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Mon, 25 Nov 2024 11:59:43 -0800
Subject: [PATCH 10/20] fix header guard

---
 .../LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H         | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index 73faaa86..61a5405c 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -1,5 +1,5 @@
-#ifndef AMREX_GMRES_MLMG_H_
-#define AMREX_GMRES_MLMG_H_
+#ifndef AMREX_GMRES_POISSON_H_
+#define AMREX_GMRES_POISSON_H_
 #include <AMReX_Config.H>
 
 #include <AMReX_GMRES.H>

From 112d2f9de5fa36ad678a88a4e200b63b51d37c47 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Mon, 25 Nov 2024 12:00:29 -0800
Subject: [PATCH 11/20] fix comments

---
 .../LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H       | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index 61a5405c..635d17b8 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -8,11 +8,7 @@
 namespace amrex {
 
 /**
- * \brief Solve using GMRES with multigrid as preconditioner
- *
- * The linear system to solve is provided by MLMG, which is also being used
- * as the preconditioner.
- *
+ * \brief Solve Poisson's equation using unpreconditioned GMRES
  */
 template <typename MF>
 class GMRESPOISSONT

From d5c09134590e97191b1954c1f4c91a5060da8ca8 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 27 Nov 2024 11:27:15 -0800
Subject: [PATCH 12/20] fix comments randomize RHS Jacobi preconditioner
 (implemented but probably buggy as convergence slows down)

---
 .../GMRES/Poisson/AMReX_GMRES_Poisson.H       | 39 ++++++++++++++++++-
 .../LinearSolvers/GMRES/Poisson/main.cpp      | 25 +++++-------
 2 files changed, 47 insertions(+), 17 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index 635d17b8..185db16b 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -73,6 +73,9 @@ private:
     DistributionMapping m_dm;
     Geometry m_geom;
     bool m_use_precond = false;
+
+    // store the original RHS needed for Jacobi preconditioner
+    MultiFab orig_rhs;
 };
 
 template <typename MF>
@@ -144,7 +147,7 @@ void GMRESPOISSONT<MF>::apply (MF& lhs, MF& rhs) const
 
     const GpuArray<Real, AMREX_SPACEDIM> dx = m_geom.CellSizeArray();
 
-    for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+    for ( MFIter mfi(lhs,TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
 
         const Box& bx = mfi.tilebox();
 
@@ -166,11 +169,42 @@ void GMRESPOISSONT<MF>::apply (MF& lhs, MF& rhs) const
 
 }
 
+// applies preconditioner to rhs. If there is no preconditioner,
+// this function should do lhs = rhs.
 template <typename MF>
 void GMRESPOISSONT<MF>::precond (MF& lhs, MF const& rhs) const
 {
     if (m_use_precond) {
 
+        MultiFab rhs_tmp(m_ba,m_dm,1,1);
+        MultiFab::Copy(rhs_tmp,rhs,0,0,1,0);
+        rhs_tmp.FillBoundary(m_geom.periodicity());
+
+        const GpuArray<Real, AMREX_SPACEDIM> dx = m_geom.CellSizeArray();
+
+        amrex::Real fac = 1.;
+        for (int d=0; d<AMREX_SPACEDIM; ++d) fac += 2./(dx[d]*dx[d]);
+
+        for ( MFIter mfi(lhs,TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+
+            const Box& bx = mfi.tilebox();
+
+            const Array4<const Real> & orig_rhs_p = orig_rhs.array(mfi);
+            const Array4<const Real> & rhs_p = rhs_tmp.array(mfi);
+            const Array4<      Real> & lhs_p = lhs.array(mfi);
+
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                lhs_p(i,j,k) = ( rhs_p(i+1,j,k) - 2.*rhs_p(i,j,k) + rhs_p(i-1,j,k) ) / (dx[0]*dx[0])
+                             + ( rhs_p(i,j+1,k) - 2.*rhs_p(i,j,k) + rhs_p(i,j-1,k) ) / (dx[1]*dx[1])
+#if (AMREX_SPACEDIM == 3)
+                             + ( rhs_p(i,j,k+1) - 2.*rhs_p(i,j,k) + rhs_p(i,j,k-1) ) / (dx[2]*dx[2])
+#endif
+                    ;
+
+                lhs_p(i,j,k) = rhs_p(i,j,k) + (orig_rhs_p(i,j,k) - lhs_p(i,j,k) ) / fac;
+            });
+        }
     } else {
         MultiFab::Copy(lhs,rhs,0,0,1,0);
     }
@@ -179,6 +213,9 @@ void GMRESPOISSONT<MF>::precond (MF& lhs, MF const& rhs) const
 template <typename MF>
 void GMRESPOISSONT<MF>::solve (MF& a_sol, MF const& a_rhs, RT a_tol_rel, RT a_tol_abs)
 {
+    orig_rhs.define(m_ba,m_dm,1,0);
+    MultiFab::Copy(orig_rhs,a_rhs,0,0,1,0);
+    
     m_gmres.solve(a_sol, a_rhs, a_tol_rel, a_tol_abs);
 }
 
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
index 595006e3..3b00201c 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
@@ -1,10 +1,7 @@
 /*
- * A simplified single file version of the HeatEquation_EX0_C exmaple.
- * This code is designed to be used with Demo_Tutorial.rst.
- *
+ * A simplified usage of the AMReX GMRES class
  */
 
-
 #include <AMReX.H>
 #include <AMReX_PlotFileUtil.H>
 #include <AMReX_ParmParse.H>
@@ -77,9 +74,6 @@ int main (int argc, char* argv[])
     // This defines a Geometry object
     geom.define(domain, real_box, amrex::CoordSys::cartesian, is_periodic);
 
-    // extract dx from the geometry object
-    amrex::GpuArray<amrex::Real,3> dx = geom.CellSizeArray();
-
     // How Boxes are distrubuted among MPI processes
     amrex::DistributionMapping dm(ba);
 
@@ -98,19 +92,18 @@ int main (int argc, char* argv[])
 
         const amrex::Array4<amrex::Real>& rhs_p = rhs.array(mfi);
 
-        // set rhs = 1 + e^(-(r-0.5)^2)
-        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
+        // fill rhs with random numbers
+        amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::RandomEngine const& engine) noexcept
         {
-            // **********************************
-            // SET VALUES FOR EACH CELL
-            // **********************************
-            amrex::Real x = (i+0.5) * dx[0];
-            amrex::Real y = (j+0.5) * dx[1];
-            amrex::Real z = (k+0.5) * dx[2];
-            rhs_p(i,j,k) = sin(2.*M_PI*x) * sin(4.*M_PI*y) * sin(8.*M_PI*z);
+            rhs_p(i,j,k) = amrex::RandomNormal(0.,1.,engine);
         });
     }
 
+    // offset data so the rhs sums to zero
+    amrex::Real sum = rhs.sum();
+    amrex::Long npts = rhs.boxArray().numPts();
+    rhs.plus(-sum/npts,0,1);
+
     WriteSingleLevelPlotfile("rhs", rhs, {"rhs"}, geom, 0., 0);
 
     amrex::GMRESPOISSON gmres_poisson(ba,dm,geom);

From 52dbb11a4c7178b1a550f99ed91962eda13a082e Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 27 Nov 2024 11:40:05 -0800
Subject: [PATCH 13/20] remove templating

---
 .../GMRES/Poisson/AMReX_GMRES_Poisson.H       | 74 ++++++++-----------
 1 file changed, 29 insertions(+), 45 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index 185db16b..b2c65665 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -10,14 +10,13 @@ namespace amrex {
 /**
  * \brief Solve Poisson's equation using unpreconditioned GMRES
  */
-template <typename MF>
-class GMRESPOISSONT
+class GMRESPOISSON
 {
 public:
-    using RT = amrex::Real; // typename MF::RT; // double or float
-    using GM = GMRES<MF,GMRESPOISSONT<MF>>;
+    using RT = amrex::Real; // typename MultiFab::RT; // double or float
+    using GM = GMRES<MultiFab,GMRESPOISSON>;
 
-    explicit GMRESPOISSONT (const BoxArray& ba, const DistributionMapping& dm, const Geometry& geom);
+    explicit GMRESPOISSON (const BoxArray& ba, const DistributionMapping& dm, const Geometry& geom);
 
     /**
      * \brief Solve the linear system
@@ -27,7 +26,7 @@ public:
      * \param a_tol_rel relative tolerance.
      * \param a_tol_abs absolute tolerance.
      */
-    void solve (MF& a_sol, MF const& a_rhs, RT a_tol_rel, RT a_tol_abs);
+    void solve (MultiFab& a_sol, MultiFab const& a_rhs, RT a_tol_rel, RT a_tol_abs);
 
     //! Sets verbosity.
     void setVerbose (int v) { m_gmres.setVerbose(v); }
@@ -36,33 +35,33 @@ public:
     GM& getGMRES () { return m_gmres; }
 
     //! Make MultiFab without ghost cells
-    MF makeVecRHS () const;
+    MultiFab makeVecRHS () const;
 
     //! Make MultiFab with ghost cells and set ghost cells to zero
-    MF makeVecLHS () const;
+    MultiFab makeVecLHS () const;
 
-    RT norm2 (MF const& mf) const;
+    RT norm2 (MultiFab const& mf) const;
 
-    static void scale (MF& mf, RT scale_factor);
+    static void scale (MultiFab& mf, RT scale_factor);
 
-    RT dotProduct (MF const& mf1, MF const& mf2) const;
+    RT dotProduct (MultiFab const& mf1, MultiFab const& mf2) const;
 
     //! lhs = 0
-    static void setToZero (MF& lhs);
+    static void setToZero (MultiFab& lhs);
 
     //! lhs = rhs
-    static void assign (MF& lhs, MF const& rhs);
+    static void assign (MultiFab& lhs, MultiFab const& rhs);
 
     //! lhs += a*rhs
-    static void increment (MF& lhs, MF const& rhs, RT a);
+    static void increment (MultiFab& lhs, MultiFab const& rhs, RT a);
 
     //! lhs = a*rhs_a + b*rhs_b
-    static void linComb (MF& lhs, RT a, MF const& rhs_a, RT b, MF const& rhs_b);
+    static void linComb (MultiFab& lhs, RT a, MultiFab const& rhs_a, RT b, MultiFab const& rhs_b);
 
     //! lhs = L(rhs)
-    void apply (MF& lhs, MF& rhs) const;
+    void apply (MultiFab& lhs, MultiFab& rhs) const;
 
-    void precond (MF& lhs, MF const& rhs) const;
+    void precond (MultiFab& lhs, MultiFab const& rhs) const;
 
     //! Control whether or not to use MLMG as preconditioner.
     bool usePrecond (bool new_flag) { return std::exchange(m_use_precond, new_flag); }
@@ -78,69 +77,58 @@ private:
     MultiFab orig_rhs;
 };
 
-template <typename MF>
-GMRESPOISSONT<MF>::GMRESPOISSONT (const BoxArray& ba, const DistributionMapping& dm, const Geometry& geom)
+GMRESPOISSON::GMRESPOISSON (const BoxArray& ba, const DistributionMapping& dm, const Geometry& geom)
     : m_ba(ba), m_dm(dm), m_geom(geom)
 {
     m_gmres.define(*this);
 }
 
-template <typename MF>
-auto GMRESPOISSONT<MF>::makeVecRHS () const -> MF
+auto GMRESPOISSON::makeVecRHS () const -> MultiFab
 {
     return MultiFab(m_ba, m_dm, 1, 0);
 }
 
-template <typename MF>
-auto GMRESPOISSONT<MF>::makeVecLHS () const -> MF
+auto GMRESPOISSON::makeVecLHS () const -> MultiFab
 {
     return MultiFab(m_ba, m_dm, 1, 1);
 }
 
-template <typename MF>
-auto GMRESPOISSONT<MF>::norm2 (MF const& mf) const -> RT
+auto GMRESPOISSON::norm2 (MultiFab const& mf) const -> RT
 {
     return mf.norm2();
 }
 
-template <typename MF>
-void GMRESPOISSONT<MF>::scale (MF& mf, RT scale_factor)
+void GMRESPOISSON::scale (MultiFab& mf, RT scale_factor)
 {
     mf.mult(scale_factor);
 }
 
-template <typename MF>
-auto GMRESPOISSONT<MF>::dotProduct (MF const& mf1, MF const& mf2) const -> RT
+auto GMRESPOISSON::dotProduct (MultiFab const& mf1, MultiFab const& mf2) const -> RT
 {
     return MultiFab::Dot(mf1,0,mf2,0,1,0);
 }
 
-template <typename MF>
-void GMRESPOISSONT<MF>::setToZero (MF& lhs)
+void GMRESPOISSON::setToZero (MultiFab& lhs)
 {
     lhs.setVal(0.);
 }
 
-template <typename MF>
-void GMRESPOISSONT<MF>::assign (MF& lhs, MF const& rhs)
+void GMRESPOISSON::assign (MultiFab& lhs, MultiFab const& rhs)
 {
     MultiFab::Copy(lhs,rhs,0,0,1,0);
 }
 
-template <typename MF>
-void GMRESPOISSONT<MF>::increment (MF& lhs, MF const& rhs, RT a)
+void GMRESPOISSON::increment (MultiFab& lhs, MultiFab const& rhs, RT a)
 {
     MultiFab::Saxpy(lhs,a,rhs,0,0,1,0);
 }
 
-template <typename MF>
-void GMRESPOISSONT<MF>::linComb (MF& lhs, RT a, MF const& rhs_a, RT b, MF const& rhs_b)
+void GMRESPOISSON::linComb (MultiFab& lhs, RT a, MultiFab const& rhs_a, RT b, MultiFab const& rhs_b)
 {
     MultiFab::LinComb(lhs,a,rhs_a,0,b,rhs_b,0,0,1,0);
 }
 
-template <typename MF>
-void GMRESPOISSONT<MF>::apply (MF& lhs, MF& rhs) const
+void GMRESPOISSON::apply (MultiFab& lhs, MultiFab& rhs) const
 {
     // apply matrix to rhs for output lhs
     rhs.FillBoundary(m_geom.periodicity());
@@ -171,8 +159,7 @@ void GMRESPOISSONT<MF>::apply (MF& lhs, MF& rhs) const
 
 // applies preconditioner to rhs. If there is no preconditioner,
 // this function should do lhs = rhs.
-template <typename MF>
-void GMRESPOISSONT<MF>::precond (MF& lhs, MF const& rhs) const
+void GMRESPOISSON::precond (MultiFab& lhs, MultiFab const& rhs) const
 {
     if (m_use_precond) {
 
@@ -210,8 +197,7 @@ void GMRESPOISSONT<MF>::precond (MF& lhs, MF const& rhs) const
     }
 }
 
-template <typename MF>
-void GMRESPOISSONT<MF>::solve (MF& a_sol, MF const& a_rhs, RT a_tol_rel, RT a_tol_abs)
+void GMRESPOISSON::solve (MultiFab& a_sol, MultiFab const& a_rhs, RT a_tol_rel, RT a_tol_abs)
 {
     orig_rhs.define(m_ba,m_dm,1,0);
     MultiFab::Copy(orig_rhs,a_rhs,0,0,1,0);
@@ -219,8 +205,6 @@ void GMRESPOISSONT<MF>::solve (MF& a_sol, MF const& a_rhs, RT a_tol_rel, RT a_to
     m_gmres.solve(a_sol, a_rhs, a_tol_rel, a_tol_abs);
 }
 
-using GMRESPOISSON = GMRESPOISSONT<MultiFab>;
-
 }
 
 #endif

From c585d3e1fa51ca22f99dd7c894c0be125937d2d2 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 27 Nov 2024 11:41:23 -0800
Subject: [PATCH 14/20] trailing whitespace

---
 ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index b2c65665..a756576c 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -201,7 +201,7 @@ void GMRESPOISSON::solve (MultiFab& a_sol, MultiFab const& a_rhs, RT a_tol_rel,
 {
     orig_rhs.define(m_ba,m_dm,1,0);
     MultiFab::Copy(orig_rhs,a_rhs,0,0,1,0);
-    
+
     m_gmres.solve(a_sol, a_rhs, a_tol_rel, a_tol_abs);
 }
 

From 4a1dcaad07c008a5d6f870a0929d4bd9dc6357f2 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 27 Nov 2024 11:43:56 -0800
Subject: [PATCH 15/20] fix comment

---
 ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index a756576c..38b065c4 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -13,7 +13,7 @@ namespace amrex {
 class GMRESPOISSON
 {
 public:
-    using RT = amrex::Real; // typename MultiFab::RT; // double or float
+    using RT = amrex::Real; // double or float
     using GM = GMRES<MultiFab,GMRESPOISSON>;
 
     explicit GMRESPOISSON (const BoxArray& ba, const DistributionMapping& dm, const Geometry& geom);

From 1dd9b6f6fd99d7d223303a2f79867ffd8b62c2c3 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 27 Nov 2024 11:47:14 -0800
Subject: [PATCH 16/20] preconditioner control

---
 .../LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H      | 4 ++--
 ExampleCodes/LinearSolvers/GMRES/Poisson/inputs            | 2 ++
 ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp          | 7 +++++++
 3 files changed, 11 insertions(+), 2 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index 38b065c4..2e77a736 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -8,7 +8,7 @@
 namespace amrex {
 
 /**
- * \brief Solve Poisson's equation using unpreconditioned GMRES
+ * \brief Solve Poisson's equation using amrex GMRES class
  */
 class GMRESPOISSON
 {
@@ -71,7 +71,7 @@ private:
     BoxArray m_ba;
     DistributionMapping m_dm;
     Geometry m_geom;
-    bool m_use_precond = false;
+    bool m_use_precond;
 
     // store the original RHS needed for Jacobi preconditioner
     MultiFab orig_rhs;
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs b/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs
index 9b5aeb18..bcf33ac7 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs
@@ -1,2 +1,4 @@
 n_cell = 32
 max_grid_size = 16
+
+use_precond = 0
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
index 3b00201c..28f3f7d3 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
@@ -22,6 +22,9 @@ int main (int argc, char* argv[])
     // size of each box (or grid)
     int max_grid_size;
 
+    // use preconditioner?
+    int use_precond;
+
     // **********************************
     // READ PARAMETER VALUES FROM INPUT DATA
     // **********************************
@@ -38,6 +41,9 @@ int main (int argc, char* argv[])
 
         // The domain is broken into boxes of size max_grid_size
         pp.get("max_grid_size",max_grid_size);
+
+        use_precond = 0;
+        pp.query("use_precond",use_precond);
     }
 
     // **********************************
@@ -111,6 +117,7 @@ int main (int argc, char* argv[])
     // initial guess
     phi.setVal(0.);
 
+    gmres_poisson.usePrecond(use_precond);
     gmres_poisson.setVerbose(2);
     gmres_poisson.solve(phi, rhs, 1.e-12, 0.);
 

From 4ed14fe4fa35a54b57448dbd5409874a4d79f1b9 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 27 Nov 2024 12:12:43 -0800
Subject: [PATCH 17/20] bugfix in Jacobi - still not quite right

---
 .../LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H         | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index 2e77a736..c8e087ad 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -169,8 +169,8 @@ void GMRESPOISSON::precond (MultiFab& lhs, MultiFab const& rhs) const
 
         const GpuArray<Real, AMREX_SPACEDIM> dx = m_geom.CellSizeArray();
 
-        amrex::Real fac = 1.;
-        for (int d=0; d<AMREX_SPACEDIM; ++d) fac += 2./(dx[d]*dx[d]);
+        amrex::Real fac = 0.;
+        for (int d=0; d<AMREX_SPACEDIM; ++d) fac -= 2./(dx[d]*dx[d]);
 
         for ( MFIter mfi(lhs,TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
 

From a0854a595f045bf5f7337146654f8ca515fbe0f4 Mon Sep 17 00:00:00 2001
From: Weiqun Zhang <weiqunzhang@lbl.gov>
Date: Mon, 2 Dec 2024 10:31:01 -0800
Subject: [PATCH 18/20] Fix Jacobi preconditioner

---
 .../GMRES/Poisson/AMReX_GMRES_Poisson.H       | 54 ++++++++++---------
 1 file changed, 29 insertions(+), 25 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index c8e087ad..b28119dd 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -162,35 +162,39 @@ void GMRESPOISSON::apply (MultiFab& lhs, MultiFab& rhs) const
 void GMRESPOISSON::precond (MultiFab& lhs, MultiFab const& rhs) const
 {
     if (m_use_precond) {
-
-        MultiFab rhs_tmp(m_ba,m_dm,1,1);
-        MultiFab::Copy(rhs_tmp,rhs,0,0,1,0);
-        rhs_tmp.FillBoundary(m_geom.periodicity());
-
         const GpuArray<Real, AMREX_SPACEDIM> dx = m_geom.CellSizeArray();
-
         amrex::Real fac = 0.;
-        for (int d=0; d<AMREX_SPACEDIM; ++d) fac -= 2./(dx[d]*dx[d]);
-
-        for ( MFIter mfi(lhs,TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
-
-            const Box& bx = mfi.tilebox();
-
-            const Array4<const Real> & orig_rhs_p = orig_rhs.array(mfi);
-            const Array4<const Real> & rhs_p = rhs_tmp.array(mfi);
-            const Array4<      Real> & lhs_p = lhs.array(mfi);
-
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
-            {
-                lhs_p(i,j,k) = ( rhs_p(i+1,j,k) - 2.*rhs_p(i,j,k) + rhs_p(i-1,j,k) ) / (dx[0]*dx[0])
-                             + ( rhs_p(i,j+1,k) - 2.*rhs_p(i,j,k) + rhs_p(i,j-1,k) ) / (dx[1]*dx[1])
+        for (int d=0; d<AMREX_SPACEDIM; ++d) { fac -= 2./(dx[d]*dx[d]); }
+
+        MultiFab tmp(lhs.boxArray(), lhs.DistributionMap(), 1, 1);
+        auto const& tmp_ma = tmp.const_arrays();
+        auto const& rh_ma = rhs.const_arrays();
+        auto const& x_ma = lhs.arrays();
+
+        const int niters = 8;
+        for (int iter = 0; iter < niters; ++iter) {
+            if (iter == 0) {
+                ParallelFor(lhs, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
+                {
+                    x_ma[b](i,j,k) = rh_ma[b](i,j,k) / fac;
+                });
+            } else {
+                MultiFab::Copy(tmp, lhs, 0, 0, 1, 0);
+                tmp.FillBoundary(m_geom.periodicity());
+                ParallelFor(lhs, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
+                {
+                    auto const& phi = tmp_ma[b];
+                    auto ax = ( phi(i+1,j,k) - 2.*phi(i,j,k) + phi(i-1,j,k) ) / (dx[0]*dx[0])
+                            + ( phi(i,j+1,k) - 2.*phi(i,j,k) + phi(i,j-1,k) ) / (dx[1]*dx[1])
 #if (AMREX_SPACEDIM == 3)
-                             + ( rhs_p(i,j,k+1) - 2.*rhs_p(i,j,k) + rhs_p(i,j,k-1) ) / (dx[2]*dx[2])
+                            + ( phi(i,j,k+1) - 2.*phi(i,j,k) + phi(i,j,k-1) ) / (dx[2]*dx[2])
 #endif
-                    ;
-
-                lhs_p(i,j,k) = rhs_p(i,j,k) + (orig_rhs_p(i,j,k) - lhs_p(i,j,k) ) / fac;
-            });
+                        ;
+                    auto res = rh_ma[b](i,j,k) - ax;
+                    x_ma[b](i,j,k) += res / fac * Real(2./3.); // 2/3: weighted jacobi
+                });
+            }
+            Gpu::streamSynchronize();
         }
     } else {
         MultiFab::Copy(lhs,rhs,0,0,1,0);

From 2751710b688b7b73b3de04dc773bf04c0907ea9e Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Mon, 2 Dec 2024 16:20:35 -0800
Subject: [PATCH 19/20] cleanup code, add comments

---
 .../GMRES/Poisson/AMReX_GMRES_Poisson.H       | 61 ++++++++++---------
 .../LinearSolvers/GMRES/Poisson/inputs        |  2 +-
 2 files changed, 32 insertions(+), 31 deletions(-)

diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index b28119dd..6d13bc93 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -9,6 +9,9 @@ namespace amrex {
 
 /**
  * \brief Solve Poisson's equation using amrex GMRES class
+ *
+ * Refer to comments in amrex/Src/LinearSolvers/AMReX_GMRES.H
+ * for details on function implementation requirements
  */
 class GMRESPOISSON
 {
@@ -72,9 +75,6 @@ private:
     DistributionMapping m_dm;
     Geometry m_geom;
     bool m_use_precond;
-
-    // store the original RHS needed for Jacobi preconditioner
-    MultiFab orig_rhs;
 };
 
 GMRESPOISSON::GMRESPOISSON (const BoxArray& ba, const DistributionMapping& dm, const Geometry& geom)
@@ -157,44 +157,48 @@ void GMRESPOISSON::apply (MultiFab& lhs, MultiFab& rhs) const
 
 }
 
-// applies preconditioner to rhs. If there is no preconditioner,
-// this function should do lhs = rhs.
 void GMRESPOISSON::precond (MultiFab& lhs, MultiFab const& rhs) const
 {
     if (m_use_precond) {
+
+        // in the preconditioner we use right-preconditioning to solve
+        // lhs = P^inv(rhs), where P^inv approximates A^inv
+        // Here we use Jacobi iterations to represent P^inv with an initial guess of lhs=0
+
         const GpuArray<Real, AMREX_SPACEDIM> dx = m_geom.CellSizeArray();
+
         amrex::Real fac = 0.;
         for (int d=0; d<AMREX_SPACEDIM; ++d) { fac -= 2./(dx[d]*dx[d]); }
 
-        MultiFab tmp(lhs.boxArray(), lhs.DistributionMap(), 1, 1);
+        MultiFab tmp(m_ba, m_dm, 1, 1);
         auto const& tmp_ma = tmp.const_arrays();
-        auto const& rh_ma = rhs.const_arrays();
-        auto const& x_ma = lhs.arrays();
+        auto const& rhs_ma = rhs.const_arrays();
+        auto const& lhs_ma = lhs.arrays();
+
+        lhs.setVal(0.);
 
         const int niters = 8;
         for (int iter = 0; iter < niters; ++iter) {
-            if (iter == 0) {
-                ParallelFor(lhs, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
-                {
-                    x_ma[b](i,j,k) = rh_ma[b](i,j,k) / fac;
-                });
-            } else {
-                MultiFab::Copy(tmp, lhs, 0, 0, 1, 0);
-                tmp.FillBoundary(m_geom.periodicity());
-                ParallelFor(lhs, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
-                {
-                    auto const& phi = tmp_ma[b];
-                    auto ax = ( phi(i+1,j,k) - 2.*phi(i,j,k) + phi(i-1,j,k) ) / (dx[0]*dx[0])
-                            + ( phi(i,j+1,k) - 2.*phi(i,j,k) + phi(i,j-1,k) ) / (dx[1]*dx[1])
+
+            MultiFab::Copy(tmp, lhs, 0, 0, 1, 0);
+            tmp.FillBoundary(m_geom.periodicity());
+
+            ParallelFor(lhs, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
+            {
+                auto const& tmp_ptr = tmp_ma[b];
+                auto ax = ( tmp_ptr(i+1,j,k) - 2.*tmp_ptr(i,j,k) + tmp_ptr(i-1,j,k) ) / (dx[0]*dx[0])
+                        + ( tmp_ptr(i,j+1,k) - 2.*tmp_ptr(i,j,k) + tmp_ptr(i,j-1,k) ) / (dx[1]*dx[1])
 #if (AMREX_SPACEDIM == 3)
-                            + ( phi(i,j,k+1) - 2.*phi(i,j,k) + phi(i,j,k-1) ) / (dx[2]*dx[2])
+                        + ( tmp_ptr(i,j,k+1) - 2.*tmp_ptr(i,j,k) + tmp_ptr(i,j,k-1) ) / (dx[2]*dx[2])
 #endif
-                        ;
-                    auto res = rh_ma[b](i,j,k) - ax;
-                    x_ma[b](i,j,k) += res / fac * Real(2./3.); // 2/3: weighted jacobi
-                });
-            }
+                    ;
+                auto res = rhs_ma[b](i,j,k) - ax;
+
+                lhs_ma[b](i,j,k) += res / fac * Real(2./3.); // 2/3: weighted jacobi
+
+            });
             Gpu::streamSynchronize();
+
         }
     } else {
         MultiFab::Copy(lhs,rhs,0,0,1,0);
@@ -203,9 +207,6 @@ void GMRESPOISSON::precond (MultiFab& lhs, MultiFab const& rhs) const
 
 void GMRESPOISSON::solve (MultiFab& a_sol, MultiFab const& a_rhs, RT a_tol_rel, RT a_tol_abs)
 {
-    orig_rhs.define(m_ba,m_dm,1,0);
-    MultiFab::Copy(orig_rhs,a_rhs,0,0,1,0);
-
     m_gmres.solve(a_sol, a_rhs, a_tol_rel, a_tol_abs);
 }
 
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs b/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs
index bcf33ac7..140443b2 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/inputs
@@ -1,4 +1,4 @@
 n_cell = 32
 max_grid_size = 16
 
-use_precond = 0
+use_precond = 1

From cc09ed1a3b13c94ee4245c0a4bdb26f945906039 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Mon, 2 Dec 2024 16:22:13 -0800
Subject: [PATCH 20/20] bare-bones documentation pointing to the existence of
 GMRES

---
 Docs/source/LinearSolvers_Tutorial.rst | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/Docs/source/LinearSolvers_Tutorial.rst b/Docs/source/LinearSolvers_Tutorial.rst
index 15819efe..e677494e 100644
--- a/Docs/source/LinearSolvers_Tutorial.rst
+++ b/Docs/source/LinearSolvers_Tutorial.rst
@@ -26,6 +26,8 @@ Step-by-step instructions for configuring AMReX to use HYPRE for this example ar
 
 ``ABecLaplacian_F`` demonstrates how to solve with cell-centered data using the Fortran interfaces.
 
+``GMRES/Poisson`` demonstrates how to solve the Poisson equation using GMRES with Jacobi preconditioning.
+
 ``NodalPoisson`` demonstrates how to set up and solve a variable coefficient Poisson equation
 with the rhs and solution data on nodes.