Skip to content

Commit

Permalink
Modified OrthoEFEOverlapOperator and its inverse classes
Browse files Browse the repository at this point in the history
  • Loading branch information
Avirup Sircar committed Aug 22, 2024
1 parent 555e3e6 commit 5749afe
Show file tree
Hide file tree
Showing 13 changed files with 465 additions and 213 deletions.
28 changes: 14 additions & 14 deletions analysis/classicalEnrichmentComparison/KSDFTClassical/TestKohnShamDft.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -393,18 +393,18 @@ int main(int argc, char** argv)
const utils::ScalarSpatialFunctionReal *externalPotentialFunction = new
utils::PointChargePotentialFunction(atomCoordinatesVec, atomChargesVec);

// // Create OperatorContext for Basisoverlap
// std::shared_ptr<const basis::CFEOverlapOperatorContext<double,
// double,
// Host,
// dim>> MContext =
// std::make_shared<basis::CFEOverlapOperatorContext<double,
// double,
// Host,
// dim>>(
// *basisManagerWaveFn,
// *feBasisData,
// 50);
// Create OperatorContext for Basisoverlap
std::shared_ptr<const basis::CFEOverlapOperatorContext<double,
double,
Host,
dim>> MContext =
std::make_shared<basis::CFEOverlapOperatorContext<double,
double,
Host,
dim>>(
*basisManagerWaveFn,
*feBasisDataElec,
50);

// Set up the quadrature rule

Expand Down Expand Up @@ -562,7 +562,7 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
50,
50,
*MContextForInv,
*MContextForInv,
*MContext,
*MInvContext);

dftefeSolve->solve();
Expand Down Expand Up @@ -611,7 +611,7 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
50,
50,
*MContextForInv,
*MContextForInv,
*MContext,
*MInvContext);

