GMRES Tutorial

ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H

@@ -0,0 +1,193 @@
+#ifndef AMREX_GMRES_MLMG_H_
+#define AMREX_GMRES_MLMG_H_
+#include <AMReX_Config.H>
+
+#include <AMReX_GMRES.H>
+#include <utility>
+
+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.
+ *
+ */
+template <typename MF>
+class GMRESPOISSONT
+{
+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, const Geometry& geom);
+
+    /**
+     * \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, RT a_tol_rel, RT a_tol_abs);
+
+    //! Sets verbosity.
+    void setVerbose (int v) { m_gmres.setVerbose(v); }
+
+    //! 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;
+
+    RT norm2 (MF const& mf) const;
+
+    static void scale (MF& mf, RT scale_factor);
+
+    RT 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, 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);
+
+    //! lhs = L(rhs)
+    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); }
+
+private:
+    GM m_gmres;
+    BoxArray m_ba;
+    DistributionMapping m_dm;
+    Geometry m_geom;
+    bool m_use_precond = false;
+};
+
+template <typename MF>
+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);
+}
+
+template <typename MF>
+auto GMRESPOISSONT<MF>::makeVecRHS () const -> MF
+{
+    return MultiFab(m_ba, m_dm, 1, 0);
+}
+
+template <typename MF>
+auto GMRESPOISSONT<MF>::makeVecLHS () const -> MF
+{
+    return MultiFab(m_ba, m_dm, 1, 1);
+}
+
+template <typename MF>
+auto GMRESPOISSONT<MF>::norm2 (MF const& mf) const -> RT
+{
+    return mf.norm2();
+}
+
+template <typename MF>
+void GMRESPOISSONT<MF>::scale (MF& mf, RT scale_factor)
+{
+    mf.mult(scale_factor);
+}
+
+template <typename MF>
+auto GMRESPOISSONT<MF>::dotProduct (MF const& mf1, MF const& mf2) const -> RT
+{
+    return MultiFab::Dot(mf1,0,mf2,0,1,0);
+}
+
+template <typename MF>
+void GMRESPOISSONT<MF>::setToZero (MF& lhs)
+{
+    lhs.setVal(0.);
+}
+
+template <typename MF>
+void GMRESPOISSONT<MF>::assign (MF& lhs, MF const& rhs)
+{
+    MultiFab::Copy(lhs,rhs,0,0,1,0);
+}
+
+template <typename MF>
+void GMRESPOISSONT<MF>::increment (MF& lhs, MF 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)
+{
+    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
+{
+    // 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>
+void GMRESPOISSONT<MF>::precond (MF& lhs, MF const& rhs) const
+{
+    if (m_use_precond) {
+
+    } else {
+        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)
+{
+    m_gmres.solve(a_sol, a_rhs, a_tol_rel, a_tol_abs);
+}
+
+using GMRESPOISSON = GMRESPOISSONT<MultiFab>;
+
+}
+
+#endif

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

ExampleCodes/LinearSolvers/GMRES/Poisson/Make.package

@@ -0,0 +1,2 @@
+CEXE_sources += main.cpp
+CEXE_headers += AMReX_GMRES_Poisson.H

ExampleCodes/LinearSolvers/GMRES/Poisson/inputs

@@ -0,0 +1,2 @@
+n_cell = 32
+max_grid_size = 16

ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp

@@ -0,0 +1,131 @@
+/*
+ * 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;
+
+    // **********************************
+    // 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);
+    }
+
+    // **********************************
+    // 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();
+
+    // 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 rhs(ba, dm, 1, 0);
+    amrex::MultiFab phi(ba, dm, 1, 1);
+
+    // **********************************
+    // INITIALIZE DATA LOOP
+    // **********************************
+
+    // loop over boxes
+    for (amrex::MFIter mfi(rhs); mfi.isValid(); ++mfi)
+    {
+        const amrex::Box& bx = mfi.validbox();
+
+        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)
+        {
+            // **********************************
+            // 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);
+        });
+    }
+
+    WriteSingleLevelPlotfile("rhs", rhs, {"rhs"}, geom, 0., 0);
+
+    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.);
+
+    WriteSingleLevelPlotfile("phi", phi, {"phi"}, geom, 0., 0);
+
+    }
+    amrex::Finalize();
+    return 0;
+}
+
+