Merge pull request #100 from dftfeDevelopers/BlockedChebyImpl

Blocked cheby impl

Author: AvirupSircar

analysis/classicalEnrichmentComparison/PSPnonlocal/KSDFTClassicalUniformQuad/TestKohnShamDft.cpp

@@ -351,13 +351,16 @@ int main(int argc, char** argv)
   const size_type dim = 3;
 
   double atomPartitionTolerance = readParameter<double>(parameterInputFileName, "atomPartitionTolerance", rootCout);
-  unsigned int num1DGaussSizeElec = readParameter<unsigned int>(parameterInputFileName, "num1DGaussSizeElec", rootCout);
+  unsigned int num1DGaussSubdividedSizeElec = readParameter<unsigned int>(parameterInputFileName, "num1DGaussSubdividedSizeElec", rootCout);
   unsigned int gaussSubdividedCopiesElec = readParameter<unsigned int>(parameterInputFileName, "gaussSubdividedCopiesElec", rootCout);
-  unsigned int num1DGaussSizeEigen = readParameter<unsigned int>(parameterInputFileName, "num1DGaussSizeEigen", rootCout);
+  unsigned int num1DGaussSubdividedSizeEigen = readParameter<unsigned int>(parameterInputFileName, "num1DGaussSubdividedSizeEigen", rootCout);
   unsigned int gaussSubdividedCopiesEigen = readParameter<unsigned int>(parameterInputFileName, "gaussSubdividedCopiesEigen", rootCout);
   bool isNumericalNuclearSolve = readParameter<bool>(parameterInputFileName, "isNumericalNuclearSolve", rootCout);
   bool isDeltaRhoPoissonSolve = readParameter<bool>(parameterInputFileName, "isDeltaRhoPoissonSolve", rootCout);
 
+  unsigned int num1DGaussSubdividedSizeGrad = readParameter<unsigned int>(parameterInputFileName, "num1DGaussSubdividedSizeGrad", rootCout);
+  unsigned int gaussSubdividedCopiesGrad = readParameter<unsigned int>(parameterInputFileName, "gaussSubdividedCopiesGrad", rootCout);
+
   unsigned int num1DGaussSubdividedSizeNonLocOperator = 14;
   unsigned int gaussSubdividedCopiesNonLocOperator = 1;
 
@@ -531,10 +534,10 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
                                                   linAlgOpContext);
 
     std::shared_ptr<quadrature::QuadratureRule> gaussSubdivQuadRuleElec =
-      std::make_shared<quadrature::QuadratureRuleGaussIterated>(dim, num1DGaussSizeElec, gaussSubdividedCopiesElec);
+      std::make_shared<quadrature::QuadratureRuleGaussIterated>(dim, num1DGaussSubdividedSizeElec, gaussSubdividedCopiesElec);
 
     std::shared_ptr<quadrature::QuadratureRule> gaussSubdivQuadRuleEigen =
-      std::make_shared<quadrature::QuadratureRuleGaussIterated>(dim, num1DGaussSizeEigen, gaussSubdividedCopiesEigen);
+      std::make_shared<quadrature::QuadratureRuleGaussIterated>(dim, num1DGaussSubdividedSizeEigen, gaussSubdividedCopiesEigen);
 
     quadrature::QuadratureRuleAttributes quadAttrGaussSubdivided(quadrature::QuadratureFamily::GAUSS_SUBDIVIDED,true);
 
@@ -625,22 +628,30 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
       (basisDofHandlerWaveFn, quadAttrGaussSubdivided, basisAttrMap, ksdft::KSDFTDefaults::CELL_BATCH_SIZE_GRAD_EVAL, *linAlgOpContext);
   feBDAtomCenterNonLocalOperator->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerGaussSubdividedAtomNonLocOp, basisAttrMap);
 
+  std::shared_ptr<quadrature::QuadratureRule> gaussSubdivQuadRuleGrad =
+    std::make_shared<quadrature::QuadratureRuleGaussIterated>(dim, num1DGaussSubdividedSizeGrad, gaussSubdividedCopiesGrad);
+
+  std::shared_ptr<quadrature::QuadratureRuleContainer> quadRuleContainerAdaptiveGrad =  
+    std::make_shared<quadrature::QuadratureRuleContainer>
+    (quadAttrGaussSubdivided, 
+    gaussSubdivQuadRuleGrad, 
+    triangulationBase, 
+    *cellMapping); 
+
   basisAttrMap[basis::BasisStorageAttributes::StoreValues] = false;
   basisAttrMap[basis::BasisStorageAttributes::StoreGradient] = true;
   basisAttrMap[basis::BasisStorageAttributes::StoreHessian] = false;
   basisAttrMap[basis::BasisStorageAttributes::StoreOverlap] = false;
   basisAttrMap[basis::BasisStorageAttributes::StoreGradNiGradNj] = false;
   basisAttrMap[basis::BasisStorageAttributes::StoreJxW] = true;
 