dftefeSolve->solve();
Expand Down
18 changes: 9 additions & 9 deletions analysis/classicalEnrichmentComparison/KSDFTClassical/param.in
View file
Original file line number Diff line number Diff line change
@@ -1,18 +1,18 @@
xmax = 26.001
ymax = 26.001
zmax = 26.001
xmax = 24
ymax = 24
zmax = 24
radiusAtAtom = 0
meshSizeAtAtom = 2
radiusAroundAtom = 4
meshSizeAroundAtom = 2
meshSizeAtAtom = 0.02
radiusAroundAtom = 2
meshSizeAroundAtom = 0.1
rc = 0.6
feOrderElectrostatics = 6
feOrderEigenSolve = 3
feOrderElectrostatics = 5
feOrderEigenSolve = 5
smearingTemperature = 500.
fermiEnergyTolerance = 1e-8
fracOccupancyTolerance = 1e-3
eigenSolveResidualTolerance = 1e-3
chebyshevPolynomialDegree = 50
chebyshevPolynomialDegree = 1250
maxChebyshevFilterPass = 40
numWantedEigenvalues = 20
scfDensityResidualNormTolerance = 1e-5
Expand Down
4 changes: 2 additions & 2 deletions analysis/classicalEnrichmentComparison/KSDFTOrthoEFE/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,6 @@ add_subdirectory(${MAIN_PROJECT_DIR}/src/ksdft ${MAIN_PROJECT_DIR}/test/ksdft/li

if(ENABLE_MPI)
add_compile_definitions(DFTEFE_WITH_MPI)
add_executable(TestKohnShamDft ${MAIN_PROJECT_DIR}/analysis/classicalEnrichmentComparison/KSDFTOrthoEFE/TestKohnShamDft.cpp)
target_link_libraries(TestKohnShamDft PUBLIC dft-efe-basis dft-efe-electrostatics dft-efe-ksdft dft-efe-utils dft-efe-quadrature dft-efe-linalg dft-efe-atoms)
add_executable(TestKohnShamDft_ideal ${MAIN_PROJECT_DIR}/analysis/classicalEnrichmentComparison/KSDFTOrthoEFE/TestKohnShamDft.cpp)
target_link_libraries(TestKohnShamDft_ideal PUBLIC dft-efe-basis dft-efe-electrostatics dft-efe-ksdft dft-efe-utils dft-efe-quadrature dft-efe-linalg dft-efe-atoms)
endif()
44 changes: 25 additions & 19 deletions analysis/classicalEnrichmentComparison/KSDFTOrthoEFE/TestKohnShamDft.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -754,7 +754,7 @@ int main(int argc, char** argv)
<double, Host, dim>>
enrichClassIntfceOrbital = std::make_shared<basis::EnrichmentClassicalInterfaceSpherical
<double, Host, dim>>
(cfeBasisDataStorageAdaptiveOrbital,
(cfeBasisDataStorageGLLEigen,
cfeBasisDataStorageAdaptiveOrbital,
atomSphericalDataContainer,
atomPartitionTolerance,
Expand Down Expand Up @@ -933,21 +933,22 @@ int main(int argc, char** argv)
// utils::mpi::MPIBarrier(comm);
// start = std::chrono::high_resolution_clock::now();

// // Create OperatorContext for Basisoverlap

// std::shared_ptr<const basis::OrthoEFEOverlapOperatorContext<double,
// double,
// Host,
// dim>> MContext =
// std::make_shared<basis::OrthoEFEOverlapOperatorContext<double,
// double,
// Host,
// dim>>(
// *basisManagerWaveFn,
// *cfeBasisDataStorageAdaptiveOrbital,
// *efeBasisDataAdaptiveOrbital,
// /**cfeBasisDataStorageGLLEigen,*/
// 50);
// Create OperatorContext for Basisoverlap

std::shared_ptr<const basis::OrthoEFEOverlapOperatorContext<double,
double,
Host,
dim>> MContext =
std::make_shared<basis::OrthoEFEOverlapOperatorContext<double,
double,
Host,
dim>>(
*basisManagerWaveFn,
*cfeBasisDataStorageAdaptiveOrbital,
*efeBasisDataAdaptiveOrbital,
*cfeBasisDataStorageAdaptiveOrbital,
50,
true);

// // add device synchronize for gpu
// utils::mpi::MPIBarrier(comm);
Expand Down Expand Up @@ -994,6 +995,7 @@ int main(int argc, char** argv)
*basisManagerWaveFn,
*cfeBasisDataStorageGLLEigen,
*efeBasisDataAdaptiveOrbital,
*cfeBasisDataStorageGLLEigen,
linAlgOpContext);

// add device synchronize for gpu
Expand All @@ -1016,8 +1018,10 @@ int main(int argc, char** argv)
Host,
dim>>
(*basisManagerWaveFn,
*cfeBasisDataStorageGLLEigen,
*MContext,
/**cfeBasisDataStorageGLLEigen,
*efeBasisDataAdaptiveOrbital,
*cfeBasisDataStorageGLLEigen,*/
linAlgOpContext);

// add device synchronize for gpu
Expand Down Expand Up @@ -1114,7 +1118,8 @@ int main(int argc, char** argv)
50,
50,
*MContextForInv,
*MContextForInv,
*MContext,
/**MContextForInv,*/
*MInvContext);

// add device synchronize for gpu
Expand Down Expand Up @@ -1176,7 +1181,8 @@ int main(int argc, char** argv)
50,
50,
*MContextForInv,
*MContextForInv,
*MContext,
/**MContextForInv,*/
*MInvContext);

