From 19d813a6c8644a5b7ef828a91ec4c68512ca10b6 Mon Sep 17 00:00:00 2001 From: Avirup Sircar Date: Sun, 16 Feb 2025 07:23:37 -0800 Subject: [PATCH] fixed compilation error in DealiiMatrixFree PoissonSolve --- src/basis/DealiiFEEvaluationWrapper.cpp | 78 +++++------- src/basis/DealiiFEEvaluationWrapper.h | 27 ++-- .../PoissonSolverDealiiMatrixFreeFE.h | 14 ++- .../PoissonSolverDealiiMatrixFreeFE.t.cpp | 118 +++++++++++------- 4 files changed, 134 insertions(+), 103 deletions(-) diff --git a/src/basis/DealiiFEEvaluationWrapper.cpp b/src/basis/DealiiFEEvaluationWrapper.cpp index 54fc0774..414007cf 100644 --- a/src/basis/DealiiFEEvaluationWrapper.cpp +++ b/src/basis/DealiiFEEvaluationWrapper.cpp @@ -20,15 +20,13 @@ ******************************************************************************/ /* - * @author DFTFE + * @author dftefe */ -#include +#include - - -namespace dftfe +namespace dftefe { - namespace dftUtils + namespace basis { template - // void FEEvaluationWrapperDerived:: - // submitInterpolatedGradientsAndMultiply( - // dealii::AlignedVector< - // dealii::Tensor<1, 3, dealii::VectorizedArray>> &alpha) - // { - // AssertThrow(false, dftUtils::ExcNotImplementedYet()); - // // if (positiveFlag) - // // for (unsigned int q = 0; q < d_dealiiFEEvaluation->n_q_points; - // ++q) - // // { - // // d_dealiiFEEvaluation->submit_gradient( - // // alpha[q] * d_dealiiFEEvaluation->get_gradient(q), q); - // // } - // } - - template @@ -463,14 +441,16 @@ namespace dftfe { d_dealiiFEEvaluation->distribute_local_to_global(tempvec); } -#ifdef DFTFE_MINIMAL_COMPILE -# define RANGE_FEORDER ((1)(2)(3)(4)(5)(6)(7)) -# define RANGE_QUADRATURE ((2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)) -#else -# define RANGE_FEORDER ((1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)) -# define RANGE_QUADRATURE \ - ((2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(18)) -#endif +// #ifdef DFTEFE_MINIMAL_COMPILE +// # define RANGE_FEORDER ((1)(2)(3)(4)(5)(6)(7)) +// # define RANGE_QUADRATURE ((2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)) +// #else +#define RANGE_FEORDER \ + ((3)(4)(5)(6)(7)) //(1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14) +#define RANGE_QUADRATURE \ + ((10)(12)(14)(16)(18)( \ + 20)) //(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(18)(19)(20) +// #endif #define MACRO(r, p) \ template class FEEvaluationWrapperDerived + const FEEvaluationWrapperBase & + DealiiFEEvaluationWrapper::getFEEvaluationWrapperBase() + const + { + return *d_feEvaluationBase; + } + template class DealiiFEEvaluationWrapper<1>; // template class DealiiFEEvaluationWrapper<3>; - } // End of namespace dftUtils -} // End of namespace dftfe + } // End of namespace basis +} // End of namespace dftefe diff --git a/src/basis/DealiiFEEvaluationWrapper.h b/src/basis/DealiiFEEvaluationWrapper.h index fcc3154d..2bf97b9f 100644 --- a/src/basis/DealiiFEEvaluationWrapper.h +++ b/src/basis/DealiiFEEvaluationWrapper.h @@ -20,28 +20,31 @@ ******************************************************************************/ /* - * @author DFTFE + * @author dftefe */ #ifndef FEEvaluationWrapper_h #define FEEvaluationWrapper_h -#include -#include #include #include #include -#include -#include -#include -namespace dftfe +#include +#include + +namespace dftefe { - namespace dftUtils + namespace basis { class FEEvaluationWrapperBase { public: + template + using distributedCPUVec = + dealii::LinearAlgebra::distributed::Vector; + /** * @brief Returns the total number of quadrature points in all 3 directions */ @@ -285,6 +288,10 @@ namespace dftfe const unsigned int matrixFreeVectorComponent, const unsigned int matrixFreeQuadratureComponent); + const FEEvaluationWrapperBase & + getFEEvaluationWrapperBase() const; + + private: unsigned int d_feDegree; unsigned int d_num1dQuad; unsigned int d_matrixFreeVectorComponent; @@ -297,9 +304,9 @@ namespace dftfe }; // end of DealiiFEEvaluationWrapper - } // end of namespace dftUtils + } // end of namespace basis -} // end of namespace dftfe +} // end of namespace dftefe #endif // FEEvaluationWrapper_h diff --git a/src/electrostatics/PoissonSolverDealiiMatrixFreeFE.h b/src/electrostatics/PoissonSolverDealiiMatrixFreeFE.h index 821cfd71..3522bfca 100644 --- a/src/electrostatics/PoissonSolverDealiiMatrixFreeFE.h +++ b/src/electrostatics/PoissonSolverDealiiMatrixFreeFE.h @@ -32,6 +32,10 @@ #include #include #include +#include +#include +#include +#include #include #include #include @@ -77,8 +81,7 @@ namespace dftefe template using distributedCPUVec = - dealii::LinearAlgebra::distributed::Vector; + basis::FEEvaluationWrapperBase::distributedCPUVec; using ValueType = linearAlgebra::blasLapack::scalar_type &solution); @@ -216,9 +219,10 @@ namespace dftefe d_dealiiMatrixFree; std::shared_ptr> d_dealiiDofHandler; dealii::AffineConstraints - * d_dealiiAffineConstraintMatrix; + *d_dealiiAffineConstraintMatrix; + dealii::AffineConstraints *d_constraintsInfo; unsigned int d_num1DQuadPointsStiffnessMatrix; - std::map d_num1DPointsRhs; + std::map d_num1DQuadPointsRhs; size_type d_feOrder; unsigned int d_dofHandlerIndex; diff --git a/src/electrostatics/PoissonSolverDealiiMatrixFreeFE.t.cpp b/src/electrostatics/PoissonSolverDealiiMatrixFreeFE.t.cpp index 7d0b828a..523c1381 100644 --- a/src/electrostatics/PoissonSolverDealiiMatrixFreeFE.t.cpp +++ b/src/electrostatics/PoissonSolverDealiiMatrixFreeFE.t.cpp @@ -31,30 +31,41 @@ namespace dftefe namespace PoissonSolverDealiiMatrixFreeFEInternal { // solve + template void - CGsolve(dealiiLinearSolverProblem &problem, - const double absTolerance, - const unsigned int maxNumberIterations, - bool distributeFlag) + CGsolve(PoissonSolverDealiiMatrixFreeFE &problem, + const double absTolerance, + const unsigned int maxNumberIterations, + bool distributeFlag) { - int this_process; - MPI_Comm_rank(mpi_communicator, &this_process); - MPI_Barrier(mpi_communicator); - double start_time = MPI_Wtime(); + utils::mpi::MPIBarrier(problem.getMPIComm()); + double start_time = utils::mpi::MPIWtime(); double time; // compute RHS - distributedCPUVec rhs; + basis::FEEvaluationWrapperBase::distributedCPUVec rhs, gvec, + dvec, hvec; problem.getRhs(rhs); - MPI_Barrier(mpi_communicator); - time = MPI_Wtime(); + MPI_Barrier(problem.getMPIComm()); + time = utils::mpi::MPIWtime(); + + int rank; + utils::mpi::MPICommRank(problem.getMPIComm(), &rank); + utils::ConditionalOStream pcout(std::cout, rank == 0); pcout << "Time for compute rhs: " << time - start_time << std::endl; bool conv = false; // false : converged; true : converged - distributedCPUVec x = problem.getInitialGuess(); + basis::FEEvaluationWrapperBase::distributedCPUVec x = + problem.getInitialGuess(); double res = 0.0, initial_res = 0.0; int it = 0; @@ -164,14 +175,18 @@ namespace dftefe << " , current abs. residual: " << res << " , nsteps: " << it << " , abs. tolerance criterion: " << absTolerance << "\n\n"; - MPI_Barrier(mpi_communicator); - time = MPI_Wtime() - time; + utils::mpi::MPIBarrier(problem.getMPIComm()); + time = utils::mpi::MPIWtime() - time; pcout << "Time for Poisson/Helmholtz problem CG iterations: " << time << std::endl; } - + template + void getDealiiQuadratureRule( std::shared_ptr< const basis::FEBasisDataStorage> @@ -267,11 +282,12 @@ namespace dftefe // Check wether the dofhandler and constrints come from classical basis or // not - const basis::FEBasisDofHandler - &cfeBasisDofHandlerDealii = std::dynamic_cast< - const basis:: - CFEBasisDofHandlerDealii &>( - feBasisManagerField->getBasisDofHandler()); + const basis::CFEBasisDofHandlerDealii + &cfeBasisDofHandlerDealii = + dynamic_cast &>( + feBasisManagerField->getBasisDofHandler()); utils::throwException( cfeBasisDofHandlerDealii != nullptr, "Could not cast BasisDofHandler to CFEBasisDofHandlerDealii " @@ -293,9 +309,9 @@ namespace dftefe const basis::CFEConstraintsLocalDealii &cfeConstraintsLocalDealii = - std::dynamic_cast< - const basis:: - CFEConstraintsLocalDealii &>( + dynamic_cast &>( feBasisManagerHomo->getConstraints()); utils::throwException( cfeConstraintsLocalDealii != nullptr, @@ -445,16 +461,16 @@ namespace dftefe const basis::CFEConstraintsLocalDealii &cfeConstraintsLocalDealii = - std::dynamic_cast< - const basis:: - CFEConstraintsLocalDealii &>( + dynamic_cast &>( d_feBasisManagerField->getConstraints()); utils::throwException( cfeConstraintsLocalDealii != nullptr, "Could not cast ConstraintsLocal to CFEConstraintsLocalDealii " "in PoissonSolverDealiiMatrixFreeFE."); - *d_constraintsInfo = &cfeConstraintsLocalDealii.getAffineConstraints(); + d_constraintsInfo = &cfeConstraintsLocalDealii.getAffineConstraints(); // Compute RHS computeRhs(d_rhs, inpRhs); @@ -498,7 +514,7 @@ namespace dftefe ValueTypeOperator, ValueTypeOperand, memorySpace, - dim>::setSolution(distributedCPUVec &x) + dim>::setSolution(const distributedCPUVec &x) { d_x = x; } @@ -521,7 +537,7 @@ namespace dftefe for (size_type i = 0; i < solution.locallyOwnedSize(); i++) { - solution.data()[j] = d_x.data()[j]; + solution.data()[i] = d_x.data()[i]; } solution.updateGhostValues(); @@ -538,7 +554,7 @@ namespace dftefe typename ValueTypeOperand, utils::MemorySpace memorySpace, size_type dim> - const distributedCPUVec & + const basis::FEEvaluationWrapperBase::distributedCPUVec & PoissonSolverDealiiMatrixFreeFE - const distributedCPUVec & + const basis::FEEvaluationWrapperBase::distributedCPUVec< + linearAlgebra::blasLapack::scalar_type> & PoissonSolverDealiiMatrixFreeFE( + *this, absTolerance, maxNumberIterations, true); } // Ax @@ -602,13 +619,16 @@ namespace dftefe // d_dofHandlerIndex, // d_matrixFreeQuadCompStiffnessMatrix); - basis::DealiiFEEvaluationWrapper<1> fe_eval( + basis::DealiiFEEvaluationWrapper<1> fe_eval_wrap( d_feOrder, d_num1DQuadPointsStiffnessMatrix, *d_dealiiMatrixFree, d_dofHandlerIndex, d_matrixFreeQuadCompStiffnessMatrix); + basis::FEEvaluationWrapperBase fe_eval = + fe_eval_wrap.getFEEvaluationWrapperBase(); + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) { @@ -753,7 +773,7 @@ namespace dftefe ValueTypeOperand>, memorySpace> &> &inpRhs) { - dealii::DoFHandler::active_cell_iterator subCellPtr; + typename dealii::DoFHandler::active_cell_iterator subCellPtr; rhs.reinit(d_x); rhs = 0; // if (d_isStoreSmearedChargeRhs) @@ -783,19 +803,23 @@ namespace dftefe // d_dofHandlerIndex, // d_matrixFreeQuadCompStiffnessMatrix); - basis::DealiiFEEvaluationWrapper<1> fe_eval( + basis::DealiiFEEvaluationWrapper<1> fe_eval_wrap( d_feOrder, d_num1DQuadPointsStiffnessMatrix, *d_dealiiMatrixFree, d_dofHandlerIndex, d_matrixFreeQuadCompStiffnessMatrix); + basis::FEEvaluationWrapperBase fe_eval = + fe_eval_wrap.getFEEvaluationWrapperBase(); + const dealii::Quadrature &quadratureRuleAxTemp = d_dealiiMatrixFree->get_quadrature(d_matrixFreeQuadCompStiffnessMatrix); int isPerformStaticCondensation = (tempvec.linfty_norm() > 1e-10) ? 1 : 0; - MPI_Bcast(&isPerformStaticCondensation, 1, MPI_INT, 0, mpi_communicator); + utils::mpi::MPIBcast( + &isPerformStaticCondensation, 1, utils::mpi::MPIInt, 0, getMPIComm()); if (isPerformStaticCondensation == 1) { @@ -831,15 +855,19 @@ namespace dftefe // d_dofHandlerIndex, // matrixFreeQuadratureComponentRhs); - basis::DealiiFEEvaluationWrapper<1> fe_eval_density( + basis::DealiiFEEvaluationWrapper<1> fe_eval_density_wrap( d_feOrder, - d_num1DPointsRhs[iter->first], + d_num1DQuadPointsRhs[iter->first], *d_dealiiMatrixFree, d_dofHandlerIndex, matrixFreeQuadratureComponentRhs); + basis::FEEvaluationWrapperBase fe_eval_density = + fe_eval_density_wrap.getFEEvaluationWrapperBase(); + dealii::AlignedVector> rhoQuads( - fe_eval_density.n_q_points, dealii::make_vectorized_array(0.0)); + fe_eval_density.totalNumberofQuadraturePoints(), + dealii::make_vectorized_array(0.0)); for (unsigned int macrocell = 0; macrocell < d_dealiiMatrixFree->n_cell_batches(); ++macrocell) @@ -863,9 +891,11 @@ namespace dftefe d_basisOperationsPtr->cellIndex(subCellId); const double *tempVec = inpRhs[iter->first].data() + - cellIndex * fe_eval_density.n_q_points; + cellIndex * fe_eval_density.totalNumberofQuadraturePoints(); - for (unsigned int q = 0; q < fe_eval_density.n_q_points; ++q) + for (unsigned int q = 0; + q < fe_eval_density.totalNumberofQuadraturePoints(); + ++q) rhoQuads[q][iSubCell] = tempVec[q]; }