PoissonSolverDealiiMatrixFreeFE.h

Go to the documentation of this file.

/******************************************************************************
 * Copyright (c) 2021.                                                        *
 * The Regents of the University of Michigan and DFT-EFE developers.          *
 *                                                                            *
 * This file is part of the DFT-EFE code.                                     *
 *                                                                            *
 * DFT-EFE is free software: you can redistribute it and/or modify            *
 *   it under the terms of the Lesser GNU General Public License as           *
 *   published by the Free Software Foundation, either version 3 of           *
 *   the License, or (at your option) any later version.                      *
 *                                                                            *
 * DFT-EFE is distributed in the hope that it will be useful, but             *
 *   WITHOUT ANY WARRANTY; without even the implied warranty                  *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.                     *
 *   See the Lesser GNU General Public License for more details.              *
 *                                                                            *
 * You should have received a copy of the GNU Lesser General Public           *
 *   License at the top level of DFT-EFE distribution.  If not, see           *
 *   <https://www.gnu.org/licenses/>.                                         *
 ******************************************************************************/
 
/*
 * @author Avirup Sircar
 */
 
#ifndef dftefePoissonSolverDealiiMatrixFreeFE_h
#define dftefePoissonSolverDealiiMatrixFreeFE_h
 
#include <utils/TypeConfig.h>
#include <utils/MemorySpaceType.h>
#include <basis/FEBasisManager.h>
#include <basis/FEBasisOperations.h>
#include <basis/FEBasisDataStorage.h>
#include <quadrature/QuadratureValuesContainer.h>
#include <basis/DealiiFEEvaluationWrapper.h>
#include <basis/CFEBasisDataStorageDealii.h>
#include <basis/CFEConstraintsLocalDealii.h>
#include <basis/CFEBDSOnTheFlyComputeDealii.h>
#include <vector>
#include <memory>
#include <utils/Profiler.h>
 
namespace dftefe
{
  namespace electrostatics
  {
    template <typename ValueTypeOperator,
              typename ValueTypeOperand,
              utils::MemorySpace memorySpace,
              size_type          dim>
    class PoissonSolverDealiiMatrixFreeFE
    {
    public:
      template <typename T>
      using distributedCPUVec =
        basis::FEEvaluationWrapperBase::distributedCPUVec<T>;
 
      using ValueType =
        linearAlgebra::blasLapack::scalar_type<ValueTypeOperator,
                                               ValueTypeOperand>;
 
    public:
      PoissonSolverDealiiMatrixFreeFE(
        std::shared_ptr<const basis::FEBasisManager<ValueTypeOperand,
                                                    ValueTypeOperator,
                                                    memorySpace,
                                                    dim>> feBasisManagerField,
        std::shared_ptr<
          const basis::FEBasisDataStorage<ValueTypeOperator, memorySpace>>
          feBasisDataStorageStiffnessMatrix,
        const std::map<
          std::string,
          std::shared_ptr<
            const basis::FEBasisDataStorage<ValueTypeOperator, memorySpace>>>
          &feBasisDataStorageRhs,
        const std::map<
          std::string,
          const quadrature::QuadratureValuesContainer<ValueType, memorySpace> &>
          &                                     inpRhs,
        const linearAlgebra::PreconditionerType pcType,
        std::shared_ptr<linearAlgebra::LinAlgOpContext<memorySpace>>
          linAlgOpContext);
 
      PoissonSolverDealiiMatrixFreeFE(
        std::shared_ptr<const basis::FEBasisManager<ValueTypeOperand,
                                                    ValueTypeOperator,
                                                    memorySpace,
                                                    dim>> feBasisManagerField,
        std::shared_ptr<
          const basis::FEBasisDataStorage<ValueTypeOperator, memorySpace>>
          feBasisDataStorageStiffnessMatrix,
        std::shared_ptr<
          const basis::FEBasisDataStorage<ValueTypeOperator, memorySpace>>
          feBasisDataStorageRhs,
        const quadrature::QuadratureValuesContainer<ValueType, memorySpace>
          &                                     inpRhs,
        const linearAlgebra::PreconditionerType pcType,
        std::shared_ptr<linearAlgebra::LinAlgOpContext<memorySpace>>
          linAlgOpContext);
 
      void
      reinit(
        std::shared_ptr<const basis::FEBasisManager<ValueTypeOperand,
                                                    ValueTypeOperator,
                                                    memorySpace,
                                                    dim>> feBasisManagerField,
        const std::map<
          std::string,
          const quadrature::QuadratureValuesContainer<ValueType, memorySpace> &>
          &inpRhs);
 
