Solver design refactoring

alpine/ChargedParticles.hpp

@@ -9,10 +9,10 @@
 
 #include "Utility/TypeUtils.h"
 
-#include "Solver/ElectrostaticsCG.h"
 #include "Solver/FFTPeriodicPoissonSolver.h"
 #include "Solver/FFTPoissonSolver.h"
 #include "Solver/P3MSolver.h"
+#include "Solver/PoissonCG.h"
 
 unsigned LoggingPeriod = 1;
 
@@ -54,7 +54,7 @@ using VField_t = Field<Vector_t<T, Dim>, Dim, ViewArgs...>;
 
 // heFFTe does not support 1D FFTs, so we switch to CG in the 1D case
 template <typename T = double, unsigned Dim = 3>
-using CGSolver_t = ippl::ElectrostaticsCG<Field<T, Dim>, Field_t<Dim>>;
+using CGSolver_t = ippl::PoissonCG<Field<T, Dim>, Field_t<Dim>>;
 
 using ippl::detail::ConditionalType, ippl::detail::VariantFromConditionalTypes;
 

src/LinearSolvers/CMakeLists.txt

@@ -0,0 +1,16 @@
+set (_SRCS
+    )
+
+set (_HDRS
+    SolverAlgorithm.h
+    PCG.h
+)
+
+include_DIRECTORIES (
+    ${CMAKE_CURRENT_SOURCE_DIR}
+    )
+
+add_ippl_sources (${_SRCS})
+add_ippl_headers (${_HDRS})
+
+install (FILES ${_HDRS} DESTINATION include/LinearSolvers)

src/Solver/CMakeLists.txt

@@ -2,11 +2,8 @@ set (_SRCS
     )
 
 set (_HDRS
-    SolverAlgorithm.h
-    ElectrostaticsCG.h
-    Electrostatics.h
-    PCG.h
-    Solver.h
+    PoissonCG.h
+    Poisson.h
 )
 
 if (ENABLE_FFT)

src/Solver/FFTPeriodicPoissonSolver.h

@@ -1,6 +1,7 @@
 //
 // Class FFTPeriodicPoissonSolver
-//   Solves periodic electrostatics problems using Fourier transforms
+//   Solves the periodic Poisson problem using Fourier transforms
+//   cf. https://math.mit.edu/~stevenj/fft-deriv.pdf Algorithm 5
 //
 //
 
@@ -11,15 +12,15 @@
 
 #include "Types/ViewTypes.h"
 
-#include "Electrostatics.h"
 #include "FFT/FFT.h"
 #include "FieldLayout/FieldLayout.h"
 #include "Index/NDIndex.h"
