FEMMaxwellDiffusionSolver.h

// Class FEMMaxwellDifussionSolver
//   Solves the electric diffusion probelm given by curl(curl(E)) + E = f in
//   the domain and n x E = 0 on the boundary.
 
#ifndef IPPL_FEM_MAXWELL_DIFFUSION_SOLVER_H
#define IPPL_FEM_MAXWELL_DIFFUSION_SOLVER_H
 
#include <fstream>
#include <iomanip>
#include <iostream>
 
#include "FEM/FEMQuadratureData.h"
#include "LinearSolvers/PCG.h"
#include "Maxwell.h"
 
namespace ippl {
 
    template <typename T, unsigned Dim, unsigned numElementDOFs>
    struct EvalFunctor {
        const Vector<T, Dim> DPhiInvT;
 
        const T absDetDPhi;
 
        EvalFunctor(Vector<T, Dim> DPhiInvT, T absDetDPhi)
            : DPhiInvT(DPhiInvT)
            , absDetDPhi(absDetDPhi) {}
 
        KOKKOS_FUNCTION auto operator()(
            size_t i, size_t j,
            const QuadratureData<Vector<T, Dim>, Vector<T, Dim>, numElementDOFs>& qd) const {
            T curlTerm = dot(DPhiInvT * qd.deriv_q[j], DPhiInvT * qd.deriv_q[i]).apply();
            T massTerm = dot(qd.val_q[j], qd.val_q[i]).apply();
            return (curlTerm + massTerm) * absDetDPhi;
        }
    };
 
    template <typename FieldType>
    class FEMMaxwellDiffusionSolver : public Maxwell<FieldType, FieldType> {
        constexpr static unsigned Dim = FieldType::dim;
 
        // we call value_type twice as in theory we expect the field to store
        // vector data represented by an ippl::Vector
        using T = typename FieldType::value_type::value_type;
 
        typedef Vector<T, Dim> point_t;
 
    public:
        using Base     = Maxwell<FieldType, FieldType>;
        using MeshType = typename FieldType::Mesh_t;
 
        // PCG (Preconditioned Conjugate Gradient) is the solver algorithm used
        using PCGSolverAlgorithm_t = CG<FEMVector<T>, FEMVector<T>, FEMVector<T>, FEMVector<T>,
                                        FEMVector<T>, FEMVector<T>, FEMVector<T>>;
 
        // FEM Space types
        using ElementType =
            std::conditional_t<Dim == 2, ippl::QuadrilateralElement<T>, ippl::HexahedralElement<T>>;
 
        using QuadratureType = GaussJacobiQuadrature<T, 5, ElementType>;
 
        using NedelecType = NedelecSpace<T, Dim, 1, ElementType, QuadratureType, FieldType>;
 
        // default constructor (compatibility with Alpine)
        FEMMaxwellDiffusionSolver()
            : Base()
            , rhsVector_m(nullptr)
            , refElement_m()
            , quadrature_m(refElement_m, 0.0, 0.0)
            , nedelecSpace_m(*(new MeshType(NDIndex<Dim>(Vector<unsigned, Dim>(0)),
                                            Vector<T, Dim>(0), Vector<T, Dim>(0))),
                             refElement_m, quadrature_m) {}
 
        FEMMaxwellDiffusionSolver(FieldType& lhs, FieldType& rhs,
                                  const FEMVector<point_t>& rhsVector)
            : Base(lhs, lhs, rhs)
            , rhsVector_m(nullptr)
            , refElement_m()
            , quadrature_m(refElement_m, 0.0, 0.0)
            , nedelecSpace_m(rhs.get_mesh(), refElement_m, quadrature_m, rhs.getLayout()) {
            static_assert(std::is_floating_point<T>::value, "Not a floating point type");
            setDefaultParameters();
 
            // Calcualte the rhs, using the Nedelec space
            rhsVector_m =
                std::make_unique<FEMVector<T>>(nedelecSpace_m.evaluateLoadVector(rhsVector));
 
            rhsVector_m->accumulateHalo();
            rhsVector_m->fillHalo();
        }
 