-  quadrature::QuadratureRuleAttributes quadAttrEigen(quadrature::QuadratureFamily::GAUSS,true,feOrderEigen+1);
-
   // Set up the FE Basis Data Storage
   std::shared_ptr<basis::FEBasisDataStorage<double, Host>> feBDHamStiffnessMatrix =
     std::make_shared<basis::CFEBDSOnTheFlyComputeDealii<double, double, Host,dim>>
-    (basisDofHandlerWaveFn, quadAttrEigen, basisAttrMap, ksdft::KSDFTDefaults::CELL_BATCH_SIZE_GRAD_EVAL, *linAlgOpContext);
+    (basisDofHandlerWaveFn, quadAttrGaussSubdivided, basisAttrMap, ksdft::KSDFTDefaults::CELL_BATCH_SIZE_GRAD_EVAL, *linAlgOpContext);
 
   // evaluate basis data
-  feBDHamStiffnessMatrix->evaluateBasisData(quadAttrEigen, basisAttrMap);
+  feBDHamStiffnessMatrix->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerAdaptiveGrad, basisAttrMap);
 
   std::shared_ptr<const basis::FEBasisDataStorage<double,Host>> feBDKineticHamiltonian =  feBDHamStiffnessMatrix;
   std::shared_ptr<const basis::FEBasisDataStorage<double, Host>> feBDEXCHamiltonian = feBDElectrostaticsHamiltonian;

analysis/classicalEnrichmentComparison/PSPnonlocal/KSDFTOrthoEFEUniformQuad/TestKohnShamDft.cpp

@@ -490,6 +490,8 @@ int main(int argc, char** argv)
   adaptiveMesh.createMesh(*triangulationBase); 
     p.registerEnd("Create Mesh");
 
+  utils::printCurrentMemoryUsage(comm, "Create Mesh");
+
     p.registerStart("Quadrature Rule Creation");
     //quadrature::QuadratureRuleAttributes quadAttrAdaptive(quadrature::QuadratureFamily::ADAPTIVE,false);
 
@@ -580,6 +582,7 @@ int main(int argc, char** argv)
   rootCout << "Number of quadrature points in wave function adaptive quadrature: "<<nQuad<<"\n";
 
   p.registerEnd("Quadrature Rule Creation");
+    utils::printCurrentMemoryUsage(comm, "Quadrature Rule Creation");
   p.registerStart("Ortho EFE basis manager creation");
 
   std::vector<double> smearedChargeRadiusVec(atomCoordinatesVec.size(),rc);
@@ -658,19 +661,6 @@ int main(int argc, char** argv)
   std::shared_ptr<basis::EnrichmentClassicalInterfaceSpherical
                           <double, Host, dim>>
         enrichClassIntfceTotalPot = nullptr;
-  if (!isDeltaRhoPoissonSolve)
-    // Create the enrichmentClassicalInterface object for vtotal
-    enrichClassIntfceTotalPot = std::make_shared<basis::EnrichmentClassicalInterfaceSpherical
-                        <double, Host, dim>>
-                        (cfeBasisDataStorageGLLElec,
-                        cfeBasisDataStorageGaussSubdividedElec,
-                        atomSphericalDataContainer,
-                        atomPartitionTolerance,
-                        atomSymbolVec,
-                        atomCoordinatesVec,
-                        "vtotal",
-                        linAlgOpContext,
-                        comm);
 
     // Create the enrichmentClassicalInterface object for wavefn
   std::shared_ptr<basis::EnrichmentClassicalInterfaceSpherical
@@ -691,9 +681,22 @@ int main(int argc, char** argv)
 
   std::shared_ptr<basis::FEBasisDofHandler<double, Host,dim>> basisDofHandlerTotalPot = nullptr;
   if (!isDeltaRhoPoissonSolve)
+  {
+    enrichClassIntfceTotalPot = std::make_shared<basis::EnrichmentClassicalInterfaceSpherical
+                        <double, Host, dim>>
+                        (cfeBasisDataStorageGLLElec,
+                        cfeBasisDataStorageGaussSubdividedElec,
+                        atomSphericalDataContainer,
+                        atomPartitionTolerance,
+                        atomSymbolVec,
+                        atomCoordinatesVec,
+                        "vtotal",
+                        linAlgOpContext,
+                        comm);
      basisDofHandlerTotalPot =  
       std::make_shared<basis::EFEBasisDofHandlerDealii<double, double,Host,dim>>(
         enrichClassIntfceTotalPot, comm);
+  }
   else
     basisDofHandlerTotalPot = cfeBasisDofHandlerElec;
 
@@ -719,6 +722,7 @@ int main(int argc, char** argv)
   */
 
   p.registerEnd("Ortho EFE basis manager creation");