+#include "Poisson.h"
 
 namespace ippl {
 
     template <typename FieldLHS, typename FieldRHS>
-    class FFTPeriodicPoissonSolver : public Electrostatics<FieldLHS, FieldRHS> {
+    class FFTPeriodicPoissonSolver : public Poisson<FieldLHS, FieldRHS> {
         constexpr static unsigned Dim = FieldLHS::dim;
         using Trhs                    = typename FieldRHS::value_type;
         using mesh_type               = typename FieldRHS::Mesh_t;
@@ -32,7 +33,7 @@ namespace ippl {
         using Layout_t  = FieldLayout<Dim>;
         using Vector_t  = Vector<Trhs, Dim>;
 
-        using Base = Electrostatics<FieldLHS, FieldRHS>;
+        using Base = Poisson<FieldLHS, FieldRHS>;
         using typename Base::lhs_type, typename Base::rhs_type;
         using scalar_type = typename FieldLHS::Mesh_t::value_type;
         using vector_type = typename FieldLHS::Mesh_t::vector_type;

src/Solver/FFTPeriodicPoissonSolver.hpp

@@ -1,6 +1,7 @@
 //
 // Class FFTPeriodicPoissonSolver
-//   Solves periodic electrostatics problems using Fourier transforms
+//   Solves the periodic Poisson problem using Fourier transforms
+//   cf. https://math.mit.edu/~stevenj/fft-deriv.pdf Algorithm 5
 //
 //
 

src/Solver/FFTPoissonSolver.h

@@ -19,11 +19,11 @@
 #include "Field/Field.h"
 
 #include "Communicate/Archive.h"
-#include "Electrostatics.h"
 #include "FFT/FFT.h"
 #include "Field/HaloCells.h"
 #include "FieldLayout/FieldLayout.h"
 #include "Meshes/UniformCartesian.h"
+#include "Poisson.h"
 
 namespace ippl {
 
@@ -70,15 +70,15 @@ namespace ippl {
     }  // namespace detail
 
     template <typename FieldLHS, typename FieldRHS>
-    class FFTPoissonSolver : public Electrostatics<FieldLHS, FieldRHS> {
+    class FFTPoissonSolver : public Poisson<FieldLHS, FieldRHS> {
         constexpr static unsigned Dim = FieldLHS::dim;
         using Trhs                    = typename FieldRHS::value_type;
         using mesh_type               = typename FieldLHS::Mesh_t;
         using Tg                      = typename FieldLHS::value_type::value_type;
 
     public:
         // type of output
-        using Base = Electrostatics<FieldLHS, FieldRHS>;
+        using Base = Poisson<FieldLHS, FieldRHS>;
 
         // types for LHS and RHS
         using typename Base::lhs_type, typename Base::rhs_type;

src/Solver/Maxwell.h

@@ -0,0 +1,85 @@
+//
+// Class Maxwell
+//   Base class for solvers for Maxwell's equations
+//
+
+#ifndef IPPL_MAXWELL_H
+#define IPPL_MAXWELL_H
+
+#include "Types/Vector.h"
+
+#include "Field/Field.h"
+
+#include "FieldLayout/FieldLayout.h"
+#include "Meshes/UniformCartesian.h"
+
+namespace ippl {
+
+    template <typename EMField, typename SourceField>
+    class Maxwell {
+    public:
+        // typedefs for the different fields and vector fields
+        using typeR = typename SourceField::value_type;
+
+        constexpr static unsigned Dim = EMField::dim;
+        typedef typename EMField::Mesh_t Mesh_t;
+        typedef typename EMField::Centering_t Centering_t;
+        typedef Vector<typeR, Dim> Vector_t;
+        typedef Field<Vector_t, Dim, Mesh_t, Centering_t> VectorSourceField_t;
+
+        // define type for field layout
+        typedef FieldLayout<Dim> FieldLayout_t;
+
+        // type for communication buffers
+        using memory_space = typename SourceField::memory_space;
+        using buffer_type  = Communicate::buffer_type<memory_space>;
+
+        /*!
+         * Default constructor for Maxwell solvers;
+         */
+        Maxwell() {}
+
+        Maxwell(SourceField& charge, VectorSourceField_t& current, EMField& E, EMField& B) {
+            setSources(charge, current);
+            setEMFields(E, B);
+        }
+
+        /*!
+         * Set the problem RHS (charge & current densities)
+         * @param charge Reference to rho field
+         * @param current Reference to J field
+         */
+        virtual void setSources(SourceField& charge, VectorSourceField_t& current) {
+            rhoN_mp = &charge;
+            JN_mp   = &current;
+        }
+
+        /*!
+         * Set the problem LHS (electromagnetic fields)
+         * @param E Reference to electric field
+         * @param B Reference to magnetic field
+         */
+        void setEMFields(EMField& E, EMField& B) {
+            En_mp = &E;
+            Bn_mp = &B;
+        }
+
+        /*!
+         * Solve the electromagnetic problem (Maxwell's eqs)
+         */
+        virtual void solve() = 0;
+
+        virtual ~Maxwell() {}
+
+    protected:
+        // fields containing reference to charge and current
+        SourceField* rhoN_mp       = nullptr;
+        VectorSourceField_t* JN_mp = nullptr;
+
+        // E and B fields
+        EMField* En_mp = nullptr;
+        EMField* Bn_mp = nullptr;
+    };
+}  // namespace ippl
+
+#endif

src/Solver/P3MSolver.h

@@ -17,21 +17,21 @@
 
 #include "Field/Field.h"
 
-#include "Electrostatics.h"
 #include "FFT/FFT.h"
 #include "FieldLayout/FieldLayout.h"
 #include "Meshes/UniformCartesian.h"
+#include "Poisson.h"
 
 namespace ippl {
     template <typename FieldLHS, typename FieldRHS>
-    class P3MSolver : public Electrostatics<FieldLHS, FieldRHS> {
+    class P3MSolver : public Poisson<FieldLHS, FieldRHS> {
         constexpr static unsigned Dim = FieldLHS::dim;
         using Trhs                    = typename FieldRHS::value_type;
         using mesh_type               = typename FieldRHS::Mesh_t;
 
     public:
         // type of output
-        using Base = Electrostatics<FieldLHS, FieldRHS>;
+        using Base = Poisson<FieldLHS, FieldRHS>;
 
         // types for LHS and RHS
         using typename Base::lhs_type, typename Base::rhs_type;