diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index cff526fc27..c8dfe22706 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -39,6 +39,13 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ - The maximum number of vertices allowed in polygon primitives can now be passed as a template argument to `axom::bump::TopologyMapper`, `axom::bump::PrimalAdaptor`, and `axom::mir::ElviraAlgorithm`. +- Material views in `axom::bump::views` were enhanced with `const_iterator` classes that + enable traversal of material data for zones so kernels do not need to use large fixed size + buffers to gather that data inside kernels. +- Material views in `axom::bump::views` were enhanced with an overloaded `zoneMaterials()` + method that allows data to be gathered into `axom::ArrayView` objects. +- A new `heavily_mixed` example program was added in `axom::mir` to demonstrate running MIR on + meshes with heavily mixed zones. - `saveDocument()` now has a `AUTODETECT` protocol for the file type - Sina fortran can now handle multiple records rather than a single record per application - Most Sina Fortran call can now pass the record for which the call is desired (`sina_add`, `sina_add_file`, `sina_add_curveset`, `sina_add_curve`) diff --git a/src/axom/bump/MatsetSlicer.hpp b/src/axom/bump/MatsetSlicer.hpp index e44145b544..6916333237 100644 --- a/src/axom/bump/MatsetSlicer.hpp +++ b/src/axom/bump/MatsetSlicer.hpp @@ -112,15 +112,12 @@ class MatsetSlicer const auto size = static_cast(sizesView[index]); const auto offset = offsetsView[index]; - typename MatsetView::IDList ids; - typename MatsetView::VFList vfs; - deviceMatsetView.zoneMaterials(deviceSelectedZonesView[index], ids, vfs); - - for(int i = 0; i < size; i++) + auto zoneMat = deviceMatsetView.beginZone(deviceSelectedZonesView[index]); + for(int i = 0; i < size; i++, zoneMat++) { const auto destIndex = offset + i; - materialIdsView[destIndex] = ids[i]; - volumeFractionsView[destIndex] = vfs[i]; + materialIdsView[destIndex] = zoneMat.material_id(); + volumeFractionsView[destIndex] = zoneMat.volume_fraction(); indicesView[destIndex] = destIndex; } }); diff --git a/src/axom/bump/MergeMeshes.hpp b/src/axom/bump/MergeMeshes.hpp index eb0ca81e3c..9d5ce9f08d 100644 --- a/src/axom/bump/MergeMeshes.hpp +++ b/src/axom/bump/MergeMeshes.hpp @@ -1811,8 +1811,6 @@ class MergeMeshesAndMatsets : public MergeMeshes axom::IndexType nzones, axom::IndexType zOffset) const { - using IDList = typename decltype(matsetView)::IDList; - using VFList = typename decltype(matsetView)::VFList; using MatID = typename decltype(matsetView)::IndexType; // Make some maps for renumbering material numbers. @@ -1844,19 +1842,18 @@ class MergeMeshesAndMatsets : public MergeMeshes nzones, AXOM_LAMBDA(axom::IndexType zoneIndex) { // Get this zone's materials. - IDList ids; - VFList vfs; - matsetView.zoneMaterials(zoneIndex, ids, vfs); + auto zoneMat = matsetView.beginZone(zoneIndex); + const auto nmats = zoneMat.size(); // Store the materials in the new material. const auto zoneStart = offsetsView[zOffset + zoneIndex]; - for(axom::IndexType mi = 0; mi < ids.size(); mi++) + for(axom::IndexType mi = 0; mi < nmats; mi++, zoneMat++) { const auto destIndex = zoneStart + mi; - volumeFractionsView[destIndex] = vfs[mi]; + volumeFractionsView[destIndex] = zoneMat.volume_fraction(); // Get the index of the material number in the local map. - const auto mapIndex = axom::utilities::binary_search(localView, ids[mi]); + const auto mapIndex = axom::utilities::binary_search(localView, zoneMat.material_id()); SLIC_ASSERT(mapIndex != -1); // We'll store the all materials number. const auto allMatno = allView[mapIndex]; diff --git a/src/axom/bump/TopologyMapper.hpp b/src/axom/bump/TopologyMapper.hpp index aa2f5a966f..b1342766be 100644 --- a/src/axom/bump/TopologyMapper.hpp +++ b/src/axom/bump/TopologyMapper.hpp @@ -656,11 +656,9 @@ class TopologyMapper // Get the src material - there should just be one because we assume // that a clean matset is being mapped. - typename SrcMatsetView::IDList zoneMatIds; - typename SrcMatsetView::VFList zoneMatVFs; - srcMatsetView.zoneMaterials(srcZone, zoneMatIds, zoneMatVFs); - SLIC_ASSERT(zoneMatIds.size() == 1); - const auto mat = zoneMatIds[0]; + auto zoneMat = srcMatsetView.beginZone(srcZone); + SLIC_ASSERT(zoneMat.size() == 1); + const auto mat = zoneMat.material_id(); #if defined(AXOM_DEBUG_TOPOLOGY_MAPPER) && !defined(AXOM_DEVICE_CODE) std::cout << "\tintersection:" << std::endl diff --git a/src/axom/bump/tests/bump_views.cpp b/src/axom/bump/tests/bump_views.cpp index 6d84b8bc75..07471d4587 100644 --- a/src/axom/bump/tests/bump_views.cpp +++ b/src/axom/bump/tests/bump_views.cpp @@ -658,6 +658,78 @@ struct test_braid2d_mat { EXPECT_EQ(results[i], resultsHost[i]); } + + // Test iterators. + test_matsetview_iterators(matsetView, allocatorID); + } + + template + static void test_matsetview_iterators(MatsetView matsetView, int allocatorID) + { + using ZoneIndex = typename MatsetView::ZoneIndex; + // Allocate results array on device. + const int nResults = matsetView.numberOfZones(); + axom::Array resultsArrayDevice(nResults, nResults, allocatorID); + auto resultsView = resultsArrayDevice.view(); + + axom::for_all( + matsetView.numberOfZones(), + AXOM_LAMBDA(axom::IndexType index) { + typename MatsetView::IDList ids {}; + typename MatsetView::VFList vfs {}; + matsetView.zoneMaterials(index, ids, vfs); + + // Get the end iterator for the zone. + const auto end = matsetView.endZone(index); + + int eq_count = 0; + int count = 0; + + // Make sure the iterator is for the right zone. + eq_count += (end.zoneIndex() == static_cast(index)) ? 1 : 0; + count++; + + // Make sure incrementing the last iterator has no effect. + auto end2 = end; + end2++; + eq_count += (end == end2) ? 1 : 0; + count++; + + // Make sure the iterator order is the same as for the values we got from zoneMaterials(). + int i = 0; + for(auto it = matsetView.beginZone(index); it != end; it++, i++) + { + eq_count += (vfs[i] == it.volume_fraction() && ids[i] == it.material_id()) ? 1 : 0; + count++; + } + + // Test ArrayView version of zoneMaterials(). + using IndexType = typename MatsetView::IndexType; + using FloatType = typename MatsetView::FloatType; + constexpr int ARRAY_SIZE = 10; + IndexType idStorage[ARRAY_SIZE]; + FloatType vfStorage[ARRAY_SIZE]; + axom::ArrayView idView(idStorage, ARRAY_SIZE); + axom::ArrayView vfView(vfStorage, ARRAY_SIZE); + const auto nmats = matsetView.zoneMaterials(index, idView, vfView); + eq_count += (nmats == ids.size()) ? 1 : 0; + count++; + for(axom::IndexType j = 0; j < nmats; j++) + { + eq_count += (vfs[j] == vfView[j] && ids[j] == idView[j]) ? 1 : 0; + count++; + } + + resultsView[index] = (eq_count == count) ? 1 : 0; + }); + + // Get containsView data to the host and compare results + std::vector resultsHost(nResults); + axom::copy(resultsHost.data(), resultsView.data(), sizeof(int) * nResults); + for(int i = 0; i < nResults; i++) + { + EXPECT_EQ(resultsHost[i], 1); + } } }; diff --git a/src/axom/bump/views/MaterialView.hpp b/src/axom/bump/views/MaterialView.hpp index d97549ff6f..5a58658f11 100644 --- a/src/axom/bump/views/MaterialView.hpp +++ b/src/axom/bump/views/MaterialView.hpp @@ -75,7 +75,7 @@ MaterialInformation materials(const conduit::Node &matset); \endverbatim */ -template +template class UnibufferMaterialView { public: @@ -136,6 +136,26 @@ class UnibufferMaterialView } } + AXOM_HOST_DEVICE + axom::IndexType zoneMaterials(ZoneIndex zi, + axom::ArrayView &ids, + axom::ArrayView &vfs) const + { + SLIC_ASSERT(zi < static_cast(numberOfZones())); + + const auto sz = numberOfMaterials(zi); + SLIC_ASSERT(sz <= ids.size()); + const auto offset = m_offsets[zi]; + for(axom::IndexType i = 0; i < sz; i++) + { + const auto idx = m_indices[offset + i]; + + ids[i] = m_material_ids[idx]; + vfs[i] = m_volume_fractions[idx]; + } + return sz; + } + AXOM_HOST_DEVICE bool zoneContainsMaterial(ZoneIndex zi, MaterialID mat) const { @@ -163,6 +183,106 @@ class UnibufferMaterialView return false; } + /*! + * \brief An iterator class for iterating over read-only data in a zone. + * The iterator can access material ids and volume fractions for one + * material at a time in the associated zone. + */ + class const_iterator + { + // Let the material view call the const_iterator constructor. + friend class UnibufferMaterialView; + + public: + /// Get the current material id for the iterator. + MaterialID AXOM_HOST_DEVICE material_id() const + { + SLIC_ASSERT(m_currentIndex < size()); + return m_view->m_material_ids[m_index]; + } + /// Get the current volume fraction for the iterator. + FloatType AXOM_HOST_DEVICE volume_fraction() const + { + SLIC_ASSERT(m_currentIndex < size()); + return m_view->m_volume_fractions[m_index]; + } + axom::IndexType AXOM_HOST_DEVICE size() const { return m_view->m_sizes[m_zoneIndex]; } + ZoneIndex AXOM_HOST_DEVICE zoneIndex() const { return m_zoneIndex; } + + void AXOM_HOST_DEVICE operator++() { advance(true); } + void AXOM_HOST_DEVICE operator++(int) { advance(true); } + bool AXOM_HOST_DEVICE operator==(const const_iterator &rhs) const + { + return m_currentIndex == rhs.m_currentIndex && m_zoneIndex == rhs.m_zoneIndex && + m_view == rhs.m_view; + } + bool AXOM_HOST_DEVICE operator!=(const const_iterator &rhs) const + { + return m_currentIndex != rhs.m_currentIndex || m_zoneIndex != rhs.m_zoneIndex || + m_view != rhs.m_view; + } + + private: + DISABLE_DEFAULT_CTOR(const_iterator); + + /// Constructor + AXOM_HOST_DEVICE const_iterator(const UnibufferMaterialView *view, + ZoneIndex zoneIndex, + axom::IndexType currentIndex = 0) + : m_view(view) + , m_zoneIndex(zoneIndex) + , m_currentIndex(currentIndex) + , m_index(0) + { } + + void AXOM_HOST_DEVICE advance(bool doIncrement) + { + m_currentIndex += (doIncrement && m_currentIndex < size()) ? 1 : 0; + const auto idx = m_view->m_offsets[m_zoneIndex] + m_currentIndex; + if(idx < m_view->m_indices.size()) + { + m_index = m_view->m_indices[idx]; + } + } + + const UnibufferMaterialView *m_view; + ZoneIndex m_zoneIndex; + axom::IndexType m_currentIndex; + axom::IndexType m_index; // not considered in ==, != + }; + // Let the const_iterator access members. + friend class const_iterator; + + /*! + * \brief Return the iterator for the beginning of a zone's material data. + * + * \param zi The zone index being queried. + * + * \return The iterator for the beginning of a zone's material data. + */ + const_iterator AXOM_HOST_DEVICE beginZone(ZoneIndex zi) const + { + SLIC_ASSERT(zi < static_cast(numberOfZones())); + + auto it = const_iterator(this, zi, 0); + it.advance(false); + return it; + } + + /*! + * \brief Return the iterator for the end of a zone's material data. + * + * \param zi The zone index being queried. + * + * \return The iterator for the end of a zone's material data. + */ + const_iterator AXOM_HOST_DEVICE endZone(ZoneIndex zi) const + { + SLIC_ASSERT(zi < static_cast(numberOfZones())); + + return const_iterator(this, zi, m_sizes[zi]); + } + private: axom::ArrayView m_material_ids; axom::ArrayView m_volume_fractions; @@ -196,7 +316,7 @@ class UnibufferMaterialView \endverbatim */ -template +template class MultiBufferMaterialView { public: @@ -277,6 +397,31 @@ class MultiBufferMaterialView } } + AXOM_HOST_DEVICE + axom::IndexType zoneMaterials(ZoneIndex zi, + axom::ArrayView &ids, + axom::ArrayView &vfs) const + { + axom::IndexType n = 0; + for(axom::IndexType i = 0; i < m_size; i++) + { + const auto &curIndices = m_indices[i]; + const auto &curValues = m_values[i]; + + if(zi < static_cast(curIndices.size())) + { + const auto idx = curIndices[zi]; + if(curValues[idx] > 0) + { + ids[n] = m_matnos[i]; + vfs[n] = curValues[idx]; + n++; + } + } + } + return n; + } + AXOM_HOST_DEVICE bool zoneContainsMaterial(ZoneIndex zi, MaterialID mat) const { @@ -304,6 +449,123 @@ class MultiBufferMaterialView return found; } + /*! + * \brief An iterator class for iterating over read-only data in a zone. + * The iterator can access material ids and volume fractions for one + * material at a time in the associated zone. + */ + class const_iterator + { + // Let the material view call the const_iterator constructor. + friend class MultiBufferMaterialView; + + public: + /// Get the current material id for the iterator. + MaterialID AXOM_HOST_DEVICE material_id() const + { + SLIC_ASSERT(m_currentIndex < m_view->m_size); + return m_view->m_matnos[m_currentIndex]; + } + /// Get the current volume fraction for the iterator. + FloatType AXOM_HOST_DEVICE volume_fraction() const + { + SLIC_ASSERT(m_currentIndex < m_view->m_size); + const auto &curIndices = m_view->m_indices[m_currentIndex]; + const auto &curValues = m_view->m_values[m_currentIndex]; + const auto idx = curIndices[m_zoneIndex]; + return curValues[idx]; + } + axom::IndexType AXOM_HOST_DEVICE size() const { return m_view->numberOfMaterials(m_zoneIndex); } + ZoneIndex AXOM_HOST_DEVICE zoneIndex() const { return m_zoneIndex; } + void AXOM_HOST_DEVICE operator++() + { + m_currentIndex += (m_currentIndex < m_view->m_size) ? 1 : 0; + advance(); + } + void AXOM_HOST_DEVICE operator++(int) + { + m_currentIndex += (m_currentIndex < m_view->m_size) ? 1 : 0; + advance(); + } + bool AXOM_HOST_DEVICE operator==(const const_iterator &rhs) const + { + return m_currentIndex == rhs.m_currentIndex && m_zoneIndex == rhs.m_zoneIndex && + m_view == rhs.m_view; + } + bool AXOM_HOST_DEVICE operator!=(const const_iterator &rhs) const + { + return m_currentIndex != rhs.m_currentIndex || m_zoneIndex != rhs.m_zoneIndex || + m_view != rhs.m_view; + } + + private: + DISABLE_DEFAULT_CTOR(const_iterator); + + /// Constructor + AXOM_HOST_DEVICE const_iterator(const MultiBufferMaterialView *view, + ZoneIndex zoneIndex, + axom::IndexType currentIndex = 0) + : m_view(view) + , m_zoneIndex(zoneIndex) + , m_currentIndex(currentIndex) + { } + + /// Advance to the next valid material slot for the zone. + void AXOM_HOST_DEVICE advance() + { + while(m_currentIndex < m_view->m_size) + { + const auto &curIndices = m_view->m_indices[m_currentIndex]; + const auto &curValues = m_view->m_values[m_currentIndex]; + + if(m_zoneIndex < static_cast(curIndices.size())) + { + const auto idx = curIndices[m_zoneIndex]; + if(curValues[idx] > 0) + { + break; + } + } + m_currentIndex++; + } + } + + const MultiBufferMaterialView *m_view; + ZoneIndex m_zoneIndex; + axom::IndexType m_currentIndex; + }; + // Let the const_iterator access members. + friend class const_iterator; + + /*! + * \brief Return the iterator for the beginning of a zone's material data. + * + * \param zi The zone index being queried. + * + * \return The iterator for the beginning of a zone's material data. + */ + const_iterator AXOM_HOST_DEVICE beginZone(ZoneIndex zi) const + { + SLIC_ASSERT(zi < static_cast(numberOfZones())); + + auto it = const_iterator(this, zi, 0); + it.advance(); + return it; + } + + /*! + * \brief Return the iterator for the end of a zone's material data. + * + * \param zi The zone index being queried. + * + * \return The iterator for the end of a zone's material data. + */ + const_iterator AXOM_HOST_DEVICE endZone(ZoneIndex zi) const + { + SLIC_ASSERT(zi < static_cast(numberOfZones())); + return const_iterator(this, zi, m_size); + } + private: AXOM_HOST_DEVICE axom::IndexType indexOfMaterialID(MaterialID mat) const @@ -349,7 +611,7 @@ class MultiBufferMaterialView \endverbatim */ -template +template class ElementDominantMaterialView { public: @@ -422,6 +684,26 @@ class ElementDominantMaterialView } } + AXOM_HOST_DEVICE + axom::IndexType zoneMaterials(ZoneIndex zi, + axom::ArrayView &ids, + axom::ArrayView &vfs) const + { + axom::IndexType n = 0; + for(axom::IndexType i = 0; i < m_volume_fractions.size(); i++) + { + const auto ¤tVF = m_volume_fractions[i]; + SLIC_ASSERT(zi < currentVF.size()); + if(currentVF[zi] > 0) + { + ids[n] = m_matnos[i]; + vfs[n] = currentVF[zi]; + n++; + } + } + return n; + } + AXOM_HOST_DEVICE bool zoneContainsMaterial(ZoneIndex zi, MaterialID mat) const { @@ -445,6 +727,113 @@ class ElementDominantMaterialView return found; } + /*! + * \brief An iterator class for iterating over read-only data in a zone. + * The iterator can access material ids and volume fractions for one + * material at a time in the associated zone. + */ + class const_iterator + { + // Let the material view call the const_iterator constructor. + friend class ElementDominantMaterialView; + + public: + /// Get the current material id for the iterator. + MaterialID AXOM_HOST_DEVICE material_id() const + { + SLIC_ASSERT(m_currentIndex < m_view->m_volume_fractions.size()); + return m_view->m_matnos[m_currentIndex]; + } + /// Get the current volume fraction for the iterator. + FloatType AXOM_HOST_DEVICE volume_fraction() const + { + SLIC_ASSERT(m_currentIndex < m_view->m_volume_fractions.size()); + return m_view->m_volume_fractions[m_currentIndex][m_zoneIndex]; + } + axom::IndexType AXOM_HOST_DEVICE size() const { return m_view->numberOfMaterials(m_zoneIndex); } + ZoneIndex AXOM_HOST_DEVICE zoneIndex() const { return m_zoneIndex; } + void AXOM_HOST_DEVICE operator++() + { + m_currentIndex += (m_currentIndex < m_view->m_volume_fractions.size()) ? 1 : 0; + advance(); + } + void AXOM_HOST_DEVICE operator++(int) + { + m_currentIndex += (m_currentIndex < m_view->m_volume_fractions.size()) ? 1 : 0; + advance(); + } + bool AXOM_HOST_DEVICE operator==(const const_iterator &rhs) const + { + return m_currentIndex == rhs.m_currentIndex && m_zoneIndex == rhs.m_zoneIndex && + m_view == rhs.m_view; + } + bool AXOM_HOST_DEVICE operator!=(const const_iterator &rhs) const + { + return m_currentIndex != rhs.m_currentIndex || m_zoneIndex != rhs.m_zoneIndex || + m_view != rhs.m_view; + } + + private: + DISABLE_DEFAULT_CTOR(const_iterator); + + /// Constructor + AXOM_HOST_DEVICE const_iterator(const ElementDominantMaterialView *view, + ZoneIndex zoneIndex, + axom::IndexType currentIndex = 0) + : m_view(view) + , m_zoneIndex(zoneIndex) + , m_currentIndex(currentIndex) + { } + + /// Advance to the next valid material slot for the zone. + void AXOM_HOST_DEVICE advance() + { + while(m_currentIndex < m_view->m_volume_fractions.size()) + { + if(m_view->m_volume_fractions[m_currentIndex][m_zoneIndex] > 0) + { + break; + } + m_currentIndex++; + } + } + + const ElementDominantMaterialView *m_view; + ZoneIndex m_zoneIndex; + axom::IndexType m_currentIndex; + }; + // Let the const_iterator access members. + friend class const_iterator; + + /*! + * \brief Return the iterator for the beginning of a zone's material data. + * + * \param zi The zone index being queried. + * + * \return The iterator for the beginning of a zone's material data. + */ + const_iterator AXOM_HOST_DEVICE beginZone(ZoneIndex zi) const + { + SLIC_ASSERT(zi < static_cast(numberOfZones())); + + auto it = const_iterator(this, zi, 0); + it.advance(); + return it; + } + + /*! + * \brief Return the iterator for the end of a zone's material data. + * + * \param zi The zone index being queried. + * + * \return The iterator for the end of a zone's material data. + */ + const_iterator AXOM_HOST_DEVICE endZone(ZoneIndex zi) const + { + SLIC_ASSERT(zi < static_cast(numberOfZones())); + return const_iterator(this, zi, m_volume_fractions.size()); + } + private: #if !defined(AXOM_DEVICE_CODE) void checkSizes() const @@ -508,7 +897,7 @@ class ElementDominantMaterialView \note This matset type does not seem so GPU friendly since there is some work to do for some of the queries. */ -template +template class MaterialDominantMaterialView { public: @@ -602,6 +991,31 @@ class MaterialDominantMaterialView } } + AXOM_HOST_DEVICE + axom::IndexType zoneMaterials(ZoneIndex zi, + axom::ArrayView &ids, + axom::ArrayView &vfs) const + { + axom::IndexType n = 0; + for(axom::IndexType mi = 0; mi < m_size; mi++) + { + const auto &element_ids = m_element_ids[mi]; + const auto &volume_fractions = m_volume_fractions[mi]; + const auto sz = element_ids.size(); + for(axom::IndexType i = 0; i < sz; i++) + { + if(element_ids[i] == zi) + { + ids[n] = m_matnos[mi]; + vfs[n] = volume_fractions[i]; + n++; + break; + } + } + } + return n; + } + AXOM_HOST_DEVICE bool zoneContainsMaterial(ZoneIndex zi, MaterialID mat) const { @@ -633,6 +1047,129 @@ class MaterialDominantMaterialView return found; } + /*! + * \brief An iterator class for iterating over read-only data in a zone. + * The iterator can access material ids and volume fractions for one + * material at a time in the associated zone. + */ + class const_iterator + { + // Let the material view call the const_iterator constructor. + friend class MaterialDominantMaterialView; + + public: + /// Get the current material id for the iterator. + MaterialID AXOM_HOST_DEVICE material_id() const + { + SLIC_ASSERT(m_miIndex < m_view->m_size); + return m_view->m_matnos[m_miIndex]; + } + /// Get the current volume fraction for the iterator. + FloatType AXOM_HOST_DEVICE volume_fraction() const + { + SLIC_ASSERT(m_miIndex < m_view->m_size && m_index < m_view->m_element_ids[m_miIndex].size()); + return m_view->m_volume_fractions[m_miIndex][m_index]; + } + ZoneIndex AXOM_HOST_DEVICE zoneIndex() const { return m_zoneIndex; } + axom::IndexType AXOM_HOST_DEVICE size() const { return m_view->numberOfMaterials(m_zoneIndex); } + void AXOM_HOST_DEVICE operator++() { advance(true); } + void AXOM_HOST_DEVICE operator++(int) { advance(true); } + bool AXOM_HOST_DEVICE operator==(const const_iterator &rhs) const + { + return m_miIndex == rhs.m_miIndex && m_index == rhs.m_index && + m_zoneIndex == rhs.m_zoneIndex && m_view == rhs.m_view; + } + bool AXOM_HOST_DEVICE operator!=(const const_iterator &rhs) const + { + return m_miIndex != rhs.m_miIndex || m_index != rhs.m_index || + m_zoneIndex != rhs.m_zoneIndex || m_view != rhs.m_view; + } + + private: + DISABLE_DEFAULT_CTOR(const_iterator); + + /// Constructor + AXOM_HOST_DEVICE const_iterator(const MaterialDominantMaterialView *view, + ZoneIndex zoneIndex, + axom::IndexType miIndex, + axom::IndexType index) + : m_view(view) + , m_zoneIndex(zoneIndex) + , m_miIndex(miIndex) + , m_index(index) + { } + + /// Advance to the next valid material slot for the zone. + void AXOM_HOST_DEVICE advance(bool doIncrement) + { + if(doIncrement) + { + if(m_miIndex < m_view->m_size) + { + m_index = 0; + m_miIndex++; + } + if(m_miIndex == m_view->m_size) + { + m_index = m_view->m_element_ids[m_view->m_size - 1].size(); + } + } + + // Look for the next m_miIndex,m_index pair that contains material for the selected zone index. + for(; m_miIndex < m_view->m_size; m_miIndex++) + { + const auto &element_ids = m_view->m_element_ids[m_miIndex]; + const auto sz = element_ids.size(); + for(; m_index < sz; m_index++) + { + if(element_ids[m_index] == m_zoneIndex) + { + return; + } + } + m_index = (m_miIndex + 1 == m_view->m_size) ? m_index : 0; + } + } + + const MaterialDominantMaterialView *m_view; + ZoneIndex m_zoneIndex; + axom::IndexType m_miIndex; + axom::IndexType m_index; + }; + // Let the const_iterator access members. + friend class const_iterator; + + /*! + * \brief Return the iterator for the beginning of a zone's material data. + * + * \param zi The zone index being queried. + * + * \return The iterator for the beginning of a zone's material data. + */ + const_iterator AXOM_HOST_DEVICE beginZone(ZoneIndex zi) const + { + SLIC_ASSERT(zi < static_cast(numberOfZones())); + + auto it = const_iterator(this, zi, 0, 0); + it.advance(false); + return it; + } + + /*! + * \brief Return the iterator for the end of a zone's material data. + * + * \param zi The zone index being queried. + * + * \return The iterator for the end of a zone's material data. + */ + const_iterator AXOM_HOST_DEVICE endZone(ZoneIndex zi) const + { + SLIC_ASSERT(zi < static_cast(numberOfZones())); + const axom::IndexType miIndex = m_size; + const axom::IndexType index = (m_size > 0) ? (m_volume_fractions[m_size - 1].size()) : 0; + return const_iterator(this, zi, miIndex, index); + } + private: AXOM_HOST_DEVICE axom::IndexType indexOfMaterialID(MaterialID mat) const diff --git a/src/axom/mir/ElviraAlgorithm.hpp b/src/axom/mir/ElviraAlgorithm.hpp index a278a12e25..ea08ce7ed1 100644 --- a/src/axom/mir/ElviraAlgorithm.hpp +++ b/src/axom/mir/ElviraAlgorithm.hpp @@ -35,13 +35,13 @@ #include // Uncomment to save inputs and outputs. -//#define AXOM_ELVIRA_DEBUG +// #define AXOM_ELVIRA_DEBUG // Uncomment to debug make fragments. -//#define AXOM_ELVIRA_DEBUG_MAKE_FRAGMENTS +// #define AXOM_ELVIRA_DEBUG_MAKE_FRAGMENTS // Uncomment to gather ELVIRA data and save to YAML file. -//#define AXOM_ELVIRA_GATHER_INFO +// #define AXOM_ELVIRA_GATHER_INFO #if defined(AXOM_ELVIRA_DEBUG) #include @@ -650,11 +650,15 @@ class ElviraAlgorithm : public axom::mir::MIRAlgorithm axom::Array fragmentVFStencil(numFragmentsStencil, numFragmentsStencil, allocatorID); auto fragmentVFStencilView = fragmentVFStencil.view(); - // Sorted material ids for each zone. + // Sorted material ids / vfs for each zone. axom::Array sortedMaterialIds(numFragments, numFragments, allocatorID); + axom::Array sortedMaterialVfs(numFragments, + numFragments, + allocatorID); auto sortedMaterialIdsView = sortedMaterialIds.view(); + auto sortedMaterialVfsView = sortedMaterialVfs.view(); // Coordinate stencil data for each zone. const auto nzonesStencil = nzones * StencilSize; @@ -677,22 +681,18 @@ class ElviraAlgorithm : public axom::mir::MIRAlgorithm // Where to begin writing this zone's fragment data. const auto offset = matOffsetView[szIndex]; - // Get materials for this zone from the matset. - typename MatsetView::IDList ids; - typename MatsetView::VFList vfs; - deviceMatsetView.zoneMaterials(matZoneIndex, ids, vfs); + // Determine the views for this zone's material data. + axom::ArrayView ids(sortedMaterialIdsView.data() + offset, + matCount); + axom::ArrayView vfs(sortedMaterialVfsView.data() + offset, + matCount); + // Get materials for this zone from the matset, directly into the "sorted" views. + [[maybe_unused]] auto ids_size = deviceMatsetView.zoneMaterials(matZoneIndex, ids, vfs); // Reverse sort the materials by the volume fraction so the larger VFs are first. - SLIC_ASSERT(ids.size() == matCount); + SLIC_ASSERT(ids_size == matCount); axom::utilities::reverse_sort_multiple(vfs.data(), ids.data(), matCount); - // Save sorted ids in sortedMaterialIdsView. - for(axom::IndexType m = 0; m < matCount; m++) - { - const auto fragmentIndex = offset + m; - sortedMaterialIdsView[fragmentIndex] = ids[m]; - } - // Retrieve the stencil data from neighbor zones. auto logical = deviceTopologyView.indexing().IndexToLogicalIndex(zoneIndex); for(int si = 0; si < StencilSize; si++) @@ -920,6 +920,11 @@ class ElviraAlgorithm : public axom::mir::MIRAlgorithm // Get the starting shape. const auto inputShape = deviceShapeView.getShape(zoneIndex); +#if defined(AXOM_ELVIRA_DEBUG_MAKE_FRAGMENTS) && !defined(AXOM_DEVICE_CODE) + // Get the shape's bounding box and enlarge it a little. + auto inputShapeBBox = axom::primal::compute_bounding_box(inputShape); + inputShapeBBox.scale(1.05); +#endif // Get the zone's actual volume. const double zoneVol = utils::ComputeShapeAmount::execute(inputShape); @@ -1002,27 +1007,38 @@ class ElviraAlgorithm : public axom::mir::MIRAlgorithm // Emit clippedShape as material matId buildView.addShape(zoneIndex, fragmentIndex, clippedShape, matId, pt, planeOffset, normalPtr); +#if defined(AXOM_ELVIRA_DEBUG_MAKE_FRAGMENTS) && !defined(AXOM_DEVICE_CODE) + // Examine clippedShape's bounding box. It should NEVER be larger than the + // original inputShape's bounding box. If so, there was probably an error + // in clipping. + const auto clippedShapeBBox = axom::primal::compute_bounding_box(clippedShape); + if(!inputShapeBBox.contains(clippedShapeBBox)) + { + SLIC_ERROR("\tclip: BAD CLIPPED SHAPE IN ZONE " + << zoneIndex << "\n\t\tinputShape=" << inputShape << "\n\t\tinputShapeBBox=" + << inputShapeBBox << "\n\t\tclippedShape=" << clippedShape + << "\n\t\tclippedShapeBBox=" << clippedShapeBBox); + } +#endif + // Clip in the other direction to get the remaining fragment for the next material. if(m == 0) { #if defined(AXOM_ELVIRA_DEBUG_MAKE_FRAGMENTS) && !defined(AXOM_DEVICE_CODE) - SLIC_DEBUG("\tclip: P=" << P << ", before=" << inputShape); -#endif - remaining = axom::primal::clip(inputShape, P); -#if defined(AXOM_ELVIRA_DEBUG_MAKE_FRAGMENTS) && !defined(AXOM_DEVICE_CODE) - SLIC_DEBUG("\tclip: after=" << clippedShape); + SLIC_DEBUG("\tclip: before=" << inputShape << ", P=" << P << ", pt=" << pt); #endif + remaining = axom::primal::clip(inputShape, P, detail::clip_precision::eps); } else { #if defined(AXOM_ELVIRA_DEBUG_MAKE_FRAGMENTS) && !defined(AXOM_DEVICE_CODE) - SLIC_DEBUG("\tclip: P=" << P << ", before=" << remaining); + SLIC_DEBUG("\tclip: before=" << remaining << ", P=" << P << ", pt=" << pt); #endif - remaining = axom::primal::clip(remaining, P); + remaining = axom::primal::clip(remaining, P, detail::clip_precision::eps); + } #if defined(AXOM_ELVIRA_DEBUG_MAKE_FRAGMENTS) && !defined(AXOM_DEVICE_CODE) - SLIC_DEBUG("\tclip: after=" << remaining); + SLIC_DEBUG("\tclip: after=" << remaining); #endif - } } // Emit the last leftover fragment. @@ -1037,6 +1053,20 @@ class ElviraAlgorithm : public axom::mir::MIRAlgorithm lastNormal[d] = -normal[d]; } buildView.addShape(zoneIndex, fragmentIndex, remaining, matId, pt, -planeOffset, lastNormal); + +#if defined(AXOM_ELVIRA_DEBUG_MAKE_FRAGMENTS) && !defined(AXOM_DEVICE_CODE) + // Examine remaining's bounding box. It should NEVER be larger than the + // original inputShape's bounding box. If so, there was probably an error + // in clipping. + const auto remainingBBox = axom::primal::compute_bounding_box(remaining); + if(!inputShapeBBox.contains(remainingBBox)) + { + SLIC_ERROR("\tclip: BAD CLIPPED SHAPE IN ZONE " + << zoneIndex << "\n\t\tinputShape=" << inputShape + << "\n\t\tinputShapeBBox=" << inputShapeBBox << "\n\t\tremaining=" << remaining + << "\n\t\tremainingBBox=" << remainingBBox); + } +#endif }); reportErrors(__LINE__); } diff --git a/src/axom/mir/detail/elvira_detail.hpp b/src/axom/mir/detail/elvira_detail.hpp index 41c5cf0307..2b3cba7fb6 100644 --- a/src/axom/mir/detail/elvira_detail.hpp +++ b/src/axom/mir/detail/elvira_detail.hpp @@ -72,6 +72,21 @@ class ELVIRAOptions : public axom::bump::Options namespace detail { +//------------------------------------------------------------------------------ +template +struct clip_precision +{ + // The default precision used in axom::primal::clip + static constexpr double eps = 1.e-10; +}; + +template <> +struct clip_precision +{ + // Make the axom::primal::clip more precise for double. + static constexpr double eps = 1.e-15; +}; + //------------------------------------------------------------------------------ /*! * \brief Compute a range given by 2 points given a shape and a normal. The @@ -162,7 +177,7 @@ AXOM_HOST_DEVICE inline ClipResultType clipToVolume(const ShapeType &shape, const auto P = axom::primal::Plane(clipNormal, pt, false); // Clip the shape at the current plane. - clippedShape = axom::primal::clip(shape, P); + clippedShape = axom::primal::clip(shape, P, clip_precision::eps); // Find the volume of the clipped shape. const double fragmentVolume = utils::ComputeShapeAmount::execute(clippedShape); diff --git a/src/axom/mir/examples/CMakeLists.txt b/src/axom/mir/examples/CMakeLists.txt index 538202e59d..fc7d11078d 100644 --- a/src/axom/mir/examples/CMakeLists.txt +++ b/src/axom/mir/examples/CMakeLists.txt @@ -4,4 +4,5 @@ # SPDX-License-Identifier: (BSD-3-Clause) add_subdirectory(concentric_circles) +add_subdirectory(heavily_mixed) add_subdirectory(tutorial_simple) diff --git a/src/axom/mir/examples/concentric_circles/MIRApplication.cpp b/src/axom/mir/examples/concentric_circles/MIRApplication.cpp index 92e86f4dcc..94ff2dafd8 100644 --- a/src/axom/mir/examples/concentric_circles/MIRApplication.cpp +++ b/src/axom/mir/examples/concentric_circles/MIRApplication.cpp @@ -42,7 +42,8 @@ int MIRApplication::initialize(int argc, char **argv) app.add_option("--gridsize", gridSize) ->check(axom::CLI::PositiveNumber) ->description("The number of zones along an axis."); - app.add_option("--method", method)->description("The MIR method name (equiz, elvira)"); + app.add_option("--method", method) + ->description("The MIR method (or operation) name (equiz, elvira, traversal)"); app.add_option("--numcircles", numCircles) ->check(axom::CLI::PositiveNumber) ->description("The number of circles to use for material creation."); @@ -108,6 +109,7 @@ int MIRApplication::initialize(int argc, char **argv) int MIRApplication::execute() { axom::slic::SimpleLogger logger(axom::slic::message::Info); + axom::slic::setLoggingMsgLevel(axom::slic::message::Debug); if(handler) { @@ -172,6 +174,7 @@ int MIRApplication::runMIR() } // Begin material interface reconstruction + timer.reset(); timer.start(); conduit::Node options, resultMesh; options["matset"] = "mat"; diff --git a/src/axom/mir/examples/concentric_circles/runMIR.hpp b/src/axom/mir/examples/concentric_circles/runMIR.hpp index bcf033338c..59406a32a9 100644 --- a/src/axom/mir/examples/concentric_circles/runMIR.hpp +++ b/src/axom/mir/examples/concentric_circles/runMIR.hpp @@ -10,6 +10,54 @@ #include "axom/bump.hpp" #include "axom/mir.hpp" // for Mir classes & functions +template +void test_matset_traversal(MatsetView matsetView) +{ + AXOM_ANNOTATE_SCOPE("test_matset_traversal"); + double vf1, vf2; + { + AXOM_ANNOTATE_SCOPE("zoneMaterials"); + axom::ReduceSum vfSum(0.); + axom::for_all( + matsetView.numberOfZones(), + AXOM_LAMBDA(axom::IndexType zoneIndex) { + typename MatsetView::IDList ids; + typename MatsetView::VFList vfs; + matsetView.zoneMaterials(zoneIndex, ids, vfs); + double sum = 0.; + for(axom::IndexType i = 0; i < vfs.size(); i++) + { + sum += vfs[i]; + } + vfSum += sum; + }); + vf1 = vfSum.get(); + } + { + AXOM_ANNOTATE_SCOPE("iterators"); + axom::ReduceSum vfSum(0.); + axom::for_all( + matsetView.numberOfZones(), + AXOM_LAMBDA(axom::IndexType zoneIndex) { + const auto end = matsetView.endZone(zoneIndex); + double sum = 0.; + for(auto it = matsetView.beginZone(zoneIndex); it != end; it++) + { + sum += it.volume_fraction(); + } + vfSum += sum; + }); + vf2 = vfSum.get(); + } + + const double eps = 1.e-10; + SLIC_INFO(axom::fmt::format("test_matset_traversal: vf1={}, vf2={}, nzones={}, result={}", + vf1, + vf2, + matsetView.numberOfZones(), + (axom::utilities::abs(vf1 - vf2) < eps) ? "pass" : "fail")); +} + template int runMIR(const conduit::Node &hostMesh, const conduit::Node &options, conduit::Node &hostResult) { @@ -98,6 +146,12 @@ int runMIR(const conduit::Node &hostMesh, const conduit::Node &options, conduit: MIR m(topologyView, coordsetView, matsetView); m.execute(deviceMesh, options, deviceResult); } + else if(method == "traversal") + { + using namespace axom::bump::views; + auto matsetView = make_unibuffer_matset::view(n_matset); + test_matset_traversal(matsetView); + } else { SLIC_ERROR(axom::fmt::format("Unsupported MIR method {}", method)); diff --git a/src/axom/mir/examples/heavily_mixed/CMakeLists.txt b/src/axom/mir/examples/heavily_mixed/CMakeLists.txt new file mode 100644 index 0000000000..4729454873 --- /dev/null +++ b/src/axom/mir/examples/heavily_mixed/CMakeLists.txt @@ -0,0 +1,59 @@ +# Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +# other Axom Project Developers. See the top-level COPYRIGHT file for details. +# +# SPDX-License-Identifier: (BSD-3-Clause) + +set( mir_example_dependencies + core + slic + mir + bump + ) + +# Speed up compilation by separating out the execspaces into different runMIR_xxx.cpp files. +axom_add_executable( + NAME mir_heavily_mixed + SOURCES mir_heavily_mixed.cpp HMApplication.cpp runMIR_seq.cpp runMIR_omp.cpp runMIR_cuda.cpp runMIR_hip.cpp + OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} + DEPENDS_ON ${mir_example_dependencies} + FOLDER axom/mir/examples + ) + +if(AXOM_ENABLE_TESTS) + set (_policies "seq") + if(RAJA_FOUND AND UMPIRE_FOUND) + blt_list_append(TO _policies ELEMENTS "omp" IF AXOM_ENABLE_OPENMP) + blt_list_append(TO _policies ELEMENTS "cuda" IF AXOM_ENABLE_CUDA) + blt_list_append(TO _policies ELEMENTS "hip" IF AXOM_ENABLE_HIP) + endif() + + foreach(_policy ${_policies}) + set(_testname "mir_heavily_mixed_2D_${_policy}") + axom_add_test( + NAME ${_testname} + COMMAND mir_heavily_mixed + --dims 100 100 + --refinement 20 + --materials 10 + --policy ${_policy} + --disable-write + ) + + set_tests_properties(${_testname} PROPERTIES + PASS_REGULAR_EXPRESSION "Material interface reconstruction time:") + + set(_testname "mir_heavily_mixed_3D_${_policy}") + axom_add_test( + NAME ${_testname} + COMMAND mir_heavily_mixed + --dims 50 50 10 + --refinement 10 + --materials 10 + --policy ${_policy} + --disable-write + ) + + set_tests_properties(${_testname} PROPERTIES + PASS_REGULAR_EXPRESSION "Material interface reconstruction time:") + endforeach() +endif() diff --git a/src/axom/mir/examples/heavily_mixed/HMApplication.cpp b/src/axom/mir/examples/heavily_mixed/HMApplication.cpp new file mode 100644 index 0000000000..b0a14eecfd --- /dev/null +++ b/src/axom/mir/examples/heavily_mixed/HMApplication.cpp @@ -0,0 +1,439 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include "axom/config.hpp" +#include "axom/core.hpp" // for axom macros +#include "axom/slic.hpp" +#include "axom/bump.hpp" +#include "axom/mir.hpp" // for Mir classes & functions +#include "runMIR.hpp" +#include "HMApplication.hpp" + +#include +#include + +#include + +using RuntimePolicy = axom::runtime_policy::Policy; + +namespace detail +{ +/*! + * \brief Turn a field on a fine mesh into a matset on the coarse mesh. + * + * \param topoName The name of the topology. + * \param dims The number of cells in each logical dimension. + * \param refinement The amount of refinement between coarse/fine meshes. + * \param n_coarse The coarse mesh (it gets the matset). + * \param n_field The field used for matset creation. + * \param nmats The number of materials to make. + */ +void heavily_mixed_matset(const std::string &topoName, + int dims[3], + int refinement, + conduit::Node &n_coarse, + const conduit::Node &n_field, + int nmats) +{ + const auto fine = n_field.as_float64_accessor(); + int nzones = dims[0] * dims[1] * dims[2]; + int nslots = nzones * nmats; + std::vector vfs(nslots, 0.); + + // break the data range into nmats parts. + const double matSize = 1000. / nmats; //fine.max() / nmats; + + const int rdims[] = {dims[0] * refinement, dims[1] * refinement, dims[2] * refinement}; + + for(int k = 0; k < dims[2]; k++) + { + for(int j = 0; j < dims[1]; j++) + { + for(int i = 0; i < dims[0]; i++) + { + const int zoneIndex = k * dims[1] * dims[0] + j * dims[0] + i; + const int kr = k * refinement; + const int jr = j * refinement; + const int ir = i * refinement; + for(int jj = 0; jj < refinement; jj++) + for(int ii = 0; ii < refinement; ii++) + { + const int fine_index = (kr * rdims[0] * rdims[1]) + ((jr + jj) * rdims[0]) + (ir + ii); + const int matid = + axom::utilities::clampVal(static_cast(fine[fine_index] / matSize), 0, nmats - 1); + const int matslot = zoneIndex * nmats + matid; + vfs[matslot] += 1. / (refinement * refinement); + } + } + } + } + + std::vector material_ids; + std::vector volume_fractions; + std::vector indices; + std::vector sizes; + std::vector offsets; + for(int k = 0; k < dims[2]; k++) + { + for(int j = 0; j < dims[1]; j++) + { + for(int i = 0; i < dims[0]; i++) + { + const int zoneIndex = k * dims[0] * dims[1] + j * dims[0] + i; + + int size = 0; + offsets.push_back(indices.size()); + for(int m = 0; m < nmats; m++) + { + int matslot = zoneIndex * nmats + m; + if(vfs[matslot] > 0) + { + indices.push_back(material_ids.size()); + material_ids.push_back(m + 1); + volume_fractions.push_back(vfs[matslot]); + size++; + } + } + sizes.push_back(size); + } + } + } + conduit::Node &n_matset = n_coarse["matsets/mat"]; + n_matset["topology"] = topoName; + conduit::Node &n_material_map = n_matset["material_map"]; + for(int i = 0; i < nmats; i++) + { + int matno = i + 1; + const std::string name = axom::fmt::format("mat{:02d}", matno); + n_material_map[name] = matno; + } + n_matset["material_ids"].set(material_ids); + n_matset["volume_fractions"].set(volume_fractions); + n_matset["indices"].set(indices); + n_matset["sizes"].set(sizes); + n_matset["offsets"].set(offsets); + + n_coarse["fields/nmats/association"] = "element"; + n_coarse["fields/nmats/topology"] = topoName; + n_coarse["fields/nmats/values"].set(sizes); +} + +template +void heavily_mixed(conduit::Node &n_mesh, int dims[3], int refinement, int nmats) +{ + const int rdims[] = {refinement * dims[0], refinement * dims[1], refinement * dims[2]}; + + // Default window + const conduit::float64 x_min = -0.6; + const conduit::float64 x_max = 0.6; + const conduit::float64 y_min = -0.5; + const conduit::float64 y_max = 0.5; + const conduit::float64 c_re = -0.5125; + const conduit::float64 c_im = 0.5213; + + conduit::blueprint::mesh::examples::julia(dims[0], dims[1], x_min, x_max, y_min, y_max, c_re, c_im, n_mesh); + if(dims[2] > 1) + { + // Add another dimension to the coordset. + const conduit::float64 z_min = 0.; + const conduit::float64 z_max = x_max - x_min; + std::vector z; + for(int i = 0; i <= dims[2]; i++) + { + const auto t = static_cast(i) / dims[2]; + const auto zc = axom::utilities::lerp(z_min, z_max, t); + z.push_back(zc); + } + n_mesh["coordsets/coords/values/z"].set(z); + + // Destination window + const conduit::float64 s = 0.9; + const conduit::float64 x1_min = x_min * s; + const conduit::float64 x1_max = x_max * s; + const conduit::float64 y1_min = y_min * s; + const conduit::float64 y1_max = y_max * s; + + conduit::Node n_field; + n_field.set(conduit::DataType::int32(rdims[0] * rdims[1] * rdims[2])); + conduit::int32 *destPtr = n_field.as_int32_ptr(); + axom::for_all( + rdims[2], + AXOM_LAMBDA(int k) { + const auto t = static_cast(k) / (dims[2] - 1); + // Interpolate the window + const conduit::float64 x0 = axom::utilities::lerp(x_min, x1_min, t); + const conduit::float64 x1 = axom::utilities::lerp(x_max, x1_max, t); + const conduit::float64 y0 = axom::utilities::lerp(y_min, y1_min, t); + const conduit::float64 y1 = axom::utilities::lerp(y_max, y1_max, t); + conduit::Node n_rmesh; + conduit::blueprint::mesh::examples::julia(rdims[0], rdims[1], x0, x1, y0, y1, c_re, c_im, n_rmesh); + const conduit::Node &n_src_field = n_rmesh["fields/iters/values"]; + const conduit::int32 *srcPtr = n_src_field.as_int32_ptr(); + conduit::int32 *currentDestPtr = destPtr + k * rdims[0] * rdims[1]; + axom::copy(currentDestPtr, srcPtr, rdims[0] * rdims[1] * sizeof(conduit::int32)); + SLIC_INFO(axom::fmt::format("Made slice {}/{}", k + 1, rdims[2])); + }); + + // Make a matset based on the higher resolution julia field. + heavily_mixed_matset("topo", dims, refinement, n_mesh, n_field, nmats); + } + else + { + // Generate the same julia set at higher resolution to use as materials. + conduit::Node n_rmesh; + conduit::blueprint::mesh::examples::julia(rdims[0], + rdims[1], + x_min, + x_max, + y_min, + y_max, + c_re, + c_im, + n_rmesh); + + // Make a matset based on the higher resolution julia field. + const conduit::Node &n_field = n_rmesh["fields/iters/values"]; + heavily_mixed_matset("topo", dims, refinement, n_mesh, n_field, nmats); + } +} + +} // end namespace detail + +//-------------------------------------------------------------------------------- +HMApplication::HMApplication() + : m_handler(true) + , m_dims {100, 100, 1} + , m_numMaterials(40) + , m_refinement(40) + , m_numTrials(1) + , m_writeFiles(true) + , m_outputFilePath("output") + , m_method("elvira") + , m_policy(RuntimePolicy::seq) + , m_annotationMode("report") +{ } + +//-------------------------------------------------------------------------------- +int HMApplication::initialize(int argc, char **argv) +{ + axom::CLI::App app; + app.add_flag("--handler", m_handler) + ->description("Install a custom error handler that loops forever.") + ->capture_default_str(); + + std::vector dims; + app.add_option("--dims", dims, "Dimensions in x,y,z (z optional)") + ->expected(2, 3) + ->check(axom::CLI::PositiveNumber); + + app.add_option("--method", m_method) + ->description("The MIR method (or operation) name (equiz, elvira, traversal)"); + app.add_option("--materials", m_numMaterials) + ->check(axom::CLI::PositiveNumber) + ->description("The number of materials to create."); + app.add_option("--refinement", m_refinement) + ->check(axom::CLI::PositiveNumber) + ->description("The refinement within a zone when creating materials."); + app.add_option("--output", m_outputFilePath) + ->description("The file path for HDF5/YAML output files"); + bool disable_write = !m_writeFiles; + app.add_flag("--disable-write", disable_write)->description("Disable writing data files"); + app.add_option("--trials", m_numTrials) + ->check(axom::CLI::PositiveNumber) + ->description("The number of MIR trials to run on the mesh."); + +#if defined(AXOM_USE_CALIPER) + app.add_option("--caliper", m_annotationMode) + ->description( + "caliper annotation mode. Valid options include 'none' and 'report'. " + "Use 'help' to see full list.") + ->capture_default_str() + ->check(axom::utilities::ValidCaliperMode); +#endif + + std::stringstream pol_sstr; + pol_sstr << "Set MIR runtime policy method."; +#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) + pol_sstr << "\nSet to 'seq' or 0 to use the RAJA sequential policy."; + #ifdef AXOM_USE_OPENMP + pol_sstr << "\nSet to 'omp' or 1 to use the RAJA OpenMP policy."; + #endif + #ifdef AXOM_USE_CUDA + pol_sstr << "\nSet to 'cuda' or 2 to use the RAJA CUDA policy."; + #endif + #ifdef AXOM_USE_HIP + pol_sstr << "\nSet to 'hip' or 3 to use the RAJA HIP policy."; + #endif +#endif + app.add_option("-p, --policy", m_policy, pol_sstr.str()) + ->capture_default_str() + ->transform(axom::CLI::CheckedTransformer(axom::runtime_policy::s_nameToPolicy)); + + // Parse command line options. + int retval = 0; + try + { + app.parse(argc, argv); + m_writeFiles = !disable_write; + } + catch(axom::CLI::CallForHelp &e) + { + std::cout << app.help() << std::endl; + retval = -1; + } + catch(axom::CLI::ParseError &e) + { + // Handle other parsing errors + std::cerr << e.what() << std::endl; + retval = -2; + } + + for(size_t i = 0; i < axom::utilities::min(dims.size(), static_cast(3)); i++) + { + m_dims[i] = dims[i]; + } + + return retval; +} + +//-------------------------------------------------------------------------------- +int HMApplication::execute() +{ + axom::slic::SimpleLogger logger(axom::slic::message::Info); + axom::slic::setLoggingMsgLevel(axom::slic::message::Debug); + + if(m_handler) + { + conduit::utils::set_error_handler(conduit_debug_err_handler); + } +#if defined(AXOM_USE_CALIPER) + axom::utilities::raii::AnnotationsWrapper annotations_raii_wrapper(m_annotationMode); +#endif + int retval = 0; + try + { + retval = runMIR(); + } + catch(std::invalid_argument const &e) + { + SLIC_WARNING("Bad input. " << e.what()); + retval = -2; + } + catch(std::out_of_range const &e) + { + SLIC_WARNING("Integer overflow. " << e.what()); + retval = -3; + } + return retval; +} + +//-------------------------------------------------------------------------------- +int HMApplication::runMIR() +{ + // Initialize a mesh for testing MIR + auto timer = axom::utilities::Timer(true); + conduit::Node mesh; + { + AXOM_ANNOTATE_SCOPE("generate"); + int dims[3]; + // NOTE: Use axom::copy to copy m_dims into dims. This way, the Umpire resource + // manager gets created now so by the time we have multiple threads using + // axom::copy, there is no race condition. + axom::copy(dims, m_dims, sizeof(int) * 3); + + SLIC_INFO(axom::fmt::format("dims: {},{},{}", dims[0], dims[1], dims[2])); + SLIC_INFO(axom::fmt::format("refinement: {}", m_refinement)); + SLIC_INFO(axom::fmt::format("numMaterials: {}", m_numMaterials)); +#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) && defined(AXOM_USE_OPENMP) + using CPUExecSpace = axom::OMP_EXEC; +#else + using CPUExecSpace = axom::SEQ_EXEC; +#endif + detail::heavily_mixed(mesh, dims, m_refinement, m_numMaterials); + } + timer.stop(); + SLIC_INFO("Mesh init time: " << timer.elapsedTimeInMilliSec() << " ms."); + + // Output initial mesh. + if(m_writeFiles) + { + saveMesh(mesh, "heavily_mixed"); + } + + // Begin material interface reconstruction + timer.reset(); + timer.start(); + conduit::Node options, resultMesh; + options["matset"] = "mat"; + options["method"] = m_method; // pass method via options. + options["trials"] = m_numTrials; + + const int dimension = (m_dims[2] > 1) ? 3 : 2; + int retval = 0; + if(m_policy == RuntimePolicy::seq) + { + retval = runMIR_seq(dimension, mesh, options, resultMesh); + } +#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) + #if defined(AXOM_USE_OPENMP) + else if(m_policy == RuntimePolicy::omp) + { + retval = runMIR_omp(dimension, mesh, options, resultMesh); + } + #endif + #if defined(AXOM_USE_CUDA) + else if(m_policy == RuntimePolicy::cuda) + { + constexpr int CUDA_BLOCK_SIZE = 256; + using cuda_exec = axom::CUDA_EXEC; + retval = runMIR_cuda(dimension, mesh, options, resultMesh); + } + #endif + #if defined(AXOM_USE_HIP) + else if(m_policy == RuntimePolicy::hip) + { + retval = runMIR_hip(dimension, mesh, options, resultMesh); + } + #endif +#endif + else + { + retval = -1; + SLIC_ERROR("Unhandled policy."); + } + timer.stop(); + SLIC_INFO("Material interface reconstruction time: " << timer.elapsedTimeInMilliSec() << " ms."); + + // Output results + if(m_writeFiles) + { + AXOM_ANNOTATE_SCOPE("save_output"); + saveMesh(resultMesh, m_outputFilePath); + } + + return retval; +} + +//-------------------------------------------------------------------------------- +void HMApplication::adjustMesh(conduit::Node &) { } + +//-------------------------------------------------------------------------------- +void HMApplication::saveMesh(const conduit::Node &n_mesh, const std::string &path) +{ +#if defined(CONDUIT_RELAY_IO_HDF5_ENABLED) + std::string protocol("hdf5"); +#else + std::string protocol("yaml"); +#endif + conduit::relay::io::blueprint::save_mesh(n_mesh, path, protocol); +} + +//-------------------------------------------------------------------------------- +void HMApplication::conduit_debug_err_handler(const std::string &s1, const std::string &s2, int i1) +{ + SLIC_ERROR(axom::fmt::format("Error from Conduit: s1={}, s2={}, i1={}", s1, s2, i1)); + // This is on purpose. + while(1); +} diff --git a/src/axom/mir/examples/heavily_mixed/HMApplication.hpp b/src/axom/mir/examples/heavily_mixed/HMApplication.hpp new file mode 100644 index 0000000000..2f8fc7787c --- /dev/null +++ b/src/axom/mir/examples/heavily_mixed/HMApplication.hpp @@ -0,0 +1,73 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) +#ifndef AXOM_MIR_EXAMPLES_HEAVILY_MIXED_APPLICATION_HPP +#define AXOM_MIR_EXAMPLES_HEAVILY_MIXED_APPLICATION_HPP +#include "axom/config.hpp" +#include "axom/core.hpp" + +#include + +#include + +/*! + * \brief The heavily mixed materials application. + */ +class HMApplication +{ +public: + /// Constructor + HMApplication(); + + /*! + * \brief Initialize the application from command line args. + * \return 0 on success; less than zero otherwise. + */ + int initialize(int argc, char **argv); + + /*! + * \brief Execute the main application logic. + * \return 0 on success; less than zero otherwise. + */ + int execute(); + +protected: + /*! + * \brief Invoke the MIR appropriate for the selected runtime policy. + * \return 0 on success; less than zero otherwise. + */ + int runMIR(); + + /*! + * \brief Make any adjustments to the mesh. + */ + virtual void adjustMesh(conduit::Node &); + + /*! + * \brief Save the mesh to a file. + * + * \param path The filepath where the file will be saved. + * \param n_mesh The mesh to be saved. + */ + virtual void saveMesh(const conduit::Node &n_mesh, const std::string &path); + + /*! + * \brief A static error handler for Conduit. + */ + static void conduit_debug_err_handler(const std::string &s1, const std::string &s2, int i1); + + bool m_handler; + axom::StackArray m_dims; + int m_numMaterials; + int m_refinement; + int m_numTrials; + bool m_writeFiles; + std::string m_outputFilePath; + std::string m_method; + axom::runtime_policy::Policy m_policy; + std::string m_annotationMode; + std::string m_protocol; +}; + +#endif diff --git a/src/axom/mir/examples/heavily_mixed/mir_heavily_mixed.cpp b/src/axom/mir/examples/heavily_mixed/mir_heavily_mixed.cpp new file mode 100644 index 0000000000..3df9133de0 --- /dev/null +++ b/src/axom/mir/examples/heavily_mixed/mir_heavily_mixed.cpp @@ -0,0 +1,18 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include "HMApplication.hpp" + +int main(int argc, char **argv) +{ + HMApplication app; + int retval = app.initialize(argc, argv); + if(retval == 0) + { + retval = app.execute(); + } + + return retval; +} diff --git a/src/axom/mir/examples/heavily_mixed/runMIR.hpp b/src/axom/mir/examples/heavily_mixed/runMIR.hpp new file mode 100644 index 0000000000..b87527ba5c --- /dev/null +++ b/src/axom/mir/examples/heavily_mixed/runMIR.hpp @@ -0,0 +1,120 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) +#ifndef AXOM_MIR_EXAMPLES_HEAVILY_MIXED_RUNMIR_HPP +#define AXOM_MIR_EXAMPLES_HEAVILY_MIXED_RUNMIR_HPP +#include "axom/config.hpp" +#include "axom/core.hpp" +#include "axom/slic.hpp" +#include "axom/bump.hpp" +#include "axom/mir.hpp" + +template +int runMIR(const conduit::Node &hostMesh, const conduit::Node &options, conduit::Node &hostResult) +{ + AXOM_ANNOTATE_SCOPE("runMIR"); + + namespace utils = axom::bump::utilities; + + // Pick the method out of the options. + std::string method("equiz"); + if(options.has_child("method")) + { + method = options["method"].as_string(); + } + SLIC_INFO(axom::fmt::format("Using policy {} for {} {}D", + axom::execution_space::name(), + method, + NDIMS)); + + // Get the number of times we want to run MIR. + int trials = 1; + if(options.has_child("trials")) + { + trials = std::max(1, options["trials"].to_int()); + } + + // Check materials. + constexpr int MAXMATERIALS = 100; + auto materialInfo = axom::bump::views::materials(hostMesh["matsets/mat"]); + if(materialInfo.size() >= MAXMATERIALS) + { + SLIC_WARNING( + axom::fmt::format("To use more than {} materials, recompile with " + "larger MAXMATERIALS value.", + MAXMATERIALS)); + return -4; + } + + conduit::Node deviceMesh; + { + AXOM_ANNOTATE_SCOPE("host->device"); + utils::copy(deviceMesh, hostMesh); + } + + const conduit::Node &n_coordset = deviceMesh["coordsets/coords"]; + const conduit::Node &n_topology = deviceMesh["topologies/topo"]; + const conduit::Node &n_matset = deviceMesh["matsets/mat"]; + conduit::Node deviceResult; + for(int trial = 0; trial < trials; trial++) + { + deviceResult.reset(); + + // Make views + using namespace axom::bump::views; + auto coordsetView = make_rectilinear_coordset::view(n_coordset); + using CoordsetView = decltype(coordsetView); + + auto topologyView = make_rectilinear_topology::view(n_topology); + using TopologyView = decltype(topologyView); + + auto matsetView = make_unibuffer_matset::view(n_matset); + using MatsetView = decltype(matsetView); + + if(method == "equiz") + { + using MIR = axom::mir::EquiZAlgorithm; + MIR m(topologyView, coordsetView, matsetView); + m.execute(deviceMesh, options, deviceResult); + } + else if(method == "elvira") + { + using IndexingPolicy = typename TopologyView::IndexingPolicy; + using MIR = axom::mir::ElviraAlgorithm; + MIR m(topologyView, coordsetView, matsetView); + m.execute(deviceMesh, options, deviceResult); + } + else + { + SLIC_ERROR(axom::fmt::format("Unsupported MIR method {}", method)); + } + } + + { + AXOM_ANNOTATE_SCOPE("device->host"); + utils::copy(hostResult, deviceResult); + } + + return 0; +} + +// Prototypes. +int runMIR_seq(int dimension, + const conduit::Node &mesh, + const conduit::Node &options, + conduit::Node &result); +int runMIR_omp(int dimension, + const conduit::Node &mesh, + const conduit::Node &options, + conduit::Node &result); +int runMIR_cuda(int dimension, + const conduit::Node &mesh, + const conduit::Node &options, + conduit::Node &result); +int runMIR_hip(int dimension, + const conduit::Node &mesh, + const conduit::Node &options, + conduit::Node &result); + +#endif diff --git a/src/axom/mir/examples/heavily_mixed/runMIR_cuda.cpp b/src/axom/mir/examples/heavily_mixed/runMIR_cuda.cpp new file mode 100644 index 0000000000..e7ef3f0188 --- /dev/null +++ b/src/axom/mir/examples/heavily_mixed/runMIR_cuda.cpp @@ -0,0 +1,34 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) +#include "runMIR.hpp" + +#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) && defined(AXOM_USE_CUDA) +int runMIR_cuda(int dimension, + const conduit::Node &mesh, + const conduit::Node &options, + conduit::Node &result) +{ + constexpr int CUDA_BLOCK_SIZE = 256; + using cuda_exec = axom::CUDA_EXEC; + int retval = 0; + if(dimension == 3) + { + retval = runMIR(mesh, options, result); + } + else + { + retval = runMIR(mesh, options, result); + } + return retval; +} +#else +int runMIR_cuda(int AXOM_UNUSED_PARAM(dimension), + const conduit::Node &AXOM_UNUSED_PARAM(mesh), + const conduit::Node &AXOM_UNUSED_PARAM(options), + conduit::Node &AXOM_UNUSED_PARAM(result)) +{ + return 0; +} +#endif diff --git a/src/axom/mir/examples/heavily_mixed/runMIR_hip.cpp b/src/axom/mir/examples/heavily_mixed/runMIR_hip.cpp new file mode 100644 index 0000000000..b03bff1b8f --- /dev/null +++ b/src/axom/mir/examples/heavily_mixed/runMIR_hip.cpp @@ -0,0 +1,34 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) +#include "runMIR.hpp" + +#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) && defined(AXOM_USE_HIP) +int runMIR_hip(int dimension, + const conduit::Node &mesh, + const conduit::Node &options, + conduit::Node &result) +{ + constexpr int HIP_BLOCK_SIZE = 64; + using hip_exec = axom::HIP_EXEC; + int retval = 0; + if(dimension == 3) + { + retval = runMIR(mesh, options, result); + } + else + { + retval = runMIR(mesh, options, result); + } + return retval; +} +#else +int runMIR_hip(int AXOM_UNUSED_PARAM(dimension), + const conduit::Node &AXOM_UNUSED_PARAM(mesh), + const conduit::Node &AXOM_UNUSED_PARAM(options), + conduit::Node &AXOM_UNUSED_PARAM(result)) +{ + return 0; +} +#endif diff --git a/src/axom/mir/examples/heavily_mixed/runMIR_omp.cpp b/src/axom/mir/examples/heavily_mixed/runMIR_omp.cpp new file mode 100644 index 0000000000..6b57981575 --- /dev/null +++ b/src/axom/mir/examples/heavily_mixed/runMIR_omp.cpp @@ -0,0 +1,32 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) +#include "runMIR.hpp" + +#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) && defined(AXOM_USE_OPENMP) +int runMIR_omp(int dimension, + const conduit::Node &mesh, + const conduit::Node &options, + conduit::Node &result) +{ + int retval = 0; + if(dimension == 3) + { + retval = runMIR(mesh, options, result); + } + else + { + retval = runMIR(mesh, options, result); + } + return retval; +} +#else +int runMIR_omp(int AXOM_UNUSED_PARAM(dimension), + const conduit::Node &AXOM_UNUSED_PARAM(mesh), + const conduit::Node &AXOM_UNUSED_PARAM(options), + conduit::Node &AXOM_UNUSED_PARAM(result)) +{ + return 0; +} +#endif diff --git a/src/axom/mir/examples/heavily_mixed/runMIR_seq.cpp b/src/axom/mir/examples/heavily_mixed/runMIR_seq.cpp new file mode 100644 index 0000000000..af38756a9c --- /dev/null +++ b/src/axom/mir/examples/heavily_mixed/runMIR_seq.cpp @@ -0,0 +1,22 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) +#include "runMIR.hpp" + +int runMIR_seq(int dimension, + const conduit::Node &mesh, + const conduit::Node &options, + conduit::Node &result) +{ + int retval = 0; + if(dimension == 3) + { + retval = runMIR(mesh, options, result); + } + else + { + retval = runMIR(mesh, options, result); + } + return retval; +} diff --git a/src/axom/mir/tests/mir_elvira3d.cpp b/src/axom/mir/tests/mir_elvira3d.cpp index 52707f9824..ec46a03774 100644 --- a/src/axom/mir/tests/mir_elvira3d.cpp +++ b/src/axom/mir/tests/mir_elvira3d.cpp @@ -300,20 +300,17 @@ struct test_Elvira3D axom::for_all( matsetView.numberOfZones(), AXOM_LAMBDA(axom::IndexType zi) { - // Get the materials for this zone. - typename MatsetView::IDList ids; - typename MatsetView::VFList vfs; - matsetView.zoneMaterials(zi, ids, vfs); - // Add the material volumes to the total volumes. - for(axom::IndexType i = 0; i < ids.size(); i++) + const auto end = matsetView.endZone(zi); + for(auto zoneMat = matsetView.beginZone(zi); zoneMat != end; zoneMat++) { - auto index = axom::utilities::binary_search(sortedIdsView, ids[i]); + auto index = axom::utilities::binary_search(sortedIdsView, zoneMat.material_id()); // RelWithDebInfo workaround - "sortedIdsView.size()" substitutes lambda capture device failure for "nmats" SLIC_ASSERT(index >= 0 && index < sortedIdsView.size()); // Use an atomic to sum the value. - axom::atomicAdd(totalVolumeView.data() + index, zoneVolumes[zi] * vfs[i]); + axom::atomicAdd(totalVolumeView.data() + index, + zoneVolumes[zi] * zoneMat.volume_fraction()); } });