Skip to content

Commit

Permalink
Minor modifications in KSDFT/Defaults.h
Browse files Browse the repository at this point in the history
  • Loading branch information
Avirup Sircar committed Dec 7, 2024
1 parent fdf5a86 commit 9c0b1c3
Show file tree
Hide file tree
Showing 6 changed files with 71 additions and 21 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -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<basis::FEBasisDataStorage<double, Host>> feBasisDataElec =
std::shared_ptr<basis::FEBasisDataStorage<double, Host>> feBDTotalChargeStiffnessMatrix =
std::make_shared<basis::CFEBDSOnTheFlyComputeDealii<double, double, Host,dim>>
(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<const utils::ScalarSpatialFunctionReal>
zeroFunction = std::make_shared
Expand Down Expand Up @@ -804,8 +804,6 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
(*basisManagerWaveFn,
*feBasisDataGLLEigen,
linAlgOpContext);

std::shared_ptr<const basis::FEBasisDataStorage<double, Host>> feBDTotalChargeStiffnessMatrix = feBasisDataElec;

std::shared_ptr<quadrature::QuadratureRule> gaussSubdivQuadRuleElec =
std::make_shared<quadrature::QuadratureRuleGaussIterated>(dim, num1DGaussSizeElec, gaussSubdividedCopiesElec);
Expand Down Expand Up @@ -870,7 +868,7 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
rootCout << "Number of quadrature points in gauss subdivided quadrature eigen: "<< nQuad<<"\n";

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;
Expand All @@ -881,14 +879,36 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
(basisDofHandlerWaveFn, quadAttrGaussSubdivided, basisAttrMap, ksdft::KSDFTDefaults::CELL_BATCH_SIZE_GRAD_EVAL, *linAlgOpContext);
feBDElectrostaticsHamiltonian->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerGaussSubdividedEigen, basisAttrMap);

std::shared_ptr<basis::FEBasisDataStorage<double, Host>> feBDElecChargeRhs =
std::make_shared<basis::CFEBDSOnTheFlyComputeDealii<double, double, Host,dim>>
(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<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);

// evaluate basis data
feBDHamStiffnessMatrix->evaluateBasisData(quadAttrEigen, basisAttrMap);

std::shared_ptr<const quadrature::QuadratureRuleContainer> quadRuleContainerRho =
feBDElectrostaticsHamiltonian->getQuadratureRuleContainer();

// scale the electronic charges
quadrature::QuadratureValuesContainer<double, Host>
electronChargeDensity(quadRuleContainerRho, 1, 0.0);

std::shared_ptr<const basis::FEBasisDataStorage<double,Host>> feBDKineticHamiltonian = feBDElectrostaticsHamiltonian;
std::shared_ptr<const basis::FEBasisDataStorage<double,Host>> feBDKineticHamiltonian = feBDHamStiffnessMatrix;
std::shared_ptr<const basis::FEBasisDataStorage<double, Host>> feBDEXCHamiltonian = feBDElectrostaticsHamiltonian;

// Create OperatorContext for Basisoverlap
Expand Down Expand Up @@ -918,8 +938,6 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
(basisDofHandlerTotalPot, quadAttrGaussSubdivided, basisAttrMap, ksdft::KSDFTDefaults::CELL_BATCH_SIZE_GRAD_EVAL, *linAlgOpContext);
feBDNucChargeRhs->evaluateBasisData(quadAttrGaussSubdivided, quadRuleContainerGaussSubdividedElec, basisAttrMap);

std::shared_ptr<basis::FEBasisDataStorage<double, Host>> feBDElecChargeRhs = feBDElectrostaticsHamiltonian;

