From ae34bdd2a2f36de22ac599a3ad9ee27ccb79d36b Mon Sep 17 00:00:00 2001 From: Avirup Sircar Date: Sun, 3 Nov 2024 10:35:55 -0500 Subject: [PATCH] Modified CFBasisDataStorageOnTheFlyComputeDealii --- src/basis/CFEBDSOnTheFlyComputeDealii.h | 2 +- src/basis/CFEBDSOnTheFlyComputeDealii.t.cpp | 42 +- src/basis/EFEBasisDataStorageDealii.t.cpp | 375 ++++++------------ .../OrthoEFEOverlapInverseOpContextGLL.t.cpp | 14 +- .../OrthoEFEOverlapOperatorContext.t.cpp | 10 +- test/basis/src/TestBasisDataStorageMemOpt.cpp | 303 ++++++++------ 6 files changed, 343 insertions(+), 403 deletions(-) diff --git a/src/basis/CFEBDSOnTheFlyComputeDealii.h b/src/basis/CFEBDSOnTheFlyComputeDealii.h index 147a0822..bfcb9301 100644 --- a/src/basis/CFEBDSOnTheFlyComputeDealii.h +++ b/src/basis/CFEBDSOnTheFlyComputeDealii.h @@ -232,7 +232,7 @@ namespace dftefe std::vector d_cellStartIdsBasisHessianQuadStorage; std::vector d_cellStartIdsBasisQuadStorage; size_type d_maxCellBlock; - Storage d_tmpGradientBlock; + std::shared_ptr d_tmpGradientBlock; linearAlgebra::LinAlgOpContext &d_linAlgOpContext; }; // end of CFEBDSOnTheFlyComputeDealii diff --git a/src/basis/CFEBDSOnTheFlyComputeDealii.t.cpp b/src/basis/CFEBDSOnTheFlyComputeDealii.t.cpp index 8a6667cc..1ccb1043 100644 --- a/src/basis/CFEBDSOnTheFlyComputeDealii.t.cpp +++ b/src/basis/CFEBDSOnTheFlyComputeDealii.t.cpp @@ -268,11 +268,11 @@ namespace dftefe if (basisStorageAttributesBoolMap .find(BasisStorageAttributes::StoreGradient) ->second) - dealiiUpdateFlagsPara |= dealii::update_gradients; + dealiiUpdateFlagsPara = dealii::update_gradients; // This is for getting the gradient in parametric cell dealii::FEValues dealiiFEValuesPara(feBDH->getReferenceFE(cellId), dealiiQuadratureRule, - dealiiUpdateFlags); + dealiiUpdateFlagsPara); if (basisStorageAttributesBoolMap .find(BasisStorageAttributes::StoreGradient) @@ -546,8 +546,7 @@ namespace dftefe { d_dofsInCell[iCell] = d_feBDH->nCellDofs(iCell); } - d_tmpGradientBlock = - typename BasisDataStorage::Storage(0); + d_tmpGradientBlock = nullptr; } template ( + d_dofsInCell[0] * nQuadPointsInCell[0] * dim * d_maxCellBlock); size_type gradientParaCellSize = basisGradientParaCellQuadStorage->size(); for (size_type iCell = 0; iCell < d_maxCellBlock; ++iCell) { - d_tmpGradientBlock.template copyFrom( + d_tmpGradientBlock->template copyFrom( basisGradientParaCellQuadStorage->data(), gradientParaCellSize, 0, @@ -802,13 +801,13 @@ namespace dftefe d_basisJacobianInvQuadStorage = basisJacobianInvQuadStorage; d_cellStartIdsBasisJacobianInvQuadStorage = cellStartIdsBasisJacobianInvQuadStorage; - d_tmpGradientBlock.resize(d_dofsInCell[0] * nQuadPointsInCell[0] * - dim * d_maxCellBlock); + d_tmpGradientBlock = std::make_shared( + d_dofsInCell[0] * nQuadPointsInCell[0] * dim * d_maxCellBlock); size_type gradientParaCellSize = basisGradientParaCellQuadStorage->size(); for (size_type iCell = 0; iCell < d_maxCellBlock; ++iCell) { - d_tmpGradientBlock.template copyFrom( + d_tmpGradientBlock->template copyFrom( basisGradientParaCellQuadStorage->data(), gradientParaCellSize, 0, @@ -945,7 +944,7 @@ namespace dftefe // typename BasisDataStorage::Storage // dummy( // 0); - return d_tmpGradientBlock; + return *d_tmpGradientBlock; } template second, "Basis gradient values are not evaluated for the given QuadratureRuleAttributes"); + std::shared_ptr tmpGradientBlock = nullptr; if ((cellRange.second - cellRange.first) > d_maxCellBlock) { - d_tmpGradientBlock.resize(d_dofsInCell[0] * d_nQuadPointsIncell[0] * - dim * (cellRange.second - cellRange.first)); - d_tmpGradientBlock.setValue(0); + tmpGradientBlock = std::make_shared( + d_dofsInCell[0] * d_nQuadPointsIncell[0] * dim * + (cellRange.second - cellRange.first)); size_type gradientParaCellSize = d_basisGradientParaCellQuadStorage->size(); for (size_type iCell = 0; iCell < (cellRange.second - cellRange.first); ++iCell) { - d_tmpGradientBlock.template copyFrom( + tmpGradientBlock->template copyFrom( d_basisGradientParaCellQuadStorage->data(), gradientParaCellSize, 0, gradientParaCellSize * iCell); } } + else + { + tmpGradientBlock = d_tmpGradientBlock; + } CFEBDSOnTheFlyComputeDealiiInternal:: computeJacobianInvTimesGradPara( cellRange, @@ -1155,7 +1159,7 @@ namespace dftefe d_nQuadPointsIncell, d_basisJacobianInvQuadStorage, d_cellStartIdsBasisJacobianInvQuadStorage, - d_tmpGradientBlock, + *tmpGradientBlock, d_linAlgOpContext, basisGradientData); } @@ -1375,7 +1379,7 @@ namespace dftefe // typename BasisDataStorage::Storage // dummy( // 0); - return d_tmpGradientBlock; + return *d_tmpGradientBlock; } template resize(0); } template ::Storage // dummy( // 0); - return d_tmpGradientBlock; + return *d_tmpGradientBlock; } template > efeBDH, - std::shared_ptr< - FEBasisDataStorage> + std::shared_ptr> cfeBasisDataStorage, std::vector &classicalComponentInQuadValues) { @@ -165,10 +164,10 @@ namespace dftefe size_type numEnrichmentIdsInCell = efeBDH->nCellDofs(cellIndex) - classicalDofsPerCell; - dftefe::utils::MemoryStorage + dftefe::utils::MemoryStorage basisValInCell = cfeBasisDataStorage->getBasisDataInCell(cellIndex); + ValueTypeBasisData *B = basisValInCell.data(); // Do a gemm (\Sigma c_i N_i^classical) // and get the quad values in std::vector @@ -184,7 +183,7 @@ namespace dftefe (ValueTypeBasisData)1.0, coeffsInCell.data(), numEnrichmentIdsInCell, - basisValInCell.data(), + B, classicalDofsPerCell, (ValueTypeBasisData)0.0, classicalComponentInQuadValues.data(), @@ -205,8 +204,7 @@ namespace dftefe ValueTypeBasisData, memorySpace, dim>> efeBDH, - std::shared_ptr< - FEBasisDataStorage> + std::shared_ptr> cfeBasisDataStorage, std::vector &classicalComponentInQuadGradients) { @@ -217,8 +215,7 @@ namespace dftefe efeBDH->nCellDofs(cellIndex) - classicalDofsPerCell; // saved as cell->quad->dim->node - dftefe::utils::MemoryStorage + dftefe::utils::MemoryStorage basisGradInCell = cfeBasisDataStorage->getBasisGradientDataInCell(cellIndex); @@ -282,8 +279,7 @@ namespace dftefe std::vector &cellStartIdsBasisGradientQuadStorage, std::vector &cellStartIdsBasisHessianQuadStorage, const BasisStorageAttributesBoolMap basisStorageAttributesBoolMap, - std::shared_ptr< - FEBasisDataStorage> + std::shared_ptr> cfeBasisDataStorage = nullptr) { bool numCellsZero = efeBDH->nLocallyOwnedCells() == 0 ? true : false; @@ -574,138 +570,6 @@ namespace dftefe cfeBasisDataStorage, classicalComponentInQuadGradients); } - - /** - cfeBasisManager->getCellDofsLocalIds(cellIndex, - vecClassicalLocalNodeId); - - std::vector coeffsInCell( - classicalDofsPerCell * numEnrichmentIdsInCell, 0); - - for (size_type cellEnrichId = 0; - cellEnrichId < numEnrichmentIdsInCell; - cellEnrichId++) - { - // get the enrichmentIds - global_size_type enrichmentId = - efeBDH->getEnrichmentClassicalInterface() - ->getEnrichmentId(cellIndex, cellEnrichId); - - // get the vectors of non-zero localIds and coeffs - auto iter = - enrichmentIdToInterfaceCoeffMap->find(enrichmentId); - auto it = - enrichmentIdToClassicalLocalIdMap->find(enrichmentId); - if (iter != enrichmentIdToInterfaceCoeffMap->end() && - it != enrichmentIdToClassicalLocalIdMap->end()) - { - const std::vector - &coeffsInLocalIdsMap = iter->second; - - for (size_type i = 0; i < classicalDofsPerCell; i++) - { - size_type pos = 0; - bool found = false; - it->second.getPosition(vecClassicalLocalNodeId[i], - pos, - found); - if (found) - { - coeffsInCell[numEnrichmentIdsInCell * i + - cellEnrichId] = - coeffsInLocalIdsMap[pos]; - } - } - } - } - - if (basisStorageAttributesBoolMap - .find(BasisStorageAttributes::StoreValues) - ->second || - basisStorageAttributesBoolMap - .find(BasisStorageAttributes::StoreOverlap) - ->second) - { - dftefe::utils::MemoryStorage - basisValInCell = - cfeBasisDataStorage->getBasisDataInCell(cellIndex); - - // Do a gemm (\Sigma c_i N_i^classical) - // and get the quad values in std::vector - - linearAlgebra::blasLapack::gemm( - linearAlgebra::blasLapack::Layout::ColMajor, - linearAlgebra::blasLapack::Op::NoTrans, - linearAlgebra::blasLapack::Op::Trans, - numEnrichmentIdsInCell, - nQuadPointInCell, - classicalDofsPerCell, - (ValueTypeBasisData)1.0, - coeffsInCell.data(), - numEnrichmentIdsInCell, - basisValInCell.data(), - nQuadPointInCell, - (ValueTypeBasisData)0.0, - classicalComponentInQuadValues.data(), - numEnrichmentIdsInCell, - *efeBDH->getEnrichmentClassicalInterface() - ->getLinAlgOpContext()); - } - - if (basisStorageAttributesBoolMap - .find(BasisStorageAttributes::StoreGradient) - ->second) - { - // saved as cell->dim->node->quad - dftefe::utils::MemoryStorage - basisGradInCell = - cfeBasisDataStorage->getBasisGradientDataInCell( - cellIndex); - - // Do a gemm (\Sigma c_i N_i^classical) - // and get the quad values in std::vector - - for (size_type iDim = 0; iDim < dim; iDim++) - { - ValueTypeBasisData *B = - basisGradInCell.data() + - iDim * nQuadPointInCell * classicalDofsPerCell; - - std::vector - classicalComponentInQuadDim(nQuadPointInCell * - numEnrichmentIdsInCell, - (ValueTypeBasisData)0); - - linearAlgebra::blasLapack::gemm< - ValueTypeBasisData, - ValueTypeBasisData, - utils::MemorySpace::HOST>( - linearAlgebra::blasLapack::Layout::ColMajor, - linearAlgebra::blasLapack::Op::NoTrans, - linearAlgebra::blasLapack::Op::Trans, - numEnrichmentIdsInCell, - nQuadPointInCell, - classicalDofsPerCell, - (ValueTypeBasisData)1.0, - coeffsInCell.data(), - numEnrichmentIdsInCell, - B, - nQuadPointInCell, - (ValueTypeBasisData)0.0, - classicalComponentInQuadDim.data(), - numEnrichmentIdsInCell, - *efeBDH->getEnrichmentClassicalInterface() - ->getLinAlgOpContext()); - - classicalComponentInQuadGradients[iDim] = - classicalComponentInQuadDim; - } - } - */ } if (basisStorageAttributesBoolMap @@ -1031,8 +895,7 @@ namespace dftefe const quadrature::QuadratureRuleAttributes &quadratureRuleAttributes, std::shared_ptr quadratureRuleContainer, - std::shared_ptr< - FEBasisDataStorage> + std::shared_ptr> cfeBasisDataStorage = nullptr) { @@ -1480,8 +1343,7 @@ namespace dftefe std::vector &cellStartIdsBasisGradientQuadStorage, std::vector &cellStartIdsBasisHessianQuadStorage, const BasisStorageAttributesBoolMap basisStorageAttributesBoolMap, - std::shared_ptr< - FEBasisDataStorage> + std::shared_ptr> cfeBasisDataStorage = nullptr) { // Assert if QuadFamily is not having variable quadpoints in each cell @@ -2267,8 +2129,7 @@ namespace dftefe const quadrature::QuadratureRuleAttributes &quadratureRuleAttributes, std::shared_ptr quadratureRuleContainer, - std::shared_ptr< - FEBasisDataStorage> + std::shared_ptr> cfeBasisDataStorage = nullptr) { const quadrature::QuadratureFamily quadratureFamily = @@ -2806,15 +2667,15 @@ namespace dftefe // Set up the FE Basis Data Storage // In HOST !! - cfeBasisDataStorage = std::make_shared< - CFEBasisDataStorageDealii>( - d_efeBDH->getEnrichmentClassicalInterface() - ->getCFEBasisDofHandler(), - quadratureRuleAttributes, - basisAttrMap); + cfeBasisDataStorage = + std::make_shared>( + d_efeBDH->getEnrichmentClassicalInterface() + ->getCFEBasisDofHandler(), + quadratureRuleAttributes, + basisAttrMap); cfeBasisDataStorage->evaluateBasisData(quadratureRuleAttributes, basisAttrMap); @@ -2836,15 +2697,15 @@ namespace dftefe // Set up the FE Basis Data Storage - cfeBasisDataStorage = std::make_shared< - CFEBasisDataStorageDealii>( - d_efeBDH->getEnrichmentClassicalInterface() - ->getCFEBasisDofHandler(), - quadratureRuleAttributes, - basisAttrMap); + cfeBasisDataStorage = + std::make_shared>( + d_efeBDH->getEnrichmentClassicalInterface() + ->getCFEBasisDofHandler(), + quadratureRuleAttributes, + basisAttrMap); cfeBasisDataStorage->evaluateBasisData(quadratureRuleAttributes, basisAttrMap); @@ -2861,15 +2722,15 @@ namespace dftefe // Set up the FE Basis Data Storage - cfeBasisDataStorage = std::make_shared< - CFEBasisDataStorageDealii>( - d_efeBDH->getEnrichmentClassicalInterface() - ->getCFEBasisDofHandler(), - quadratureRuleAttributes, - basisAttrMap); + cfeBasisDataStorage = + std::make_shared>( + d_efeBDH->getEnrichmentClassicalInterface() + ->getCFEBasisDofHandler(), + quadratureRuleAttributes, + basisAttrMap); cfeBasisDataStorage->evaluateBasisData(quadratureRuleAttributes, basisAttrMap); @@ -3049,15 +2910,15 @@ namespace dftefe // Set up the FE Basis Data Storage - cfeBasisDataStorage = std::make_shared< - CFEBasisDataStorageDealii>( - d_efeBDH->getEnrichmentClassicalInterface() - ->getCFEBasisDofHandler(), - quadratureRuleAttributes, - basisAttrMap); + cfeBasisDataStorage = + std::make_shared>( + d_efeBDH->getEnrichmentClassicalInterface() + ->getCFEBasisDofHandler(), + quadratureRuleAttributes, + basisAttrMap); cfeBasisDataStorage->evaluateBasisData(quadratureRuleAttributes, d_quadratureRuleContainer, @@ -3080,15 +2941,15 @@ namespace dftefe // Set up the FE Basis Data Storage - cfeBasisDataStorage = std::make_shared< - CFEBasisDataStorageDealii>( - d_efeBDH->getEnrichmentClassicalInterface() - ->getCFEBasisDofHandler(), - quadratureRuleAttributes, - basisAttrMap); + cfeBasisDataStorage = + std::make_shared>( + d_efeBDH->getEnrichmentClassicalInterface() + ->getCFEBasisDofHandler(), + quadratureRuleAttributes, + basisAttrMap); cfeBasisDataStorage->evaluateBasisData(quadratureRuleAttributes, d_quadratureRuleContainer, @@ -3106,15 +2967,15 @@ namespace dftefe // Set up the FE Basis Data Storage - cfeBasisDataStorage = std::make_shared< - CFEBasisDataStorageDealii>( - d_efeBDH->getEnrichmentClassicalInterface() - ->getCFEBasisDofHandler(), - quadratureRuleAttributes, - basisAttrMap); + cfeBasisDataStorage = + std::make_shared>( + d_efeBDH->getEnrichmentClassicalInterface() + ->getCFEBasisDofHandler(), + quadratureRuleAttributes, + basisAttrMap); cfeBasisDataStorage->evaluateBasisData(quadratureRuleAttributes, d_quadratureRuleContainer, @@ -3334,15 +3195,15 @@ namespace dftefe // Set up the FE Basis Data Storage - cfeBasisDataStorage = std::make_shared< - CFEBasisDataStorageDealii>( - d_efeBDH->getEnrichmentClassicalInterface() - ->getCFEBasisDofHandler(), - quadratureRuleAttributes, - basisAttrMap); + cfeBasisDataStorage = + std::make_shared>( + d_efeBDH->getEnrichmentClassicalInterface() + ->getCFEBasisDofHandler(), + quadratureRuleAttributes, + basisAttrMap); cfeBasisDataStorage->evaluateBasisData(quadratureRuleAttributes, d_quadratureRuleContainer, @@ -3365,15 +3226,15 @@ namespace dftefe // Set up the FE Basis Data Storage - cfeBasisDataStorage = std::make_shared< - CFEBasisDataStorageDealii>( - d_efeBDH->getEnrichmentClassicalInterface() - ->getCFEBasisDofHandler(), - quadratureRuleAttributes, - basisAttrMap); + cfeBasisDataStorage = + std::make_shared>( + d_efeBDH->getEnrichmentClassicalInterface() + ->getCFEBasisDofHandler(), + quadratureRuleAttributes, + basisAttrMap); cfeBasisDataStorage->evaluateBasisData(quadratureRuleAttributes, d_quadratureRuleContainer, @@ -3391,15 +3252,15 @@ namespace dftefe // Set up the FE Basis Data Storage - cfeBasisDataStorage = std::make_shared< - CFEBasisDataStorageDealii>( - d_efeBDH->getEnrichmentClassicalInterface() - ->getCFEBasisDofHandler(), - quadratureRuleAttributes, - basisAttrMap); + cfeBasisDataStorage = + std::make_shared>( + d_efeBDH->getEnrichmentClassicalInterface() + ->getCFEBasisDofHandler(), + quadratureRuleAttributes, + basisAttrMap); cfeBasisDataStorage->evaluateBasisData(quadratureRuleAttributes, d_quadratureRuleContainer, @@ -3601,15 +3462,15 @@ namespace dftefe // Set up the FE Basis Data Storage - cfeBasisDataStorage = std::make_shared< - CFEBasisDataStorageDealii>( - d_efeBDH->getEnrichmentClassicalInterface() - ->getCFEBasisDofHandler(), - quadratureRuleAttributes, - basisAttrMap); + cfeBasisDataStorage = + std::make_shared>( + d_efeBDH->getEnrichmentClassicalInterface() + ->getCFEBasisDofHandler(), + quadratureRuleAttributes, + basisAttrMap); cfeBasisDataStorage->evaluateBasisData(quadratureRuleAttributes, d_quadratureRuleContainer, @@ -3632,15 +3493,15 @@ namespace dftefe // Set up the FE Basis Data Storage - cfeBasisDataStorage = std::make_shared< - CFEBasisDataStorageDealii>( - d_efeBDH->getEnrichmentClassicalInterface() - ->getCFEBasisDofHandler(), - quadratureRuleAttributes, - basisAttrMap); + cfeBasisDataStorage = + std::make_shared>( + d_efeBDH->getEnrichmentClassicalInterface() + ->getCFEBasisDofHandler(), + quadratureRuleAttributes, + basisAttrMap); cfeBasisDataStorage->evaluateBasisData(quadratureRuleAttributes, d_quadratureRuleContainer, @@ -3658,15 +3519,15 @@ namespace dftefe // Set up the FE Basis Data Storage - cfeBasisDataStorage = std::make_shared< - CFEBasisDataStorageDealii>( - d_efeBDH->getEnrichmentClassicalInterface() - ->getCFEBasisDofHandler(), - quadratureRuleAttributes, - basisAttrMap); + cfeBasisDataStorage = + std::make_shared>( + d_efeBDH->getEnrichmentClassicalInterface() + ->getCFEBasisDofHandler(), + quadratureRuleAttributes, + basisAttrMap); cfeBasisDataStorage->evaluateBasisData(quadratureRuleAttributes, d_quadratureRuleContainer, diff --git a/src/basis/OrthoEFEOverlapInverseOpContextGLL.t.cpp b/src/basis/OrthoEFEOverlapInverseOpContextGLL.t.cpp index cca2639b..c7f063ba 100644 --- a/src/basis/OrthoEFEOverlapInverseOpContextGLL.t.cpp +++ b/src/basis/OrthoEFEOverlapInverseOpContextGLL.t.cpp @@ -417,15 +417,14 @@ namespace dftefe } } - utils::MemoryStorage + utils::MemoryStorage basisValInCellEC = enrichmentBlockClassicalBasisDataStorage.getBasisDataInCell( cellIndex); // Do a gemm (\Sigma c_i N_i^classical) // and get the quad values in std::vector - + ValueTypeOperator *B = basisValInCellEC.data(); linearAlgebra::blasLapack::gemm( @@ -438,21 +437,20 @@ namespace dftefe (ValueTypeOperator)1.0, coeffsInCell.data(), numEnrichmentIdsInCell, - basisValInCellEC.data(), + B, dofsPerCellCFE, (ValueTypeOperator)0.0, classicalComponentInQuadValuesEC.data(), numEnrichmentIdsInCell, *eci->getLinAlgOpContext()); - utils::MemoryStorage + utils::MemoryStorage basisValInCellEE = enrichmentBlockEnrichmentBasisDataStorage .getBasisDataInCell(cellIndex); // Do a gemm (\Sigma c_i N_i^classical) // and get the quad values in std::vector - + B = basisValInCellEE.data(); linearAlgebra::blasLapack::gemm( @@ -465,7 +463,7 @@ namespace dftefe (ValueTypeOperator)1.0, coeffsInCell.data(), numEnrichmentIdsInCell, - basisValInCellEE.data(), + B, dofsPerCell, (ValueTypeOperator)0.0, classicalComponentInQuadValuesEE.data(), diff --git a/src/basis/OrthoEFEOverlapOperatorContext.t.cpp b/src/basis/OrthoEFEOverlapOperatorContext.t.cpp index 5bffa6a9..ce1d8770 100644 --- a/src/basis/OrthoEFEOverlapOperatorContext.t.cpp +++ b/src/basis/OrthoEFEOverlapOperatorContext.t.cpp @@ -517,12 +517,12 @@ namespace dftefe } } - dftefe::utils::MemoryStorage + dftefe::utils::MemoryStorage basisValInCellEC = enrichmentBlockClassicalBasisDataStorage.getBasisDataInCell( cellIndex); + ValueTypeOperator *B = basisValInCellEC.data(); // Do a gemm (\Sigma c_i N_i^classical) // and get the quad values in std::vector @@ -538,7 +538,7 @@ namespace dftefe (ValueTypeOperator)1.0, coeffsInCell.data(), numEnrichmentIdsInCell, - basisValInCellEC.data(), + B, dofsPerCellCFE, (ValueTypeOperator)0.0, classicalComponentInQuadValuesEC.data(), @@ -546,13 +546,13 @@ namespace dftefe *eci->getLinAlgOpContext()); - dftefe::utils::MemoryStorage + dftefe::utils::MemoryStorage basisValInCellEE = enrichmentBlockEnrichmentBasisDataStorage .getBasisDataInCell(cellIndex); // Do a gemm (\Sigma c_i N_i^classical) // and get the quad values in std::vector + B = basisValInCellEE.data(); linearAlgebra::blasLapack::gemm #include #include - +using namespace dftefe; int main() { @@ -41,45 +41,45 @@ int main() // NOTE : The test case only works for orthogonalized EFE basis int mpiInitFlag = 0; - dftefe::utils::mpi::MPIInitialized(&mpiInitFlag); + utils::mpi::MPIInitialized(&mpiInitFlag); if(!mpiInitFlag) { - dftefe::utils::mpi::MPIInit(NULL, NULL); + utils::mpi::MPIInit(NULL, NULL); } - dftefe::utils::mpi::MPIComm comm = dftefe::utils::mpi::MPICommWorld; + utils::mpi::MPIComm comm = utils::mpi::MPICommWorld; // Get the rank of the process int rank; - dftefe::utils::mpi::MPICommRank(comm, &rank); + utils::mpi::MPICommRank(comm, &rank); // Get nProcs int numProcs; - dftefe::utils::mpi::MPICommSize(comm, &numProcs); + utils::mpi::MPICommSize(comm, &numProcs); int blasQueue = 0; int lapackQueue = 0; - std::shared_ptr> blasQueuePtr = std::make_shared - >(blasQueue); - std::shared_ptr> lapackQueuePtr = std::make_shared - >(lapackQueue); - std::shared_ptr> linAlgOpContext = - std::make_shared>(blasQueuePtr, lapackQueuePtr); + std::shared_ptr> blasQueuePtr = std::make_shared + >(blasQueue); + std::shared_ptr> lapackQueuePtr = std::make_shared + >(lapackQueue); + std::shared_ptr> linAlgOpContext = + std::make_shared>(blasQueuePtr, lapackQueuePtr); // Set up Triangulation const unsigned int dim = 3; - std::shared_ptr triangulationBase = - std::make_shared>(comm); - std::vector subdivisions = {10, 10, 10}; + std::shared_ptr triangulationBase = + std::make_shared>(comm); + std::vector subdivisions = {5, 5, 5}; std::vector isPeriodicFlags(dim, false); - std::vector domainVectors(dim, - dftefe::utils::Point(dim, 0.0)); + std::vector domainVectors(dim, + utils::Point(dim, 0.0)); double xmax = 10.; double ymax = 10.; @@ -107,7 +107,7 @@ int main() fstream.open(inputFileName, std::fstream::in); // read the input file and create atomSymbolVec vector and atom coordinates vector. - std::vector atomCoordinatesVec; + std::vector atomCoordinatesVec; std::vector coordinates; coordinates.resize(dim,0.); std::vector atomSymbolVec; @@ -123,7 +123,7 @@ int main() atomCoordinatesVec.push_back(coordinates); atomSymbolVec.push_back(symbol); } - dftefe::utils::mpi::MPIBarrier(comm); + utils::mpi::MPIBarrier(comm); std::map atomSymbolToFilename; for (auto i:atomSymbolVec ) @@ -133,8 +133,8 @@ int main() std::vector fieldNames{"vnuclear"}; std::vector metadataNames{ "symbol", "Z", "charge", "NR", "r" }; - std::shared_ptr atomSphericalDataContainer = - std::make_shared(atomSymbolToFilename, + std::shared_ptr atomSphericalDataContainer = + std::make_shared(atomSymbolToFilename, fieldNames, metadataNames); @@ -153,7 +153,7 @@ int main() for( ; triaCellIter != triangulationBase->endLocal(); triaCellIter++) { refineFlag = false; - dftefe::utils::Point centerPoint(dim, 0.0); + utils::Point centerPoint(dim, 0.0); (*triaCellIter)->center(centerPoint); double dist = (centerPoint[0] - 5)* (centerPoint[0] - 5); dist += (centerPoint[1] - 5)* (centerPoint[1] - 5); @@ -170,16 +170,16 @@ int main() triangulationBase->executeCoarseningAndRefinement(); triangulationBase->finalizeTriangulationConstruction(); // Mpi_allreduce that all the flags are 1 (mpi_max) - int err = dftefe::utils::mpi::MPIAllreduce( + int err = utils::mpi::MPIAllreduce( &flag, &mpiReducedFlag, 1, - dftefe::utils::mpi::MPIInt, - dftefe::utils::mpi::MPIMin, + utils::mpi::MPIInt, + utils::mpi::MPIMin, comm); std::pair mpiIsSuccessAndMsg = - dftefe::utils::mpi::MPIErrIsSuccessAndMsg(err); - dftefe::utils::throwException(mpiIsSuccessAndMsg.first, + utils::mpi::MPIErrIsSuccessAndMsg(err); + utils::throwException(mpiIsSuccessAndMsg.first, "MPI Error:" + mpiIsSuccessAndMsg.second); } } @@ -188,14 +188,14 @@ int main() unsigned int feOrder = 3; // Set up the vector of scalarSpatialRealFunctions for adaptive quadrature - std::vector> functionsVec(0); + std::vector> functionsVec(0); unsigned int numfun = 2; functionsVec.resize(numfun); // Enrichment Functions std::vector tolerances(numfun); std::vector integralThresholds(numfun); for ( unsigned int i=0 ;i < functionsVec.size() ; i++ ) { - functionsVec[i] = std::make_shared>( + functionsVec[i] = std::make_shared>( atomSphericalDataContainer, atomSymbolVec, atomCoordinatesVec, @@ -210,21 +210,21 @@ int main() //Set up quadAttr for Rhs and OverlapMatrix - dftefe::quadrature::QuadratureRuleAttributes quadAttrAdaptive(dftefe::quadrature::QuadratureFamily::ADAPTIVE,false); + quadrature::QuadratureRuleAttributes quadAttrAdaptive(quadrature::QuadratureFamily::ADAPTIVE,false); // Set up base quadrature rule for adaptive quadrature unsigned int num1DGaussSize = feOrder + 1; - std::shared_ptr baseQuadRule = - std::make_shared(dim, num1DGaussSize); + std::shared_ptr baseQuadRule = + std::make_shared(dim, num1DGaussSize); - dftefe::quadrature::QuadratureRuleAttributes quadAttrGauss(dftefe::quadrature::QuadratureFamily::GAUSS,true,num1DGaussSize); + quadrature::QuadratureRuleAttributes quadAttrGauss(quadrature::QuadratureFamily::GAUSS,true,num1DGaussSize); - std::shared_ptr cellMapping = std::make_shared>(); - std::shared_ptr parentToChildCellsManager = std::make_shared>(); + std::shared_ptr cellMapping = std::make_shared>(); + std::shared_ptr parentToChildCellsManager = std::make_shared>(); - std::shared_ptr quadRuleContainerAdaptive = - std::make_shared + std::shared_ptr quadRuleContainerAdaptive = + std::make_shared (quadAttrAdaptive, baseQuadRule, triangulationBase, @@ -238,20 +238,20 @@ int main() maxRecursion); // Set the CFE basis manager and handler for bassiInterfaceCoeffcient distributed vector - std::shared_ptr> cfeBasisDofHandler = - std::make_shared>(triangulationBase, feOrder, comm); + std::shared_ptr> cfeBasisDofHandler = + std::make_shared>(triangulationBase, feOrder, comm); - dftefe::basis::BasisStorageAttributesBoolMap basisAttrMap; - basisAttrMap[dftefe::basis::BasisStorageAttributes::StoreValues] = true; - basisAttrMap[dftefe::basis::BasisStorageAttributes::StoreGradient] = true; - basisAttrMap[dftefe::basis::BasisStorageAttributes::StoreHessian] = false; - basisAttrMap[dftefe::basis::BasisStorageAttributes::StoreOverlap] = false; - basisAttrMap[dftefe::basis::BasisStorageAttributes::StoreGradNiGradNj] = false; - basisAttrMap[dftefe::basis::BasisStorageAttributes::StoreJxW] = true; + basis::BasisStorageAttributesBoolMap basisAttrMap; + basisAttrMap[basis::BasisStorageAttributes::StoreValues] = true; + basisAttrMap[basis::BasisStorageAttributes::StoreGradient] = true; + basisAttrMap[basis::BasisStorageAttributes::StoreHessian] = false; + basisAttrMap[basis::BasisStorageAttributes::StoreOverlap] = false; + basisAttrMap[basis::BasisStorageAttributes::StoreGradNiGradNj] = false; + basisAttrMap[basis::BasisStorageAttributes::StoreJxW] = false; // Set up the CFE Basis Data Storage for Rhs - std::shared_ptr> cfeBasisDataStorageAdaptUniformQuad = - std::make_shared> + std::shared_ptr> cfeBasisDataStorageAdaptUniformQuad = + std::make_shared> (cfeBasisDofHandler, quadAttrAdaptive, basisAttrMap); // evaluate basis data cfeBasisDataStorageAdaptUniformQuad->evaluateBasisData(quadAttrAdaptive, quadRuleContainerAdaptive, basisAttrMap); @@ -259,78 +259,155 @@ int main() unsigned int cellBlockSize = 10; // Set up the CFE Basis Data Storage for Rhs - std::shared_ptr> cfeBasisDataStorageUniformQuad = - std::make_shared> - (cfeBasisDofHandler, quadAttrGauss, basisAttrMap, cellBlockSize ,*linAlgOpContext); - // evaluate basis data + std::shared_ptr> cfeBasisDataStorageUniformQuad = + std::make_shared> + (cfeBasisDofHandler, quadAttrGauss, basisAttrMap, cellBlockSize, *linAlgOpContext); + // check basis data cfeBasisDataStorageUniformQuad->evaluateBasisData(quadAttrGauss, basisAttrMap); size_type numLocallyOwnedCells = cfeBasisDofHandler->nLocallyOwnedCells(); + size_type nDofs = cfeBasisDofHandler->nCellDofs(0); + size_type nQuadPtInCell = cfeBasisDataStorageUniformQuad->getQuadratureRuleContainer()->nCellQuadraturePoints(0); + for(int i = 0 ; i < numLocallyOwnedCells ; i += cellBlockSize) + { + size_type cellStartId = i; + size_type cellEndId = std::min(i+cellBlockSize, numLocallyOwnedCells); + utils::MemoryStorage + basisDataStorageInCellRange((cellEndId - cellStartId)*nDofs*nQuadPtInCell); + cfeBasisDataStorageUniformQuad->getBasisDataInCellRange(std::make_pair(cellStartId, cellEndId), basisDataStorageInCellRange); + for(int j = 0 ; j < cellEndId - cellStartId ; j++) + { + auto basisDataStorageInCellAdaptUniform = cfeBasisDataStorageAdaptUniformQuad->getBasisDataInCell(j); + for(int k = 0 ; k < basisDataStorageInCellAdaptUniform.size() ; k++) + { + double val1 = *(basisDataStorageInCellRange.data() + k); + double val2 = *(basisDataStorageInCellAdaptUniform.data() + k); + if(std::abs(val1 - val2) > 1e-12) + std::cout << val1 << "\t" << val2 << "\n"; + } + } + } + // check basis gradient data + for(int i = 0 ; i < numLocallyOwnedCells ; i += cellBlockSize) + { + size_type cellStartId = i; + size_type cellEndId = std::min(i+cellBlockSize, numLocallyOwnedCells); + utils::MemoryStorage + basisGradientDataStorageInCellRange((cellEndId - cellStartId)*nDofs*nQuadPtInCell*dim); + cfeBasisDataStorageUniformQuad->getBasisGradientDataInCellRange(std::make_pair(cellStartId, cellEndId), basisGradientDataStorageInCellRange); + for(int j = 0 ; j < cellEndId - cellStartId ; j++) + { + auto basisGradientDataStorageInCellAdaptUniform = cfeBasisDataStorageAdaptUniformQuad->getBasisGradientDataInCell(j); + for(int k = 0 ; k < basisGradientDataStorageInCellAdaptUniform.size() ; k++) + { + double val1 = *(basisGradientDataStorageInCellRange.data() + k); + double val2 = *(basisGradientDataStorageInCellAdaptUniform.data() + k); + if(std::abs(val1 - val2) > 1e-12) + std::cout << val1 << "\t" << val2 << "\n"; + } + } + } + + // Create the enrichmentClassicalInterface object + std::shared_ptr> + enrichClassIntfce = std::make_shared> + (cfeBasisDataStorageUniformQuad, + cfeBasisDataStorageUniformQuad, + atomSphericalDataContainer, + atomPartitionTolerance, + atomSymbolVec, + atomCoordinatesVec, + fieldName, + linAlgOpContext, + comm); + + // initialize the basis + std::shared_ptr> basisDofHandler = + std::make_shared>( + enrichClassIntfce, comm); + + std::map dofCoords; + basisDofHandler->getBasisCenters(dofCoords); + + std::cout << "Locally owned cells : " <nLocallyOwnedCells() << "\n"; + std::cout << "Total Number of dofs : " << basisDofHandler->nGlobalNodes() << "\n"; + + basisAttrMap[basis::BasisStorageAttributes::StoreValues] = true; + basisAttrMap[basis::BasisStorageAttributes::StoreGradient] = true; + basisAttrMap[basis::BasisStorageAttributes::StoreHessian] = false; + basisAttrMap[basis::BasisStorageAttributes::StoreOverlap] = false; + basisAttrMap[basis::BasisStorageAttributes::StoreGradNiGradNj] = false; + basisAttrMap[basis::BasisStorageAttributes::StoreJxW] = false; + + // Set up Adaptive quadrature for EFE Basis Data Storage + std::shared_ptr> efeBasisDataStorageAdaptUniformQuad = + std::make_shared> + (basisDofHandler, quadAttrAdaptive, basisAttrMap); + + // evaluate basis data + efeBasisDataStorageAdaptUniformQuad->evaluateBasisData(quadAttrAdaptive, quadRuleContainerAdaptive, basisAttrMap); + + // Set up the CFE Basis Data Storage for Rhs + std::shared_ptr> efeBasisDataStorageUniformQuad = + std::make_shared> + (basisDofHandler, quadAttrGauss, basisAttrMap, cellBlockSize, *linAlgOpContext); + // check basis data + efeBasisDataStorageUniformQuad->evaluateBasisData(quadAttrGauss, basisAttrMap); + numLocallyOwnedCells = basisDofHandler->nLocallyOwnedCells(); + nDofs = basisDofHandler->nCellDofs(0); + nQuadPtInCell = efeBasisDataStorageUniformQuad->getQuadratureRuleContainer()->nCellQuadraturePoints(0); for(int i = 0 ; i < numLocallyOwnedCells ; i += cellBlockSize) { - size_type cellStartId = i; - size_type cellEndId = std::min(i+cellBlockSize, numLocallyOwnedCells); - auto basisDataStorageInCellRange = - cfeBasisDataStorageUniformQuad->getBasisDataStorageInCellRange(std::make_pair(cellStartId, cellEndId)); - for(int j = 0 ; j < cellEndId - cellStartId ; j++) + size_type cellStartId = i; + size_type cellEndId = std::min(i+cellBlockSize, numLocallyOwnedCells); + utils::MemoryStorage + basisDataStorageInCellRange((cellEndId - cellStartId)*nDofs*nQuadPtInCell); + efeBasisDataStorageUniformQuad->getBasisDataInCellRange(std::make_pair(cellStartId, cellEndId), basisDataStorageInCellRange); + for(int j = 0 ; j < cellEndId - cellStartId ; j++) + { + auto basisDataStorageInCellAdaptUniform = efeBasisDataStorageAdaptUniformQuad->getBasisDataInCell(j); + for(int k = 0 ; k < basisDataStorageInCellAdaptUniform.size() ; k++) { - auto basisDataStorageInCellAdaptUniform = cfeBasisDataStorageAdaptUniformQuad->getBasisDataInCell(j); - for(int k = 0 ; k < basisDataStorageInCellAdaptUniform.size() ; k++) - { - std::cout << *(basisDataStorageInCellRange.data() + dofsPerCell * nQuadPerCell * j + k) << - "\t" << *(basisDataStorageInCellAdaptUniform.data() + k) << "\n"; - } + double val1 = *(basisDataStorageInCellRange.data() + k); + double val2 = *(basisDataStorageInCellAdaptUniform.data() + k); + if(std::abs(val1 - val2) > 1e-12) + std::cout << val1 << "\t" << val2 << "\n"; } + } + } + // check basis gradient data + for(int i = 0 ; i < numLocallyOwnedCells ; i += cellBlockSize) + { + size_type cellStartId = i; + size_type cellEndId = std::min(i+cellBlockSize, numLocallyOwnedCells); + utils::MemoryStorage + basisGradientDataStorageInCellRange((cellEndId - cellStartId)*nDofs*nQuadPtInCell*dim); + efeBasisDataStorageUniformQuad->getBasisGradientDataInCellRange(std::make_pair(cellStartId, cellEndId), basisGradientDataStorageInCellRange); + for(int j = 0 ; j < cellEndId - cellStartId ; j++) + { + auto basisGradientDataStorageInCellAdaptUniform = efeBasisDataStorageAdaptUniformQuad->getBasisGradientDataInCell(j); + for(int k = 0 ; k < basisGradientDataStorageInCellAdaptUniform.size() ; k++) + { + double val1 = *(basisGradientDataStorageInCellRange.data() + k); + double val2 = *(basisGradientDataStorageInCellAdaptUniform.data() + k); + if(std::abs(val1 - val2) > 1e-12) + std::cout << val1 << "\t" << val2 << "\n"; + } + } } -// // Create the enrichmentClassicalInterface object -// std::shared_ptr> -// (cfeBasisDataStorageOverlapMatrix, -// cfeBasisDataStorageRhs, -// atomSphericalDataContainer, -// atomPartitionTolerance, -// atomSymbolVec, -// atomCoordinatesVec, -// fieldName, -// linAlgOpContext, -// comm); - -// // initialize the basis -// std::shared_ptr> basisDofHandler = -// std::make_shared>( -// enrichClassIntfce, feOrder, comm); - -// std::map dofCoords; -// basisDofHandler->getBasisCenters(dofCoords); - -// std::cout << "Locally owned cells : " <nLocallyOwnedCells() << "\n"; -// std::cout << "Total Number of dofs : " << basisDofHandler->nGlobalNodes() << "\n"; - -// basisAttrMap[dftefe::basis::BasisStorageAttributes::StoreValues] = true; -// basisAttrMap[dftefe::basis::BasisStorageAttributes::StoreGradient] = false; -// basisAttrMap[dftefe::basis::BasisStorageAttributes::StoreHessian] = false; -// basisAttrMap[dftefe::basis::BasisStorageAttributes::StoreOverlap] = false; -// basisAttrMap[dftefe::basis::BasisStorageAttributes::StoreGradNiGradNj] = false; -// basisAttrMap[dftefe::basis::BasisStorageAttributes::StoreJxW] = false; - -// // Set up Adaptive quadrature for EFE Basis Data Storage -// std::shared_ptr> feBasisData = -// std::make_shared> -// (basisDofHandler, quadAttrAdaptive, basisAttrMap); - -// // evaluate basis data -// feBasisData->evaluateBasisData(quadAttrAdaptive, quadRuleContainerAdaptive, basisAttrMap); - - - dftefe::utils::mpi::MPIBarrier(comm); + utils::mpi::MPIBarrier(comm); //gracefully end MPI int mpiFinalFlag = 0; - dftefe::utils::mpi::MPIFinalized(&mpiFinalFlag); + utils::mpi::MPIFinalized(&mpiFinalFlag); if(!mpiFinalFlag) { - dftefe::utils::mpi::MPIFinalize(); + utils::mpi::MPIFinalize(); } }