Merge branch '111-merge-solver-repo-into-ippl' into 'master'

Resolve "Merge solver repo into Ippl"

Closes #111

See merge request OPAL/Libraries/ippl!119

CMakeLists.txt

@@ -62,6 +62,8 @@ if (ENABLE_FFT)
     message (STATUS "Found Heffte_DIR: ${Heffte_DIR}")
 endif ()
 
+option (ENABLE_SOLVERS "Enable IPPL solvers" OFF)
+
 add_subdirectory (src)
 
 option (ENABLE_IPPLTESTS "IPPL Tests" OFF)
@@ -71,8 +73,11 @@ endif ()
 
 option (ENABLE_ALPINE, "Alpine" OFF)
 if (ENABLE_ALPINE)
-    if (NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/src/Solver/CMakeLists.txt")
-        message (FATAL_ERROR "Cannot enable Alpine since Solver not available!")
+    if (NOT ENABLE_SOLVERS)
+        message (FATAL_ERROR "Cannot enable Alpine since Solver not enabled (-DENABLE_SOLVERS=ON)!")
+    endif ()
+    if (NOT ENABLE_FFT)
+        message (FATAL_ERROR "Cannot enable Alpine since FFT not enabled (-DENABLE_FFT=ON)!")
     endif ()
     message (STATUS "Enabling Alpine")
     add_subdirectory (alpine)

src/CMakeLists.txt

@@ -79,10 +79,10 @@ add_subdirectory (Expression)
 add_subdirectory (Types)
 add_subdirectory (Partition)
 
-# if (EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Solver/CMakeLists.txt")
-#     message (STATUS "Adding Solver subdirectory")
-#     add_subdirectory (Solver)
-# endif ()
+if (ENABLE_SOLVERS)
+    message (STATUS "Adding Solver subdirectory")
+    add_subdirectory (Solver)
+endif ()
 
 if (ENABLE_AMR)
     add_subdirectory(AmrParticle)

src/Solver/CMakeLists.txt

@@ -0,0 +1,30 @@
+set (_SRCS
+    )
+
+set (_HDRS
+    SolverAlgorithm.h
+    ElectrostaticsCG.h
+    Electrostatics.h
+    PCG.h
+    Solver.h
+)
+
+if (ENABLE_FFT)
+    list (APPEND _HDRS
+                 FFTPoissonSolver.h
+                 FFTPoissonSolver.hpp
+                 FFTPeriodicPoissonSolver.h
+                 FFTPeriodicPoissonSolver.hpp
+    )
+endif ()
+
+include_DIRECTORIES (
+    ${CMAKE_CURRENT_SOURCE_DIR}
+    )
+
+add_ippl_sources (${_SRCS})
+add_ippl_headers (${_HDRS})
+
+add_subdirectory(test)
+
+install (FILES ${_HDRS} DESTINATION include/Solver)

src/Solver/Electrostatics.h

@@ -0,0 +1,88 @@
+//
+// Class Electrostatics
+//   Base class for solvers for electrostatics problems
+//
+// Copyright (c) 2021 Alessandro Vinciguerra, ETH Zürich, Zurich, Switzerland
+// All rights reserved
+//
+// This file is part of IPPL.
+//
+// IPPL is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// You should have received a copy of the GNU General Public License
+// along with IPPL. If not, see <https://www.gnu.org/licenses/>.
+//
+
+#ifndef IPPL_ELECTROSTATICS_H
+#define IPPL_ELECTROSTATICS_H
+
+#include "Solver/Solver.h"
+
+namespace ippl {
+
+    template <typename Tlhs, typename Trhs, unsigned Dim,
+              class M=UniformCartesian<double, Dim>,
+              class C=typename M::DefaultCentering>
+    class Electrostatics : public Solver<Tlhs, Trhs, Dim, M, C>
+    {
+    public:
+        using grad_type = Field<Vector<Tlhs, Dim>, Dim, M, C>;
+        using lhs_type = typename Solver<Tlhs, Trhs, Dim, M, C>::lhs_type;
+        using rhs_type = typename Solver<Tlhs, Trhs, Dim, M, C>::rhs_type;
+
+        /*!
+         * Represents the types of fields that should
+         * be output by the solver
+         */
+        enum OutputType {
+            SOL             = 0b01,
+            GRAD            = 0b10,
+            SOL_AND_GRAD    = 0b11
+        };
+
+        /*!
+         * Default constructor for electrostatic solvers;
+         * desired output type defaults to solution only
+         */
+        Electrostatics()
+            : Solver<Tlhs, Trhs, Dim, M, C>()
+            , grad_mp(nullptr)
+        {
+            setDefaultParameters();
+        }
+
+        Electrostatics(lhs_type& lhs, rhs_type& rhs)
+            : Solver<Tlhs, Trhs, Dim, M, C>(lhs, rhs)
+            , grad_mp(nullptr)
+        {
+            setDefaultParameters();
+        }
+
+        /*!
+         * Set the field in which the gradient of the computed potential
+         * should be stored
+         * @param grad Reference to field in which to store the gradient
+         */
+        void setGradient(grad_type& grad) { grad_mp = &grad; }
+
+        /*!
+         * Solve the electrostatics problem described by
+         * -laplace(lhs) = rhs
+         */
+        virtual void solve() = 0;
+
+        virtual ~Electrostatics() { }
+
+    protected:
+        grad_type* grad_mp;
+
+        virtual void setDefaultParameters() override {
+            this->params_m.add("output_type", SOL);
+        }
+    };
+}
+
+#endif

src/Solver/ElectrostaticsCG.h

@@ -0,0 +1,89 @@
+//
+// Class ElectrostaticsCG
+//   Solves electrostatics problems with the CG algorithm
+//
+// Copyright (c) 2021
+// All rights reserved
+//
+// This file is part of IPPL.
+//
+// IPPL is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// You should have received a copy of the GNU General Public License
+// along with IPPL. If not, see <https://www.gnu.org/licenses/>.
+//
+
+#ifndef IPPL_ELECTROSTATICS_CG_H
+#define IPPL_ELECTROSTATICS_CG_H
+
+#include "Electrostatics.h"
+#include "PCG.h"
+
+namespace ippl {
+
+    // Expands to a lambda that acts as a wrapper for a differential operator
+    // fun: the function for which to create the wrapper, such as ippl::laplace
+    // type: the argument type, which should match the LHS type for the solver
+    #define IPPL_SOLVER_OPERATOR_WRAPPER(fun, type)    \
+    [] (type arg) {                                    \
+        return fun(arg);                               \
+    }
+
+    template <typename Tlhs, typename Trhs, unsigned Dim,
+              class M=UniformCartesian<double, Dim>,
+              class C=typename M::DefaultCentering>
+    class ElectrostaticsCG : public Electrostatics<Tlhs, Trhs, Dim, M, C>
+    {
+    public:
+        using lhs_type = typename Solver<Tlhs, Trhs, Dim, M, C>::lhs_type;
+        using rhs_type = typename Solver<Tlhs, Trhs, Dim, M, C>::rhs_type;
+        using OpRet = UnaryMinus<detail::meta_laplace<lhs_type>>;
+        using algo = PCG<Tlhs, Trhs, Dim, OpRet, M, C>;
+        using Base = Electrostatics<Tlhs, Trhs, Dim, M, C>;
+
+        ElectrostaticsCG()
+            : Base()
+        {
+            setDefaultParameters();
+        }
+
+        ElectrostaticsCG(lhs_type& lhs, rhs_type& rhs)
+            : Base(lhs, rhs)
+        {
+            setDefaultParameters();
+        }
+
+        void solve() override {
+            algo_m.setOperator(IPPL_SOLVER_OPERATOR_WRAPPER(-laplace, lhs_type));
+            algo_m(*(this->lhs_mp), *(this->rhs_mp), this->params_m);
+
+            int output = this->params_m.template get<int>("output_type");
+            if (output & Base::GRAD) {
+                *(this->grad_mp) = grad(*(this->lhs_mp));
+            }
+        }
+
+        /*!
+         * Query how many iterations were required to obtain the solution
+         * the last time this solver was used
+         * @return Iteration count of last solve
+         */
+        int getIterationCount() {
+            return algo_m.getIterationCount();
+        }
+
+    protected:
+        algo algo_m = algo();
+
+        virtual void setDefaultParameters() override {
+            this->params_m.add("max_iterations", 1000);
+            this->params_m.add("tolerance", (Tlhs)1e-13);
+        }
+    };
+
+}
+
+#endif

src/Solver/FFTPeriodicPoissonSolver.h

@@ -0,0 +1,111 @@
+//
+// Class FFTPeriodicPoissonSolver
+//   Solves periodic electrostatics problems using Fourier transforms
+//
+// Copyright (c) 2021, Sriramkrishnan Muralikrishnan,
+// Paul Scherrer Institut, Villigen, Switzerland
+// All rights reserved
+//
+// This file is part of IPPL.
+//
+// IPPL is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// You should have received a copy of the GNU General Public License
+// along with IPPL. If not, see <https://www.gnu.org/licenses/>.
+//
+
+#ifndef IPPL_FFT_PERIODIC_POISSON_SOLVER_H
+#define IPPL_FFT_PERIODIC_POISSON_SOLVER_H
+
+#include "FieldLayout/FieldLayout.h"
+#include "Index/NDIndex.h"
+#include "Types/ViewTypes.h"
+#include "FFT/FFT.h"
+#include "Electrostatics.h"
+
+namespace ippl {
+
+    template <typename Tlhs, typename Trhs, unsigned Dim,
+              class M=UniformCartesian<double, Dim>,
+              class C=typename M::DefaultCentering>
+    class FFTPeriodicPoissonSolver : public Electrostatics<Tlhs, Trhs, Dim, M, C>
+    {
+    public:
+        using Field_t = Field<Trhs, Dim>;
+        using FFT_t = FFT<RCTransform, Dim, Trhs>;
+        using Complex_t = Kokkos::complex<Trhs>;
+        using CxField_t = Field<Complex_t, Dim>;
+        using Layout_t = FieldLayout<Dim>;
+        using Vector_t = Vector<Trhs, Dim>;
+
+        using Base = Electrostatics<Tlhs, Trhs, Dim, M, C>;
+        using lhs_type = typename Solver<Tlhs, Trhs, Dim, M, C>::lhs_type;
+        using rhs_type = typename Solver<Tlhs, Trhs, Dim, M, C>::rhs_type;
+
+        FFTPeriodicPoissonSolver()
+            : Base()
+        {
+            setDefaultParameters();
+        }
+
+        FFTPeriodicPoissonSolver(lhs_type& lhs,
+                                 rhs_type& rhs)
+            : Base(lhs, rhs)
+        {
+            setDefaultParameters();
+        }
+
+        //~FFTPeriodicPoissonSolver() {}
+
+        void setRhs(rhs_type& rhs);
+
+        void solve() override;
+
+    private:
+
+        void initialize();
+
+        std::shared_ptr<FFT_t> fft_mp;
+        CxField_t fieldComplex_m;
+        CxField_t tempFieldComplex_m;
+        NDIndex<Dim> domain_m;
+        std::shared_ptr<Layout_t> layoutComplex_mp;
+
+    protected:
+
+        virtual void setDefaultParameters() override {
+            using heffteBackend = typename FFT_t::heffteBackend;
+            heffte::plan_options opts = heffte::default_options<heffteBackend>();
+            this->params_m.add("use_pencils", opts.use_pencils);
+            this->params_m.add("use_reorder", opts.use_reorder);
+            this->params_m.add("use_gpu_aware", opts.use_gpu_aware);
+            this->params_m.add("r2c_direction", 0);
+
+            switch (opts.algorithm) {
+                case heffte::reshape_algorithm::alltoall :
+                    this->params_m.add("comm", a2a);
+                    break;
+                case heffte::reshape_algorithm::alltoallv :
+                    this->params_m.add("comm", a2av);
+                    break;
+                case heffte::reshape_algorithm::p2p :
+                    this->params_m.add("comm", p2p);
+                    break;
+                case heffte::reshape_algorithm::p2p_plined :
+                    this->params_m.add("comm", p2p_pl);
+                    break;
+                default:
+                    throw IpplException("FFTPeriodicPoissonSolver::setDefaultParameters",
+                                       "Unrecognized heffte communication type");
+            }
+
+        }
+
+    };
+}
+
+#include "Solver/FFTPeriodicPoissonSolver.hpp"
+#endif