// add device synchronize for gpu
Expand Down
18 changes: 9 additions & 9 deletions analysis/classicalEnrichmentComparison/KSDFTOrthoEFE/param.in
View file
Original file line number Diff line number Diff line change
@@ -1,24 +1,24 @@
xmax = 26.001
xmax = 26.001
ymax = 26.001
zmax = 26.001
radiusAtAtom = 0
meshSizeAtAtom = 0.5
radiusAroundAtom = 2
meshSizeAroundAtom = 0.5
meshSizeAtAtom = 2.8
radiusAroundAtom = 4
meshSizeAroundAtom = 2.8
rc = 0.6
feOrderElectrostatics = 5
feOrderElectrostatics = 6
feOrderEigenSolve = 2
atomPartitionTolerance = 1e-14
smallestCellVolume = 1e-14
maxRecursion = 100
adaptiveQuadAbsTolerance = 1e-7
adaptiveQuadRelTolerance = 1e-7
integralThreshold = 1e-7
adaptiveQuadAbsTolerance = 1e-6
adaptiveQuadRelTolerance = 1e-6
integralThreshold = 1e-6
smearingTemperature = 500.
fermiEnergyTolerance = 1e-8
fracOccupancyTolerance = 1e-3
eigenSolveResidualTolerance = 1e-3
chebyshevPolynomialDegree = 50
chebyshevPolynomialDegree = 20
maxChebyshevFilterPass = 40
numWantedEigenvalues = 20
scfDensityResidualNormTolerance = 1e-5
Expand Down
18 changes: 18 additions & 0 deletions src/atoms/SphericalHarmonicFunctions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,7 @@ namespace dftefe
return 1.0 / sqrt(M_PI);
}

/*
double
Clm(const int l, const int m)
{
Expand All @@ -116,6 +117,23 @@ namespace dftefe
return sqrt(((2.0 * l + 1) * boost::math::factorial<double>(l - abs(m))) /
(2.0 * boost::math::factorial<double>(l + abs(m))));
}
*/

// Implement this instead of above function to remove underflow/overflow issues in factorial
double Blm(const int l, const int m)
{
if (m==0)
return sqrt((2.0*l+1)/2.0);
else
return Blm(l, m-1)/sqrt((l-m+1.0)*(l+m));
}

double Clm(const int l, const int m)
{
// assert(m >= 0);
assert(std::abs(m) <= l);
return Blm(l, abs(m));
}

double
Qm(const int m, const double phi)
Expand Down
2 changes: 1 addition & 1 deletion src/basis/GenerateMesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -201,7 +201,7 @@ namespace dftefe
{
for (auto &i : baseMeshSize)
{
i = std::pow(2, round(log2((4.0) / largestMeshSizeAroundAtom))) *
i = std::pow(2, round(log2((std::max(4.0, largestMeshSizeAroundAtom)/*4.0*/) / largestMeshSizeAroundAtom))) *
largestMeshSizeAroundAtom;
}
}
Expand Down
4 changes: 2 additions & 2 deletions src/basis/OrthoEFEOverlapInverseOpContextGLL.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,8 +70,8 @@ namespace dftefe
&classicalBlockGLLBasisDataStorage,
const FEBasisDataStorage<ValueTypeOperator, memorySpace>
&enrichmentBlockEnrichmentBasisDataStorage,
/*const FEBasisDataStorage<ValueTypeOperator, memorySpace>
&enrichmentBlockClassicalBasisDataStorage,*/
const FEBasisDataStorage<ValueTypeOperator, memorySpace>
&enrichmentBlockClassicalBasisDataStorage,
std::shared_ptr<linearAlgebra::LinAlgOpContext<memorySpace>>
linAlgOpContext);

Expand Down
29 changes: 15 additions & 14 deletions src/basis/OrthoEFEOverlapInverseOpContextGLL.t.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -223,8 +223,8 @@ namespace dftefe
&classicalBlockGLLBasisDataStorage,
const FEBasisDataStorage<ValueTypeOperator, memorySpace>
&enrichmentBlockEnrichmentBasisDataStorage,
/*const FEBasisDataStorage<ValueTypeOperator, memorySpace>
&enrichmentBlockClassicalBasisDataStorage, */
const FEBasisDataStorage<ValueTypeOperator, memorySpace>
&enrichmentBlockClassicalBasisDataStorage,
std::shared_ptr<
const FEBasisDofHandler<ValueTypeOperand, memorySpace, dim>> cfeBDH,
std::shared_ptr<const EFEBasisDofHandler<ValueTypeOperand,
Expand Down Expand Up @@ -285,7 +285,7 @@ namespace dftefe