for (size_type iCell = 0; iCell < electronChargeDensity.nCells(); iCell++)
{
size_type quadId = 0;
Expand Down Expand Up @@ -957,7 +975,7 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
// Set up the FE Basis Data Storage
std::shared_ptr<basis::FEBasisDataStorage<double, Host>> feBDNuclearChargeRhs = feBDNucChargeRhs;

std::shared_ptr<const basis::FEBasisDataStorage<double,Host>> feBDNuclearChargeStiffnessMatrix = feBasisDataElec;
std::shared_ptr<const basis::FEBasisDataStorage<double,Host>> feBDNuclearChargeStiffnessMatrix = feBDTotalChargeStiffnessMatrix;

std::shared_ptr<ksdft::KohnShamDFT<double,
double,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -616,6 +616,9 @@ int main(int argc, char** argv)
unsigned int num1DGaussSubdividedSizeEigen = readParameter<unsigned int>(parameterInputFileName, "num1DGaussSubdividedSizeEigen", rootCout);
unsigned int gaussSubdividedCopiesEigen = readParameter<unsigned int>(parameterInputFileName, "gaussSubdividedCopiesEigen", rootCout);

unsigned int num1DGaussSubdividedSizeGrad = readParameter<unsigned int>(parameterInputFileName, "num1DGaussSubdividedSizeGrad", rootCout);
unsigned int gaussSubdividedCopiesGrad = readParameter<unsigned int>(parameterInputFileName, "gaussSubdividedCopiesGrad", rootCout);

bool isNumericalNuclearSolve = readParameter<bool>(parameterInputFileName, "isNumericalNuclearSolve", rootCout);

// Set up Triangulation
Expand Down Expand Up @@ -967,6 +970,16 @@ int main(int argc, char** argv)
triangulationBase,
*cellMapping);

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);

// add device synchronize for gpu
utils::mpi::MPIBarrier(comm);
stop = std::chrono::high_resolution_clock::now();
Expand Down Expand Up @@ -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<const basis::FEBasisDataStorage<double, Host>> feBDTotalChargeStiffnessMatrix = efeBasisDataAdaptiveTotPot;

Expand Down Expand Up @@ -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;
Expand All @@ -1165,14 +1178,27 @@ 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);

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<basis::FEBasisDataStorage<double, Host>> efeBasisDataAdaptiveGrad =
std::make_shared<basis::EFEBDSOnTheFlyComputeDealii<double, double, Host,dim>>
(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<const basis::FEBasisDataStorage<double, Host>> feBDElectrostaticsHamiltonian = efeBasisDataAdaptiveOrbital;
std::shared_ptr<const basis::FEBasisDataStorage<double,Host>> feBDKineticHamiltonian = efeBasisDataAdaptiveOrbital;
std::shared_ptr<const basis::FEBasisDataStorage<double,Host>> feBDKineticHamiltonian = efeBasisDataAdaptiveGrad;
std::shared_ptr<const basis::FEBasisDataStorage<double, Host>> feBDEXCHamiltonian = efeBasisDataAdaptiveOrbital;

// add device synchronize for gpu
Expand Down
2 changes: 2 additions & 0 deletions src/electrostatics/PoissonLinearSolverFunctionFE.t.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -540,6 +540,8 @@ namespace dftefe

d_feBasisManagerField->getConstraints().distributeParentToChild(
solution, numComponents);

d_initial = solution;
}

template <typename ValueTypeOperator,
Expand Down
6 changes: 4 additions & 2 deletions src/ksdft/Defaults.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -34,8 +34,8 @@ namespace dftefe
const linearAlgebra::PreconditionerType PoissonProblemDefaults::PC_TYPE =
linearAlgebra::PreconditionerType::JACOBI;
const size_type PoissonProblemDefaults::MAX_ITER = 2e7;
const double PoissonProblemDefaults::ABSOLUTE_TOL = 1e-12;
const double PoissonProblemDefaults::RELATIVE_TOL = 1e-14;
const double PoissonProblemDefaults::ABSOLUTE_TOL = 1e-10;
const double PoissonProblemDefaults::RELATIVE_TOL = 1e-12;
const double PoissonProblemDefaults::DIVERGENCE_TOL = 1e10;

/**
Expand All @@ -50,6 +50,8 @@ namespace dftefe
*/
const std::map<size_type, size_type>
LinearEigenSolverDefaults::CHEBY_ORDER_LOOKUP{
{100, 15},
{300, 19},
{500, 24}, // <= 500 ~> chebyshevOrder = 24
{750, 30},
{1000, 39},
Expand Down
2 changes: 2 additions & 0 deletions src/ksdft/ElectrostaticLocalFE.h
Original file line number Diff line number Diff line change
Expand Up @@ -219,6 +219,8 @@ namespace dftefe
const std::vector<double> d_smearedChargeRadius;
RealType d_energy;
const utils::ScalarSpatialFunctionReal &d_externalPotentialFunction;

// Causing memory errors: Change these to smart pointers
quadrature::QuadratureValuesContainer<RealType, memorySpace>
*d_nuclearChargesDensity;
const quadrature::QuadratureValuesContainer<ValueType, memorySpace>
Expand Down
10 changes: 5 additions & 5 deletions src/ksdft/ElectrostaticLocalFE.t.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down

0 comments on commit 9c0b1c3

Please sign in to comment.