+  utils::printCurrentMemoryUsage(comm, "Ortho EFE basis manager creation");
 
   rootCout << "Total Number of dofs electrostatics: " << basisDofHandlerTotalPot->nGlobalNodes() << "\n";
   rootCout << "Total Number of dofs eigensolve: " << basisDofHandlerWaveFn->nGlobalNodes() << "\n";
@@ -753,6 +757,7 @@ int main(int argc, char** argv)
     feBDTotalChargeStiffnessMatrix->evaluateBasisData(quadAttrGaussElectro, basisAttrMap);
 
   p.registerEnd("Electrostatics basis grad datastorage eval");
+  utils::printCurrentMemoryUsage(comm, "Electrostatics basis grad datastorage eval");
   p.registerStart("Electrostatics basis bsmear datastorage eval");
 
   basisAttrMap[basis::BasisStorageAttributes::StoreValues] = true;
@@ -775,6 +780,7 @@ int main(int argc, char** argv)
   feBDNucChargeRhs->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerGaussSubdividedElec, basisAttrMap);
 
   p.registerEnd("Electrostatics basis bsmear datastorage eval");
+  utils::printCurrentMemoryUsage(comm, "Electrostatics basis bsmear datastorage eval");
   p.registerStart("Electrostatics basis rho datastorage eval");
 
   std::shared_ptr<basis::FEBasisDataStorage<double, Host>> feBDElecChargeRhs = nullptr;
@@ -790,7 +796,7 @@ int main(int argc, char** argv)
   feBDElecChargeRhs->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerAdaptiveOrbital, basisAttrMap);
 
   p.registerEnd("Electrostatics basis rho datastorage eval");
-
+  utils::printCurrentMemoryUsage(comm, "Electrostatics basis rho datastorage eval");
   p.registerStart("Orbital basis datastorage eval");
 
   basisAttrMap[basis::BasisStorageAttributes::StoreValues] = true;
@@ -819,6 +825,11 @@ int main(int argc, char** argv)
 
   feBDAtomCenterNonLocalOperator->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerAdaptiveAtomCenterNonLocalOperator, basisAttrMap);
 
+  p.registerEnd("Orbital basis datastorage eval");
+  utils::printCurrentMemoryUsage(comm, "Orbital basis datastorage eval");
+  
+  p.registerStart("Orbital Grad basis datastorage eval");
+
   basisAttrMap[basis::BasisStorageAttributes::StoreValues] = false;
   basisAttrMap[basis::BasisStorageAttributes::StoreGradient] = true;
   basisAttrMap[basis::BasisStorageAttributes::StoreHessian] = false;
@@ -830,16 +841,14 @@ int main(int argc, char** argv)
     std::make_shared<basis::EFEBDSOnTheFlyComputeDealii<double, double, Host,dim>>
       (basisDofHandlerWaveFn, quadAttrGaussSubdivided, basisAttrMap, ksdft::KSDFTDefaults::CELL_BATCH_SIZE_GRAD_EVAL, *linAlgOpContext);
 
-  p.registerEnd("Orbital basis datastorage eval");
-  p.registerStart("Orbital Grad basis datastorage eval");
-
    efeBasisDataAdaptiveGrad->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerAdaptiveGrad, basisAttrMap);
 
     std::shared_ptr<const basis::FEBasisDataStorage<double, Host>> feBDElectrostaticsHamiltonian = efeBasisDataAdaptiveOrbital;
     std::shared_ptr<const basis::FEBasisDataStorage<double,Host>> feBDKineticHamiltonian =  efeBasisDataAdaptiveGrad;
     std::shared_ptr<const basis::FEBasisDataStorage<double, Host>> feBDEXCHamiltonian = efeBasisDataAdaptiveOrbital;
 
     p.registerEnd("Orbital Grad basis datastorage eval");
+    utils::printCurrentMemoryUsage(comm, "Orbital Grad basis datastorage eval");
     p.print();
 
     std::shared_ptr<const utils::ScalarSpatialFunctionReal>
@@ -934,6 +943,7 @@ int main(int argc, char** argv)
                                                     linAlgOpContext);    
 
   p.registerEnd("Hamiltonian Basis overlap inverse eval");
+  utils::printCurrentMemoryUsage(comm, "Hamiltonian Basis overlap and inv");
   p.registerStart("Kohn Sham DFT Class Init");
   std::shared_ptr<ksdft::KohnShamDFT<double,
                                         double,
@@ -1048,6 +1058,7 @@ int main(int argc, char** argv)
     utils::throwException(false, "Option not there for KohnShamDFT class creation.");
   }
   p.registerEnd("Kohn Sham DFT Class Init"); 
+  utils::printCurrentMemoryUsage(comm, "Kohn Sham DFT Class Init");
   p.print();
 
   pTot.registerEnd("Initilization");