locallyOwnedCellIter = efeBDH->beginLocallyOwnedCells();

/**

const std::unordered_map<global_size_type,
utils::OptimizedIndexSet<size_type>>
*enrichmentIdToClassicalLocalIdMap =
Expand All @@ -305,7 +305,7 @@ namespace dftefe
memorySpace,
dim>>(
eci->getCFEBasisManager());
**/


for (; locallyOwnedCellIter != efeBDH->endLocallyOwnedCells();
++locallyOwnedCellIter)
Expand All @@ -327,14 +327,14 @@ namespace dftefe
.getQuadratureRuleContainer()
->getCellJxW(cellIndex);

/**

size_type nQuadPointInCellEnrichmentBlockClassical =
enrichmentBlockClassicalBasisDataStorage.getQuadratureRuleContainer()
->nCellQuadraturePoints(cellIndex);
std::vector<double> cellJxWValuesEnrichmentBlockClassical =
enrichmentBlockClassicalBasisDataStorage.getQuadratureRuleContainer()
->getCellJxW(cellIndex);
**/


const ValueTypeOperator *cumulativeClassicalBlockDofQuadPoints =
basisDataInAllCellsClassicalBlock.data(); /*GLL Quad rule*/
Expand All @@ -344,7 +344,7 @@ namespace dftefe
basisDataInAllCellsEnrichmentBlockEnrichment.data() +
cumulativeEnrichmentBlockEnrichmentDofQuadPointsOffset;

/**

std::vector<utils::Point> quadRealPointsVec =
enrichmentBlockEnrichmentBasisDataStorage
.getQuadratureRuleContainer()
Expand Down Expand Up @@ -466,7 +466,7 @@ namespace dftefe
numEnrichmentIdsInCell,
*eci->getLinAlgOpContext());
}
**/


for (unsigned int iNode = 0; iNode < dofsPerCell; iNode++)
{
Expand All @@ -493,6 +493,7 @@ namespace dftefe

else if (iNode >= dofsPerCellCFE && jNode >= dofsPerCellCFE)
{
/**
// Ni_pristine*Ni_pristine at quadpoints
for (unsigned int qPoint = 0;
qPoint < nQuadPointInCellEnrichmentBlockEnrichment;
Expand All @@ -509,8 +510,8 @@ namespace dftefe
qPoint) *
cellJxWValuesEnrichmentBlockEnrichment[qPoint];
}

/**
**/

ValueTypeOperator NpiNpj = (ValueTypeOperator)0,
ciNciNpj = (ValueTypeOperator)0,
NpicjNcj = (ValueTypeOperator)0,
Expand Down Expand Up @@ -581,7 +582,7 @@ namespace dftefe
}
*basisOverlapTmpIter +=
NpiNpj - NpicjNcj - ciNciNpj + ciNcicjNcj;
**/

}
basisOverlapTmpIter++;
}
Expand Down Expand Up @@ -620,8 +621,8 @@ namespace dftefe
&classicalBlockGLLBasisDataStorage,
const FEBasisDataStorage<ValueTypeOperator, memorySpace>
&enrichmentBlockEnrichmentBasisDataStorage,
/*const FEBasisDataStorage<ValueTypeOperator, memorySpace>
&enrichmentBlockClassicalBasisDataStorage,*/
const FEBasisDataStorage<ValueTypeOperator, memorySpace>
&enrichmentBlockClassicalBasisDataStorage,
std::shared_ptr<linearAlgebra::LinAlgOpContext<memorySpace>>
linAlgOpContext)
: d_feBasisManager(&feBasisManager)
Expand Down Expand Up @@ -725,7 +726,7 @@ namespace dftefe
memorySpace,
dim>(classicalBlockGLLBasisDataStorage,
enrichmentBlockEnrichmentBasisDataStorage,
/*enrichmentBlockClassicalBasisDataStorage, */
enrichmentBlockClassicalBasisDataStorage,
cfeBDH,
efeBDH,
NiNjInAllCells);
Expand Down
Loading

0 comments on commit 5749afe

Please sign in to comment.