From 9c0b1c33c4ab30788c96cec9624ce69b5d875fe3 Mon Sep 17 00:00:00 2001 From: Avirup Sircar Date: Fri, 6 Dec 2024 16:52:28 -0800 Subject: [PATCH] Minor modifications in KSDFT/Defaults.h --- .../TestKohnShamDft.cpp | 38 ++++++++++++++----- .../TestKohnShamDft.cpp | 34 +++++++++++++++-- .../PoissonLinearSolverFunctionFE.t.cpp | 2 + src/ksdft/Defaults.cpp | 6 ++- src/ksdft/ElectrostaticLocalFE.h | 2 + src/ksdft/ElectrostaticLocalFE.t.cpp | 10 ++--- 6 files changed, 71 insertions(+), 21 deletions(-) diff --git a/analysis/classicalEnrichmentComparison/PSP/KSDFTClassicalUniformQuad/TestKohnShamDft.cpp b/analysis/classicalEnrichmentComparison/PSP/KSDFTClassicalUniformQuad/TestKohnShamDft.cpp index 8fe2bf24..104b77bc 100644 --- a/analysis/classicalEnrichmentComparison/PSP/KSDFTClassicalUniformQuad/TestKohnShamDft.cpp +++ b/analysis/classicalEnrichmentComparison/PSP/KSDFTClassicalUniformQuad/TestKohnShamDft.cpp @@ -733,12 +733,12 @@ int main(int argc, char** argv) basisAttrMap[basis::BasisStorageAttributes::StoreJxW] = true; // Set up the FE Basis Data Storage - std::shared_ptr> feBasisDataElec = + std::shared_ptr> feBDTotalChargeStiffnessMatrix = std::make_shared> (basisDofHandlerTotalPot, quadAttrElec, basisAttrMap, ksdft::KSDFTDefaults::CELL_BATCH_SIZE_GRAD_EVAL, *linAlgOpContext); // evaluate basis data - feBasisDataElec->evaluateBasisData(quadAttrElec, basisAttrMap); + feBDTotalChargeStiffnessMatrix->evaluateBasisData(quadAttrElec, basisAttrMap); std::shared_ptr zeroFunction = std::make_shared @@ -804,8 +804,6 @@ std::shared_ptr> feBDTotalChargeStiffnessMatrix = feBasisDataElec; std::shared_ptr gaussSubdivQuadRuleElec = std::make_shared(dim, num1DGaussSizeElec, gaussSubdividedCopiesElec); @@ -870,7 +868,7 @@ std::shared_ptrevaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerGaussSubdividedEigen, basisAttrMap); + std::shared_ptr> feBDElecChargeRhs = + std::make_shared> + (basisDofHandlerTotalPot, quadAttrGaussSubdivided, basisAttrMap, ksdft::KSDFTDefaults::CELL_BATCH_SIZE_GRAD_EVAL, *linAlgOpContext); + feBDElecChargeRhs->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerGaussSubdividedEigen, basisAttrMap); + + 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> feBDHamStiffnessMatrix = + std::make_shared> + (basisDofHandlerWaveFn, quadAttrEigen, basisAttrMap, ksdft::KSDFTDefaults::CELL_BATCH_SIZE_GRAD_EVAL, *linAlgOpContext); + + // evaluate basis data + feBDHamStiffnessMatrix->evaluateBasisData(quadAttrEigen, basisAttrMap); + std::shared_ptr quadRuleContainerRho = feBDElectrostaticsHamiltonian->getQuadratureRuleContainer(); - + // scale the electronic charges quadrature::QuadratureValuesContainer electronChargeDensity(quadRuleContainerRho, 1, 0.0); - std::shared_ptr> feBDKineticHamiltonian = feBDElectrostaticsHamiltonian; + std::shared_ptr> feBDKineticHamiltonian = feBDHamStiffnessMatrix; std::shared_ptr> feBDEXCHamiltonian = feBDElectrostaticsHamiltonian; // Create OperatorContext for Basisoverlap @@ -918,8 +938,6 @@ std::shared_ptrevaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerGaussSubdividedElec, basisAttrMap); - std::shared_ptr> feBDElecChargeRhs = feBDElectrostaticsHamiltonian; - for (size_type iCell = 0; iCell < electronChargeDensity.nCells(); iCell++) { size_type quadId = 0; @@ -957,7 +975,7 @@ std::shared_ptr> feBDNuclearChargeRhs = feBDNucChargeRhs; - std::shared_ptr> feBDNuclearChargeStiffnessMatrix = feBasisDataElec; + std::shared_ptr> feBDNuclearChargeStiffnessMatrix = feBDTotalChargeStiffnessMatrix; std::shared_ptr(parameterInputFileName, "num1DGaussSubdividedSizeEigen", rootCout); unsigned int gaussSubdividedCopiesEigen = readParameter(parameterInputFileName, "gaussSubdividedCopiesEigen", rootCout); + unsigned int num1DGaussSubdividedSizeGrad = readParameter(parameterInputFileName, "num1DGaussSubdividedSizeGrad", rootCout); + unsigned int gaussSubdividedCopiesGrad = readParameter(parameterInputFileName, "gaussSubdividedCopiesGrad", rootCout); + bool isNumericalNuclearSolve = readParameter(parameterInputFileName, "isNumericalNuclearSolve", rootCout); // Set up Triangulation @@ -967,6 +970,16 @@ int main(int argc, char** argv) triangulationBase, *cellMapping); + std::shared_ptr gaussSubdivQuadRuleGrad = + std::make_shared(dim, num1DGaussSubdividedSizeGrad, gaussSubdividedCopiesGrad); + + std::shared_ptr quadRuleContainerAdaptiveGrad = + std::make_shared + (quadAttrGaussSubdivided, + gaussSubdivQuadRuleGrad, + triangulationBase, + *cellMapping); + // add device synchronize for gpu utils::mpi::MPIBarrier(comm); stop = std::chrono::high_resolution_clock::now(); @@ -1125,7 +1138,7 @@ int main(int argc, char** argv) utils::mpi::MPIBarrier(comm); start = std::chrono::high_resolution_clock::now(); - efeBasisDataAdaptiveTotPot->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerAdaptiveOrbital, basisAttrMap); + efeBasisDataAdaptiveTotPot->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerAdaptiveGrad, basisAttrMap); std::shared_ptr> feBDTotalChargeStiffnessMatrix = efeBasisDataAdaptiveTotPot; @@ -1155,7 +1168,7 @@ int main(int argc, char** argv) rootCout << "Time for electrostatics basis datastorage evaluation is(in secs) : " << duration.count()/1e6 << std::endl; basisAttrMap[basis::BasisStorageAttributes::StoreValues] = true; - basisAttrMap[basis::BasisStorageAttributes::StoreGradient] = true; + basisAttrMap[basis::BasisStorageAttributes::StoreGradient] = false; basisAttrMap[basis::BasisStorageAttributes::StoreHessian] = false; basisAttrMap[basis::BasisStorageAttributes::StoreOverlap] = false; basisAttrMap[basis::BasisStorageAttributes::StoreGradNiGradNj] = false; @@ -1165,14 +1178,27 @@ int main(int argc, char** argv) std::make_shared> (basisDofHandlerWaveFn, quadAttrGaussSubdivided, basisAttrMap, ksdft::KSDFTDefaults::CELL_BATCH_SIZE_GRAD_EVAL, *linAlgOpContext); + efeBasisDataAdaptiveOrbital->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerAdaptiveOrbital, basisAttrMap); + + 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; + + std::shared_ptr> efeBasisDataAdaptiveGrad = + std::make_shared> + (basisDofHandlerWaveFn, quadAttrGaussSubdivided, basisAttrMap, ksdft::KSDFTDefaults::CELL_BATCH_SIZE_GRAD_EVAL, *linAlgOpContext); + // add device synchronize for gpu utils::mpi::MPIBarrier(comm); start = std::chrono::high_resolution_clock::now(); - efeBasisDataAdaptiveOrbital->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerAdaptiveOrbital, basisAttrMap); + efeBasisDataAdaptiveGrad->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerAdaptiveGrad, basisAttrMap); std::shared_ptr> feBDElectrostaticsHamiltonian = efeBasisDataAdaptiveOrbital; - std::shared_ptr> feBDKineticHamiltonian = efeBasisDataAdaptiveOrbital; + std::shared_ptr> feBDKineticHamiltonian = efeBasisDataAdaptiveGrad; std::shared_ptr> feBDEXCHamiltonian = efeBasisDataAdaptiveOrbital; // add device synchronize for gpu diff --git a/src/electrostatics/PoissonLinearSolverFunctionFE.t.cpp b/src/electrostatics/PoissonLinearSolverFunctionFE.t.cpp index d0b6a17b..72ec9cac 100644 --- a/src/electrostatics/PoissonLinearSolverFunctionFE.t.cpp +++ b/src/electrostatics/PoissonLinearSolverFunctionFE.t.cpp @@ -540,6 +540,8 @@ namespace dftefe d_feBasisManagerField->getConstraints().distributeParentToChild( solution, numComponents); + + d_initial = solution; } template LinearEigenSolverDefaults::CHEBY_ORDER_LOOKUP{ + {100, 15}, + {300, 19}, {500, 24}, // <= 500 ~> chebyshevOrder = 24 {750, 30}, {1000, 39}, diff --git a/src/ksdft/ElectrostaticLocalFE.h b/src/ksdft/ElectrostaticLocalFE.h index b9c4bb6b..7ed6eb56 100644 --- a/src/ksdft/ElectrostaticLocalFE.h +++ b/src/ksdft/ElectrostaticLocalFE.h @@ -219,6 +219,8 @@ namespace dftefe const std::vector d_smearedChargeRadius; RealType d_energy; const utils::ScalarSpatialFunctionReal &d_externalPotentialFunction; + + // Causing memory errors: Change these to smart pointers quadrature::QuadratureValuesContainer *d_nuclearChargesDensity; const quadrature::QuadratureValuesContainer diff --git a/src/ksdft/ElectrostaticLocalFE.t.cpp b/src/ksdft/ElectrostaticLocalFE.t.cpp index b64578a2..0a20259f 100644 --- a/src/ksdft/ElectrostaticLocalFE.t.cpp +++ b/src/ksdft/ElectrostaticLocalFE.t.cpp @@ -291,11 +291,11 @@ namespace dftefe delete d_correctionPotHamQuad; d_correctionPotHamQuad = nullptr; } - if (d_correctionPotRhoQuad != nullptr) - { - delete d_correctionPotRhoQuad; - d_correctionPotRhoQuad = nullptr; - } + // if (d_correctionPotRhoQuad != nullptr) + // { + // delete d_correctionPotRhoQuad; + // d_correctionPotRhoQuad = nullptr; + // } if (d_totalChargePotential != nullptr) { delete d_totalChargePotential;