        void setRhs(FieldType& rhs, const FEMVector<point_t>& rhsVector) {
            Base::setRhs(rhs);
 
            // Calcualte the rhs, using the Nedelec space
            rhsVector_m =
                std::make_unique<FEMVector<T>>(nedelecSpace_m.evaluateLoadVector(rhsVector));
 
            rhsVector_m->accumulateHalo();
            rhsVector_m->fillHalo();
        }
 
        void solve() override {
            const Vector<size_t, Dim> zeroNdIndex = Vector<size_t, Dim>(0);
 
            // We can pass the zeroNdIndex here, since the transformation
            // jacobian does not depend on translation
            const auto firstElementVertexPoints =
                nedelecSpace_m.getElementMeshVertexPoints(zeroNdIndex);
 
            // Compute Inverse Transpose Transformation Jacobian ()
            const Vector<T, Dim> DPhiInvT =
                refElement_m.getInverseTransposeTransformationJacobian(firstElementVertexPoints);
 
            // Compute absolute value of the determinant of the transformation
            // jacobian (|det D Phi_K|)
            const T absDetDPhi = Kokkos::abs(
                refElement_m.getDeterminantOfTransformationJacobian(firstElementVertexPoints));
 
            // Create the functor object which stores the function we have to
            // solve for the lhs
            EvalFunctor<T, Dim, NedelecType::numElementDOFs> maxwellDiffusionEval(DPhiInvT,
                                                                                  absDetDPhi);
 
            // The Ax operator
            const auto algoOperator = [maxwellDiffusionEval,
                                       this](FEMVector<T> vector) -> FEMVector<T> {
                vector.fillHalo();
 
                FEMVector<T> return_vector =
                    nedelecSpace_m.evaluateAx(vector, maxwellDiffusionEval);
 
                return_vector.accumulateHalo();
 
                return return_vector;
            };
 
            // setup the CG solver
            pcg_algo_m.setOperator(algoOperator);
 
            // Create the coefficient vector for the solution
            FEMVector<T> lhsVector = rhsVector_m->deepCopy();
 
            // Solve the system using CG
            try {
                pcg_algo_m(lhsVector, *rhsVector_m, this->params_m);
            } catch (IpplException& e) {
                std::string msg = e.where() + ": " + e.what() + "\n";
                Kokkos::abort(msg.c_str());
            }
 
            // store solution.
            lhsVector_m = std::make_unique<FEMVector<T>>(lhsVector);
 
            // set the boundary values to the correct values.
            lhsVector.fillHalo();
        }
 
        int getIterationCount() { return pcg_algo_m.getIterationCount(); }
 
        T getResidue() const { return pcg_algo_m.getResidue(); }
 
        Kokkos::View<point_t*> reconstructToPoints(const Kokkos::View<point_t*>& positions) {
            return this->nedelecSpace_m.reconstructToPoints(positions, *lhsVector_m);
        }
 
        template <typename F>
        T getL2Error(const F& analytic) {
            T error_norm = this->nedelecSpace_m.computeError(*lhsVector_m, analytic);
            return error_norm;
        }
 
    protected:
        PCGSolverAlgorithm_t pcg_algo_m;
 
        virtual void setDefaultParameters() {
            this->params_m.add("max_iterations", 10);
            this->params_m.add("tolerance", (T)1e-13);
        }
 
        std::unique_ptr<FEMVector<T>> rhsVector_m;
 
        std::unique_ptr<FEMVector<T>> lhsVector_m;
 
        ElementType refElement_m;
 
        QuadratureType quadrature_m;
 
        NedelecType nedelecSpace_m;
    };
 
}  // namespace ippl
 
#endif  // IPPL_FEM_MAXWELL_DIFFUSION_SOLVER_H