From cba65e6ea3a6a8bb632b815f903a8debf785f40c Mon Sep 17 00:00:00 2001 From: Avirup Sircar Date: Thu, 7 Nov 2024 10:06:28 -0500 Subject: [PATCH] Removed bugs in EFEBasisOnTheFlyComputeDealii --- src/basis/CFEBDSOnTheFlyComputeDealii.t.cpp | 3 +- src/basis/EFEBDSOnTheFlyComputeDealii.h | 14 +- src/basis/EFEBDSOnTheFlyComputeDealii.t.cpp | 313 ++++++++++-------- test/basis/src/TestBasisDataStorageMemOpt.cpp | 83 +++-- 4 files changed, 242 insertions(+), 171 deletions(-) diff --git a/src/basis/CFEBDSOnTheFlyComputeDealii.t.cpp b/src/basis/CFEBDSOnTheFlyComputeDealii.t.cpp index 9ca2549f..f648b3ee 100644 --- a/src/basis/CFEBDSOnTheFlyComputeDealii.t.cpp +++ b/src/basis/CFEBDSOnTheFlyComputeDealii.t.cpp @@ -339,7 +339,8 @@ namespace dftefe .find(BasisStorageAttributes::StoreHessian) ->second) { - std::cout << "Store Hessian is not memory optimized in CFEOnTheFlyComputeDealii.h. Contact developers for making it optimal."; + std::cout + << "Store Hessian is not memory optimized in CFEOnTheFlyComputeDealii.h. Contact developers for making it optimal."; basisHessianQuadStorage = std::make_shared::Storage>( diff --git a/src/basis/EFEBDSOnTheFlyComputeDealii.h b/src/basis/EFEBDSOnTheFlyComputeDealii.h index 0eede301..3b3776f4 100644 --- a/src/basis/EFEBDSOnTheFlyComputeDealii.h +++ b/src/basis/EFEBDSOnTheFlyComputeDealii.h @@ -214,8 +214,10 @@ namespace dftefe private: bool d_evaluateBasisData; - std::shared_ptr< - const EFEBasisDofHandlerDealii> + std::shared_ptr> d_efeBDH; std::shared_ptr d_quadratureRuleContainer; @@ -233,9 +235,11 @@ namespace dftefe size_type d_maxCellBlock; std::shared_ptr d_tmpGradientBlock; linearAlgebra::LinAlgOpContext &d_linAlgOpContext; - size_type d_classialDofsInCell; - std::unordered_map::Storage>> + size_type d_classialDofsInCell; + std::unordered_map< + size_type, + std::shared_ptr< + typename BasisDataStorage::Storage>> d_basisGradientEnrichQuadStorageMap, d_basisEnrichQuadStorageMap; }; // end of EFEBDSOnTheFlyComputeDealii diff --git a/src/basis/EFEBDSOnTheFlyComputeDealii.t.cpp b/src/basis/EFEBDSOnTheFlyComputeDealii.t.cpp index 51a13831..3a89b4e9 100644 --- a/src/basis/EFEBDSOnTheFlyComputeDealii.t.cpp +++ b/src/basis/EFEBDSOnTheFlyComputeDealii.t.cpp @@ -58,10 +58,9 @@ namespace dftefe typename BasisDataStorage::Storage &basisGradientData) { - size_type numCellsInBlock = cellRange.second - cellRange.first; - size_type numMats = 0; - for (size_type iCell = 0; iCell < numCellsInBlock; ++iCell) + for (size_type iCell = cellRange.first; iCell < cellRange.second; + ++iCell) { for (size_type iQuad = 0; iQuad < nQuadPointsInCell[iCell]; ++iQuad) { @@ -86,17 +85,19 @@ namespace dftefe std::vector strideBTmp(numMats, 0); std::vector strideCTmp(numMats, 0); - for (size_type iCell = 0; iCell < numCellsInBlock; ++iCell) + for (size_type iCell = cellRange.first; iCell < cellRange.second; + ++iCell) { for (size_type iQuad = 0; iQuad < nQuadPointsInCell[iCell]; ++iQuad) { - size_type index = iCell * nQuadPointsInCell[iCell] + iQuad; + size_type index = + (iCell - cellRange.first) * nQuadPointsInCell[iCell] + iQuad; mSizesTmp[index] = classicalDofsInCell; nSizesTmp[index] = dim; kSizesTmp[index] = dim; ldaSizesTmp[index] = mSizesTmp[index]; ldbSizesTmp[index] = kSizesTmp[index]; - ldcSizesTmp[index] = mSizesTmp[index]; + ldcSizesTmp[index] = dofsInCell[iCell]; strideATmp[index] = mSizesTmp[index] * kSizesTmp[index]; strideBTmp[index] = kSizesTmp[index] * nSizesTmp[index]; strideCTmp[index] = dofsInCell[iCell] * nSizesTmp[index]; @@ -357,21 +358,25 @@ namespace dftefe size_type dim> void storeValuesHRefinedSameQuadEveryCell( - std::shared_ptr< - const EFEBasisDofHandlerDealii> - efeBDH, + std::shared_ptr> efeBDH, std::shared_ptr< typename BasisDataStorage::Storage> &basisParaCellClassQuadStorage, - std::unordered_map::Storage>> + std::unordered_map< + size_type, + std::shared_ptr::Storage>> &basisEnrichQuadStorageMap, std::shared_ptr< typename BasisDataStorage::Storage> &basisGradientParaCellClassQuadStorage, - std::unordered_map::Storage>> + std::unordered_map< + size_type, + std::shared_ptr::Storage>> &basisGradientEnrichQuadStorageMap, std::shared_ptr< typename BasisDataStorage::Storage> @@ -458,8 +463,8 @@ namespace dftefe ->second) dealiiUpdateFlags |= dealii::update_hessians; - size_type classicalDofsPerCell = utils::mathFunctions::sizeTypePow( - (efeBDH->getFEOrder(0) + 1), dim); + size_type classicalDofsPerCell = + utils::mathFunctions::sizeTypePow((efeBDH->getFEOrder(0) + 1), dim); // NOTE: cellId 0 passed as we assume h-refine finite element mesh in // this function @@ -490,7 +495,7 @@ namespace dftefe const size_type numLocallyOwnedCells = efeBDH->nLocallyOwnedCells(); // NOTE: cellId 0 passed as we assume only H refined in this function - size_type dofsPerCell = efeBDH->nCellDofs(cellId); + size_type dofsPerCell = efeBDH->nCellDofs(cellId); const size_type nQuadPointInCell = numCellsZero ? 0 : quadratureRuleContainer->nCellQuadraturePoints(cellId); @@ -500,10 +505,11 @@ namespace dftefe nQuadPointsInCell.resize(numLocallyOwnedCells, nQuadPointInCell); std::vector basisParaCellClassQuadStorageTmp(0); std::vector basisJacobianInvQuadStorageTmp(0); - std::vector basisGradientParaCellClassQuadStorageTmp(0); + std::vector + basisGradientParaCellClassQuadStorageTmp(0); std::vector basisHessianQuadStorageTmp(0); - std::unordered_map> + std::unordered_map> basisEnrichQuadStorageMapTmp, basisGradientEnrichQuadStorageMapTmp; size_type cellIndex = 0; @@ -527,7 +533,9 @@ namespace dftefe std::make_shared::Storage>( classicalDofsPerCell * nQuadPointInCell); - basisParaCellClassQuadStorageTmp.resize(classicalDofsPerCell * nQuadPointInCell, ValueTypeBasisData(0)); + basisParaCellClassQuadStorageTmp.resize(classicalDofsPerCell * + nQuadPointInCell, + ValueTypeBasisData(0)); } if (basisStorageAttributesBoolMap @@ -544,7 +552,8 @@ namespace dftefe std::make_shared::Storage>( classicalDofsPerCell * nQuadPointInCell * dim); - basisGradientParaCellClassQuadStorageTmp.resize(classicalDofsPerCell * nQuadPointInCell * dim); + basisGradientParaCellClassQuadStorageTmp.resize( + classicalDofsPerCell * nQuadPointInCell * dim); cellStartIdsBasisJacobianInvQuadStorage.resize(numLocallyOwnedCells, 0); } @@ -552,7 +561,8 @@ namespace dftefe .find(BasisStorageAttributes::StoreHessian) ->second) { - std::cout << "Store Hessian is not memory optimized in CFEOnTheFlyComputeDealii.h. Contact developers for making it optimal."; + std::cout + << "Store Hessian is not memory optimized in CFEOnTheFlyComputeDealii.h. Contact developers for making it optimal."; basisHessianQuadStorage = std::make_shared::Storage>( @@ -574,7 +584,7 @@ namespace dftefe "Dynamic casting of FECellBase to FECellDealii not successful"); } - cellIndex = 0; + cellIndex = 0; size_type cumulativeQuadPointsxnDofs = 0; for (; locallyOwnedCellIter != efeBDH->endLocallyOwnedCells(); @@ -598,7 +608,8 @@ namespace dftefe std::vector classicalComponentInQuadValues(0); - std::vector classicalComponentInQuadGradients(0); + std::vector classicalComponentInQuadGradients( + 0); if (basisStorageAttributesBoolMap .find(BasisStorageAttributes::StoreValues) @@ -667,47 +678,47 @@ namespace dftefe .find(BasisStorageAttributes::StoreValues) ->second) { - if(locallyOwnedCellIter == efeBDH->beginLocallyOwnedCells()) - { - for (unsigned int iNode = 0; iNode < classicalDofsPerCell; iNode++) - { - for (unsigned int qPoint = 0; qPoint < nQuadPointInCell; - qPoint++) - { - auto it = - basisParaCellClassQuadStorageTmp.begin() + - qPoint * classicalDofsPerCell + - iNode; - *it = dealiiFEValues.shape_value(iNode, qPoint); - } - } - } + if (locallyOwnedCellIter == efeBDH->beginLocallyOwnedCells()) + { + for (unsigned int iNode = 0; iNode < classicalDofsPerCell; + iNode++) + { + for (unsigned int qPoint = 0; qPoint < nQuadPointInCell; + qPoint++) + { + auto it = basisParaCellClassQuadStorageTmp.begin() + + qPoint * classicalDofsPerCell + iNode; + *it = dealiiFEValues.shape_value(iNode, qPoint); + } + } + } - if(numEnrichmentIdsInCell > 0) - { - basisEnrichQuadStorageMapTmp[cellIndex].resize(nQuadPointInCell * numEnrichmentIdsInCell); - for (unsigned int iNode = 0; iNode < numEnrichmentIdsInCell; iNode++) + if (numEnrichmentIdsInCell > 0) { - for (unsigned int qPoint = 0; qPoint < nQuadPointInCell; - qPoint++) + basisEnrichQuadStorageMapTmp[cellIndex].resize( + nQuadPointInCell * numEnrichmentIdsInCell); + for (unsigned int iNode = 0; iNode < numEnrichmentIdsInCell; + iNode++) { - // std::cout << efeBDH->getEnrichmentValue( - // cellIndex, - // iNode, - // quadRealPointsVec[qPoint]) << " " << - // classicalComponentInQuadValues - // [numEnrichmentIdsInCell * qPoint + iNode] << "\n"; - *(basisEnrichQuadStorageMapTmp[cellIndex].data() - + qPoint * numEnrichmentIdsInCell + iNode) = - efeBDH->getEnrichmentValue( - cellIndex, - iNode, - quadRealPointsVec[qPoint]) - - classicalComponentInQuadValues - [numEnrichmentIdsInCell * qPoint + iNode]; + for (unsigned int qPoint = 0; qPoint < nQuadPointInCell; + qPoint++) + { + // std::cout << efeBDH->getEnrichmentValue( + // cellIndex, + // iNode, + // quadRealPointsVec[qPoint]) << " " << + // classicalComponentInQuadValues + // [numEnrichmentIdsInCell * qPoint + iNode] << + // "\n"; + *(basisEnrichQuadStorageMapTmp[cellIndex].data() + + qPoint * numEnrichmentIdsInCell + iNode) = + efeBDH->getEnrichmentValue( + cellIndex, iNode, quadRealPointsVec[qPoint]) - + classicalComponentInQuadValues + [numEnrichmentIdsInCell * qPoint + iNode]; + } } } - } } if (basisStorageAttributesBoolMap @@ -718,7 +729,8 @@ namespace dftefe cellIndex * nDimSqxNumQuad; if (locallyOwnedCellIter == efeBDH->beginLocallyOwnedCells()) { - for (unsigned int iNode = 0; iNode < classicalDofsPerCell; iNode++) + for (unsigned int iNode = 0; iNode < classicalDofsPerCell; + iNode++) { for (unsigned int qPoint = 0; qPoint < nQuadPointInCell; qPoint++) @@ -728,7 +740,8 @@ namespace dftefe for (unsigned int iDim = 0; iDim < dim; iDim++) { auto it = - basisGradientParaCellClassQuadStorageTmp.begin() + + basisGradientParaCellClassQuadStorageTmp + .begin() + qPoint * dim * classicalDofsPerCell + iDim * classicalDofsPerCell + iNode; *it = shapeGrad[iDim]; @@ -753,33 +766,35 @@ namespace dftefe } } - if(numEnrichmentIdsInCell > 0) - { - basisGradientEnrichQuadStorageMapTmp[cellIndex] - .resize(dim * nQuadPointInCell * numEnrichmentIdsInCell); - for (unsigned int iNode = 0; iNode < numEnrichmentIdsInCell; iNode++) + if (numEnrichmentIdsInCell > 0) { - for (unsigned int qPoint = 0; qPoint < nQuadPointInCell; - qPoint++) + basisGradientEnrichQuadStorageMapTmp[cellIndex].resize( + dim * nQuadPointInCell * numEnrichmentIdsInCell); + for (unsigned int iNode = 0; iNode < numEnrichmentIdsInCell; + iNode++) { - auto shapeGrad = efeBDH->getEnrichmentDerivative( - cellIndex, - iNode, - quadRealPointsVec[qPoint]); - // enriched gradient function call - for (unsigned int iDim = 0; iDim < dim; iDim++) + for (unsigned int qPoint = 0; qPoint < nQuadPointInCell; + qPoint++) { - auto it = basisGradientEnrichQuadStorageMapTmp[cellIndex].data() + - qPoint * dim * numEnrichmentIdsInCell + - iDim * numEnrichmentIdsInCell + iNode; - *it = shapeGrad[iDim] - - classicalComponentInQuadGradients - [numEnrichmentIdsInCell * dim * qPoint + - iDim * numEnrichmentIdsInCell + iNode]; + auto shapeGrad = efeBDH->getEnrichmentDerivative( + cellIndex, iNode, quadRealPointsVec[qPoint]); + // enriched gradient function call + for (unsigned int iDim = 0; iDim < dim; iDim++) + { + auto it = + basisGradientEnrichQuadStorageMapTmp + [cellIndex] + .data() + + qPoint * dim * numEnrichmentIdsInCell + + iDim * numEnrichmentIdsInCell + iNode; + *it = shapeGrad[iDim] - + classicalComponentInQuadGradients + [numEnrichmentIdsInCell * dim * qPoint + + iDim * numEnrichmentIdsInCell + iNode]; + } } } } - } } if (basisStorageAttributesBoolMap @@ -860,15 +875,16 @@ namespace dftefe basisParaCellClassQuadStorage->data(), basisParaCellClassQuadStorageTmp.data()); - for(auto &it : basisEnrichQuadStorageMapTmp) + for (auto &it : basisEnrichQuadStorageMapTmp) { - basisEnrichQuadStorageMap[it.first] = - std::make_shared::Storage>(it.second.size()); - utils::MemoryTransfer::copy( - it.second.size(), - basisEnrichQuadStorageMap[it.first]->data(), - it.second.data()); + basisEnrichQuadStorageMap[it.first] = std::make_shared< + typename BasisDataStorage::Storage>( + it.second.size()); + utils::MemoryTransfer:: + copy(it.second.size(), + basisEnrichQuadStorageMap[it.first]->data(), + it.second.data()); } } @@ -886,15 +902,16 @@ namespace dftefe basisJacobianInvQuadStorage->data(), basisJacobianInvQuadStorageTmp.data()); - for(auto &it : basisGradientEnrichQuadStorageMapTmp) + for (auto &it : basisGradientEnrichQuadStorageMapTmp) { - basisGradientEnrichQuadStorageMap[it.first] = - std::make_shared::Storage>(it.second.size()); - utils::MemoryTransfer::copy( - it.second.size(), - basisGradientEnrichQuadStorageMap[it.first]->data(), - it.second.data()); + basisGradientEnrichQuadStorageMap[it.first] = std::make_shared< + typename BasisDataStorage::Storage>( + it.second.size()); + utils::MemoryTransfer:: + copy(it.second.size(), + basisGradientEnrichQuadStorageMap[it.first]->data(), + it.second.data()); } } if (basisStorageAttributesBoolMap @@ -949,7 +966,8 @@ namespace dftefe d_dofsInCell[iCell] = d_efeBDH->nCellDofs(iCell); } d_tmpGradientBlock = nullptr; - d_classialDofsInCell = utils::mathFunctions::sizeTypePow((d_efeBDH->getFEOrder(0) + 1), dim); + d_classialDofsInCell = + utils::mathFunctions::sizeTypePow((d_efeBDH->getFEOrder(0) + 1), dim); } template (d_efeBDH, basisParaCellClassQuadStorage, - d_basisEnrichQuadStorageMap, + d_basisEnrichQuadStorageMap, basisGradientParaCellClassQuadStorage, d_basisGradientEnrichQuadStorageMap, basisJacobianInvQuadStorage, @@ -1142,15 +1160,16 @@ namespace dftefe .find(BasisStorageAttributes::StoreValues) ->second) { - d_basisParaCellClassQuadStorage = basisParaCellClassQuadStorage; + d_basisParaCellClassQuadStorage = basisParaCellClassQuadStorage; } if (basisStorageAttributesBoolMap .find(BasisStorageAttributes::StoreGradient) ->second) { - d_basisGradientParaCellClassQuadStorage = basisGradientParaCellClassQuadStorage; - d_basisJacobianInvQuadStorage = basisJacobianInvQuadStorage; + d_basisGradientParaCellClassQuadStorage = + basisGradientParaCellClassQuadStorage; + d_basisJacobianInvQuadStorage = basisJacobianInvQuadStorage; d_cellStartIdsBasisJacobianInvQuadStorage = cellStartIdsBasisJacobianInvQuadStorage; d_tmpGradientBlock = std::make_shared( @@ -1365,7 +1384,7 @@ namespace dftefe memorySpace, dim>(d_efeBDH, basisParaCellClassQuadStorage, - d_basisEnrichQuadStorageMap, + d_basisEnrichQuadStorageMap, basisGradientParaCellClassQuadStorage, d_basisGradientEnrichQuadStorageMap, basisJacobianInvQuadStorage, @@ -1382,15 +1401,16 @@ namespace dftefe .find(BasisStorageAttributes::StoreValues) ->second) { - d_basisParaCellClassQuadStorage = basisParaCellClassQuadStorage; + d_basisParaCellClassQuadStorage = basisParaCellClassQuadStorage; } if (basisStorageAttributesBoolMap .find(BasisStorageAttributes::StoreGradient) ->second) { - d_basisGradientParaCellClassQuadStorage = basisGradientParaCellClassQuadStorage; - d_basisJacobianInvQuadStorage = basisJacobianInvQuadStorage; + d_basisGradientParaCellClassQuadStorage = + basisGradientParaCellClassQuadStorage; + d_basisJacobianInvQuadStorage = basisJacobianInvQuadStorage; d_cellStartIdsBasisJacobianInvQuadStorage = cellStartIdsBasisJacobianInvQuadStorage; d_tmpGradientBlock = std::make_shared( @@ -1603,7 +1623,6 @@ namespace dftefe .find(BasisStorageAttributes::StoreValues) ->second, "Basis values are not evaluated for the given QuadratureRuleAttributes"); - } template (d_basisParaCellClassQuadStorage->data(), - d_classialDofsInCell, - d_classialDofsInCell * quadId, - cumulativeOffset + - d_dofsInCell[cellId] * quadId); - if (iter != d_basisEnrichQuadStorageMap.end()) - { - basisData.template copyFrom(iter->second->data(), - d_dofsInCell[cellId] - d_classialDofsInCell, - (d_dofsInCell[cellId] - d_classialDofsInCell) * quadId, - cumulativeOffset + - d_dofsInCell[cellId] * quadId - + d_classialDofsInCell); - } + auto iter = d_basisEnrichQuadStorageMap.find(cellId); + for (size_type quadId = 0; quadId < d_nQuadPointsIncell[cellId]; + quadId++) + { + basisData.template copyFrom( + d_basisParaCellClassQuadStorage->data(), + d_classialDofsInCell, + d_classialDofsInCell * quadId, + cumulativeOffset + d_dofsInCell[cellId] * quadId); + if (iter != d_basisEnrichQuadStorageMap.end()) + { + basisData.template copyFrom( + iter->second->data(), + d_dofsInCell[cellId] - d_classialDofsInCell, + (d_dofsInCell[cellId] - d_classialDofsInCell) * quadId, + cumulativeOffset + d_dofsInCell[cellId] * quadId + + d_classialDofsInCell); + } + } + cumulativeOffset += + d_dofsInCell[cellId] * d_nQuadPointsIncell[cellId]; } - cumulativeOffset += d_dofsInCell[cellId] * d_nQuadPointsIncell[cellId]; - } } template (iter->second->data(), - (d_dofsInCell[cellId] - d_classialDofsInCell) * dim, - (d_dofsInCell[cellId] - d_classialDofsInCell) * dim * quadId, - cumulativeOffset + - d_dofsInCell[cellId] * quadId * dim - + d_classialDofsInCell * dim); + auto iter = d_basisGradientEnrichQuadStorageMap.find(cellId); + if (iter != d_basisGradientEnrichQuadStorageMap.end()) + { + for (size_type quadId = 0; quadId < d_nQuadPointsIncell[cellId]; + quadId++) + for (size_type iDim = 0; iDim < dim; iDim++) + basisGradientData.template copyFrom( + iter->second->data(), + (d_dofsInCell[cellId] - d_classialDofsInCell), + (d_dofsInCell[cellId] - d_classialDofsInCell) * dim * + quadId + + iDim, + cumulativeOffset + d_dofsInCell[cellId] * dim * quadId + + d_dofsInCell[cellId] * iDim + d_classialDofsInCell); + } + cumulativeOffset += + d_dofsInCell[cellId] * d_nQuadPointsIncell[cellId] * dim; } - cumulativeOffset += d_dofsInCell[cellId] * d_nQuadPointsIncell[cellId] * dim; - } } template second, "Basis values are not evaluated for the given QuadraturePointAttributes"); - } template #include #include -#include +#include +using namespace std::chrono; +#include using namespace dftefe; int main() { @@ -76,7 +78,7 @@ int main() const unsigned int dim = 3; std::shared_ptr triangulationBase = std::make_shared>(comm); - std::vector subdivisions = {5, 5, 5}; + std::vector subdivisions = {10, 10, 10}; std::vector isPeriodicFlags(dim, false); std::vector domainVectors(dim, utils::Point(dim, 0.0)); @@ -185,7 +187,7 @@ int main() } // Make orthogonalized EFE basis - unsigned int feOrder = 3; + unsigned int feOrder = 5; // Set up the vector of scalarSpatialRealFunctions for adaptive quadrature std::vector> functionsVec(0); @@ -214,7 +216,7 @@ int main() // Set up base quadrature rule for adaptive quadrature - unsigned int num1DGaussSize = feOrder + 1; + unsigned int num1DGaussSize = feOrder + 5; std::shared_ptr baseQuadRule = std::make_shared(dim, num1DGaussSize); @@ -247,7 +249,7 @@ int main() basisAttrMap[basis::BasisStorageAttributes::StoreHessian] = false; basisAttrMap[basis::BasisStorageAttributes::StoreOverlap] = false; basisAttrMap[basis::BasisStorageAttributes::StoreGradNiGradNj] = false; - basisAttrMap[basis::BasisStorageAttributes::StoreJxW] = false; + basisAttrMap[basis::BasisStorageAttributes::StoreJxW] = true; // Set up the CFE Basis Data Storage for Rhs std::shared_ptr> cfeBasisDataStorageAdaptUniformQuad = @@ -256,17 +258,19 @@ int main() // evaluate basis data cfeBasisDataStorageAdaptUniformQuad->evaluateBasisData(quadAttrAdaptive, quadRuleContainerAdaptive, basisAttrMap); + size_type numLocallyOwnedCells = cfeBasisDofHandler->nLocallyOwnedCells(); 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); - // 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); + + // // check basis data for(int i = 0 ; i < numLocallyOwnedCells ; i += cellBlockSize) { size_type cellStartId = i; @@ -277,7 +281,7 @@ int main() cfeBasisDataStorageUniformQuad->getBasisDataInCellRange(std::make_pair(cellStartId, cellEndId), basisDataStorageInCellRange); for(int j = 0 ; j < cellEndId - cellStartId ; j++) { - auto basisDataStorageInCellAdaptUniform = cfeBasisDataStorageAdaptUniformQuad->getBasisDataInCell(j); + auto basisDataStorageInCellAdaptUniform = cfeBasisDataStorageAdaptUniformQuad->getBasisDataInCell(j+cellStartId); for(int k = 0 ; k < basisDataStorageInCellAdaptUniform.size() ; k++) { double val1 = *(basisDataStorageInCellRange.data() + k); @@ -288,6 +292,10 @@ int main() } } // check basis gradient data + + utils::mpi::MPIBarrier(comm); + auto startTotal = std::chrono::high_resolution_clock::now(); + for(int i = 0 ; i < numLocallyOwnedCells ; i += cellBlockSize) { size_type cellStartId = i; @@ -298,7 +306,7 @@ int main() cfeBasisDataStorageUniformQuad->getBasisGradientDataInCellRange(std::make_pair(cellStartId, cellEndId), basisGradientDataStorageInCellRange); for(int j = 0 ; j < cellEndId - cellStartId ; j++) { - auto basisGradientDataStorageInCellAdaptUniform = cfeBasisDataStorageAdaptUniformQuad->getBasisGradientDataInCell(j); + auto basisGradientDataStorageInCellAdaptUniform = cfeBasisDataStorageAdaptUniformQuad->getBasisGradientDataInCell(j+cellStartId); for(int k = 0 ; k < basisGradientDataStorageInCellAdaptUniform.size() ; k++) { double val1 = *(basisGradientDataStorageInCellRange.data() + k); @@ -308,6 +316,12 @@ int main() } } } + utils::mpi::MPIBarrier(comm); + auto stopTotal = std::chrono::high_resolution_clock::now(); + + auto durationTotal = std::chrono::duration_cast(stopTotal - startTotal); + + std::cout << "Total wall time(in secs) : " << durationTotal.count()/1e6 << std::endl; // Create the enrichmentClassicalInterface object std::shared_ptr> @@ -356,50 +370,75 @@ int main() // 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); + nDofs = 0; + for(int iCell = cellStartId ; iCell < cellEndId ; iCell++) + { + nDofs += basisDofHandler->nCellDofs(iCell); + } utils::MemoryStorage - basisDataStorageInCellRange((cellEndId - cellStartId)*nDofs*nQuadPtInCell); + basisDataStorageInCellRange(nDofs*nQuadPtInCell); efeBasisDataStorageUniformQuad->getBasisDataInCellRange(std::make_pair(cellStartId, cellEndId), basisDataStorageInCellRange); - for(int j = 0 ; j < cellEndId - cellStartId ; j++) + size_type cumulativeOffset = 0; + for(int j = cellStartId ; j < cellEndId ; j++) { auto basisDataStorageInCellAdaptUniform = efeBasisDataStorageAdaptUniformQuad->getBasisDataInCell(j); for(int k = 0 ; k < basisDataStorageInCellAdaptUniform.size() ; k++) { - double val1 = *(basisDataStorageInCellRange.data() + k); + double val1 = *(basisDataStorageInCellRange.data() + cumulativeOffset + k); double val2 = *(basisDataStorageInCellAdaptUniform.data() + k); if(std::abs(val1 - val2) > 1e-12) - std::cout << val1 << "\t" << val2 << "\n"; + std::cout << val1 << "\t" << val2 << "--- Basis " << rank << "\n"; } + cumulativeOffset += basisDofHandler->nCellDofs(j) * nQuadPtInCell; } } // check basis gradient data + + utils::mpi::MPIBarrier(comm); + startTotal = std::chrono::high_resolution_clock::now(); + for(int i = 0 ; i < numLocallyOwnedCells ; i += cellBlockSize) { size_type cellStartId = i; size_type cellEndId = std::min(i+cellBlockSize, numLocallyOwnedCells); + nDofs = 0; + for(int iCell = cellStartId ; iCell < cellEndId ; iCell++) + { + nDofs += basisDofHandler->nCellDofs(iCell); + } utils::MemoryStorage - basisGradientDataStorageInCellRange((cellEndId - cellStartId)*nDofs*nQuadPtInCell*dim); - efeBasisDataStorageUniformQuad->getBasisGradientDataInCellRange(std::make_pair(cellStartId, cellEndId), basisGradientDataStorageInCellRange); - for(int j = 0 ; j < cellEndId - cellStartId ; j++) + basisDataStorageInCellRange(nDofs*nQuadPtInCell*dim); + efeBasisDataStorageUniformQuad->getBasisGradientDataInCellRange(std::make_pair(cellStartId, cellEndId), basisDataStorageInCellRange); + size_type cumulativeOffset = 0; + for(int j = cellStartId ; j < cellEndId ; j++) { - auto basisGradientDataStorageInCellAdaptUniform = efeBasisDataStorageAdaptUniformQuad->getBasisGradientDataInCell(j); - for(int k = 0 ; k < basisGradientDataStorageInCellAdaptUniform.size() ; k++) + auto basisDataStorageInCellAdaptUniform = efeBasisDataStorageAdaptUniformQuad->getBasisGradientDataInCell(j); + for(int k = 0 ; k < basisDataStorageInCellAdaptUniform.size() ; k++) { - double val1 = *(basisGradientDataStorageInCellRange.data() + k); - double val2 = *(basisGradientDataStorageInCellAdaptUniform.data() + k); + double val1 = *(basisDataStorageInCellRange.data() + cumulativeOffset + k); + double val2 = *(basisDataStorageInCellAdaptUniform.data() + k); if(std::abs(val1 - val2) > 1e-12) - std::cout << val1 << "\t" << val2 << "\n"; + std::cout << val1 << "\t" << val2 << "--- Grad " << rank << "\n"; } + cumulativeOffset += basisDofHandler->nCellDofs(j) * nQuadPtInCell * dim; } } + utils::mpi::MPIBarrier(comm); + stopTotal = std::chrono::high_resolution_clock::now(); + + durationTotal = std::chrono::duration_cast(stopTotal - startTotal); + + std::cout << "Total wall time(in secs) : " << durationTotal.count()/1e6 << std::endl; + utils::mpi::MPIBarrier(comm); //gracefully end MPI