      void
      reinit(
        std::shared_ptr<const basis::FEBasisManager<ValueTypeOperand,
                                                    ValueTypeOperator,
                                                    memorySpace,
                                                    dim>> feBasisManagerField,
        const quadrature::QuadratureValuesContainer<ValueType, memorySpace>
          &inpRhs);
 
      ~PoissonSolverDealiiMatrixFreeFE() = default;
 
      void
      solve(const double absTolerance, const unsigned int maxNumberIterations);
 
      void
      getSolution(linearAlgebra::MultiVector<ValueType, memorySpace> &solution);
 
      const utils::mpi::MPIComm &
      getMPIComm() const;
 
    private:
      const distributedCPUVec<ValueTypeOperand> &
      getRhs() const;
 
      const distributedCPUVec<ValueType> &
      getInitialGuess() const;
 
      void
      setSolution(const distributedCPUVec<ValueType> &x);
 
      void
      computeRhs(distributedCPUVec<double> &rhs,
                 const std::map<
                   std::string,
                   const quadrature::QuadratureValuesContainer<
                     linearAlgebra::blasLapack::scalar_type<ValueTypeOperator,
                                                            ValueTypeOperand>,
                     memorySpace> &> &inpRhs);
 
      void
      vmult(distributedCPUVec<double> &Ax, distributedCPUVec<double> &x);
 
      void
      precondition_Jacobi(distributedCPUVec<double> &      dst,
                          const distributedCPUVec<double> &src) const;
 
      void
      computeDiagonalA();
 
      void
      AX(const dealii::MatrixFree<dim, double> &      matrixFreeData,
         distributedCPUVec<double> &                  dst,
         const distributedCPUVec<double> &            src,
         const std::pair<unsigned int, unsigned int> &cell_range) const;
 
      void
      CGsolve(const double       absTolerance,
              const unsigned int maxNumberIterations,
              bool               distributeFlag);
 
      size_type d_numComponents;
      std::shared_ptr<
        const basis::
          FEBasisManager<ValueTypeOperand, ValueTypeOperator, memorySpace, dim>>
                                        d_feBasisManagerField;
      linearAlgebra::PreconditionerType d_pcType;
      utils::Profiler                   d_p;
 
 
      std::shared_ptr<
        basis::
          FEBasisManager<ValueTypeOperand, ValueTypeOperator, memorySpace, dim>>
        d_feBasisManagerHomo;
 
      distributedCPUVec<ValueTypeOperator> d_x, d_rhs, d_initial, d_diagonalA;
      std::shared_ptr<dealii::MatrixFree<dim, ValueTypeOperator>>
                                                     d_dealiiMatrixFree;
      std::shared_ptr<const dealii::DoFHandler<dim>> d_dealiiDofHandler;
      const dealii::AffineConstraints<ValueTypeOperand>
        *d_dealiiAffineConstraintMatrix;
      const dealii::AffineConstraints<ValueTypeOperand> *d_constraintsInfo;
      unsigned int                        d_num1DQuadPointsStiffnessMatrix;
      std::map<std::string, unsigned int> d_num1DQuadPointsRhs;
      size_type                           d_feOrder;
      unsigned int                        d_dofHandlerIndex;
 
      std::map<
        std::string,
        std::shared_ptr<
          const basis::FEBasisDataStorage<ValueTypeOperator, memorySpace>>>
        d_feBasisDataStorageRhs;
 
      std::vector<distributedCPUVec<ValueType>> d_nonTensorSructuredQuadeRhs;
 
      unsigned int d_matrixFreeQuadCompStiffnessMatrix;
      std::map<dealii::CellId, unsigned int> d_cellIdToCellIndexMap;
 
      std::vector<dealii::Quadrature<dim>> d_dealiiQuadratureRuleVec;
      dealii::MappingQ1<dim, dim>          d_mappingDealii;
      // dealii::IndexSet d_ghostIndexSet, d_locallyOwnedIndexSet;
      utils::ConditionalOStream pcout;
 
      linearAlgebra::MultiVector<
        linearAlgebra::blasLapack::scalar_type<ValueTypeOperator,
                                               ValueTypeOperand>,
        memorySpace>
        d_scratchMultiVec;
 
    }; // end of class PoissonSolverDealiiMatrixFreeFE
  }    // namespace electrostatics
} // end of namespace dftefe
#include <electrostatics/PoissonSolverDealiiMatrixFreeFE.t.cpp>
#endif // dftefePoissonSolverDealiiMatrixFreeFE_h