From d6a18429b499cb1ce9803b941890f3e674f7db05 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Thu, 13 Aug 2026 18:15:50 -0700 Subject: [PATCH 01/16] Quest: Adds support to DCP example to read in blueprint object mesh When provided, we use that instead of an analytically generated mesh. --- ...est_distributed_distance_query_example.cpp | 68 ++++++++++++------- 1 file changed, 45 insertions(+), 23 deletions(-) diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index 8377c667f7..0f7a16c7db 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -38,6 +38,7 @@ // C/C++ includes #include #include +#include #include #include @@ -67,6 +68,8 @@ struct Input std::string distanceFile {"cp_coords"}; std::string objectFile {"object_mesh"}; + std::string objectMeshFile; + double circleRadius {1.0}; std::vector circleCenter {0.0, 0.0}; int longPointCount {100}; @@ -131,6 +134,12 @@ struct Input ->description("Name of output file containing object mesh.") ->capture_default_str(); + app.add_option("--object-mesh-file", objectMeshFile) + ->description( + "Path to a conduit blueprint point mesh root file. " + " When provided, we use this instead of an analytically generated object mesh.") + ->check(axom::CLI::ExistingFile); + app.add_flag("-v,--verbose,!--no-verbose", m_verboseOutput) ->description("Enable/disable verbose output") ->capture_default_str(); @@ -721,6 +730,16 @@ class ObjectMeshWrapper SLIC_ASSERT(group != nullptr); } + //!@brief Construct by reading a blueprint object mesh from disk. + //! Uses the topology/coordset names found in the file, exactly like the + //! query mesh reader. (The empty-topology BlueprintParticleMesh ctor lets + //! read_blueprint_mesh pick the file's actual topology.) + ObjectMeshWrapper(sidre::Group* group, const std::string& meshFilename) : m_objectMesh(group) + { + SLIC_ASSERT(group != nullptr); + m_objectMesh.read_blueprint_mesh(meshFilename); + } + BlueprintParticleMesh& getParticleMesh() { return m_objectMesh; } /// Get a pointer to the root group for this mesh @@ -1303,15 +1322,33 @@ int main(int argc, char** argv) slic::flushStreams(); const size_t spatialDim = queryMeshWrapper.getParticleMesh().dimension(); - SLIC_ASSERT(params.circleCenter.size() == spatialDim); + + const bool readObjectFromFile = !params.objectMeshFile.empty(); + + // The analytic circle/sphere generator needs a center whose dimension matches the query mesh. + // only enforce that when we are actually generating. + SLIC_ASSERT(readObjectFromFile || params.circleCenter.size() == spatialDim); //--------------------------------------------------------------------------- - // Generate object mesh + // Object (second) mesh: read from file, or generate analytically //--------------------------------------------------------------------------- - ObjectMeshWrapper objectMeshWrapper(dataStore.getRoot()->createGroup("object_mesh", true)); + std::unique_ptr objectMeshWrapperPtr; + if(readObjectFromFile) + { + objectMeshWrapperPtr = + std::make_unique(dataStore.getRoot()->createGroup("object_mesh", true), + params.objectMeshFile); + } + else + { + objectMeshWrapperPtr = + std::make_unique(dataStore.getRoot()->createGroup("object_mesh", true)); + } + ObjectMeshWrapper& objectMeshWrapper = *objectMeshWrapperPtr; objectMeshWrapper.setVerbosity(params.isVerbose()); + if(!readObjectFromFile) { SLIC_ASSERT(params.objDomainCountRange[1] >= params.objDomainCountRange[0]); const unsigned int omin = params.objDomainCountRange[0]; @@ -1336,7 +1373,11 @@ int main(int argc, char** argv) } slic::flushStreams(); - objectMeshWrapper.saveMesh(params.objectFile); + // Only re-save the generated object; a file-read object is already on disk. + if(!readObjectFromFile) + { + objectMeshWrapper.saveMesh(params.objectFile); + } slic::flushStreams(); //--------------------------------------------------------------------------- @@ -1360,25 +1401,6 @@ int main(int argc, char** argv) conduit::Node queryMeshNode; queryMeshWrapper.getBlueprintGroup()->createNativeLayout(queryMeshNode); - // To test with contiguous and interleaved coordinate storage, - // make half them contiguous. - for(int di = 0; di < objectMeshNode.number_of_children(); ++di) - { - auto& dom = objectMeshNode.child(di); - if((my_rank + di) % 2 == 1) - { - make_coords_contiguous(dom.fetch_existing("coordsets/coords/values")); - } - } - for(int di = 0; di < queryMeshNode.number_of_children(); ++di) - { - auto& dom = queryMeshNode.child(di); - if((my_rank + di) % 2 == 1) - { - make_coords_contiguous(dom.fetch_existing("coordsets/coords/values")); - } - } - // Create distributed closest point query object and set some parameters quest::DistributedClosestPoint query; query.setRuntimePolicy(params.policy); From bdd543608bd46849e722c2dde1604ba5316cfef1 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Thu, 13 Aug 2026 18:24:40 -0700 Subject: [PATCH 02/16] Quest: Increase flexibilty for blueprint mesh in DCP example Allow for topologies other than "mesh" and for unstructured topologies. --- ...est_distributed_distance_query_example.cpp | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index 0f7a16c7db..b5734cf28e 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -314,24 +314,25 @@ struct BlueprintParticleMesh if(domCount > 0) { + if(m_topologyName.empty()) + { + // No topology given. Pick the first one. + m_topologyName = mdMesh[0].fetch_existing("topologies")[0].name(); + } + auto topologyPath = axom::fmt::format("topologies/{}", m_topologyName); + + // Detect strided structured coordinates. + // Structured topologies have elements/dims, but points and unstructured topology don't. + // Guard the probe use the resolved topology name rather than a hardcoded "mesh" topology. + const std::string dimsPath = topologyPath + "/elements/dims"; m_coordsAreStrided = - mdMesh[0].fetch_existing("topologies/mesh/elements/dims").has_child("strides"); + mdMesh[0].has_path(dimsPath) && mdMesh[0].fetch_existing(dimsPath).has_child("strides"); if(m_coordsAreStrided) { SLIC_WARNING( axom::fmt::format("Mesh '{}' is strided. Stride support is under development.", meshFilename)); } - } - - if(domCount > 0) - { - if(m_topologyName.empty()) - { - // No topology given. Pick the first one. - m_topologyName = mdMesh[0].fetch_existing("topologies")[0].name(); - } - auto topologyPath = axom::fmt::format("topologies/{}", m_topologyName); m_coordsetName = mdMesh[0].fetch_existing(topologyPath + "/coordset").as_string(); const conduit::Node coordsetNode = From 1c8f7a39e6d49a599ec252414a1bafa705e5d7b2 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Thu, 13 Aug 2026 18:40:57 -0700 Subject: [PATCH 03/16] Quest: Adds some memory diagnostics to DCP example --- ...est_distributed_distance_query_example.cpp | 327 ++++++++++++++++++ 1 file changed, 327 insertions(+) diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index b5734cf28e..f7dbc7f480 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -41,6 +41,14 @@ #include #include #include +#include +#include +#include +#include +#include +#if defined(__GLIBC__) + #include // mallinfo2 / mallinfo / malloc_trim +#endif namespace quest = axom::quest; namespace slic = axom::slic; @@ -70,6 +78,11 @@ struct Input std::string objectMeshFile; + // Memory instrumentation (all off by default). + bool trackMemory {false}; // report RSS / malloc / umpire at phase boundaries + bool trimAfterQuery {false}; // malloc_trim(0) after the query, then re-report + int sampleMemoryMs {0}; // if > 0, background-sample peak RSS during the query + double circleRadius {1.0}; std::vector circleCenter {0.0, 0.0}; int longPointCount {100}; @@ -140,6 +153,26 @@ struct Input " When provided, we use this instead of an analytically generated object mesh.") ->check(axom::CLI::ExistingFile); + app.add_flag("--track-memory", trackMemory) + ->description( + "Report RSS, glibc malloc live/arena bytes (and Umpire high-water when available) " + "before/after BVH build and before/after the closest-point query, reduced across ranks.") + ->capture_default_str(); + + app.add_flag("--trim-after-query", trimAfterQuery) + ->description( + "After the query, call malloc_trim(0) and report memory again. " + "If RSS drops here (while malloc-live was already back at baseline), " + "the memory was arena retention, not a leak. Implies --track-memory.") + ->capture_default_str(); + + app.add_option("--sample-memory-ms", sampleMemoryMs) + ->description( + "If > 0, run a background thread sampling RSS at this interval (ms) " + "during the query and report the peak. Useful for observing in-flight " + "send-buffer accumulation. Implies --track-memory.") + ->capture_default_str(); + app.add_flag("-v,--verbose,!--no-verbose", m_verboseOutput) ->description("Enable/disable verbose output") ->capture_default_str(); @@ -1206,6 +1239,270 @@ void make_coords_interleaved(conduit::Node& coordValues) } } +//----------------------------------------------------------------------------- +// Optional memory instrumentation +// +// Purpose: distinguish a genuine leak (memory still referenced after the call) +// from arena retention (memory freed at the allocator but not returned to the OS). +// RSS shows what the process actually holds; glibc "malloc live" shows what is still allocated; +// "malloc arena" shows free memory kept in the arena. +// * RSS after call ~= baseline -> transient, no problem +// * RSS high, malloc-live ~= baseline -> arena retention (malloc_trim should reclaim it) +// * malloc-live high after call -> genuine leak +// Everything here is Linux/glibc-specific and degrades to "n/a" elsewhere. +// All of it is opt-in (--track-memory) and off by default. +//----------------------------------------------------------------------------- + +/// Read current (VmRSS) and peak (VmHWM) resident set size in bytes; -1 if n/a. +inline void readProcRss(long long& rssBytes, long long& peakRssBytes) +{ + rssBytes = -1; + peakRssBytes = -1; +#if defined(__linux__) + std::ifstream status("/proc/self/status"); + std::string line; + while(std::getline(status, line)) + { + long long kb = 0; + if(std::sscanf(line.c_str(), "VmRSS: %lld kB", &kb) == 1) + { + rssBytes = kb * 1024; + } + else if(std::sscanf(line.c_str(), "VmHWM: %lld kB", &kb) == 1) + { + peakRssBytes = kb * 1024; + } + } +#endif +} + +/// Reset the kernel's peak-RSS high-water mark (VmHWM) to the current RSS, so a +/// later VmHWM read reflects the peak of just the intervening phase. +/// Best-effort: requires Linux clear_refs type 5 (kernel >= 4.0). +inline void resetPeakRss() +{ +#if defined(__linux__) + std::ofstream clear("/proc/self/clear_refs"); + if(clear) + { + clear << "5\n"; + } +#endif +} + +/// Human-readable byte count. +inline std::string humanBytes(long long b) +{ + if(b < 0) + { + return "n/a"; + } + const char* units[] = {"B", "KiB", "MiB", "GiB", "TiB"}; + double v = static_cast(b); + int i = 0; + while(v >= 1024.0 && i < 4) + { + v /= 1024.0; + ++i; + } + return axom::fmt::format("{:.2f} {}", v, units[i]); +} + +/// MPI-reduce a per-rank byte count to its max (with the rank achieving it) and +/// its sum across ranks. A negative local value means "unavailable". +inline void reduceBytes(long long local, long long& maxVal, int& maxRank, long long& sumVal) +{ + struct + { + long val; + int rank; + } in {static_cast(local), my_rank}, out {0, 0}; + MPI_Allreduce(&in, &out, 1, MPI_LONG_INT, MPI_MAXLOC, MPI_COMM_WORLD); + maxVal = out.val; + maxRank = out.rank; + MPI_Allreduce(&local, &sumVal, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); +} + +/*! + * \brief Opt-in per-run memory probe for the closest-point query. + * + * Reports RSS, peak RSS, glibc live/arena bytes (and Umpire high-water when built with Umpire), + * reduced across ranks (total and hottest rank). Also offers a background sampler to capture + * the peak RSS *during* a phase (useful for observing in-flight send-buffer accumulation), + * and a resetPeak() to make VmHWM phase-local. + */ +class MemoryProbe +{ +public: + MemoryProbe(bool enabled, int sampleMs, int umpireAllocatorId = -1) + : m_enabled(enabled) + , m_sampleMs(sampleMs) + , m_umpireAllocatorId(umpireAllocatorId) + { } + + bool enabled() const { return m_enabled; } + + /// Reset VmHWM so the next report()'s peak reflects only the next phase. + void resetPeak() + { + if(m_enabled) + { + resetPeakRss(); + } + } + + /// Start a background thread sampling RSS every sampleMs; no-op if disabled or sampleMs <= 0 + /// Records the peak RSS seen until stopSampler(). + void startSampler() + { + if(!m_enabled || m_sampleMs <= 0) + { + return; + } + m_samplerPeak = 0; + m_stopSampler.store(false, std::memory_order_relaxed); + m_samplerThread = std::thread([this]() { + while(!m_stopSampler.load(std::memory_order_relaxed)) + { + long long rss = -1, peak = -1; + readProcRss(rss, peak); + if(rss > m_samplerPeak) + { + m_samplerPeak = rss; + } + std::this_thread::sleep_for(std::chrono::milliseconds(m_sampleMs)); + } + }); + } + + /// Stop the sampler and report the peak RSS during the sampled phase. + void stopSampler(const std::string& label) + { + if(!m_enabled || m_sampleMs <= 0) + { + return; + } + m_stopSampler.store(true, std::memory_order_relaxed); + if(m_samplerThread.joinable()) + { + m_samplerThread.join(); + } + long long rss = -1, peak = -1; // final read to catch the tail + readProcRss(rss, peak); + if(rss > m_samplerPeak) + { + m_samplerPeak = rss; + } + + long long maxVal, sumVal; + int maxRank; + reduceBytes(m_samplerPeak, maxVal, maxRank, sumVal); + if(my_rank == 0) + { + SLIC_INFO(axom::fmt::format( + "[mem] {}: peak RSS during phase (sampled @ {} ms): max/rank={} (rank {}), total={}", + label, + m_sampleMs, + humanBytes(maxVal), + maxRank, + humanBytes(sumVal))); + } + } + + /// Take a snapshot on every rank, reduce it, and print an aggregate (rank 0). + void report(const std::string& label) + { + if(!m_enabled) + { + return; + } + + long long rss = -1, peakRss = -1; + readProcRss(rss, peakRss); + + long long mallocLive = -1, mallocArena = -1; +#if defined(__GLIBC__) + #if defined(__GLIBC_PREREQ) && __GLIBC_PREREQ(2, 33) + struct mallinfo2 mi = mallinfo2(); // size_t fields: safe above 2 GiB + mallocLive = static_cast(mi.uordblks) + static_cast(mi.hblkhd); + mallocArena = static_cast(mi.arena); + #else + struct mallinfo mi = mallinfo(); // NOTE: int fields saturate above ~2 GiB + mallocLive = static_cast(mi.uordblks) + static_cast(mi.hblkhd); + mallocArena = static_cast(mi.arena); + #endif +#endif + + long long umpireCur = -1, umpireHwm = -1; +#if defined(AXOM_USE_UMPIRE) + if(m_umpireAllocatorId >= 0) + { + auto& rm = umpire::ResourceManager::getInstance(); + umpire::Allocator alloc = rm.getAllocator(m_umpireAllocatorId); + umpireCur = static_cast(alloc.getCurrentSize()); + umpireHwm = static_cast(alloc.getHighWatermark()); + } +#endif + + long long rssMax, rssSum, peakMax, peakSum, liveMax, liveSum, arenaMax, arenaSum; + long long umpCurMax, umpCurSum, umpHwmMax, umpHwmSum; + int rssMaxRank, peakMaxRank, liveMaxRank, arenaMaxRank, umpCurMaxRank, umpHwmMaxRank; + reduceBytes(rss, rssMax, rssMaxRank, rssSum); + reduceBytes(peakRss, peakMax, peakMaxRank, peakSum); + reduceBytes(mallocLive, liveMax, liveMaxRank, liveSum); + reduceBytes(mallocArena, arenaMax, arenaMaxRank, arenaSum); + reduceBytes(umpireCur, umpCurMax, umpCurMaxRank, umpCurSum); + reduceBytes(umpireHwm, umpHwmMax, umpHwmMaxRank, umpHwmSum); + + if(my_rank == 0) + { + std::string msg = axom::fmt::format( + "[mem] {} (total over {} ranks | max on one rank)\n" + " RSS : {:>11} | {:>11} (rank {})\n" + " RSS peak : {:>11} | {:>11} (rank {})\n" + " malloc live : {:>11} | {:>11} (rank {})\n" + " malloc arena : {:>11} | {:>11} (rank {})", + label, + num_ranks, + humanBytes(rssSum), + humanBytes(rssMax), + rssMaxRank, + humanBytes(peakSum), + humanBytes(peakMax), + peakMaxRank, + humanBytes(liveSum), + humanBytes(liveMax), + liveMaxRank, + humanBytes(arenaSum), + humanBytes(arenaMax), + arenaMaxRank); +#if defined(AXOM_USE_UMPIRE) + if(umpHwmMax >= 0) + { + msg += axom::fmt::format( + "\n umpire current: {:>11} | {:>11} (rank {})" + "\n umpire hi-water: {:>10} | {:>11} (rank {})", + humanBytes(umpCurSum), + humanBytes(umpCurMax), + umpCurMaxRank, + humanBytes(umpHwmSum), + humanBytes(umpHwmMax), + umpHwmMaxRank); + } +#endif + SLIC_INFO(msg); + } + } + +private: + bool m_enabled {false}; + int m_sampleMs {0}; + int m_umpireAllocatorId {-1}; + std::atomic m_stopSampler {false}; + std::thread m_samplerThread; + long long m_samplerPeak {0}; +}; + /// Utility function to initialize the logger void initializeLogger() { @@ -1416,6 +1713,15 @@ int main(int argc, char** argv) query.setObjectMesh(objectMeshNode.number_of_children() == 1 ? objectMeshNode[0] : objectMeshNode, objectMeshWrapper.getTopologyName()); + // Optional memory instrumentation around the index build and the query. + int memUmpireId = -1; +#if defined(AXOM_USE_UMPIRE) + memUmpireId = umpireAllocator.getId(); +#endif + const bool trackMem = params.trackMemory || params.trimAfterQuery || params.sampleMemoryMs > 0; + MemoryProbe memProbe(trackMem, params.sampleMemoryMs, memUmpireId); + memProbe.report("baseline (meshes read, before BVH)"); + // Build the spatial index over the object on each rank SLIC_INFO(init_str); slic::flushStreams(); @@ -1423,14 +1729,35 @@ int main(int argc, char** argv) query.generateBVHTree(); initTimer.stop(); + memProbe.report("after generateBVHTree (object index built)"); + // Run the distributed closest point query over the nodes of the computational mesh // To test support for single-domain format, use single-domain when possible. slic::flushStreams(); + memProbe.resetPeak(); // make the post-query VmHWM reflect only the query phase + memProbe.startSampler(); queryTimer.start(); query.computeClosestPoints( queryMeshNode.number_of_children() == 1 ? queryMeshNode[0] : queryMeshNode, queryMeshWrapper.getTopologyName()); queryTimer.stop(); + memProbe.stopSampler("computeClosestPoints"); + + memProbe.report("after computeClosestPoints"); + + // Optionally return free arena memory to the OS and re-measure. + // If RSS drops here while "malloc live" was already back at baseline, + // the retained memory was glibc arena, not a leak. + // (This is a diagnostic; the library-side fix would trim inside computeClosestPoints after releasing transfer buffers.) + if(params.trimAfterQuery) + { +#if defined(__GLIBC__) + ::malloc_trim(0); + memProbe.report("after malloc_trim(0)"); +#else + SLIC_WARNING("--trim-after-query requested but malloc_trim is glibc-only; skipping."); +#endif + } auto getDoubleMinMax = [](double inVal, double& minVal, double& maxVal, double& sumVal) { MPI_Allreduce(&inVal, &minVal, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); From 2013504f31f6655e19c31993b364015742f3b82c Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Thu, 13 Aug 2026 19:24:52 -0700 Subject: [PATCH 04/16] Quest: Adds support to DCP example for single-domain bp meshes And for input meshes without a "fields" group. --- .../detail/DistributedClosestPointImpl.hpp | 2 +- ...est_distributed_distance_query_example.cpp | 32 ++++++++++++------- 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index 5e7d427a99..a5841977b9 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -381,7 +381,7 @@ class DistributedClosestPointImpl { auto& queryDom = isMultidomain ? queryNode.child(domainNum) : queryNode; conduit::Node& xferDom = xferDoms.child(domainNum); - conduit::Node& fields = queryDom.fetch_existing("fields"); + conduit::Node& fields = queryDom.fetch("fields"); conduit::Node genericHeaders; genericHeaders["association"] = "vertex"; diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index f7dbc7f480..a0476f28ec 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -292,7 +292,8 @@ struct BlueprintParticleMesh /// Gets the parent group for the blueprint fields sidre::Group* fields_group(axom::IndexType groupIdx) const { - return domain_group(groupIdx)->getGroup("fields"); + auto* domain = domain_group(groupIdx); + return domain->hasGroup("fields") ? domain->getGroup("fields") : domain->createGroup("fields"); } const std::string& getTopologyName() const { return m_topologyName; } @@ -342,7 +343,13 @@ struct BlueprintParticleMesh conduit::Node mdMesh; conduit::relay::mpi::io::blueprint::load_mesh(meshFilename, mdMesh, MPI_COMM_WORLD); - assert(conduit::blueprint::mesh::is_multi_domain(mdMesh)); + if(!conduit::blueprint::mesh::is_multi_domain(mdMesh) && mdMesh.number_of_children() > 0) + { + conduit::Node singleDomainMesh; + singleDomainMesh.update(mdMesh); + mdMesh.reset(); + conduit::blueprint::mesh::to_multi_domain(singleDomainMesh, mdMesh); + } conduit::index_t domCount = conduit::blueprint::mesh::number_of_domains(mdMesh); if(domCount > 0) @@ -597,13 +604,13 @@ struct BlueprintParticleMesh { auto* domain = m_group->getGroup(domainIdx); auto* fields = domain->getGroup("fields"); - auto has = fields->hasGroup(fieldName); - return has; + return fields != nullptr && fields->hasGroup(fieldName); } bool hasVectorField(const std::string& fieldName, int domainIdx = 0) const { - return m_group->getGroup(domainIdx)->getGroup("fields")->hasGroup(fieldName); + auto* fields = m_group->getGroup(domainIdx)->getGroup("fields"); + return fields != nullptr && fields->hasGroup(fieldName); } template @@ -617,7 +624,7 @@ struct BlueprintParticleMesh auto* domain = m_group->getGroup(domainIdx); auto* fields = domain->getGroup("fields"); - auto* field = fields->getGroup(fieldName); + auto* field = fields != nullptr ? fields->getGroup(fieldName) : nullptr; T* data = field ? static_cast(field->getView("values")->getVoidPtr()) : nullptr; return field ? axom::ArrayView(data, numPoints(domainIdx)) : axom::ArrayView(); @@ -638,11 +645,11 @@ struct BlueprintParticleMesh // need to modify this implementation accordingly. T* data = nullptr; axom::IndexType npts = 0; - bool has = m_group->getGroup(domainIdx)->getGroup("fields")->hasGroup(fieldName); + auto* fields = m_group->getGroup(domainIdx)->getGroup("fields"); + bool has = fields != nullptr && fields->hasGroup(fieldName); if(has) { - auto xView = - m_group->getGroup(domainIdx)->getGroup("fields")->getGroup(fieldName)->getView("values/x"); + auto xView = fields->getGroup(fieldName)->getView("values/x"); data = static_cast(xView->getVoidPtr()); npts = xView->getNumElements(); } @@ -677,8 +684,8 @@ struct BlueprintParticleMesh sidre::Group* getDomain(axom::IndexType domain) { return m_group->getGroup(domain); } sidre::Group* getFields(axom::IndexType domainIdx) { - auto* fields = m_group->getGroup(domainIdx)->getGroup("fields"); - return fields; + auto* domain = m_group->getGroup(domainIdx); + return domain->hasGroup("fields") ? domain->getGroup("fields") : domain->createGroup("fields"); } /// Checks whether the blueprint is valid and prints diagnostics @@ -856,7 +863,8 @@ class QueryMeshWrapper sidre::Group& domGroup = *dstDomains->getGroup(d); const conduit::Node& domNode = isMultidomain ? node.child(d) : node; - sidre::Group& dstFieldsGroup = *domGroup.getGroup("fields"); + sidre::Group& dstFieldsGroup = + *(domGroup.hasGroup("fields") ? domGroup.getGroup("fields") : domGroup.createGroup("fields")); const conduit::Node& srcFieldsNode = domNode.fetch_existing("fields"); { if(!m_queryMesh.hasScalarField("cp_rank")) From 0d49848ee235642438963654b5ae0f1ad9ddc22f Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Thu, 13 Aug 2026 19:26:42 -0700 Subject: [PATCH 05/16] Adds Python script to generate blueprint Point meshes for DCP query --- src/tools/gen-multidom-point-mesh.py | 557 +++++++++++++++++++++++++++ 1 file changed, 557 insertions(+) create mode 100755 src/tools/gen-multidom-point-mesh.py diff --git a/src/tools/gen-multidom-point-mesh.py b/src/tools/gen-multidom-point-mesh.py new file mode 100755 index 0000000000..7778f2c75d --- /dev/null +++ b/src/tools/gen-multidom-point-mesh.py @@ -0,0 +1,557 @@ +#!/usr/bin/env python3 + +# gen-multidom-point-mesh.py +# Write a point mesh following the Conduit mesh blueprint. +# +# The generated Blueprint hierarchy is: +# (single-domain output) +# / (multidomain output) +# |-- state +# | `-- domain_id == +# |-- topologies +# | `-- +# | |-- coordset == +# | |-- type == "points" or "unstructured" +# | `-- elements (only for "unstructured") +# | |-- shape == "point" +# | |-- connectivity (int32, point_count) +# | |-- sizes (int32, point_count, all 1) +# | `-- offsets (int32, point_count) +# |-- coordsets +# | `-- +# | |-- type == "explicit" +# | `-- values +# | |-- x (float64, point_count) +# | |-- y (float64, point_count) +# | `-- [z] (float64, point_count, present in 3D) +# `-- fields (optional; omitted when there are no fields) +# `-- global_id (only with --id-field) +# |-- topology == +# |-- association == "vertex" +# `-- values (int64, point_count) + +# This script requires a conduit installation configured with python3 and hdf5. +# Make sure PYTHONPATH includes /path/to/conduit/install/python-modules, +# or use Axom's convenience script /path/to/axom_build_dir/bin/run_python_with_axom.sh +# that includes Conduit in PYTHONPATH. + +import itertools +import math +import sys +from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser, ArgumentTypeError + +try: + import conduit + import conduit.blueprint + import conduit.relay +except ModuleNotFoundError as e: + print( + f'{e}\nMake sure your PYTHONPATH includes /path/to/conduit/install/python-modules\n' + 'Conduit must be configured with python and hdf5.\n' + 'Alternatively, you can use the convenience script\n' + '/path/to/axom_build_dir/bin/run_python_with_axom.sh\n' + 'that includes Conduit in PYTHONPATH.' + ) + sys.exit(-1) + +try: + import numpy as np +except ModuleNotFoundError as e: + print(f'{e}\nThis script requires numpy.') + sys.exit(-1) + + +AXES = 'xyz' + + +def get_mpi_context(): + '''Return MPI context when mpi4py and Conduit MPI relay are available.''' + try: + from mpi4py import MPI + import conduit.blueprint.mpi.mesh + import conduit.relay.mpi + import conduit.relay.mpi.io + import conduit.relay.mpi.io.blueprint + except ModuleNotFoundError: + return {'enabled': False, 'comm': None, 'rank': 0, 'size': 1} + + comm = MPI.COMM_WORLD.py2f() + return { + 'enabled': True, + 'comm': comm, + 'rank': conduit.relay.mpi.rank(comm), + 'size': conduit.relay.mpi.size(comm), + } + + +def csv_values(s, converter, value_name): + '''Convert a comma-separated string to a list.''' + try: + values = [converter(token) for token in s.split(',')] + except ValueError as e: + raise ArgumentTypeError(f'{value_name} must be a comma-separated list') from e + + if len(values) == 0: + raise ArgumentTypeError(f'{value_name} must not be empty') + + return values + + +def i_c(s): + '''Convert comma-separated string to list of integers.''' + return csv_values(s, int, 'integer values') + + +def f_c(s): + '''Convert comma-separated string to list of floating point numbers.''' + return csv_values(s, float, 'floating point values') + + +def positive_int(s): + '''Convert a string to a positive integer.''' + try: + value = int(s) + except ValueError as e: + raise ArgumentTypeError('value must be an integer') from e + + if value < 1: + raise ArgumentTypeError('value must be positive') + + return value + + +def parse_args(): + parser = ArgumentParser( + description='Write a single-domain or multidomain blueprint point mesh.', + formatter_class=ArgumentDefaultsHelpFormatter) + parser.add_argument( + 'shape', + choices=('circle', 'sphere', 'torus', 'grid', 'gaussian', 'uniform'), + help='Point distribution to generate') + parser.add_argument( + '-d', + '--dimension', + type=int, + choices=(2, 3), + help='Spatial dimension for dimension-flexible shapes') + parser.add_argument( + '-n', + '--point-count', + type=positive_int, + default=1024, + help='Number of points, or target number for grid when --grid-size is omitted') + parser.add_argument( + '--grid-size', + type=i_c, + help='Point counts in each grid direction, e.g. 101,101 or 65,65,65') + parser.add_argument( + '-dc', + '--domain-counts', + type=i_c, + help='Domain counts in one chunk direction or each coordinate direction') + parser.add_argument( + '--domains', + type=positive_int, + help='Total domain count for a one-dimensional chunk partition') + parser.add_argument( + '--single-domain', + action='store_true', + help='Write a single-domain point mesh instead of multidomain output') + parser.add_argument('-ml', type=f_c, help='Lower coordinates for grid and uniform shapes') + parser.add_argument('-mu', type=f_c, help='Upper coordinates for grid and uniform shapes') + parser.add_argument('--center', type=f_c, help='Center for circle, sphere, torus, and gaussian') + parser.add_argument('--radius', type=float, default=1.0, help='Radius for circle or sphere') + parser.add_argument('--major-radius', type=float, default=1.0, help='Major radius for torus') + parser.add_argument('--minor-radius', type=float, default=0.25, help='Minor radius for torus') + parser.add_argument('--stddev', + type=f_c, + help='Gaussian standard deviation; scalar or one value per dimension') + parser.add_argument('--seed', type=int, default=0, help='Random seed for stochastic shapes') + parser.add_argument('--use-list', + action='store_true', + help='Put domains in a list instead of a map') + parser.add_argument('--topology-name', default='mesh', help='Blueprint topology name') + parser.add_argument('--coordset-name', default='coords', help='Blueprint coordset name') + parser.add_argument( + '--topology-type', + choices=('unstructured', 'points'), + default='unstructured', + help='Blueprint point topology representation') + parser.add_argument( + '--id-field', + action='store_true', + help='Add a vertex-associated int64 global_id field') + parser.add_argument('--protocol', default='hdf5', help='Conduit relay protocol for save_mesh') + parser.add_argument('-o', '--output', default='mdpoint', help='Output file base name') + parser.add_argument('-v', '--verbose', action='store_true', help='Print additional info') + + opts, unkn = parser.parse_known_args() + if opts.verbose: + print(opts, unkn) + if unkn: + print('Unrecognized arguments:', *unkn) + sys.exit(1) + + return opts + + +def infer_dimension(opts): + '''Infer and validate the spatial dimension.''' + fixed_dims = {'circle': 2, 'sphere': 3, 'torus': 3} + if opts.shape in fixed_dims: + dim = fixed_dims[opts.shape] + if opts.dimension is not None and opts.dimension != dim: + raise RuntimeError(f'{opts.shape} requires dimension {dim}') + return dim + + dim = opts.dimension + inferred = [] + for values in (opts.grid_size, opts.ml, opts.mu, opts.center, opts.domain_counts): + if values is not None: + inferred.append(len(values)) + + if dim is None and inferred: + dim = inferred[0] + if dim is None: + dim = 2 + + if dim not in (2, 3): + raise RuntimeError('dimension must be 2 or 3') + + for value_dim in inferred: + if value_dim not in (1, dim): + raise RuntimeError( + 'dimensioned options must have one value or one value per spatial dimension') + + return dim + + +def vector_option(values, dim, default, name): + '''Return a dimension-sized numpy vector from an optional scalar/list argument.''' + if values is None: + return np.full(dim, default, dtype=float) + if len(values) == 1: + return np.full(dim, values[0], dtype=float) + if len(values) != dim: + raise RuntimeError(f'{name} must have one value or {dim} values') + return np.array(values, dtype=float) + + +def domain_counts_from_options(opts, dim, mpi_size): + '''Return domain counts as an integer vector.''' + if opts.single_domain: + return np.array([1], dtype=int) + + if opts.domain_counts is not None and opts.domains is not None: + raise RuntimeError('Use --domain-counts or --domains, not both') + + if opts.domain_counts is not None: + counts = np.array(opts.domain_counts, dtype=int) + if len(counts) not in (1, dim): + raise RuntimeError( + f'--domain-counts must have one value or {dim} values; use --domains for chunks') + elif opts.domains is not None: + counts = np.array([opts.domains], dtype=int) + else: + counts = np.array([mpi_size], dtype=int) + + if np.any(counts < 1): + raise RuntimeError('domain counts must be positive') + + return counts + + +def grid_counts_from_options(opts, dim): + '''Return point counts for the regular grid shape.''' + if opts.grid_size is not None: + counts = np.array(opts.grid_size, dtype=int) + if len(counts) != dim: + raise RuntimeError(f'--grid-size must have {dim} values') + else: + base = max(1, int(round(opts.point_count**(1.0 / dim)))) + counts = np.full(dim, base, dtype=int) + while int(np.prod(counts)) < opts.point_count: + counts[np.argmin(counts)] += 1 + + if np.any(counts < 1): + raise RuntimeError('grid sizes must be positive') + + return counts + + +def split_range(total, parts, part): + '''Return the half-open item range for one chunk in an even partition.''' + base = total // parts + rem = total % parts + begin = part * base + min(part, rem) + end = begin + base + (1 if part < rem else 0) + return begin, end + + +def spatially_sort_points(points): + '''Return points ordered so contiguous chunks have reasonable spatial coherence.''' + if points.shape[0] == 0: + return points + + spans = np.ptp(points, axis=0) + primary = int(np.argmax(spans)) + secondary = [d for d in range(points.shape[1]) if d != primary] + keys = [points[:, d] for d in reversed(secondary)] + keys.append(points[:, primary]) + return points[np.lexsort(tuple(keys))] + + +def split_points(points, domain_counts): + '''Split points into nearly equal contiguous chunks.''' + points = spatially_sort_points(points) + point_count = points.shape[0] + domain_count = int(np.prod(domain_counts)) + if point_count < domain_count: + raise RuntimeError( + f'point count ({point_count}) must be at least the domain count ({domain_count})') + + chunks = [] + domain_indices = list(itertools.product(*[range(c) for c in domain_counts])) + for domain_id, domain_index in enumerate(domain_indices): + begin, end = split_range(point_count, domain_count, domain_id) + chunks.append((domain_id, domain_index, points[begin:end], begin)) + return chunks + + +def domain_name(domain_index): + '''Create a stable map key for a domain.''' + if len(domain_index) == 1: + return f'domain_{domain_index[0]:06d}' + return 'domain_' + '_'.join(f'{idx:03d}' for idx in domain_index) + + +def generate_circle(point_count, center, radius): + theta = np.linspace(0.0, 2.0 * math.pi, point_count, endpoint=False) + points = np.empty((point_count, 2), dtype=np.float64) + points[:, 0] = center[0] + radius * np.cos(theta) + points[:, 1] = center[1] + radius * np.sin(theta) + return points + + +def generate_sphere(point_count, center, radius): + # Fibonacci sphere points give a deterministic, nearly uniform surface sampling. + idx = np.arange(point_count, dtype=np.float64) + 0.5 + golden_angle = math.pi * (3.0 - math.sqrt(5.0)) + z = 1.0 - 2.0 * idx / point_count + rxy = np.sqrt(np.maximum(0.0, 1.0 - z * z)) + theta = golden_angle * idx + + points = np.empty((point_count, 3), dtype=np.float64) + points[:, 0] = center[0] + radius * rxy * np.cos(theta) + points[:, 1] = center[1] + radius * rxy * np.sin(theta) + points[:, 2] = center[2] + radius * z + return points + + +def generate_torus(point_count, center, major_radius, minor_radius): + # Use a deterministic irrational sequence to avoid requiring two surface counts. + idx = np.arange(point_count, dtype=np.float64) + golden_ratio_conj = (math.sqrt(5.0) - 1.0) / 2.0 + theta = 2.0 * math.pi * idx / point_count + phi = 2.0 * math.pi * np.mod(idx * golden_ratio_conj, 1.0) + ring_radius = major_radius + minor_radius * np.cos(phi) + + points = np.empty((point_count, 3), dtype=np.float64) + points[:, 0] = center[0] + ring_radius * np.cos(theta) + points[:, 1] = center[1] + ring_radius * np.sin(theta) + points[:, 2] = center[2] + minor_radius * np.sin(phi) + return points + + +def generate_gaussian(point_count, center, stddev, seed): + rng = np.random.default_rng(seed) + return rng.normal(loc=center, scale=stddev, size=(point_count, len(center))) + + +def generate_uniform(point_count, lower, upper, seed): + rng = np.random.default_rng(seed) + return rng.uniform(low=lower, high=upper, size=(point_count, len(lower))) + + +def generate_grid_points(grid_counts, lower, upper): + axes = [np.linspace(lower[d], upper[d], grid_counts[d]) for d in range(len(grid_counts))] + coords = np.meshgrid(*axes, indexing='ij') + return np.column_stack([coord.ravel() for coord in coords]) + + +def generate_grid_domains(grid_counts, lower, upper, domain_counts): + '''Generate grid points one spatial domain at a time.''' + for d in range(len(grid_counts)): + if domain_counts[d] > grid_counts[d]: + raise RuntimeError( + f'domain count {domain_counts[d]} exceeds grid size {grid_counts[d]} ' + f'in {AXES[d]} direction') + + axes = [np.linspace(lower[d], upper[d], grid_counts[d]) for d in range(len(grid_counts))] + chunks = [] + global_start = 0 + for domain_id, domain_index in enumerate(itertools.product(*[range(c) for c in domain_counts])): + local_axes = [] + for d, idx in enumerate(domain_index): + begin, end = split_range(grid_counts[d], domain_counts[d], idx) + local_axes.append(axes[d][begin:end]) + coords = np.meshgrid(*local_axes, indexing='ij') + points = np.column_stack([coord.ravel() for coord in coords]) + chunks.append((domain_id, domain_index, points, global_start)) + global_start += points.shape[0] + return chunks + + +def generated_domains(opts, dim, lower, upper, center, stddev, domain_counts): + '''Return generated domains as tuples: id, index, points, global_id_start.''' + if opts.shape == 'grid': + grid_counts = grid_counts_from_options(opts, dim) + if opts.verbose: + print(f'grid_size={grid_counts.tolist()} point_count={int(np.prod(grid_counts))}') + if len(domain_counts) == dim: + return generate_grid_domains(grid_counts, lower, upper, domain_counts) + + points = generate_grid_points(grid_counts, lower, upper) + return split_points(points, domain_counts) + + if opts.shape == 'circle': + points = generate_circle(opts.point_count, center, opts.radius) + elif opts.shape == 'sphere': + points = generate_sphere(opts.point_count, center, opts.radius) + elif opts.shape == 'torus': + points = generate_torus(opts.point_count, center, opts.major_radius, opts.minor_radius) + elif opts.shape == 'gaussian': + points = generate_gaussian(opts.point_count, center, stddev, opts.seed) + elif opts.shape == 'uniform': + points = generate_uniform(opts.point_count, lower, upper, opts.seed) + else: + raise RuntimeError(f'Unknown shape: {opts.shape}') + + return split_points(points, domain_counts) + + +def fill_domain(dom, opts, domain_id, points, global_id_start): + '''Fill one blueprint domain with point mesh data.''' + point_count = points.shape[0] + dim = points.shape[1] + + dom['state/domain_id'] = int(domain_id) + dom['coordsets/' + opts.coordset_name + '/type'] = 'explicit' + for d in range(dim): + values_path = f'coordsets/{opts.coordset_name}/values/{AXES[d]}' + dom[values_path].set(np.ascontiguousarray(points[:, d], dtype=np.float64)) + + topo = dom['topologies/' + opts.topology_name] + topo['type'] = opts.topology_type + topo['coordset'] = opts.coordset_name + if opts.topology_type == 'unstructured': + topo['elements/shape'] = 'point' + topo['elements/connectivity'].set(np.arange(point_count, dtype=np.int32)) + topo['elements/sizes'].set(np.ones(point_count, dtype=np.int32)) + topo['elements/offsets'].set(np.arange(point_count, dtype=np.int32)) + + if opts.id_field: + field = dom['fields/global_id'] + field['association'] = 'vertex' + field['topology'] = opts.topology_name + field['values'].set( + np.arange(global_id_start, global_id_start + point_count, dtype=np.int64)) + + +def create_domain(md_mesh, opts, domain_id, domain_index, points, global_id_start): + '''Append one domain to the multidomain mesh.''' + if opts.use_list: + dom = md_mesh.append() + else: + dom = md_mesh[domain_name(domain_index)] + + fill_domain(dom, opts, domain_id, points, global_id_start) + + +def local_chunks_for_rank(chunks, rank, size, single_domain): + '''Return the chunks that should be written by this rank.''' + if size <= 1: + return chunks + if single_domain: + return chunks if rank == 0 else [] + + begin, end = split_range(len(chunks), size, rank) + return chunks[begin:end] + + +def verify_mesh(mesh, mpi): + '''Verify a mesh with the serial or MPI Blueprint verifier.''' + info = conduit.Node() + if mpi['enabled'] and mpi['size'] > 1: + valid = conduit.blueprint.mpi.mesh.verify(mesh, info, mpi['comm']) + else: + valid = conduit.blueprint.mesh.verify(mesh, info) + + return valid, info + + +def save_mesh(mesh, opts, mpi): + '''Save a mesh with the serial or MPI Blueprint relay.''' + if mpi['enabled'] and mpi['size'] > 1: + try: + conduit.relay.mpi.io.blueprint.save_mesh(mesh, opts.output, opts.protocol, mpi['comm']) + except TypeError: + conduit.relay.mpi.io.blueprint.save_mesh(mesh, opts.output, mpi['comm'], opts.protocol) + else: + conduit.relay.io.blueprint.save_mesh(mesh, opts.output, opts.protocol) + + +def main(): + mpi = get_mpi_context() + opts = parse_args() + dim = infer_dimension(opts) + lower = vector_option(opts.ml, dim, -1.0, '-ml') + upper = vector_option(opts.mu, dim, 1.0, '-mu') + center = vector_option(opts.center, dim, 0.0, '--center') + stddev = vector_option(opts.stddev, dim, 1.0, '--stddev') + domain_counts = domain_counts_from_options(opts, dim, mpi['size']) + + if np.any(upper <= lower): + raise RuntimeError('-mu values must be greater than -ml values') + if opts.radius <= 0.0: + raise RuntimeError('--radius must be positive') + if opts.major_radius <= 0.0 or opts.minor_radius <= 0.0: + raise RuntimeError('--major-radius and --minor-radius must be positive') + if np.any(stddev <= 0.0): + raise RuntimeError('--stddev values must be positive') + + chunks = generated_domains(opts, dim, lower, upper, center, stddev, domain_counts) + local_chunks = local_chunks_for_rank(chunks, mpi['rank'], mpi['size'], opts.single_domain) + + mesh = conduit.Node() + if opts.single_domain: + if local_chunks: + domain_id, _, points, global_id_start = local_chunks[0] + fill_domain(mesh, opts, domain_id, points, global_id_start) + else: + for domain_id, domain_index, points, global_id_start in local_chunks: + create_domain(mesh, opts, domain_id, domain_index, points, global_id_start) + + if opts.verbose: + print(f'rank {mpi["rank"]} mesh:') + print(mesh) + + valid, info = verify_mesh(mesh, mpi) + if not valid: + print('Mesh failed blueprint verification. Info:') + print(info) + return 2 + + save_mesh(mesh, opts, mpi) + if mpi['rank'] == 0: + total_points = sum(points.shape[0] for _, _, points, _ in chunks) + mesh_kind = 'single-domain' if opts.single_domain else 'multidomain' + print( + f'Wrote {total_points} {dim}D {opts.shape} points as a {mesh_kind} mesh ' + f'with {len(chunks)} domain(s) to {opts.output} using {opts.protocol}') + + return 0 + + +if __name__ == '__main__': + sys.exit(main()) From 8bc5d405a0aef64527e25d8a6b7b4ab923fe6ad9 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Thu, 13 Aug 2026 20:25:27 -0700 Subject: [PATCH 06/16] Quest: Improves python script to generate DCP inputs ... and remove analytical generation from DCP example. --- src/axom/quest/examples/CMakeLists.txt | 150 +++-- ...est_distributed_distance_query_example.cpp | 456 ++------------- src/tools/gen-multidom-point-mesh.py | 517 +++++++++++++----- 3 files changed, 511 insertions(+), 612 deletions(-) diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt index 5261e39afd..f5760c3c7d 100644 --- a/src/axom/quest/examples/CMakeLists.txt +++ b/src/axom/quest/examples/CMakeLists.txt @@ -378,83 +378,141 @@ if(AXOM_ENABLE_MPI AND AXOM_ENABLE_SIDRE AND HDF5_FOUND) FOLDER axom/quest/examples ) - if(AXOM_ENABLE_TESTS AND AXOM_DATA_DIR) + if(AXOM_ENABLE_TESTS AND AXOM_ENABLE_PYTHON_TESTS) set(_nranks 3) + set(_dcp_mesh_generator ${PROJECT_SOURCE_DIR}/tools/gen-multidom-point-mesh.py) + + set(_dcp_object_mesh ${CMAKE_CURRENT_BINARY_DIR}/dcp_object_circle_empty_domains) + set(_dcp_gen_object_test quest_distributed_closest_point_gen_object_circle_empty_domains) + axom_add_python_test( + NAME ${_dcp_gen_object_test} + COMMAND ${Python_EXECUTABLE} + ${_dcp_mesh_generator} + circle + --point-count 4 + --center 0.7,0.9 + --radius 0.9 + --domains 6 + --output ${_dcp_object_mesh} + NUM_MPI_TASKS ${_nranks}) - # Run the distributed closest point example on N ranks for each enabled policy - # Non-zero empty-rank probability tests domain underloading case - set(_meshes "mdmesh.2x1" "mdmesh.2x3" "mdmesh.2x2x1" "mdmeshg.2x2x1") - # The mdmesh.* files were generated by these commands: - # src/tools/gen-multidom-structured-mesh.py -ml=0,0 -mu=2,2 -ms=100,100 -dc=2,1 -o mdmesh.2x1 - # src/tools/gen-multidom-structured-mesh.py -ml=0,0 -mu=2,2 -ms=100,100 -dc=2,3 -o mdmesh.2x3 - # src/tools/gen-multidom-structured-mesh.py -ml=0,0,0 -mu=2,2,2 -ms=20,20,15 -dc=2,2,1 -o mdmeshg.2x2x1 --strided - foreach(_pol ${AXOM_EXECUTION_POLICIES}) - # sets the number of threads for the omp policy; leave empty for non-omp policies - set(_num_threads) - if(_pol STREQUAL "omp") - set(_num_threads ${AXOM_TEST_NUM_OMP_THREADS}) + set(_dcp_object_mesh_3d ${CMAKE_CURRENT_BINARY_DIR}/dcp_object_sphere) + set(_dcp_gen_object_test_3d quest_distributed_closest_point_gen_object_sphere) + axom_add_python_test( + NAME ${_dcp_gen_object_test_3d} + COMMAND ${Python_EXECUTABLE} + ${_dcp_mesh_generator} + sphere + --long-point-count 12 + --lat-point-count 6 + --center 0.7,0.9,0.5 + --radius 0.9 + --domains 6 + --output ${_dcp_object_mesh_3d} + NUM_MPI_TASKS ${_nranks}) + + # Run the distributed closest point example on N ranks for each enabled policy. + # Query meshes are generated here instead of checked into the data repository. + set(_query_meshes structured_quads unstructured_quads structured_hexes unstructured_hexes) + foreach(_mesh ${_query_meshes}) + set(_query_mesh ${CMAKE_CURRENT_BINARY_DIR}/dcp_query_${_mesh}) + set(_gen_query_test quest_distributed_closest_point_gen_query_${_mesh}) + if(_mesh MATCHES "quads$") + set(_dim 2) + set(_grid_size 30,30) + set(_min 0,0) + set(_max 2,2) + set(_domain_counts 2,1) + set(_object_mesh ${_dcp_object_mesh}) + set(_object_test ${_dcp_gen_object_test}) + else() + set(_dim 3) + set(_grid_size 8,8,6) + set(_min 0,0,0) + set(_max 2,2,2) + set(_domain_counts 2,2,1) + set(_object_mesh ${_dcp_object_mesh_3d}) + set(_object_test ${_dcp_gen_object_test_3d}) + endif() + if(_mesh MATCHES "^structured") + set(_grid_topology structured) + else() + set(_grid_topology unstructured) endif() - foreach(_mesh ${_meshes}) - # Determine problem dimension by mesh filename. - # and set dimension-dependent arguments. - string(REGEX MATCH "\\.[0-9]+(x[0-9]+)+$" _sizes "${_mesh}") - string(REGEX MATCHALL "[0-9]+" _sizes ${_sizes}) - list(LENGTH _sizes _ndim) + axom_add_python_test( + NAME ${_gen_query_test} + COMMAND ${Python_EXECUTABLE} + ${_dcp_mesh_generator} + grid + --grid-size ${_grid_size} + --min ${_min} + --max ${_max} + --domain-counts ${_domain_counts} + --grid-topology ${_grid_topology} + --output ${_query_mesh} + NUM_MPI_TASKS ${_nranks}) - if(_ndim EQUAL 2) - set(_center 0.7 0.9) - elseif(_ndim EQUAL 3) - set(_center 0.7 0.9 0.5) + foreach(_pol ${AXOM_EXECUTION_POLICIES}) + # sets the number of threads for the omp policy; leave empty for non-omp policies + set(_num_threads) + if(_pol STREQUAL "omp") + set(_num_threads ${AXOM_TEST_NUM_OMP_THREADS}) endif() - set(_test "quest_distributed_closest_point_run_${_ndim}D_${_pol}_${_mesh}") + set(_test "quest_distributed_closest_point_run_${_dim}D_${_pol}_${_mesh}") axom_add_test( NAME ${_test} COMMAND quest_distributed_distance_query_ex - --mesh-file ${quest_data_dir}/${_mesh}.root - --long-point-count 60 - --center ${_center} - --radius 0.9 - --lat-point-count 30 - --obj-domain-count-range 0 2 + --mesh-file ${_query_mesh}.root + --object-mesh-file ${_object_mesh}.root --dist-threshold .3 - --no-random-spacing - --check-results --dynamic-distance-filtering --policy ${_pol} - --object-file dcp_object_mesh_${_ndim}d_${_pol}_${_mesh} - --distance-file dcp_closest_point_2d_${_pol}_${_mesh} + --distance-file dcp_closest_point_${_dim}d_${_pol}_${_mesh} NUM_MPI_TASKS ${_nranks} NUM_OMP_THREADS ${_num_threads}) + set_tests_properties(${_test} + PROPERTIES DEPENDS "${_object_test};${_gen_query_test}") - if(_pol STREQUAL "seq" AND (_mesh STREQUAL "mdmesh.2x1" OR _mesh STREQUAL "mdmesh.2x2x1")) + if(_pol STREQUAL "seq" AND + (_mesh STREQUAL "structured_quads" OR _mesh STREQUAL "structured_hexes")) # Keep fallback coverage for the non-dynamic filtering path # with one 2D and one 3D seq case. set(_static_test "${_test}_no_dynamic_distance_filtering") axom_add_test( NAME ${_static_test} COMMAND quest_distributed_distance_query_ex - --mesh-file ${quest_data_dir}/${_mesh}.root - --long-point-count 60 - --center ${_center} - --radius 0.9 - --lat-point-count 30 - --obj-domain-count-range 0 2 + --mesh-file ${_query_mesh}.root + --object-mesh-file ${_object_mesh}.root --dist-threshold .3 - --no-random-spacing - --check-results --no-dynamic-distance-filtering --policy ${_pol} - --object-file dcp_object_mesh_${_ndim}d_${_pol}_${_mesh}_no_dynamic - --distance-file dcp_closest_point_2d_${_pol}_${_mesh}_no_dynamic + --distance-file dcp_closest_point_${_dim}d_${_pol}_${_mesh}_no_dynamic NUM_MPI_TASKS ${_nranks} NUM_OMP_THREADS ${_num_threads}) + set_tests_properties(${_static_test} + PROPERTIES DEPENDS "${_object_test};${_gen_query_test}") endif() endforeach() endforeach() - unset(optional_dependency) + unset(_dcp_gen_object_test) + unset(_dcp_gen_object_test_3d) + unset(_dcp_mesh_generator) + unset(_dcp_object_mesh) + unset(_dcp_object_mesh_3d) + unset(_dim) + unset(_domain_counts) + unset(_gen_query_test) + unset(_grid_size) + unset(_grid_topology) + unset(_max) + unset(_min) + unset(_object_mesh) + unset(_object_test) + unset(_query_mesh) + unset(_query_meshes) unset(_nranks) unset(_test) unset(_static_test) diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index a0476f28ec..0b9bbcd316 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -37,7 +37,6 @@ // C/C++ includes #include -#include #include #include #include @@ -54,10 +53,7 @@ namespace quest = axom::quest; namespace slic = axom::slic; namespace sidre = axom::sidre; namespace slam = axom::slam; -namespace spin = axom::spin; namespace primal = axom::primal; -namespace mint = axom::mint; -namespace numerics = axom::numerics; using RuntimePolicy = axom::runtime_policy::Policy; @@ -74,8 +70,6 @@ struct Input public: std::string meshFile; std::string distanceFile {"cp_coords"}; - std::string objectFile {"object_mesh"}; - std::string objectMeshFile; // Memory instrumentation (all off by default). @@ -83,26 +77,12 @@ struct Input bool trimAfterQuery {false}; // malloc_trim(0) after the query, then re-report int sampleMemoryMs {0}; // if > 0, background-sample peak RSS during the query - double circleRadius {1.0}; - std::vector circleCenter {0.0, 0.0}; - int longPointCount {100}; - - // Latitudinal direction of object mesh - std::vector latRange {-90.0, 90.0}; - int latPointCount {20}; - RuntimePolicy policy {RuntimePolicy::seq}; double distThreshold {axom::numeric_limits::max()}; bool dynamicDistanceFiltering {true}; - bool checkResults {false}; - - bool randomSpacing {true}; - - std::vector objDomainCountRange {1, 1}; - private: bool m_verboseOutput {false}; @@ -143,15 +123,12 @@ struct Input ->description("Name of output mesh file containing closest distance.") ->capture_default_str(); - app.add_option("-o,--object-file", objectFile) - ->description("Name of output file containing object mesh.") - ->capture_default_str(); - app.add_option("--object-mesh-file", objectMeshFile) ->description( "Path to a conduit blueprint point mesh root file. " - " When provided, we use this instead of an analytically generated object mesh.") - ->check(axom::CLI::ExistingFile); + "Generate this mesh with src/tools/gen-multidom-point-mesh.py.") + ->check(axom::CLI::ExistingFile) + ->required(); app.add_flag("--track-memory", trackMemory) ->description( @@ -177,34 +154,6 @@ struct Input ->description("Enable/disable verbose output") ->capture_default_str(); - app.add_option("-r,--radius", circleRadius)->description("Radius for sphere")->capture_default_str(); - - auto* object_options = - app.add_option_group("sphere", "Options for setting up object points on the sphere"); - object_options->add_option("--center", circleCenter) - ->description("Center for object (x,y[,z])") - ->expected(2, 3); - - object_options->add_option("--obj-domain-count-range", objDomainCountRange) - ->description("Range of object domain counts per rank (min, max)") - ->expected(2); - - object_options->add_flag("--random-spacing,!--no-random-spacing", randomSpacing) - ->description("Enable/disable random spacing of object points") - ->capture_default_str(); - - object_options->add_option("-n,--long-point-count", longPointCount) - ->description("Number of points around the longitudinal direction") - ->capture_default_str(); - - object_options->add_option("--lat-range", latRange) - ->description("Latitude range in degrees from the equator (3D only)") - ->expected(2); - - object_options->add_option("--lat-point-count", latPointCount) - ->description("Number of points in the latitudinal direction (3D only)") - ->capture_default_str(); - app.add_option("-d,--dist-threshold", distThreshold) ->check(axom::CLI::NonNegativeNumber) ->description("Distance threshold to search") @@ -221,10 +170,6 @@ struct Input ->capture_default_str() ->transform(axom::CLI::CheckedTransformer(axom::runtime_policy::s_nameToPolicy)); - app.add_flag("-c,--check-results,!--no-check-results", checkResults) - ->description("Enable/disable checking results against analytical solution") - ->capture_default_str(); - app.get_formatter()->column_width(60); // could throw an exception @@ -577,6 +522,20 @@ struct BlueprintParticleMesh fld->copyView(strides); } + if(SZ == 0) + { + fld->createViewAndAllocate("values/x", sidre::detail::SidreTT::id, 0); + if(DIM > 1) + { + fld->createViewAndAllocate("values/y", sidre::detail::SidreTT::id, 0); + } + if(DIM > 2) + { + fld->createViewAndAllocate("values/z", sidre::detail::SidreTT::id, 0); + } + continue; + } + // create views into a shared buffer for the coordinates, with stride DIM auto* buf = domain_group(dIdx) ->getDataStore() @@ -766,11 +725,6 @@ struct BlueprintParticleMesh class ObjectMeshWrapper { public: - ObjectMeshWrapper(sidre::Group* group) : m_objectMesh(group, "mesh", "coords") - { - SLIC_ASSERT(group != nullptr); - } - //!@brief Construct by reading a blueprint object mesh from disk. //! Uses the topology/coordset names found in the file, exactly like the //! query mesh reader. (The empty-topology BlueprintParticleMesh ctor lets @@ -789,19 +743,8 @@ class ObjectMeshWrapper std::string getTopologyName() const { return m_objectMesh.getTopologyName(); } std::string getCoordsetName() const { return m_objectMesh.getCoordsetName(); } - void setVerbosity(bool verbose) { m_verbose = verbose; } - - /// Outputs the object mesh to disk - void saveMesh(const std::string& filename = "object_mesh") - { - SLIC_INFO(banner(axom::fmt::format("Saving object mesh '{}' to disk", filename))); - - m_objectMesh.saveMesh(filename); - } - private: BlueprintParticleMesh m_objectMesh; - bool m_verbose {false}; }; class QueryMeshWrapper @@ -910,281 +853,30 @@ class QueryMeshWrapper int dim = srcNode.fetch_existing("values").number_of_children(); for(int d = 0; d < dim; ++d) { - conduit::float64_array dst = dstGroup->getGroup("values")->getView(d)->getArray(); - const conduit::float64_array src = srcNode.fetch_existing("values").child(d).value(); - SLIC_ASSERT(src.number_of_elements() == dst.number_of_elements()); - int nPts = src.number_of_elements(); - for(int i = 0; i < nPts; ++i) - { - dst[i] = src[i]; - } - } - } - } - } - - /** - * Check for error in the search. - * - check that points within threshold have a closest point - * on the object. - * - check that found closest-point is near its corresponding - * closest point on the sphere (within tolerance) - * - * Return number of errors found on the local mesh partition. - * Populate "error_flag" field with the number of errors, for - * visualization. - * - * Randomizing points (--random-spacing switch) can cause false - * positives, so when it's on, distance inaccuracy is a warning (not - * an error) for the purpose of checking. - */ - template - int checkClosestPoints(const axom::primal::Sphere& sphere, const Input& params) - { - using PointType = axom::primal::Point; - - m_queryMesh.registerNodalScalarField("error_flag"); - - int sumErrCount = 0; - int sumWarningCount = 0; - for(axom::IndexType dIdx = 0; dIdx < m_queryMesh.domain_count(); ++dIdx) - { - auto queryPts = m_queryMesh.getPoints(dIdx); - - axom::ArrayView cpCoords = - m_queryMesh.getNodalVectorField("cp_coords", dIdx); - SLIC_INFO(axom::fmt::format("Closest points ({}):", cpCoords.size())); - - axom::ArrayView cpIndices = - m_queryMesh.getNodalScalarField("cp_index", dIdx); - - axom::ArrayView errorFlag = - m_queryMesh.getNodalScalarField("error_flag", dIdx); - - SLIC_ASSERT(queryPts.size() == cpCoords.size()); - SLIC_ASSERT(queryPts.size() == cpIndices.size()); - - if(params.isVerbose()) - { - SLIC_INFO(axom::fmt::format("Closest points ({}):", cpCoords.size())); - } - - /* - Allowable slack is half the arclength between 2 adjacent object - points. A query point on the object can correctly have that - closest-distance, even though the analytical distance is zero. - If spacing is random, distance between adjacent points is not - predictable, leading to false positives. We don't claim errors - for this in when using random. - */ - const double longSpacing = 2 * M_PI * params.circleRadius / params.longPointCount; - const double latSpacing = params.circleRadius * M_PI / 180 * - (params.latRange[1] - params.latRange[0]) / params.latPointCount; - const double avgObjectRes = - DIM == 2 ? longSpacing : std::sqrt(longSpacing * longSpacing + latSpacing * latSpacing); - const double allowableSlack = avgObjectRes / 2; - - using IndexSet = slam::PositionSet<>; - for(auto i : IndexSet(queryPts.size())) - { - bool errf = false; - - const auto& qPt = queryPts[i]; - const auto& cpCoord = cpCoords[i]; - double analyticalDist = std::fabs(sphere.computeSignedDistance(qPt)); - const bool closestPointFound = (cpIndices[i] == -1); - if(closestPointFound) - { - if(analyticalDist < params.distThreshold - allowableSlack) - { - errf = true; - SLIC_INFO( - axom::fmt::format("***Error: Query point {} ({}) is within " - "threshold by {} but lacks closest point.", - i, - qPt, - params.distThreshold - analyticalDist)); - } - } - else - { - if(analyticalDist >= params.distThreshold + allowableSlack) + auto* dstView = dstGroup->getGroup("values")->getView(d); + const auto& srcComponent = srcNode.fetch_existing("values").child(d); + int nPts = srcComponent.dtype().number_of_elements(); + SLIC_ASSERT(nPts == dstView->getNumElements()); + if(nPts == 0) { - errf = true; - SLIC_INFO( - axom::fmt::format("***Error: Query point {} ({}) is outside " - "threshold by {} but has closest point at {}.", - i, - qPt, - analyticalDist - params.distThreshold, - cpCoord)); + continue; } - if(!axom::utilities::isNearlyEqual(sphere.computeSignedDistance(cpCoord), 0.0)) - { - errf = true; - SLIC_INFO( - axom::fmt::format("***Error: Closest point ({}) for index {} " - "({}) is not on the sphere.", - cpCoords[i], - i, - qPt)); - } - - double dist = sqrt(primal::squared_distance(qPt, cpCoord)); - if(!axom::utilities::isNearlyEqual(dist, analyticalDist, allowableSlack)) + conduit::float64_array dst = dstView->getArray(); + const conduit::float64_array src = srcComponent.value(); + for(int i = 0; i < nPts; ++i) { - if(params.randomSpacing) - { - ++sumWarningCount; - SLIC_INFO( - axom::fmt::format("***Warning: Closest distance for {} (index " - "{}, cp {}) is {}, off by {}.", - qPt, - i, - cpCoords[i], - dist, - dist - analyticalDist)); - } - else - { - errf = true; - SLIC_INFO( - axom::fmt::format("***Warning: Closest distance for {} (index " - "{}, cp {}) is {}, off by {}.", - qPt, - i, - cpCoords[i], - dist, - dist - analyticalDist)); - } + dst[i] = src[i]; } } - errorFlag[i] = errf; - sumErrCount += errf; } } - - SLIC_INFO( - axom::fmt::format("Local partition has {} errors, {} warnings in closest distance results.", - sumErrCount, - sumWarningCount)); - - return sumErrCount; } private: BlueprintParticleMesh m_queryMesh; }; -/** - * Generates points on a sphere, partitioned into multiple domains. - * Point spacing in the longitudinal direction can be random (default) or uniform. - * 3D points cover the given latitude range. - */ -void generateObjectPoints(BlueprintParticleMesh& particleMesh, - int spatialDimension, - const std::vector& center, - double radius, - int longPointCount, - int localDomainCount, - bool randomSpacing = true, - double minLatitude = 0.0, - double maxLatitude = 0.0, - int latPointCount = 1) -{ - using axom::utilities::random_real; - - int rank = particleMesh.getRank(); - int nranks = particleMesh.getNumRanks(); - - // rank scan to sum longPointCount and determine local range of longitudinal angles. - axom::Array sums(nranks, nranks); - { - axom::Array indivDomainCounts(nranks, nranks); - indivDomainCounts.fill(-1); - MPI_Allgather(&localDomainCount, 1, MPI_INT, indivDomainCounts.data(), 1, MPI_INT, MPI_COMM_WORLD); - - SLIC_DEBUG_IF( - params.isVerbose(), - axom::fmt::format("After all gather: [{}]", axom::fmt::join(indivDomainCounts, ","))); - - sums[0] = indivDomainCounts[0]; - for(int i = 1; i < nranks; ++i) - { - sums[i] = sums[i - 1] + indivDomainCounts[i]; - } - // If no rank has any domains, force last one to have 1 domain. - if(sums[nranks - 1] == 0) - { - sums[nranks - 1] = 1; - if(rank == nranks - 1) - { - localDomainCount = 1; - } - } - } - - SLIC_DEBUG_IF(params.isVerbose(), - axom::fmt::format("After scan: [{}]", axom::fmt::join(sums, ","))); - - int globalDomainCount = sums[nranks - 1]; - longPointCount = std::max(longPointCount, globalDomainCount); - int longPtsPerDomain = longPointCount / globalDomainCount; - int domainsWithExtraPt = longPointCount % globalDomainCount; - - int myDomainBegin = rank == 0 ? 0 : sums[rank - 1]; - int myDomainEnd = sums[rank]; - SLIC_ASSERT(myDomainEnd - myDomainBegin == localDomainCount); - - if(spatialDimension == 2) - { - minLatitude = 0.0; - maxLatitude = 0.0; - latPointCount = 1; - } - minLatitude *= M_PI / 180; - maxLatitude *= M_PI / 180; - const double longSpacing = 2. * M_PI / longPointCount; - const double latSpacing = - latPointCount < 2 || latPointCount == 1 ? 0 : (maxLatitude - minLatitude) / (latPointCount - 1); - - for(int di = myDomainBegin; di < myDomainEnd; ++di) - { - int pBegin = di * longPtsPerDomain + std::min(di, domainsWithExtraPt); - int pEnd = (di + 1) * longPtsPerDomain + std::min((di + 1), domainsWithExtraPt); - int domainPointCount = pEnd - pBegin; - domainPointCount *= latPointCount; - axom::Array pts(domainPointCount, spatialDimension); - axom::IndexType iPts = 0; - - for(int li = 0; li < latPointCount; ++li) - { - double latAngle = minLatitude + li * latSpacing; - double xyRadius = radius * std::cos(latAngle); // Project radius onto x-y plane. - double z = spatialDimension == 2 ? 0 : center[2] + radius * std::sin(latAngle); - for(int pi = pBegin; pi < pEnd; ++pi) - { - const double ang = - randomSpacing ? random_real(longSpacing * pBegin, longSpacing * pEnd) : pi * longSpacing; - const double rsinT = center[1] + xyRadius * std::sin(ang); - const double rcosT = center[0] + xyRadius * std::cos(ang); - pts[iPts][0] = rcosT; - pts[iPts][1] = rsinT; - if(spatialDimension > 2) - { - pts[iPts][2] = z; - } - ++iPts; - } - } - particleMesh.setPoints(di, pts); - } - - axom::slic::flushStreams(); - SLIC_ASSERT(particleMesh.isValid()); -} - //--------------------------------------------------------------------------- // Transform closest points to distances and directions //--------------------------------------------------------------------------- @@ -1578,15 +1270,6 @@ int main(int argc, char** argv) exit(retval); } - // Issue warning about result-checking requiring good resolution. - if(params.checkResults && params.randomSpacing) - { - SLIC_INFO( - axom::fmt::format("***Warning: Result-checking may yield false positive (warnings) when " - "sphere points have random spacing. High resolution helps limit this." - "We recommend at least 500 points for each radius length unit.")); - } - #if defined(AXOM_USE_UMPIRE) //--------------------------------------------------------------------------- // Memory resource. For testing, choose device memory if appropriate. @@ -1629,49 +1312,12 @@ int main(int argc, char** argv) const size_t spatialDim = queryMeshWrapper.getParticleMesh().dimension(); - const bool readObjectFromFile = !params.objectMeshFile.empty(); - - // The analytic circle/sphere generator needs a center whose dimension matches the query mesh. - // only enforce that when we are actually generating. - SLIC_ASSERT(readObjectFromFile || params.circleCenter.size() == spatialDim); - //--------------------------------------------------------------------------- - // Object (second) mesh: read from file, or generate analytically + // Object (second) mesh //--------------------------------------------------------------------------- - std::unique_ptr objectMeshWrapperPtr; - if(readObjectFromFile) - { - objectMeshWrapperPtr = - std::make_unique(dataStore.getRoot()->createGroup("object_mesh", true), - params.objectMeshFile); - } - else - { - objectMeshWrapperPtr = - std::make_unique(dataStore.getRoot()->createGroup("object_mesh", true)); - } - ObjectMeshWrapper& objectMeshWrapper = *objectMeshWrapperPtr; - objectMeshWrapper.setVerbosity(params.isVerbose()); - - if(!readObjectFromFile) - { - SLIC_ASSERT(params.objDomainCountRange[1] >= params.objDomainCountRange[0]); - const unsigned int omin = params.objDomainCountRange[0]; - const unsigned int omax = params.objDomainCountRange[1]; - const double prob = axom::utilities::random_real(0., 1.); - int localDomainCount = omin + int(0.5 + prob * (omax - omin)); - generateObjectPoints(objectMeshWrapper.getParticleMesh(), - spatialDim, - params.circleCenter, - params.circleRadius, - params.longPointCount, - localDomainCount, - params.randomSpacing, - params.latRange.size() > 0 ? params.latRange[0] : 0.0, - params.latRange.size() > 1 ? params.latRange[1] : 0.0, - params.latPointCount); - } + ObjectMeshWrapper objectMeshWrapper(dataStore.getRoot()->createGroup("object_mesh", true), + params.objectMeshFile); if(params.isVerbose()) { @@ -1679,19 +1325,15 @@ int main(int argc, char** argv) } slic::flushStreams(); - // Only re-save the generated object; a file-read object is already on disk. - if(!readObjectFromFile) - { - objectMeshWrapper.saveMesh(params.objectFile); - } - slic::flushStreams(); - //--------------------------------------------------------------------------- // Initialize spatial index for object points, and run query //--------------------------------------------------------------------------- + int globalObjectPointCount = objectMeshWrapper.getParticleMesh().numPoints(); + MPI_Allreduce(MPI_IN_PLACE, &globalObjectPointCount, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + auto init_str = - banner(axom::fmt::format("Initializing BVH tree over {} points", params.longPointCount)); + banner(axom::fmt::format("Initializing BVH tree over {} object points", globalObjectPointCount)); axom::utilities::Timer initTimer(false); axom::utilities::Timer queryTimer(false); @@ -1796,25 +1438,6 @@ int main(int argc, char** argv) slic::flushStreams(); queryMeshWrapper.update_closest_points(queryMeshNode); - int errCount = 0; - int localErrCount = 0; - if(params.checkResults) - { - if(spatialDim == 2) - { - primal::Point center(params.circleCenter.data()); - primal::Sphere sphere(center, params.circleRadius); - localErrCount = queryMeshWrapper.checkClosestPoints(sphere, params); - } - else if(spatialDim == 3) - { - primal::Point center(params.circleCenter.data()); - primal::Sphere sphere(center, params.circleRadius); - localErrCount = queryMeshWrapper.checkClosestPoints(sphere, params); - } - } - MPI_Allreduce(&localErrCount, &errCount, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); - if(spatialDim == 2) { computeDistancesAndDirections<2>(queryMeshWrapper.getParticleMesh(), @@ -1837,17 +1460,10 @@ int main(int argc, char** argv) queryMeshWrapper.saveMesh(params.distanceFile); - if(errCount) - { - SLIC_INFO(axom::fmt::format(" Error exit: {} errors found.", errCount)); - } - else - { - SLIC_INFO("Normal exit."); - } + SLIC_INFO("Normal exit."); finalizeLogger(); MPI_Finalize(); - return errCount != 0; + return 0; } diff --git a/src/tools/gen-multidom-point-mesh.py b/src/tools/gen-multidom-point-mesh.py index 7778f2c75d..5c80f3e83d 100755 --- a/src/tools/gen-multidom-point-mesh.py +++ b/src/tools/gen-multidom-point-mesh.py @@ -1,7 +1,7 @@ #!/usr/bin/env python3 # gen-multidom-point-mesh.py -# Write a point mesh following the Conduit mesh blueprint. +# Write a mesh following the Conduit mesh blueprint for distributed closest point (DCP) testing. # # The generated Blueprint hierarchy is: # (single-domain output) @@ -11,16 +11,17 @@ # |-- topologies # | `-- # | |-- coordset == -# | |-- type == "points" or "unstructured" -# | `-- elements (only for "unstructured") -# | |-- shape == "point" -# | |-- connectivity (int32, point_count) -# | |-- sizes (int32, point_count, all 1) -# | `-- offsets (int32, point_count) +# | |-- type == "points", "structured", or "unstructured" +# | `-- elements (present for structured/unstructured) +# | |-- dims/{i,j,[k]} (structured cell dimensions) +# | |-- shape (unstructured "point", "quad", or "hex") +# | |-- connectivity (unstructured int32 connectivity) +# | |-- sizes (unstructured point meshes only) +# | `-- offsets (unstructured point meshes only) # |-- coordsets # | `-- # | |-- type == "explicit" -# | `-- values +# | `-- values (i-fastest node ordering for grids) # | |-- x (float64, point_count) # | |-- y (float64, point_count) # | `-- [z] (float64, point_count, present in 3D) @@ -37,6 +38,7 @@ import itertools import math +import os import sys from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser, ArgumentTypeError @@ -45,13 +47,11 @@ import conduit.blueprint import conduit.relay except ModuleNotFoundError as e: - print( - f'{e}\nMake sure your PYTHONPATH includes /path/to/conduit/install/python-modules\n' - 'Conduit must be configured with python and hdf5.\n' - 'Alternatively, you can use the convenience script\n' - '/path/to/axom_build_dir/bin/run_python_with_axom.sh\n' - 'that includes Conduit in PYTHONPATH.' - ) + print(f'{e}\nMake sure your PYTHONPATH includes /path/to/conduit/install/python-modules\n' + 'Conduit must be configured with python and hdf5.\n' + 'Alternatively, you can use the convenience script\n' + '/path/to/axom_build_dir/bin/run_python_with_axom.sh\n' + 'that includes Conduit in PYTHONPATH.') sys.exit(-1) try: @@ -60,12 +60,22 @@ print(f'{e}\nThis script requires numpy.') sys.exit(-1) - AXES = 'xyz' +LOGICAL_AXES = 'ijk' def get_mpi_context(): '''Return MPI context when mpi4py and Conduit MPI relay are available.''' + mpi_size_env = 1 + for var_name in ('OMPI_COMM_WORLD_SIZE', 'PMI_SIZE', 'PMIX_SIZE', 'SLURM_STEP_NUM_TASKS'): + try: + mpi_size_env = max(mpi_size_env, int(os.environ.get(var_name, '1'))) + except ValueError: + pass + + if mpi_size_env <= 1: + return {'enabled': False, 'comm': None, 'rank': 0, 'size': 1} + try: from mpi4py import MPI import conduit.blueprint.mpi.mesh @@ -120,70 +130,171 @@ def positive_int(s): return value +def add_common_options(parser): + '''Add output and distribution options shared by all mesh generators.''' + dist = parser.add_argument_group('distribution') + dist.add_argument('-dc', + '--domain-counts', + type=i_c, + help='Domain counts in one chunk direction or each coordinate direction') + dist.add_argument('--domains', + type=positive_int, + help='Total domain count for a one-dimensional chunk partition') + dist.add_argument('--single-domain', + action='store_true', + help='Write a single-domain mesh instead of multidomain output') + dist.add_argument('--use-list', + action='store_true', + help='Put multidomain domains in a list instead of a map') + + bp = parser.add_argument_group('blueprint names and output') + bp.add_argument('--topology-name', default='mesh', help='Blueprint topology name') + bp.add_argument('--coordset-name', default='coords', help='Blueprint coordset name') + bp.add_argument('--id-field', + action='store_true', + help='Add a vertex-associated int64 global_id field') + bp.add_argument('--protocol', default='hdf5', help='Conduit relay protocol for save_mesh') + bp.add_argument('-o', '--output', default='mdmesh', help='Output file base name') + bp.add_argument('-v', '--verbose', action='store_true', help='Print additional info') + + +def add_point_topology_options(parser): + parser.add_argument('--topology-type', + choices=('unstructured', 'points'), + default='unstructured', + help='Point topology representation') + + +def add_center_radius_options(parser, fixed_dim): + parser.add_argument('--center', + type=f_c, + help=f'Center coordinates; scalar or {fixed_dim} comma-separated values') + parser.add_argument('--radius', type=float, default=1.0, help='Circle/sphere radius') + + +def add_min_max_options(parser): + parser.add_argument('--min', + dest='lower', + type=f_c, + help='Lower coordinate bounds; scalar or one value per dimension') + parser.add_argument('--max', + dest='upper', + type=f_c, + help='Upper coordinate bounds; scalar or one value per dimension') + + def parse_args(): - parser = ArgumentParser( - description='Write a single-domain or multidomain blueprint point mesh.', - formatter_class=ArgumentDefaultsHelpFormatter) - parser.add_argument( - 'shape', - choices=('circle', 'sphere', 'torus', 'grid', 'gaussian', 'uniform'), - help='Point distribution to generate') - parser.add_argument( - '-d', - '--dimension', - type=int, - choices=(2, 3), - help='Spatial dimension for dimension-flexible shapes') - parser.add_argument( - '-n', - '--point-count', - type=positive_int, - default=1024, - help='Number of points, or target number for grid when --grid-size is omitted') - parser.add_argument( - '--grid-size', - type=i_c, - help='Point counts in each grid direction, e.g. 101,101 or 65,65,65') - parser.add_argument( - '-dc', - '--domain-counts', - type=i_c, - help='Domain counts in one chunk direction or each coordinate direction') - parser.add_argument( - '--domains', - type=positive_int, - help='Total domain count for a one-dimensional chunk partition') - parser.add_argument( - '--single-domain', - action='store_true', - help='Write a single-domain point mesh instead of multidomain output') - parser.add_argument('-ml', type=f_c, help='Lower coordinates for grid and uniform shapes') - parser.add_argument('-mu', type=f_c, help='Upper coordinates for grid and uniform shapes') - parser.add_argument('--center', type=f_c, help='Center for circle, sphere, torus, and gaussian') - parser.add_argument('--radius', type=float, default=1.0, help='Radius for circle or sphere') - parser.add_argument('--major-radius', type=float, default=1.0, help='Major radius for torus') - parser.add_argument('--minor-radius', type=float, default=0.25, help='Minor radius for torus') - parser.add_argument('--stddev', + parser = ArgumentParser(description='Write blueprint meshes for DCP tests.', + formatter_class=ArgumentDefaultsHelpFormatter) + subparsers = parser.add_subparsers(dest='shape', required=True) + + circle = subparsers.add_parser('circle', + formatter_class=ArgumentDefaultsHelpFormatter, + help='2D analytic circle point mesh') + add_common_options(circle) + add_point_topology_options(circle) + add_center_radius_options(circle, 2) + circle.add_argument('-n', + '--point-count', + '--long-point-count', + dest='point_count', + type=positive_int, + default=1024, + help='Number of points around the circle') + circle.add_argument('--random-spacing', + action='store_true', + help='Use random angular spacing instead of uniform spacing') + circle.add_argument('--seed', type=int, default=0, help='Random seed for random spacing') + + sphere = subparsers.add_parser('sphere', + formatter_class=ArgumentDefaultsHelpFormatter, + help='3D analytic sphere point mesh') + add_common_options(sphere) + add_point_topology_options(sphere) + add_center_radius_options(sphere, 3) + sphere.add_argument('-n', + '--point-count', + type=positive_int, + default=1024, + help='Number of Fibonacci sphere points when not using latitude sampling') + sphere.add_argument('--long-point-count', + type=positive_int, + help='Number of longitudinal samples for latitude sampling') + sphere.add_argument('--lat-point-count', + type=positive_int, + help='Number of latitudinal samples for latitude sampling') + sphere.add_argument('--lat-range', type=f_c, - help='Gaussian standard deviation; scalar or one value per dimension') - parser.add_argument('--seed', type=int, default=0, help='Random seed for stochastic shapes') - parser.add_argument('--use-list', + default=(-90.0, 90.0), + help='Latitude range in degrees for latitude sampling') + sphere.add_argument('--random-spacing', action='store_true', - help='Put domains in a list instead of a map') - parser.add_argument('--topology-name', default='mesh', help='Blueprint topology name') - parser.add_argument('--coordset-name', default='coords', help='Blueprint coordset name') - parser.add_argument( - '--topology-type', - choices=('unstructured', 'points'), - default='unstructured', - help='Blueprint point topology representation') - parser.add_argument( - '--id-field', - action='store_true', - help='Add a vertex-associated int64 global_id field') - parser.add_argument('--protocol', default='hdf5', help='Conduit relay protocol for save_mesh') - parser.add_argument('-o', '--output', default='mdpoint', help='Output file base name') - parser.add_argument('-v', '--verbose', action='store_true', help='Print additional info') + help='Use random longitudinal spacing with latitude sampling') + sphere.add_argument('--seed', type=int, default=0, help='Random seed for random spacing') + + torus = subparsers.add_parser('torus', + formatter_class=ArgumentDefaultsHelpFormatter, + help='3D analytic torus point mesh') + add_common_options(torus) + add_point_topology_options(torus) + torus.add_argument('-n', + '--point-count', + type=positive_int, + default=1024, + help='Number of torus surface points') + torus.add_argument('--center', type=f_c, help='Center coordinates; scalar or 3 values') + torus.add_argument('--major-radius', type=float, default=1.0, help='Torus major radius') + torus.add_argument('--minor-radius', type=float, default=0.25, help='Torus minor radius') + + grid = subparsers.add_parser('grid', + formatter_class=ArgumentDefaultsHelpFormatter, + help='Regular grid mesh over --min/--max') + add_common_options(grid) + add_min_max_options(grid) + grid.add_argument('-d', '--dimension', type=int, choices=(2, 3), help='Spatial dimension') + grid.add_argument('-n', + '--point-count', + type=positive_int, + default=1024, + help='Target point count when --grid-size is omitted') + grid.add_argument('--grid-size', + type=i_c, + help='Point counts in each grid direction, e.g. 101,101') + grid.add_argument('--grid-topology', + choices=('points', 'unstructured-points', 'structured', 'unstructured'), + default='unstructured-points', + help='Topology to generate over the regular grid coordinates') + + gaussian = subparsers.add_parser('gaussian', + formatter_class=ArgumentDefaultsHelpFormatter, + help='Random Gaussian point mesh') + add_common_options(gaussian) + add_point_topology_options(gaussian) + gaussian.add_argument('-d', '--dimension', type=int, choices=(2, 3), help='Spatial dimension') + gaussian.add_argument('-n', + '--point-count', + type=positive_int, + default=1024, + help='Number of points') + gaussian.add_argument('--center', type=f_c, help='Distribution center; scalar or per dimension') + gaussian.add_argument('--stddev', + type=f_c, + help='Gaussian standard deviation; scalar or per dimension') + gaussian.add_argument('--seed', type=int, default=0, help='Random seed') + + uniform = subparsers.add_parser('uniform', + formatter_class=ArgumentDefaultsHelpFormatter, + help='Random uniform point mesh over --min/--max') + add_common_options(uniform) + add_point_topology_options(uniform) + add_min_max_options(uniform) + uniform.add_argument('-d', '--dimension', type=int, choices=(2, 3), help='Spatial dimension') + uniform.add_argument('-n', + '--point-count', + type=positive_int, + default=1024, + help='Number of points') + uniform.add_argument('--seed', type=int, default=0, help='Random seed') opts, unkn = parser.parse_known_args() if opts.verbose: @@ -199,14 +310,12 @@ def infer_dimension(opts): '''Infer and validate the spatial dimension.''' fixed_dims = {'circle': 2, 'sphere': 3, 'torus': 3} if opts.shape in fixed_dims: - dim = fixed_dims[opts.shape] - if opts.dimension is not None and opts.dimension != dim: - raise RuntimeError(f'{opts.shape} requires dimension {dim}') - return dim + return fixed_dims[opts.shape] dim = opts.dimension inferred = [] - for values in (opts.grid_size, opts.ml, opts.mu, opts.center, opts.domain_counts): + for name in ('grid_size', 'lower', 'upper', 'center', 'stddev', 'domain_counts'): + values = getattr(opts, name, None) if values is not None: inferred.append(len(values)) @@ -240,7 +349,7 @@ def vector_option(values, dim, default, name): def domain_counts_from_options(opts, dim, mpi_size): '''Return domain counts as an integer vector.''' if opts.single_domain: - return np.array([1], dtype=int) + return np.ones(dim, dtype=int) if opts.domain_counts is not None and opts.domains is not None: raise RuntimeError('Use --domain-counts or --domains, not both') @@ -268,7 +377,7 @@ def grid_counts_from_options(opts, dim): if len(counts) != dim: raise RuntimeError(f'--grid-size must have {dim} values') else: - base = max(1, int(round(opts.point_count**(1.0 / dim)))) + base = max(1, int(round(opts.point_count ** (1.0 / dim)))) counts = np.full(dim, base, dtype=int) while int(np.prod(counts)) < opts.point_count: counts[np.argmin(counts)] += 1 @@ -306,15 +415,18 @@ def split_points(points, domain_counts): points = spatially_sort_points(points) point_count = points.shape[0] domain_count = int(np.prod(domain_counts)) - if point_count < domain_count: - raise RuntimeError( - f'point count ({point_count}) must be at least the domain count ({domain_count})') chunks = [] domain_indices = list(itertools.product(*[range(c) for c in domain_counts])) for domain_id, domain_index in enumerate(domain_indices): begin, end = split_range(point_count, domain_count, domain_id) - chunks.append((domain_id, domain_index, points[begin:end], begin)) + chunks.append({ + 'domain_id': domain_id, + 'domain_index': domain_index, + 'points': points[begin:end], + 'global_id_start': begin, + 'grid_counts': None, + }) return chunks @@ -325,31 +437,65 @@ def domain_name(domain_index): return 'domain_' + '_'.join(f'{idx:03d}' for idx in domain_index) -def generate_circle(point_count, center, radius): - theta = np.linspace(0.0, 2.0 * math.pi, point_count, endpoint=False) +def generate_circle(point_count, center, radius, random_spacing, seed): + if random_spacing: + rng = np.random.default_rng(seed) + theta = np.sort(rng.uniform(0.0, 2.0 * math.pi, point_count)) + else: + theta = np.linspace(0.0, 2.0 * math.pi, point_count, endpoint=False) + points = np.empty((point_count, 2), dtype=np.float64) points[:, 0] = center[0] + radius * np.cos(theta) points[:, 1] = center[1] + radius * np.sin(theta) return points -def generate_sphere(point_count, center, radius): - # Fibonacci sphere points give a deterministic, nearly uniform surface sampling. - idx = np.arange(point_count, dtype=np.float64) + 0.5 - golden_angle = math.pi * (3.0 - math.sqrt(5.0)) - z = 1.0 - 2.0 * idx / point_count - rxy = np.sqrt(np.maximum(0.0, 1.0 - z * z)) - theta = golden_angle * idx +def generate_sphere(opts, center, radius): + if opts.long_point_count is None and opts.lat_point_count is None: + point_count = opts.point_count + idx = np.arange(point_count, dtype=np.float64) + 0.5 + golden_angle = math.pi * (3.0 - math.sqrt(5.0)) + z = 1.0 - 2.0 * idx / point_count + rxy = np.sqrt(np.maximum(0.0, 1.0 - z * z)) + theta = golden_angle * idx + + points = np.empty((point_count, 3), dtype=np.float64) + points[:, 0] = center[0] + radius * rxy * np.cos(theta) + points[:, 1] = center[1] + radius * rxy * np.sin(theta) + points[:, 2] = center[2] + radius * z + return points + + long_count = opts.long_point_count if opts.long_point_count is not None else opts.point_count + lat_count = opts.lat_point_count if opts.lat_point_count is not None else 1 + if len(opts.lat_range) != 2: + raise RuntimeError('--lat-range must have two values') + + min_lat = math.radians(opts.lat_range[0]) + max_lat = math.radians(opts.lat_range[1]) + lat_spacing = 0.0 if lat_count == 1 else (max_lat - min_lat) / (lat_count - 1) + long_spacing = 2.0 * math.pi / long_count + rng = np.random.default_rng(opts.seed) + + points = np.empty((long_count * lat_count, 3), dtype=np.float64) + idx = 0 + for li in range(lat_count): + lat = min_lat + li * lat_spacing + xy_radius = radius * math.cos(lat) + z = center[2] + radius * math.sin(lat) + if opts.random_spacing: + theta_values = np.sort(rng.uniform(0.0, 2.0 * math.pi, long_count)) + else: + theta_values = np.arange(long_count, dtype=np.float64) * long_spacing + for theta in theta_values: + points[idx, 0] = center[0] + xy_radius * math.cos(theta) + points[idx, 1] = center[1] + xy_radius * math.sin(theta) + points[idx, 2] = z + idx += 1 - points = np.empty((point_count, 3), dtype=np.float64) - points[:, 0] = center[0] + radius * rxy * np.cos(theta) - points[:, 1] = center[1] + radius * rxy * np.sin(theta) - points[:, 2] = center[2] + radius * z return points def generate_torus(point_count, center, major_radius, minor_radius): - # Use a deterministic irrational sequence to avoid requiring two surface counts. idx = np.arange(point_count, dtype=np.float64) golden_ratio_conj = (math.sqrt(5.0) - 1.0) / 2.0 theta = 2.0 * math.pi * idx / point_count @@ -376,17 +522,11 @@ def generate_uniform(point_count, lower, upper, seed): def generate_grid_points(grid_counts, lower, upper): axes = [np.linspace(lower[d], upper[d], grid_counts[d]) for d in range(len(grid_counts))] coords = np.meshgrid(*axes, indexing='ij') - return np.column_stack([coord.ravel() for coord in coords]) + return np.column_stack([coord.ravel(order='F') for coord in coords]) def generate_grid_domains(grid_counts, lower, upper, domain_counts): - '''Generate grid points one spatial domain at a time.''' - for d in range(len(grid_counts)): - if domain_counts[d] > grid_counts[d]: - raise RuntimeError( - f'domain count {domain_counts[d]} exceeds grid size {grid_counts[d]} ' - f'in {AXES[d]} direction') - + '''Generate regular grid points one spatial domain at a time.''' axes = [np.linspace(lower[d], upper[d], grid_counts[d]) for d in range(len(grid_counts))] chunks = [] global_start = 0 @@ -395,19 +535,38 @@ def generate_grid_domains(grid_counts, lower, upper, domain_counts): for d, idx in enumerate(domain_index): begin, end = split_range(grid_counts[d], domain_counts[d], idx) local_axes.append(axes[d][begin:end]) - coords = np.meshgrid(*local_axes, indexing='ij') - points = np.column_stack([coord.ravel() for coord in coords]) - chunks.append((domain_id, domain_index, points, global_start)) + + local_counts = np.array([len(axis) for axis in local_axes], dtype=int) + if np.any(local_counts == 0): + points = np.empty((0, len(grid_counts)), dtype=np.float64) + else: + coords = np.meshgrid(*local_axes, indexing='ij') + points = np.column_stack([coord.ravel(order='F') for coord in coords]) + + chunks.append({ + 'domain_id': domain_id, + 'domain_index': domain_index, + 'points': points, + 'global_id_start': global_start, + 'grid_counts': local_counts, + }) global_start += points.shape[0] + return chunks def generated_domains(opts, dim, lower, upper, center, stddev, domain_counts): - '''Return generated domains as tuples: id, index, points, global_id_start.''' + '''Return generated domain chunk dictionaries.''' if opts.shape == 'grid': grid_counts = grid_counts_from_options(opts, dim) if opts.verbose: print(f'grid_size={grid_counts.tolist()} point_count={int(np.prod(grid_counts))}') + + if opts.grid_topology in ('structured', 'unstructured') and len(domain_counts) != dim: + raise RuntimeError( + 'structured and unstructured grid topologies require --domain-counts with ' + f'{dim} values') + if len(domain_counts) == dim: return generate_grid_domains(grid_counts, lower, upper, domain_counts) @@ -415,9 +574,10 @@ def generated_domains(opts, dim, lower, upper, center, stddev, domain_counts): return split_points(points, domain_counts) if opts.shape == 'circle': - points = generate_circle(opts.point_count, center, opts.radius) + points = generate_circle(opts.point_count, center, opts.radius, opts.random_spacing, + opts.seed) elif opts.shape == 'sphere': - points = generate_sphere(opts.point_count, center, opts.radius) + points = generate_sphere(opts, center, opts.radius) elif opts.shape == 'torus': points = generate_torus(opts.point_count, center, opts.major_radius, opts.minor_radius) elif opts.shape == 'gaussian': @@ -430,42 +590,109 @@ def generated_domains(opts, dim, lower, upper, center, stddev, domain_counts): return split_points(points, domain_counts) -def fill_domain(dom, opts, domain_id, points, global_id_start): - '''Fill one blueprint domain with point mesh data.''' +def add_unstructured_point_topology(topo, point_count): + topo['type'] = 'unstructured' + topo['elements/shape'] = 'point' + topo['elements/connectivity'].set(np.arange(point_count, dtype=np.int32)) + topo['elements/sizes'].set(np.ones(point_count, dtype=np.int32)) + topo['elements/offsets'].set(np.arange(point_count, dtype=np.int32)) + + +def point_index(i, j, k, counts): + return i + counts[0] * (j + counts[1] * k) + + +def grid_connectivity(counts): + '''Return unstructured quad/hex connectivity for a regular grid domain.''' + dim = len(counts) + if np.any(counts < 2): + return np.empty(0, dtype=np.int32), 'quad' if dim == 2 else 'hex' + + conn = [] + if dim == 2: + for j in range(counts[1] - 1): + for i in range(counts[0] - 1): + conn.extend([ + point_index(i, j, 0, counts), + point_index(i + 1, j, 0, counts), + point_index(i + 1, j + 1, 0, counts), + point_index(i, j + 1, 0, counts), + ]) + return np.array(conn, dtype=np.int32), 'quad' + + for k in range(counts[2] - 1): + for j in range(counts[1] - 1): + for i in range(counts[0] - 1): + conn.extend([ + point_index(i, j, k, counts), + point_index(i + 1, j, k, counts), + point_index(i + 1, j + 1, k, counts), + point_index(i, j + 1, k, counts), + point_index(i, j, k + 1, counts), + point_index(i + 1, j, k + 1, counts), + point_index(i + 1, j + 1, k + 1, counts), + point_index(i, j + 1, k + 1, counts), + ]) + return np.array(conn, dtype=np.int32), 'hex' + + +def fill_domain(dom, opts, chunk): + '''Fill one blueprint domain with mesh data.''' + points = chunk['points'] point_count = points.shape[0] dim = points.shape[1] - dom['state/domain_id'] = int(domain_id) - dom['coordsets/' + opts.coordset_name + '/type'] = 'explicit' + dom['state/domain_id'] = int(chunk['domain_id']) + dom[f'coordsets/{opts.coordset_name}/type'] = 'explicit' for d in range(dim): values_path = f'coordsets/{opts.coordset_name}/values/{AXES[d]}' dom[values_path].set(np.ascontiguousarray(points[:, d], dtype=np.float64)) - topo = dom['topologies/' + opts.topology_name] - topo['type'] = opts.topology_type + topo = dom[f'topologies/{opts.topology_name}'] topo['coordset'] = opts.coordset_name - if opts.topology_type == 'unstructured': - topo['elements/shape'] = 'point' - topo['elements/connectivity'].set(np.arange(point_count, dtype=np.int32)) - topo['elements/sizes'].set(np.ones(point_count, dtype=np.int32)) - topo['elements/offsets'].set(np.arange(point_count, dtype=np.int32)) + + if opts.shape == 'grid': + grid_topology = opts.grid_topology + if grid_topology == 'points': + topo['type'] = 'points' + elif grid_topology == 'unstructured-points': + add_unstructured_point_topology(topo, point_count) + elif grid_topology == 'structured': + topo['type'] = 'structured' + counts = chunk['grid_counts'] + for d in range(dim): + topo[f'elements/dims/{LOGICAL_AXES[d]}'] = max(int(counts[d]) - 1, 0) + elif grid_topology == 'unstructured': + topo['type'] = 'unstructured' + counts = chunk['grid_counts'] + conn, shape = grid_connectivity(counts) + topo['elements/shape'] = shape + topo['elements/connectivity'].set(conn) + else: + raise RuntimeError(f'Unsupported grid topology: {grid_topology}') + elif opts.topology_type == 'points': + topo['type'] = 'points' + else: + add_unstructured_point_topology(topo, point_count) if opts.id_field: field = dom['fields/global_id'] field['association'] = 'vertex' field['topology'] = opts.topology_name field['values'].set( - np.arange(global_id_start, global_id_start + point_count, dtype=np.int64)) + np.arange(chunk['global_id_start'], + chunk['global_id_start'] + point_count, + dtype=np.int64)) -def create_domain(md_mesh, opts, domain_id, domain_index, points, global_id_start): +def create_domain(md_mesh, opts, chunk): '''Append one domain to the multidomain mesh.''' if opts.use_list: dom = md_mesh.append() else: - dom = md_mesh[domain_name(domain_index)] + dom = md_mesh[domain_name(chunk['domain_index'])] - fill_domain(dom, opts, domain_id, points, global_id_start) + fill_domain(dom, opts, chunk) def local_chunks_for_rank(chunks, rank, size, single_domain): @@ -505,17 +732,17 @@ def main(): mpi = get_mpi_context() opts = parse_args() dim = infer_dimension(opts) - lower = vector_option(opts.ml, dim, -1.0, '-ml') - upper = vector_option(opts.mu, dim, 1.0, '-mu') - center = vector_option(opts.center, dim, 0.0, '--center') - stddev = vector_option(opts.stddev, dim, 1.0, '--stddev') + lower = vector_option(getattr(opts, 'lower', None), dim, -1.0, '--min') + upper = vector_option(getattr(opts, 'upper', None), dim, 1.0, '--max') + center = vector_option(getattr(opts, 'center', None), dim, 0.0, '--center') + stddev = vector_option(getattr(opts, 'stddev', None), dim, 1.0, '--stddev') domain_counts = domain_counts_from_options(opts, dim, mpi['size']) if np.any(upper <= lower): - raise RuntimeError('-mu values must be greater than -ml values') - if opts.radius <= 0.0: + raise RuntimeError('--max values must be greater than --min values') + if hasattr(opts, 'radius') and opts.radius <= 0.0: raise RuntimeError('--radius must be positive') - if opts.major_radius <= 0.0 or opts.minor_radius <= 0.0: + if hasattr(opts, 'major_radius') and (opts.major_radius <= 0.0 or opts.minor_radius <= 0.0): raise RuntimeError('--major-radius and --minor-radius must be positive') if np.any(stddev <= 0.0): raise RuntimeError('--stddev values must be positive') @@ -526,11 +753,10 @@ def main(): mesh = conduit.Node() if opts.single_domain: if local_chunks: - domain_id, _, points, global_id_start = local_chunks[0] - fill_domain(mesh, opts, domain_id, points, global_id_start) + fill_domain(mesh, opts, local_chunks[0]) else: - for domain_id, domain_index, points, global_id_start in local_chunks: - create_domain(mesh, opts, domain_id, domain_index, points, global_id_start) + for chunk in local_chunks: + create_domain(mesh, opts, chunk) if opts.verbose: print(f'rank {mpi["rank"]} mesh:') @@ -544,11 +770,10 @@ def main(): save_mesh(mesh, opts, mpi) if mpi['rank'] == 0: - total_points = sum(points.shape[0] for _, _, points, _ in chunks) + total_points = sum(chunk['points'].shape[0] for chunk in chunks) mesh_kind = 'single-domain' if opts.single_domain else 'multidomain' - print( - f'Wrote {total_points} {dim}D {opts.shape} points as a {mesh_kind} mesh ' - f'with {len(chunks)} domain(s) to {opts.output} using {opts.protocol}') + print(f'Wrote {total_points} {dim}D {opts.shape} points as a {mesh_kind} mesh ' + f'with {len(chunks)} domain(s) to {opts.output} using {opts.protocol}') return 0 From aa0908c65c1979ed22114ebc5d256fa95963e8cf Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Thu, 13 Aug 2026 21:01:56 -0700 Subject: [PATCH 07/16] Adds back verification for analytic test cases in DCP example --- ...est_distributed_distance_query_example.cpp | 428 +++++++++++++++++- src/tools/gen-multidom-point-mesh.py | 347 +++++++++++++- 2 files changed, 765 insertions(+), 10 deletions(-) diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index 0b9bbcd316..b74898f458 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -36,6 +36,7 @@ #include "mpi.h" // C/C++ includes +#include #include #include #include @@ -45,6 +46,7 @@ #include #include #include +#include #if defined(__GLIBC__) #include // mallinfo2 / mallinfo / malloc_trim #endif @@ -64,6 +66,17 @@ int my_rank = -1, num_ranks = -1; // padded on both sides with '=' symbols std::string banner(const std::string& str) { return axom::fmt::format("{:=^80}", str); } +void broadcastString(std::string& value, int root, MPI_Comm comm) +{ + int length = static_cast(value.size()); + MPI_Bcast(&length, 1, MPI_INT, root, comm); + value.resize(length); + if(length > 0) + { + MPI_Bcast(value.data(), length, MPI_CHAR, root, comm); + } +} + /// Struct to parse and store the input parameters struct Input { @@ -274,6 +287,9 @@ struct BlueprintParticleMesh } int dimension() const { return m_dimension; } + bool hasVerification() const { return m_hasVerification; } + const conduit::Node& verification() const { return m_verification; } + const std::string& description() const { return m_description; } /*! @brief Read a blueprint mesh and store it internally in m_group. @@ -297,6 +313,40 @@ struct BlueprintParticleMesh } conduit::index_t domCount = conduit::blueprint::mesh::number_of_domains(mdMesh); + m_hasVerification = false; + m_verification.reset(); + m_description.clear(); + for(conduit::index_t d = 0; d < domCount; ++d) + { + const conduit::Node& domain = mdMesh.child(d); + if(m_description.empty() && domain.has_path("state/description")) + { + m_description = domain.fetch_existing("state/description").as_string(); + } + if(!m_hasVerification && domain.has_path("state/verification")) + { + m_verification.update(domain.fetch_existing("state/verification")); + m_hasVerification = true; + } + } + + int verificationRank = m_hasVerification ? m_rank : m_nranks; + MPI_Allreduce(MPI_IN_PLACE, &verificationRank, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); + if(verificationRank < m_nranks) + { + std::string description = m_rank == verificationRank ? m_description : std::string {}; + std::string verificationYaml = + m_rank == verificationRank ? m_verification.to_yaml() : std::string {}; + + broadcastString(description, verificationRank, MPI_COMM_WORLD); + broadcastString(verificationYaml, verificationRank, MPI_COMM_WORLD); + + m_description = description; + m_verification.reset(); + m_verification.parse(verificationYaml, "yaml"); + m_hasVerification = true; + } + if(domCount > 0) { if(m_topologyName.empty()) @@ -708,8 +758,11 @@ struct BlueprintParticleMesh private: //!@brief Whether stride/offsets are given for blueprint mesh coordinates data. bool m_coordsAreStrided = false; + bool m_hasVerification = false; std::string m_topologyName; std::string m_coordsetName; + std::string m_description; + conduit::Node m_verification; /// Parent group for the entire mesh sidre::Group* m_group; @@ -742,6 +795,9 @@ class ObjectMeshWrapper std::string getTopologyName() const { return m_objectMesh.getTopologyName(); } std::string getCoordsetName() const { return m_objectMesh.getCoordsetName(); } + bool hasVerification() const { return m_objectMesh.hasVerification(); } + const conduit::Node& verification() const { return m_objectMesh.verification(); } + const std::string& description() const { return m_objectMesh.description(); } private: BlueprintParticleMesh m_objectMesh; @@ -919,6 +975,345 @@ void computeDistancesAndDirections(BlueprintParticleMesh& queryMesh, } } +std::vector conduitVector(const conduit::Node& node) +{ + conduit::Node tmp; + const conduit::Node* src = &node; + if(!node.dtype().is_float64()) + { + node.to_float64_array(tmp); + src = &tmp; + } + + conduit::float64_array values = src->as_float64_array(); + std::vector result(values.number_of_elements()); + for(conduit::index_t i = 0; i < values.number_of_elements(); ++i) + { + result[i] = values[i]; + } + return result; +} + +double conduitDouble(const conduit::Node& node, const std::string& path, double defaultValue) +{ + return node.has_path(path) ? node.fetch_existing(path).to_double() : defaultValue; +} + +template +primal::Point pointFromVector(const std::vector& values) +{ + primal::Point result(0.0); + for(int d = 0; d < DIM && d < static_cast(values.size()); ++d) + { + result[d] = values[d]; + } + return result; +} + +template +primal::Vector vectorFromVector(const std::vector& values) +{ + primal::Vector result(0.0); + for(int d = 0; d < DIM && d < static_cast(values.size()); ++d) + { + result[d] = values[d]; + } + return result; +} + +template +struct AnalyticTorus +{ + using PointType = primal::Point; + + PointType center; + double majorRadius {0.0}; + double minorRadius {0.0}; + + double computeSignedDistance(const PointType& pt) const + { + if constexpr(DIM == 2) + { + const double radialDistance = std::sqrt(primal::squared_distance(pt, center)); + return std::fabs(radialDistance - majorRadius) - minorRadius; + } + else + { + primal::Point pt_xy({pt[0], pt[1]}); + primal::Point center_xy({center[0], center[1]}); + const double radialDistance = + std::sqrt(primal::squared_distance(pt_xy, center_xy)) - majorRadius; + const double axialDistance = pt[2] - center[2]; + return std::sqrt(radialDistance * radialDistance + axialDistance * axialDistance) - minorRadius; + } + } +}; + +template +using AnalyticPrimitive = + std::variant, primal::Plane, AnalyticTorus>; + +template +bool hasPrimitive(const AnalyticPrimitive& primitive) +{ + return !std::holds_alternative(primitive); +} + +template +bool supportsDistanceEnvelope(const AnalyticPrimitive& primitive) +{ + return hasPrimitive(primitive) && !std::holds_alternative>(primitive); +} + +template +double signedDistance(const std::monostate&, const primal::Point&) +{ + return axom::numeric_limits::max(); +} + +template +double signedDistance(const primal::Sphere& sphere, const primal::Point& pt) +{ + return sphere.computeSignedDistance(pt); +} + +template +double signedDistance(const primal::Plane& plane, const primal::Point& pt) +{ + return plane.signedDistance(pt); +} + +template +double signedDistance(const AnalyticTorus& torus, const primal::Point& pt) +{ + return torus.computeSignedDistance(pt); +} + +template +double analyticDistance(const AnalyticPrimitive& primitive, const primal::Point& pt) +{ + return std::visit([&](const auto& shape) { return std::fabs(signedDistance(shape, pt)); }, + primitive); +} + +struct AnalyticVerification +{ + std::string shapeName; + std::string description; + int dimension {-1}; + std::vector center; + std::vector normal; + double radius {0.0}; + double majorRadius {0.0}; + double minorRadius {0.0}; + double innerRadius {0.0}; + double outerRadius {0.0}; + double surfaceTolerance {1.0e-8}; + double distanceTolerance {0.0}; + + template + AnalyticPrimitive makePrimitive() const + { + if(dimension != DIM) + { + return std::monostate {}; + } + + const auto c = pointFromVector(center); + if((shapeName == "circle" || shapeName == "sphere") && radius > 0.0) + { + return primal::Sphere(c, radius); + } + + if(shapeName == "plane") + { + const auto n = vectorFromVector(normal); + return n.is_zero() ? AnalyticPrimitive {std::monostate {}} + : AnalyticPrimitive {primal::Plane(n, c)}; + } + + if constexpr(DIM == 2) + { + if(shapeName == "annulus" && outerRadius > innerRadius && innerRadius > 0.0) + { + return AnalyticTorus {c, + 0.5 * (innerRadius + outerRadius), + 0.5 * (outerRadius - innerRadius)}; + } + } + else if constexpr(DIM == 3) + { + if(shapeName == "torus" && majorRadius > 0.0 && minorRadius > 0.0) + { + return AnalyticTorus {c, majorRadius, minorRadius}; + } + } + + return std::monostate {}; + } + + static AnalyticVerification fromNode(const conduit::Node& node, const std::string& description) + { + AnalyticVerification result; + result.description = description; + if(!node.has_path("shape")) + { + return result; + } + + result.shapeName = node.fetch_existing("shape").as_string(); + result.dimension = node.has_path("dimension") ? node.fetch_existing("dimension").to_int32() : -1; + if(node.has_path("center")) + { + result.center = conduitVector(node.fetch_existing("center")); + } + if(node.has_path("normal")) + { + result.normal = conduitVector(node.fetch_existing("normal")); + } + result.radius = conduitDouble(node, "radius", result.radius); + result.majorRadius = conduitDouble(node, "major_radius", result.majorRadius); + result.minorRadius = conduitDouble(node, "minor_radius", result.minorRadius); + result.innerRadius = conduitDouble(node, "inner_radius", result.innerRadius); + result.outerRadius = conduitDouble(node, "outer_radius", result.outerRadius); + result.surfaceTolerance = conduitDouble(node, "surface_tolerance", result.surfaceTolerance); + result.distanceTolerance = conduitDouble(node, "distance_tolerance", result.distanceTolerance); + return result; + } +}; + +template +int verifyAnalyticClosestPoints(BlueprintParticleMesh& queryMesh, + const AnalyticVerification& verification, + double distThreshold) +{ + SLIC_ASSERT(queryMesh.dimension() == DIM); + + using PointType = primal::Point; + using IndexSet = slam::PositionSet<>; + + const AnalyticPrimitive primitive = verification.makePrimitive(); + if(!hasPrimitive(primitive)) + { + SLIC_WARNING( + axom::fmt::format("Skipping unsupported analytic verification '{}'", verification.shapeName)); + return 0; + } + const bool checkDistanceEnvelope = supportsDistanceEnvelope(primitive); + + queryMesh.registerNodalScalarField("verification_error"); + + int localErrCount = 0; + int localLogCount = 0; + constexpr int MAX_LOCAL_LOGS = 8; + const double distTol = std::max(verification.distanceTolerance, 1.0e-10); + const double surfaceTol = std::max(verification.surfaceTolerance, 1.0e-10); + + auto logError = [&](const std::string& message) { + if(localLogCount < MAX_LOCAL_LOGS) + { + SLIC_INFO(message); + } + ++localLogCount; + }; + + for(axom::IndexType di = 0; di < queryMesh.domain_count(); ++di) + { + auto cpCoords = queryMesh.getNodalVectorField("cp_coords", di); + auto cpIndices = queryMesh.getNodalScalarField("cp_index", di); + auto errorFlag = queryMesh.getNodalScalarField("verification_error", di); + axom::Array qPts = queryMesh.getVertexPositions(di); + + for(auto ptIdx : IndexSet(queryMesh.numPoints(di))) + { + const PointType qPt(&qPts[ptIdx][0]); + const double analyticDist = analyticDistance(primitive, qPt); + const bool hasCp = cpIndices[ptIdx] >= 0; + bool err = false; + + if(hasCp) + { + const PointType& cp = cpCoords[ptIdx]; + const double cpDist = std::sqrt(primal::squared_distance(qPt, cp)); + const double surfaceResidual = analyticDistance(primitive, cp); + const double eps = 1.0e-10 * (1.0 + std::max(cpDist, analyticDist)); + + if(surfaceResidual > surfaceTol) + { + err = true; + logError(axom::fmt::format( + "***Error: Closest point {} on domain {} has analytic residual {} for '{}'.", + cp, + di, + surfaceResidual, + verification.shapeName)); + } + + if(cpDist + eps < analyticDist) + { + err = true; + logError(axom::fmt::format( + "***Error: Discrete closest distance {} is below analytic distance {} at {}.", + cpDist, + analyticDist, + qPt)); + } + + if(checkDistanceEnvelope) + { + if(cpDist > analyticDist + distTol + eps) + { + err = true; + logError(axom::fmt::format( + "***Error: Discrete closest distance {} exceeds analytic distance {} " + "plus sampling tolerance {} at {}.", + cpDist, + analyticDist, + distTol, + qPt)); + } + + if(analyticDist > distThreshold + distTol + eps) + { + err = true; + logError( + axom::fmt::format("***Error: Query point {} is analytically outside threshold by {} " + "but has closest point {}.", + qPt, + analyticDist - distThreshold, + cp)); + } + } + } + else if(checkDistanceEnvelope && analyticDist < distThreshold - distTol) + { + err = true; + logError( + axom::fmt::format("***Error: Query point {} is analytically inside threshold by {} " + "but lacks a closest point.", + qPt, + distThreshold - analyticDist)); + } + + errorFlag[ptIdx] = err ? 1 : 0; + localErrCount += err ? 1 : 0; + } + } + + int globalErrCount = 0; + MPI_Allreduce(&localErrCount, &globalErrCount, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + int globalLogCount = 0; + MPI_Allreduce(&localLogCount, &globalLogCount, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + SLIC_INFO( + axom::fmt::format("Analytic verification for '{}' found {} errors " + "({} detailed messages{}).", + verification.shapeName, + globalErrCount, + globalLogCount, + globalLogCount > MAX_LOCAL_LOGS * num_ranks ? " shown partly" : "")); + return globalErrCount; +} + void make_coords_contiguous(conduit::Node& coordValues) { bool isInterleaved = conduit::blueprint::mcarray::is_interleaved(coordValues); @@ -1319,6 +1714,20 @@ int main(int argc, char** argv) ObjectMeshWrapper objectMeshWrapper(dataStore.getRoot()->createGroup("object_mesh", true), params.objectMeshFile); + AnalyticVerification analyticVerification; + const bool hasAnalyticVerification = objectMeshWrapper.hasVerification(); + if(hasAnalyticVerification) + { + analyticVerification = AnalyticVerification::fromNode(objectMeshWrapper.verification(), + objectMeshWrapper.description()); + if(my_rank == 0) + { + SLIC_INFO(axom::fmt::format("Loaded analytic verification metadata for '{}': {}", + analyticVerification.shapeName, + objectMeshWrapper.description())); + } + } + if(params.isVerbose()) { objectMeshWrapper.getParticleMesh().printMeshSizeStats("Object mesh"); @@ -1455,6 +1864,23 @@ int main(int argc, char** argv) "direction"); } + int globalVerificationErrors = 0; + if(hasAnalyticVerification) + { + if(spatialDim == 2) + { + globalVerificationErrors = verifyAnalyticClosestPoints<2>(queryMeshWrapper.getParticleMesh(), + analyticVerification, + params.distThreshold); + } + else if(spatialDim == 3) + { + globalVerificationErrors = verifyAnalyticClosestPoints<3>(queryMeshWrapper.getParticleMesh(), + analyticVerification, + params.distThreshold); + } + } + // queryMeshNode.print(); queryMeshNode.reset(); @@ -1465,5 +1891,5 @@ int main(int argc, char** argv) finalizeLogger(); MPI_Finalize(); - return 0; + return globalVerificationErrors == 0 ? 0 : 1; } diff --git a/src/tools/gen-multidom-point-mesh.py b/src/tools/gen-multidom-point-mesh.py index 5c80f3e83d..17c80933fb 100755 --- a/src/tools/gen-multidom-point-mesh.py +++ b/src/tools/gen-multidom-point-mesh.py @@ -8,6 +8,18 @@ # / (multidomain output) # |-- state # | `-- domain_id == +# | `-- description (analytic/input description string; stored in state +# | so Blueprint relay preserves it) +# | `-- verification (optional analytic verification metadata; stored in state) +# | |-- shape == "circle", "sphere", "torus", "plane", or "annulus" +# | |-- dimension == 2 or 3 +# | |-- center (float64, dimension) +# | |-- normal (float64, dimension; plane only) +# | |-- radius (circle/sphere) +# | |-- {major,minor}_radius (torus) +# | |-- {inner,outer}_radius (annulus) +# | |-- surface_tolerance (float64 residual tolerance for shape membership) +# | `-- distance_tolerance (float64 sampling slack for DCP checks) # |-- topologies # | `-- # | |-- coordset == @@ -30,6 +42,42 @@ # |-- topology == # |-- association == "vertex" # `-- values (int64, point_count) +# +# Usage examples: +# # 2D circle object mesh split over N MPI ranks, with some empty domains: +# srun -n N build-axom/bin/run_python_with_axom.sh \\ +# src/tools/gen-multidom-point-mesh.py circle --point-count 64 --domains 4 \\ +# --center 0.7,0.9 --radius 0.9 --output dcp_object_circle +# +# # 3D sphere object mesh split over N MPI ranks using latitude/longitude samples: +# srun -n N build-axom/bin/run_python_with_axom.sh \\ +# src/tools/gen-multidom-point-mesh.py sphere --long-point-count 60 \\ +# --lat-point-count 30 --center 0.7,0.9,0.5 --radius 0.9 \\ +# --domains 6 --output dcp_object_sphere +# +# # 3D torus object mesh: +# build-axom/bin/run_python_with_axom.sh \\ +# src/tools/gen-multidom-point-mesh.py torus --point-count 1024 \\ +# --center 0,0,0 --major-radius 0.75 --minor-radius 0.2 \\ +# --domains 4 --output dcp_object_torus +# +# # 3D codimension-1 plane object mesh: +# build-axom/bin/run_python_with_axom.sh \\ +# src/tools/gen-multidom-point-mesh.py plane --dimension 3 --grid-size 32,32 \\ +# --center 0,0,0 --normal 0,0,1 --extent 2,2 --domains 4 \\ +# --output dcp_object_plane +# +# # 2D structured quad query mesh; DCP queries are evaluated over the vertices: +# build-axom/bin/run_python_with_axom.sh \\ +# src/tools/gen-multidom-point-mesh.py grid --grid-size 100,100 \\ +# --min 0,0 --max 2,2 --domain-counts 2,1 --grid-topology structured \\ +# --output dcp_query_structured_quads +# +# # 3D unstructured hex query mesh: +# build-axom/bin/run_python_with_axom.sh \\ +# src/tools/gen-multidom-point-mesh.py grid --grid-size 20,20,15 \\ +# --min 0,0,0 --max 2,2,2 --domain-counts 2,2,1 \\ +# --grid-topology unstructured --output dcp_query_unstructured_hexes # This script requires a conduit installation configured with python3 and hdf5. # Make sure PYTHONPATH includes /path/to/conduit/install/python-modules, @@ -153,6 +201,9 @@ def add_common_options(parser): bp.add_argument('--id-field', action='store_true', help='Add a vertex-associated int64 global_id field') + bp.add_argument('--no-verification', + action='store_true', + help='Omit analytic verification metadata') bp.add_argument('--protocol', default='hdf5', help='Conduit relay protocol for save_mesh') bp.add_argument('-o', '--output', default='mdmesh', help='Output file base name') bp.add_argument('-v', '--verbose', action='store_true', help='Print additional info') @@ -246,6 +297,42 @@ def parse_args(): torus.add_argument('--major-radius', type=float, default=1.0, help='Torus major radius') torus.add_argument('--minor-radius', type=float, default=0.25, help='Torus minor radius') + annulus = subparsers.add_parser('annulus', + formatter_class=ArgumentDefaultsHelpFormatter, + help='2D analytic annulus boundary point mesh') + add_common_options(annulus) + add_point_topology_options(annulus) + annulus.add_argument('-n', + '--point-count', + type=positive_int, + default=1024, + help='Total number of points on the inner and outer circles') + annulus.add_argument('--center', type=f_c, help='Center coordinates; scalar or 2 values') + annulus.add_argument('--inner-radius', type=float, default=0.5, help='Inner radius') + annulus.add_argument('--outer-radius', type=float, default=1.0, help='Outer radius') + annulus.add_argument('--random-spacing', + action='store_true', + help='Use random angular spacing instead of uniform spacing') + annulus.add_argument('--seed', type=int, default=0, help='Random seed for random spacing') + + plane = subparsers.add_parser('plane', + formatter_class=ArgumentDefaultsHelpFormatter, + help='2D line or 3D plane point mesh') + add_common_options(plane) + add_point_topology_options(plane) + plane.add_argument('-d', '--dimension', type=int, choices=(2, 3), help='Spatial dimension') + plane.add_argument('-n', + '--point-count', + type=positive_int, + default=1024, + help='Target point count when --grid-size is omitted') + plane.add_argument('--grid-size', + type=i_c, + help='Point count along the line, or two counts on the plane') + plane.add_argument('--center', type=f_c, help='Plane center; scalar or per dimension') + plane.add_argument('--normal', type=f_c, help='Plane normal; scalar or per dimension') + plane.add_argument('--extent', type=f_c, help='Line length in 2D, or two side lengths in 3D') + grid = subparsers.add_parser('grid', formatter_class=ArgumentDefaultsHelpFormatter, help='Regular grid mesh over --min/--max') @@ -308,13 +395,17 @@ def parse_args(): def infer_dimension(opts): '''Infer and validate the spatial dimension.''' - fixed_dims = {'circle': 2, 'sphere': 3, 'torus': 3} + fixed_dims = {'circle': 2, 'sphere': 3, 'torus': 3, 'annulus': 2} if opts.shape in fixed_dims: return fixed_dims[opts.shape] dim = opts.dimension inferred = [] - for name in ('grid_size', 'lower', 'upper', 'center', 'stddev', 'domain_counts'): + names = ('grid_size', 'lower', 'upper', 'center', 'stddev', 'domain_counts') + if opts.shape == 'plane': + # Plane grid_size/extent are intrinsic coordinates, not spatial dimension. + names = ('center', 'normal', 'domain_counts') + for name in names: values = getattr(opts, name, None) if values is not None: inferred.append(len(values)) @@ -338,6 +429,11 @@ def infer_dimension(opts): def vector_option(values, dim, default, name): '''Return a dimension-sized numpy vector from an optional scalar/list argument.''' if values is None: + if hasattr(default, '__len__'): + default = np.array(default, dtype=float) + if len(default) != dim: + raise RuntimeError(f'default for {name} must have {dim} values') + return default return np.full(dim, default, dtype=float) if len(values) == 1: return np.full(dim, values[0], dtype=float) @@ -346,6 +442,20 @@ def vector_option(values, dim, default, name): return np.array(values, dtype=float) +def normalized(vec, name): + '''Return a normalized copy of a vector.''' + norm = np.linalg.norm(vec) + if norm <= 0.0: + raise RuntimeError(f'{name} must be nonzero') + return vec / norm + + +def default_normal(dim): + normal = np.zeros(dim, dtype=float) + normal[-1] = 1.0 + return normal + + def domain_counts_from_options(opts, dim, mpi_size): '''Return domain counts as an integer vector.''' if opts.single_domain: @@ -509,6 +619,81 @@ def generate_torus(point_count, center, major_radius, minor_radius): return points +def generate_annulus(point_count, center, inner_radius, outer_radius, random_spacing, seed): + if point_count < 2: + raise RuntimeError('annulus requires at least two points') + + inner_count = point_count // 2 + outer_count = point_count - inner_count + rng = np.random.default_rng(seed) + + def ring_points(count, radius): + if random_spacing: + theta = np.sort(rng.uniform(0.0, 2.0 * math.pi, count)) + else: + theta = np.linspace(0.0, 2.0 * math.pi, count, endpoint=False) + ring = np.empty((count, 2), dtype=np.float64) + ring[:, 0] = center[0] + radius * np.cos(theta) + ring[:, 1] = center[1] + radius * np.sin(theta) + return ring + + return np.vstack((ring_points(inner_count, + inner_radius), ring_points(outer_count, outer_radius))) + + +def plane_basis(normal): + '''Return one tangent in 2D or two orthonormal tangents in 3D.''' + if len(normal) == 2: + return (np.array([-normal[1], normal[0]], dtype=np.float64), ) + + ref = np.array([1.0, 0.0, 0.0], dtype=np.float64) + if abs(np.dot(ref, normal)) > 0.9: + ref = np.array([0.0, 1.0, 0.0], dtype=np.float64) + tangent0 = normalized(np.cross(normal, ref), '--normal') + tangent1 = normalized(np.cross(normal, tangent0), '--normal') + return tangent0, tangent1 + + +def plane_grid_counts(opts, dim): + if opts.grid_size is not None: + counts = np.array(opts.grid_size, dtype=int) + expected = 1 if dim == 2 else 2 + if len(counts) != expected: + raise RuntimeError(f'--grid-size for a {dim}D plane must have {expected} value(s)') + if np.any(counts < 1): + raise RuntimeError('plane grid sizes must be positive') + return counts + + if dim == 2: + return np.array([opts.point_count], dtype=int) + + base = max(1, int(round(math.sqrt(opts.point_count)))) + counts = np.array([base, base], dtype=int) + while int(np.prod(counts)) < opts.point_count: + counts[np.argmin(counts)] += 1 + return counts + + +def generate_plane(opts, dim, center, normal): + counts = plane_grid_counts(opts, dim) + extent = vector_option(getattr(opts, 'extent', None), len(counts), 2.0, '--extent') + if np.any(extent <= 0.0): + raise RuntimeError('--extent values must be positive') + + basis = plane_basis(normal) + if dim == 2: + t = np.linspace(-0.5 * extent[0], 0.5 * extent[0], counts[0]) + return center + np.outer(t, basis[0]) + + axes = [np.linspace(-0.5 * extent[d], 0.5 * extent[d], counts[d]) for d in range(2)] + coords = np.meshgrid(*axes, indexing='ij') + weights = np.column_stack([coord.ravel(order='F') for coord in coords]) + points = np.empty((weights.shape[0], 3), dtype=np.float64) + for i, weight in enumerate(weights): + points[i] = center + weight[0] * basis[0] + weight[1] * basis[1] + return points + + def generate_gaussian(point_count, center, stddev, seed): rng = np.random.default_rng(seed) return rng.normal(loc=center, scale=stddev, size=(point_count, len(center))) @@ -555,7 +740,7 @@ def generate_grid_domains(grid_counts, lower, upper, domain_counts): return chunks -def generated_domains(opts, dim, lower, upper, center, stddev, domain_counts): +def generated_domains(opts, dim, lower, upper, center, normal, stddev, domain_counts): '''Return generated domain chunk dictionaries.''' if opts.shape == 'grid': grid_counts = grid_counts_from_options(opts, dim) @@ -580,6 +765,11 @@ def generated_domains(opts, dim, lower, upper, center, stddev, domain_counts): points = generate_sphere(opts, center, opts.radius) elif opts.shape == 'torus': points = generate_torus(opts.point_count, center, opts.major_radius, opts.minor_radius) + elif opts.shape == 'annulus': + points = generate_annulus(opts.point_count, center, opts.inner_radius, opts.outer_radius, + opts.random_spacing, opts.seed) + elif opts.shape == 'plane': + points = generate_plane(opts, dim, center, normal) elif opts.shape == 'gaussian': points = generate_gaussian(opts.point_count, center, stddev, opts.seed) elif opts.shape == 'uniform': @@ -636,13 +826,141 @@ def grid_connectivity(counts): return np.array(conn, dtype=np.int32), 'hex' -def fill_domain(dom, opts, chunk): +def format_vector(values): + return '(' + ', '.join(f'{value:g}' for value in values) + ')' + + +def shape_description(opts, dim, lower, upper, center, normal, stddev): + if opts.shape == 'circle': + return f'2D circle centered at {format_vector(center)} with radius {opts.radius:g}' + if opts.shape == 'sphere': + return f'3D sphere centered at {format_vector(center)} with radius {opts.radius:g}' + if opts.shape == 'torus': + return (f'3D torus centered at {format_vector(center)} with major radius ' + f'{opts.major_radius:g} and minor radius {opts.minor_radius:g}') + if opts.shape == 'annulus': + return (f'2D annulus centered at {format_vector(center)} with inner radius ' + f'{opts.inner_radius:g} and outer radius {opts.outer_radius:g}') + if opts.shape == 'plane': + return f'{dim}D codimension-1 plane centered at {format_vector(center)}' + if opts.shape == 'grid': + return f'{dim}D regular grid over [{format_vector(lower)}, {format_vector(upper)}]' + if opts.shape == 'gaussian': + return (f'{dim}D Gaussian point cloud centered at {format_vector(center)} with stddev ' + f'{format_vector(stddev)}') + if opts.shape == 'uniform': + return f'{dim}D uniform random point cloud over [{format_vector(lower)}, {format_vector(upper)}]' + return f'{dim}D {opts.shape} point mesh' + + +def circle_spacing(radius, point_count): + return 2.0 * math.pi * radius / max(point_count, 1) + + +def sphere_spacing(opts, radius): + if opts.long_point_count is not None or opts.lat_point_count is not None: + long_count = opts.long_point_count if opts.long_point_count is not None else opts.point_count + lat_count = opts.lat_point_count if opts.lat_point_count is not None else 1 + long_spacing = 2.0 * math.pi * radius / max(long_count, 1) + lat_spacing = 0.0 + if lat_count > 1: + lat_spacing = radius * math.radians(opts.lat_range[1] - + opts.lat_range[0]) / (lat_count - 1) + return math.sqrt(long_spacing * long_spacing + lat_spacing * lat_spacing) + + return radius * math.sqrt(4.0 * math.pi / max(opts.point_count, 1)) + + +def torus_spacing(major_radius, minor_radius, point_count): + area = 4.0 * math.pi * math.pi * major_radius * minor_radius + return 2.0 * math.sqrt(area / max(point_count, 1)) + + +def annulus_spacing(opts): + inner_count = max(opts.point_count // 2, 1) + outer_count = max(opts.point_count - inner_count, 1) + return max(circle_spacing(opts.inner_radius, inner_count), + circle_spacing(opts.outer_radius, outer_count)) + + +def plane_spacing(opts, dim): + counts = plane_grid_counts(opts, dim) + extent = vector_option(getattr(opts, 'extent', None), len(counts), 2.0, '--extent') + spacing = [] + for d, count in enumerate(counts): + spacing.append(extent[d] / max(count - 1, 1)) + return max(spacing) + + +def verification_metadata(opts, dim, center, normal): + if opts.no_verification: + return None + + metadata = {'shape': opts.shape, 'dimension': dim, 'surface_tolerance': 1.0e-8} + if opts.shape == 'circle': + metadata.update({ + 'center': center, + 'radius': opts.radius, + 'distance_tolerance': circle_spacing(opts.radius, opts.point_count), + }) + elif opts.shape == 'sphere': + metadata.update({ + 'center': center, + 'radius': opts.radius, + 'distance_tolerance': sphere_spacing(opts, opts.radius), + }) + elif opts.shape == 'torus': + metadata.update({ + 'center': + center, + 'major_radius': + opts.major_radius, + 'minor_radius': + opts.minor_radius, + 'distance_tolerance': + torus_spacing(opts.major_radius, opts.minor_radius, opts.point_count), + }) + elif opts.shape == 'annulus': + metadata.update({ + 'center': center, + 'inner_radius': opts.inner_radius, + 'outer_radius': opts.outer_radius, + 'distance_tolerance': annulus_spacing(opts), + }) + elif opts.shape == 'plane': + metadata.update({ + 'center': center, + 'normal': normal, + 'distance_tolerance': plane_spacing(opts, dim), + }) + else: + return None + + return metadata + + +def add_metadata(dom, metadata): + dom['state/description'] = metadata['description'] + verification = metadata['verification'] + if verification is None: + return + + ver_node = dom['state/verification'] + for key, value in verification.items(): + if isinstance(value, np.ndarray): + ver_node[key].set(np.ascontiguousarray(value, dtype=np.float64)) + else: + ver_node[key] = value + + +def fill_domain(dom, opts, chunk, metadata): '''Fill one blueprint domain with mesh data.''' points = chunk['points'] point_count = points.shape[0] dim = points.shape[1] dom['state/domain_id'] = int(chunk['domain_id']) + add_metadata(dom, metadata) dom[f'coordsets/{opts.coordset_name}/type'] = 'explicit' for d in range(dim): values_path = f'coordsets/{opts.coordset_name}/values/{AXES[d]}' @@ -685,14 +1003,14 @@ def fill_domain(dom, opts, chunk): dtype=np.int64)) -def create_domain(md_mesh, opts, chunk): +def create_domain(md_mesh, opts, chunk, metadata): '''Append one domain to the multidomain mesh.''' if opts.use_list: dom = md_mesh.append() else: dom = md_mesh[domain_name(chunk['domain_index'])] - fill_domain(dom, opts, chunk) + fill_domain(dom, opts, chunk, metadata) def local_chunks_for_rank(chunks, rank, size, single_domain): @@ -735,6 +1053,9 @@ def main(): lower = vector_option(getattr(opts, 'lower', None), dim, -1.0, '--min') upper = vector_option(getattr(opts, 'upper', None), dim, 1.0, '--max') center = vector_option(getattr(opts, 'center', None), dim, 0.0, '--center') + normal = vector_option(getattr(opts, 'normal', None), dim, default_normal(dim), '--normal') + if opts.shape == 'plane': + normal = normalized(normal, '--normal') stddev = vector_option(getattr(opts, 'stddev', None), dim, 1.0, '--stddev') domain_counts = domain_counts_from_options(opts, dim, mpi['size']) @@ -744,19 +1065,27 @@ def main(): raise RuntimeError('--radius must be positive') if hasattr(opts, 'major_radius') and (opts.major_radius <= 0.0 or opts.minor_radius <= 0.0): raise RuntimeError('--major-radius and --minor-radius must be positive') + if hasattr(opts, 'inner_radius') and (opts.inner_radius <= 0.0 + or opts.outer_radius <= opts.inner_radius): + raise RuntimeError('--outer-radius must be greater than positive --inner-radius') if np.any(stddev <= 0.0): raise RuntimeError('--stddev values must be positive') - chunks = generated_domains(opts, dim, lower, upper, center, stddev, domain_counts) + metadata = { + 'description': shape_description(opts, dim, lower, upper, center, normal, stddev), + 'verification': verification_metadata(opts, dim, center, normal), + } + + chunks = generated_domains(opts, dim, lower, upper, center, normal, stddev, domain_counts) local_chunks = local_chunks_for_rank(chunks, mpi['rank'], mpi['size'], opts.single_domain) mesh = conduit.Node() if opts.single_domain: if local_chunks: - fill_domain(mesh, opts, local_chunks[0]) + fill_domain(mesh, opts, local_chunks[0], metadata) else: for chunk in local_chunks: - create_domain(mesh, opts, chunk) + create_domain(mesh, opts, chunk, metadata) if opts.verbose: print(f'rank {mpi["rank"]} mesh:') From ad37fcff8ea9d8ef841604fd0c5f4914f80b46ad Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Fri, 14 Aug 2026 15:48:31 -0700 Subject: [PATCH 08/16] Updates the build skill about running MPI commands Also adds troubleshooting sections about shroud and about symlinked paths. --- skills/building/SKILL.md | 41 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/skills/building/SKILL.md b/skills/building/SKILL.md index dc2f7bd575..aa6f806fd2 100644 --- a/skills/building/SKILL.md +++ b/skills/building/SKILL.md @@ -64,3 +64,44 @@ If the user does **not** specify a host-config file, determine the best match by ``` Use its output as the `-hc` argument (as shown in the examples above). If the user explicitly provides a host-config file/path, use that instead. + +### Symlinked directories in the sandbox + +If configure/build commands report that an existing path is missing, inaccessible, or outside the sandbox, check whether the source, build, install, host-config, or TPL paths include symlinks. Compare logical and physical paths with: + +```bash +pwd -P +readlink -f +``` + +Prefer physical paths resolved by `readlink -f` when invoking `config-build.py` (for `-bp`, `-ip`, `-hc`, and any explicit source/TPL paths), or restart the sandbox from the resolved workspace path. A symlinked path may appear outside the sandbox policy even when its resolved target is visible. + +### MPI and Slurm in the sandbox + +To run MPI-enabled commands from the sandbox, launch Codex with the `--mpi` command-line argument. This starts a Flux instance when needed. MPI through the sandbox is single-node only. If MPI or Slurm commands fail because no Flux allocation is available, stop and ask the user to restart the sandbox with `--mpi`. + +On Slurm-based systems, run Slurm commands through the Flux wrapper instead of the system executable, for example: + +```bash +/usr/global/tools/flux_wrappers/bin/srun -n 2 +``` + +Loading the flux wrappers before launching the agent allows you to run the srun wrapper, e.g. +```bash +module load flux_wrappers +srun -n 2 +``` + +When configuring builds whose MPI tests use `srun`, override CMake's MPI launcher: + +```bash +./config-build.py -hc "$(./skills/building/scripts/determine_host_config)" -DMPIEXEC_EXECUTABLE=/usr/global/tools/flux_wrappers/bin/srun +``` + +### Shroud in the sandbox + +If CMake configuration has trouble with Shroud in the sandbox, unset the cached Shroud executable at configure time: + +```bash +./config-build.py -hc "$(./skills/building/scripts/determine_host_config)" -USHROUD_EXECUTABLE +``` From 81c338e11b916994a4de46d28bcb93c70c5908ab Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Fri, 14 Aug 2026 17:21:40 -0700 Subject: [PATCH 09/16] Quest: DCP example runs on object meshes in test repository The test meshes match the previous analytically generated circle/sphere and are validated. --- src/axom/quest/examples/CMakeLists.txt | 143 ++++++++++--------------- src/tools/gen-multidom-point-mesh.py | 9 ++ 2 files changed, 63 insertions(+), 89 deletions(-) diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt index f5760c3c7d..c2d2e9a650 100644 --- a/src/axom/quest/examples/CMakeLists.txt +++ b/src/axom/quest/examples/CMakeLists.txt @@ -378,81 +378,28 @@ if(AXOM_ENABLE_MPI AND AXOM_ENABLE_SIDRE AND HDF5_FOUND) FOLDER axom/quest/examples ) - if(AXOM_ENABLE_TESTS AND AXOM_ENABLE_PYTHON_TESTS) + if(AXOM_ENABLE_TESTS AND AXOM_DATA_DIR) set(_nranks 3) - set(_dcp_mesh_generator ${PROJECT_SOURCE_DIR}/tools/gen-multidom-point-mesh.py) - - set(_dcp_object_mesh ${CMAKE_CURRENT_BINARY_DIR}/dcp_object_circle_empty_domains) - set(_dcp_gen_object_test quest_distributed_closest_point_gen_object_circle_empty_domains) - axom_add_python_test( - NAME ${_dcp_gen_object_test} - COMMAND ${Python_EXECUTABLE} - ${_dcp_mesh_generator} - circle - --point-count 4 - --center 0.7,0.9 - --radius 0.9 - --domains 6 - --output ${_dcp_object_mesh} - NUM_MPI_TASKS ${_nranks}) - - set(_dcp_object_mesh_3d ${CMAKE_CURRENT_BINARY_DIR}/dcp_object_sphere) - set(_dcp_gen_object_test_3d quest_distributed_closest_point_gen_object_sphere) - axom_add_python_test( - NAME ${_dcp_gen_object_test_3d} - COMMAND ${Python_EXECUTABLE} - ${_dcp_mesh_generator} - sphere - --long-point-count 12 - --lat-point-count 6 - --center 0.7,0.9,0.5 - --radius 0.9 - --domains 6 - --output ${_dcp_object_mesh_3d} - NUM_MPI_TASKS ${_nranks}) # Run the distributed closest point example on N ranks for each enabled policy. - # Query meshes are generated here instead of checked into the data repository. - set(_query_meshes structured_quads unstructured_quads structured_hexes unstructured_hexes) - foreach(_mesh ${_query_meshes}) - set(_query_mesh ${CMAKE_CURRENT_BINARY_DIR}/dcp_query_${_mesh}) - set(_gen_query_test quest_distributed_closest_point_gen_query_${_mesh}) - if(_mesh MATCHES "quads$") + # Query meshes match the original DCP example. Object meshes were + # generated with src/tools/gen-multidom-point-mesh.py. + set(_meshes "mdmesh.2x1" "mdmesh.2x3" "mdmesh.2x2x1" "mdmeshg.2x2x1") + foreach(_mesh ${_meshes}) + string(REGEX MATCH "\\.[0-9]+(x[0-9]+)+$" _sizes "${_mesh}") + string(REGEX MATCHALL "[0-9]+" _sizes ${_sizes}) + list(LENGTH _sizes _ndim) + + if(_ndim EQUAL 2) set(_dim 2) - set(_grid_size 30,30) - set(_min 0,0) - set(_max 2,2) - set(_domain_counts 2,1) - set(_object_mesh ${_dcp_object_mesh}) - set(_object_test ${_dcp_gen_object_test}) - else() + set(_shape circle) + set(_object_mesh ${quest_data_dir}/dcp_object_circle.root) + elseif(_ndim EQUAL 3) set(_dim 3) - set(_grid_size 8,8,6) - set(_min 0,0,0) - set(_max 2,2,2) - set(_domain_counts 2,2,1) - set(_object_mesh ${_dcp_object_mesh_3d}) - set(_object_test ${_dcp_gen_object_test_3d}) - endif() - if(_mesh MATCHES "^structured") - set(_grid_topology structured) - else() - set(_grid_topology unstructured) + set(_shape sphere) + set(_object_mesh ${quest_data_dir}/dcp_object_sphere.root) endif() - axom_add_python_test( - NAME ${_gen_query_test} - COMMAND ${Python_EXECUTABLE} - ${_dcp_mesh_generator} - grid - --grid-size ${_grid_size} - --min ${_min} - --max ${_max} - --domain-counts ${_domain_counts} - --grid-topology ${_grid_topology} - --output ${_query_mesh} - NUM_MPI_TASKS ${_nranks}) - foreach(_pol ${AXOM_EXECUTION_POLICIES}) # sets the number of threads for the omp policy; leave empty for non-omp policies set(_num_threads) @@ -464,8 +411,8 @@ if(AXOM_ENABLE_MPI AND AXOM_ENABLE_SIDRE AND HDF5_FOUND) axom_add_test( NAME ${_test} COMMAND quest_distributed_distance_query_ex - --mesh-file ${_query_mesh}.root - --object-mesh-file ${_object_mesh}.root + --mesh-file ${quest_data_dir}/${_mesh}.root + --object-mesh-file ${_object_mesh} --dist-threshold .3 --dynamic-distance-filtering --policy ${_pol} @@ -473,18 +420,19 @@ if(AXOM_ENABLE_MPI AND AXOM_ENABLE_SIDRE AND HDF5_FOUND) NUM_MPI_TASKS ${_nranks} NUM_OMP_THREADS ${_num_threads}) set_tests_properties(${_test} - PROPERTIES DEPENDS "${_object_test};${_gen_query_test}") + PROPERTIES + PASS_REGULAR_EXPRESSION "Analytic verification for '${_shape}' found 0 errors") if(_pol STREQUAL "seq" AND - (_mesh STREQUAL "structured_quads" OR _mesh STREQUAL "structured_hexes")) + (_mesh STREQUAL "mdmesh.2x1" OR _mesh STREQUAL "mdmesh.2x2x1")) # Keep fallback coverage for the non-dynamic filtering path # with one 2D and one 3D seq case. set(_static_test "${_test}_no_dynamic_distance_filtering") axom_add_test( NAME ${_static_test} COMMAND quest_distributed_distance_query_ex - --mesh-file ${_query_mesh}.root - --object-mesh-file ${_object_mesh}.root + --mesh-file ${quest_data_dir}/${_mesh}.root + --object-mesh-file ${_object_mesh} --dist-threshold .3 --no-dynamic-distance-filtering --policy ${_pol} @@ -492,28 +440,45 @@ if(AXOM_ENABLE_MPI AND AXOM_ENABLE_SIDRE AND HDF5_FOUND) NUM_MPI_TASKS ${_nranks} NUM_OMP_THREADS ${_num_threads}) set_tests_properties(${_static_test} - PROPERTIES DEPENDS "${_object_test};${_gen_query_test}") + PROPERTIES + PASS_REGULAR_EXPRESSION "Analytic verification for '${_shape}' found 0 errors") endif() endforeach() endforeach() - unset(_dcp_gen_object_test) - unset(_dcp_gen_object_test_3d) - unset(_dcp_mesh_generator) - unset(_dcp_object_mesh) - unset(_dcp_object_mesh_3d) + foreach(_pol ${AXOM_EXECUTION_POLICIES}) + # sets the number of threads for the omp policy; leave empty for non-omp policies + set(_num_threads) + if(_pol STREQUAL "omp") + set(_num_threads ${AXOM_TEST_NUM_OMP_THREADS}) + endif() + + set(_test "quest_distributed_closest_point_run_2D_${_pol}_single_domain_object") + axom_add_test( + NAME ${_test} + COMMAND quest_distributed_distance_query_ex + --mesh-file ${quest_data_dir}/mdmesh.2x1.root + --object-mesh-file ${quest_data_dir}/dcp_object_circle_single_domain.root + --dist-threshold .3 + --policy ${_pol} + --distance-file dcp_closest_point_2d_${_pol}_single_domain_object + NUM_MPI_TASKS ${_nranks} + NUM_OMP_THREADS ${_num_threads}) + set_tests_properties(${_test} + PROPERTIES + PASS_REGULAR_EXPRESSION "Analytic verification for 'circle' found 0 errors") + endforeach() + unset(_dim) - unset(_domain_counts) - unset(_gen_query_test) - unset(_grid_size) - unset(_grid_topology) - unset(_max) - unset(_min) + unset(_mesh) + unset(_meshes) + unset(_ndim) + unset(_num_threads) unset(_object_mesh) - unset(_object_test) - unset(_query_mesh) - unset(_query_meshes) + unset(_pol) unset(_nranks) + unset(_shape) + unset(_sizes) unset(_test) unset(_static_test) endif() diff --git a/src/tools/gen-multidom-point-mesh.py b/src/tools/gen-multidom-point-mesh.py index 17c80933fb..73f0b250de 100755 --- a/src/tools/gen-multidom-point-mesh.py +++ b/src/tools/gen-multidom-point-mesh.py @@ -10,6 +10,7 @@ # | `-- domain_id == # | `-- description (analytic/input description string; stored in state # | so Blueprint relay preserves it) +# | `-- command_line (sanitized command line used to generate this mesh) # | `-- verification (optional analytic verification metadata; stored in state) # | |-- shape == "circle", "sphere", "torus", "plane", or "annulus" # | |-- dimension == 2 or 3 @@ -87,6 +88,7 @@ import itertools import math import os +import shlex import sys from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser, ArgumentTypeError @@ -939,8 +941,14 @@ def verification_metadata(opts, dim, center, normal): return metadata +def sanitized_command_line(): + '''Return the generator command line with the script path reduced to its basename.''' + return shlex.join([os.path.basename(sys.argv[0])] + sys.argv[1:]) + + def add_metadata(dom, metadata): dom['state/description'] = metadata['description'] + dom['state/command_line'] = metadata['command_line'] verification = metadata['verification'] if verification is None: return @@ -1073,6 +1081,7 @@ def main(): metadata = { 'description': shape_description(opts, dim, lower, upper, center, normal, stddev), + 'command_line': sanitized_command_line(), 'verification': verification_metadata(opts, dim, center, normal), } From 575ab5fc2f38137bf397d1061dbf8b8e1f53789c Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Fri, 14 Aug 2026 17:50:07 -0700 Subject: [PATCH 10/16] Extents convert_sidre_protocol script to support reading conduit blueprint protocols --- src/tools/CMakeLists.txt | 19 +++- src/tools/convert_sidre_protocol.py | 160 +++++++++++++++++++++++----- 2 files changed, 152 insertions(+), 27 deletions(-) diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt index e8782b9992..5c327cdb4c 100644 --- a/src/tools/CMakeLists.txt +++ b/src/tools/CMakeLists.txt @@ -263,6 +263,23 @@ if(NANOBIND_FOUND) set_tests_properties(${_testname} PROPERTIES PASS_REGULAR_EXPRESSION "Writing out datastore") + + set(_testname "convert_sidre_protocol_py_blueprint") + set(mesh_dir "${AXOM_DATA_DIR}/quest/mdmesh.2x1.root") + + axom_add_python_test( + NAME ${_testname} + COMMAND ${Python_EXECUTABLE} + ${PROJECT_BINARY_DIR}/bin/convert_sidre_protocol.py + --input ${mesh_dir} + --input-type blueprint + --output csp_blueprint_output + --protocol json + NUM_MPI_TASKS 2 + ) + + set_tests_properties(${_testname} PROPERTIES + PASS_REGULAR_EXPRESSION "Writing out Blueprint mesh") else() axom_add_python_test( NAME ${_testname} @@ -272,7 +289,7 @@ if(NANOBIND_FOUND) ) set_tests_properties(${_testname} PROPERTIES - PASS_REGULAR_EXPRESSION "Sidre protocol converter") + PASS_REGULAR_EXPRESSION "Sidre/Blueprint protocol converter") endif() endif() endif() diff --git a/src/tools/convert_sidre_protocol.py b/src/tools/convert_sidre_protocol.py index 4ec3c82dfa..116d6a49f7 100644 --- a/src/tools/convert_sidre_protocol.py +++ b/src/tools/convert_sidre_protocol.py @@ -4,11 +4,12 @@ # # SPDX-License-Identifier: (BSD-3-Clause) """ -Convert a Sidre datastore from the sidre_hdf5 protocol to another protocol. +Convert a Sidre datastore or Conduit Blueprint mesh to another protocol. -Users must supply a path to a sidre_hdf5 rootfile and base name for -the output datastores. Optional command line arguments include -a ``--protocol`` option (the default is ``json``) +Users must supply a path to an input rootfile and base name for the output. +Sidre datastore conversion remains the default. Use ``--input-type blueprint`` +to convert Conduit Blueprint meshes, including HDF5 Blueprint root files. +Optional command line arguments include a ``--protocol`` option (the default is ``json``) and a ``--strip`` option to truncate the array data to at most N elements. The strip option also prepends each array with its original size, the new size and a filler entry of 0 for integer arrays or nan for floating point @@ -31,7 +32,7 @@ import numpy as np import axom.sidre as sidre -VALID_PROTOCOLS = ( +SIDRE_PROTOCOLS = ( "json", "sidre_hdf5", "sidre_conduit_json", @@ -41,20 +42,37 @@ "conduit_json", ) +BLUEPRINT_PROTOCOLS = ( + "hdf5", + "json", + "yaml", + "conduit_bin", + "conduit_hdf5", + "conduit_json", +) + +VALID_PROTOCOLS = tuple(dict.fromkeys(SIDRE_PROTOCOLS + BLUEPRINT_PROTOCOLS)) + def parse_args() -> argparse.Namespace: - parser = argparse.ArgumentParser(description="Sidre protocol converter") + parser = argparse.ArgumentParser(description="Sidre/Blueprint protocol converter") parser.add_argument( "-i", "--input", required=True, - help="Filename of input sidre-hdf5 datastore", + help="Filename of input sidre-hdf5 datastore or Blueprint mesh root file", ) parser.add_argument( "-o", "--output", required=True, - help="Filename of output datastore (without extension)", + help="Filename of output datastore/mesh (without extension)", + ) + parser.add_argument( + "--input-type", + choices=("sidre", "blueprint"), + default="sidre", + help="Type of input file to convert", ) parser.add_argument( "-p", @@ -79,6 +97,25 @@ def parse_args() -> argparse.Namespace: return parser.parse_args() +def initialize_mpi() -> tuple[object | None, bool, int, int]: + if not sidre.AXOM_ENABLE_MPI: + return None, False, 1, 0 + + try: + from mpi4py import MPI + except ImportError as exc: + raise RuntimeError( + "convert_sidre_protocol.py requires mpi4py when Axom is built with MPI support", + ) from exc + + initialized_mpi = False + if not MPI.Is_initialized(): + MPI.Init() + initialized_mpi = True + + return MPI, initialized_mpi, MPI.COMM_WORLD.Get_size(), MPI.COMM_WORLD.Get_rank() + + # # Allocate storage for external data of the input datastore. # @@ -191,25 +228,86 @@ def truncate_bulk_data(group: sidre.Group, max_size: int, verbose: bool) -> None truncate_bulk_data(child, max_size, verbose) -def main() -> int: - args = parse_args() +def blueprint_protocol(protocol: str) -> str: + if protocol == "conduit_hdf5": + return "hdf5" + if protocol in BLUEPRINT_PROTOCOLS: + return protocol - if not sidre.AXOM_ENABLE_MPI: - raise RuntimeError("sidre.IOManager bindings require an MPI-enabled Axom build") + valid = ", ".join(BLUEPRINT_PROTOCOLS) + raise RuntimeError( + f"Protocol '{protocol}' is not valid for Blueprint mesh output. " + f"Use one of: {valid}", ) - try: - from mpi4py import MPI - except ImportError as exc: - raise RuntimeError( - "convert_sidre_protocol.py requires mpi4py when Axom is built with MPI support", - ) from exc - initialized_mpi = False - if not MPI.Is_initialized(): - MPI.Init() - initialized_mpi = True +def blueprint_domain_count(mesh) -> int: + if mesh.has_path("coordsets") and mesh.has_path("topologies"): + return 1 + return mesh.number_of_children() + + +def convert_blueprint_mesh(args: argparse.Namespace, MPI: object | None, comm_size: int, + rank: int) -> int: + if args.strip is not None: + raise RuntimeError("--strip is only supported for Sidre datastore conversion") + + import conduit + import conduit.blueprint + import conduit.relay.io.blueprint + + input_path = Path(args.input) + protocol = blueprint_protocol(args.protocol) + mesh = conduit.Node() + + if MPI is not None and comm_size > 1: + import conduit.blueprint.mpi.mesh + import conduit.relay.mpi.io.blueprint + + comm = MPI.COMM_WORLD.py2f() + if rank == 0: + print( + f"Loading Blueprint mesh from {input_path} on {comm_size} MPI rank(s)", + ) + conduit.relay.mpi.io.blueprint.load_mesh(mesh, str(input_path), comm) + + info = conduit.Node() + valid = conduit.blueprint.mpi.mesh.verify(mesh, info, comm) + local_domains = blueprint_domain_count(mesh) + total_domains = MPI.COMM_WORLD.allreduce(local_domains) + if rank == 0: + print(f"Input Blueprint mesh layout: {total_domains} domain(s)") + if not valid: + raise RuntimeError(f"Input Blueprint mesh failed verification:\n{info.to_yaml()}") + + if rank == 0: + print( + f"Writing out Blueprint mesh in '{protocol}' protocol to file(s) " + f"with base name {args.output}", + ) + conduit.relay.mpi.io.blueprint.save_mesh(mesh, args.output, comm, protocol) + else: + print(f"Loading Blueprint mesh from {input_path}") + conduit.relay.io.blueprint.load_mesh(mesh, str(input_path)) + + info = conduit.Node() + valid = conduit.blueprint.mesh.verify(mesh, info) + domains = blueprint_domain_count(mesh) + print(f"Input Blueprint mesh layout: {domains} domain(s)") + if not valid: + raise RuntimeError(f"Input Blueprint mesh failed verification:\n{info.to_yaml()}") + + print( + f"Writing out Blueprint mesh in '{protocol}' protocol to file(s) " + f"with base name {args.output}", + ) + conduit.relay.io.blueprint.save_mesh(mesh, args.output, protocol) - comm_size = MPI.COMM_WORLD.Get_size() + return 0 + + +def convert_sidre_datastore(args: argparse.Namespace, comm_size: int) -> int: + if not sidre.AXOM_ENABLE_MPI: + raise RuntimeError("sidre.IOManager bindings require an MPI-enabled Axom build") input_path = Path(args.input) manager = sidre.IOManager() @@ -267,11 +365,21 @@ def main() -> int: ) manager.write(root, num_files, args.output, args.protocol) - if initialized_mpi and not MPI.Is_finalized(): - MPI.Finalize() - return 0 +def main() -> int: + args = parse_args() + + MPI, initialized_mpi, comm_size, rank = initialize_mpi() + try: + if args.input_type == "blueprint": + return convert_blueprint_mesh(args, MPI, comm_size, rank) + return convert_sidre_datastore(args, comm_size) + finally: + if MPI is not None and initialized_mpi and not MPI.Is_finalized(): + MPI.Finalize() + + if __name__ == "__main__": sys.exit(main()) From 47174c5bf1d7f0b6d8756a25bb69e83245b833b2 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Fri, 14 Aug 2026 18:15:16 -0700 Subject: [PATCH 11/16] Adds support for truncating data using --strip --- src/tools/CMakeLists.txt | 3 +- src/tools/convert_sidre_protocol.py | 101 ++++++++++++++++++++++------ 2 files changed, 83 insertions(+), 21 deletions(-) diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt index 5c327cdb4c..cbd6a462f3 100644 --- a/src/tools/CMakeLists.txt +++ b/src/tools/CMakeLists.txt @@ -275,11 +275,12 @@ if(NANOBIND_FOUND) --input-type blueprint --output csp_blueprint_output --protocol json + --strip 4 NUM_MPI_TASKS 2 ) set_tests_properties(${_testname} PROPERTIES - PASS_REGULAR_EXPRESSION "Writing out Blueprint mesh") + PASS_REGULAR_EXPRESSION "Truncated [0-9]+ numeric Blueprint array") else() axom_add_python_test( NAME ${_testname} diff --git a/src/tools/convert_sidre_protocol.py b/src/tools/convert_sidre_protocol.py index 116d6a49f7..044eae2e8d 100644 --- a/src/tools/convert_sidre_protocol.py +++ b/src/tools/convert_sidre_protocol.py @@ -17,9 +17,13 @@ and the user passed in ``--strip 3``, the array would be converted to [6, 3, nan, 1.01, 2.02, 3.03]. -The strip option is intended as a temporary solution to truncating -a dataset to allow easier debugging. In the future, the conversion and -truncation/display functionality may be separated into distinct utilities. +For Blueprint meshes, ``--strip`` applies to numeric array leaves with more than +one element and leaves scalar/string metadata intact. +A ``state/Note`` node is added to each output domain to record that the mesh was stripped. + +The resulting stripped output is intended for debugging, not for use as a valid mesh. +In the future, the conversion and truncation/display functionality may be +separated into distinct utilities. """ from __future__ import annotations @@ -246,11 +250,61 @@ def blueprint_domain_count(mesh) -> int: return mesh.number_of_children() +def strip_note(data_kind: str, max_size: int) -> str: + return (f"This {data_kind} was created by axom's 'convert_sidre_protocol' utility " + f"with option '--strip {max_size}'. To simplify debugging, the bulk " + f"data in this {data_kind} has been truncated to have at most {max_size} " + "original values per array. Three values have been prepended to each array: " + "the size of the original array, the number of retained elements and a zero/Nan.") + + +def add_blueprint_strip_note(mesh, note: str) -> None: + if mesh.has_path("coordsets") and mesh.has_path("topologies"): + mesh["state/Note"] = note + return + + for child_idx in range(mesh.number_of_children()): + domain = mesh.child(child_idx) + if domain.has_path("coordsets") and domain.has_path("topologies"): + domain["state/Note"] = note + + +def truncate_conduit_numeric_array(node, max_size: int, verbose: bool) -> int: + dtype = node.dtype() + original_size = dtype.number_of_elements() + if not dtype.is_number() or original_size <= 1: + return 0 + + retained_size = min(max_size, original_size) + values = np.asarray(node.value()).reshape(-1) + retained = values[:retained_size].copy() + + new_values = np.empty(retained_size + 3, dtype=values.dtype) + np.copyto(new_values[:2], np.asarray([original_size, retained_size]), casting="unsafe") + new_values[2] = math.nan if np.issubdtype(new_values.dtype, np.floating) else 0 + if retained_size > 0: + np.copyto(new_values[3:], retained, casting="unsafe") + + if verbose: + print(f"Truncating node {node.path()} from {original_size} to {retained_size}") + + node.reset() + node.set(new_values) + return 1 + + +def truncate_conduit_bulk_data(node, max_size: int, verbose: bool) -> int: + if node.number_of_children() == 0: + return truncate_conduit_numeric_array(node, max_size, verbose) + + truncated_count = 0 + for child_idx in range(node.number_of_children()): + truncated_count += truncate_conduit_bulk_data(node.child(child_idx), max_size, verbose) + return truncated_count + + def convert_blueprint_mesh(args: argparse.Namespace, MPI: object | None, comm_size: int, rank: int) -> int: - if args.strip is not None: - raise RuntimeError("--strip is only supported for Sidre datastore conversion") - import conduit import conduit.blueprint import conduit.relay.io.blueprint @@ -265,9 +319,7 @@ def convert_blueprint_mesh(args: argparse.Namespace, MPI: object | None, comm_si comm = MPI.COMM_WORLD.py2f() if rank == 0: - print( - f"Loading Blueprint mesh from {input_path} on {comm_size} MPI rank(s)", - ) + print(f"Loading Blueprint mesh from {input_path} on {comm_size} MPI rank(s)", ) conduit.relay.mpi.io.blueprint.load_mesh(mesh, str(input_path), comm) info = conduit.Node() @@ -279,11 +331,19 @@ def convert_blueprint_mesh(args: argparse.Namespace, MPI: object | None, comm_si if not valid: raise RuntimeError(f"Input Blueprint mesh failed verification:\n{info.to_yaml()}") + if args.strip is not None: + if rank == 0: + print(f"Truncating numeric Blueprint arrays to at most {args.strip} elements.") + local_truncated = truncate_conduit_bulk_data(mesh, args.strip, args.verbose) + total_truncated = MPI.COMM_WORLD.allreduce(local_truncated) + if rank == 0: + print(f"Truncated {total_truncated} numeric Blueprint array(s).") + add_blueprint_strip_note(mesh, strip_note("Blueprint mesh", args.strip)) + if rank == 0: print( f"Writing out Blueprint mesh in '{protocol}' protocol to file(s) " - f"with base name {args.output}", - ) + f"with base name {args.output}", ) conduit.relay.mpi.io.blueprint.save_mesh(mesh, args.output, comm, protocol) else: print(f"Loading Blueprint mesh from {input_path}") @@ -296,10 +356,15 @@ def convert_blueprint_mesh(args: argparse.Namespace, MPI: object | None, comm_si if not valid: raise RuntimeError(f"Input Blueprint mesh failed verification:\n{info.to_yaml()}") + if args.strip is not None: + print(f"Truncating numeric Blueprint arrays to at most {args.strip} elements.") + truncated_count = truncate_conduit_bulk_data(mesh, args.strip, args.verbose) + print(f"Truncated {truncated_count} numeric Blueprint array(s).") + add_blueprint_strip_note(mesh, strip_note("Blueprint mesh", args.strip)) + print( f"Writing out Blueprint mesh in '{protocol}' protocol to file(s) " - f"with base name {args.output}", - ) + f"with base name {args.output}", ) conduit.relay.io.blueprint.save_mesh(mesh, args.output, protocol) return 0 @@ -352,13 +417,7 @@ def convert_sidre_datastore(args: argparse.Namespace, comm_size: int) -> int: if args.strip is not None: print(f"Truncating views to at most {args.strip} elements.") truncate_bulk_data(root, args.strip, args.verbose) - note = ("This datastore was created by axom's 'convert_sidre_protocol' utility " - f"with option '--strip {args.strip}'. To simplify debugging, the bulk " - f"data in this datastore has been truncated to have at most {args.strip} " - "original values per array. Three values have been prepended to each " - "array: the size of the original array, the number of retained elements " - "and a zero/Nan.") - root.createViewString("Note", note) + root.createViewString("Note", strip_note("datastore", args.strip)) print( f"Writing out datastore in '{args.protocol}' protocol to file(s) with base name {args.output}", @@ -370,6 +429,8 @@ def convert_sidre_datastore(args: argparse.Namespace, comm_size: int) -> int: def main() -> int: args = parse_args() + if args.strip is not None and args.strip < 0: + raise RuntimeError("--strip must be nonnegative") MPI, initialized_mpi, comm_size, rank = initialize_mpi() try: From 26b92bbeb3d3902612231aea433d95d857c56e25 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sun, 16 Aug 2026 16:30:04 -0700 Subject: [PATCH 12/16] Quest: Refactors and improves memory diagnostics --- src/axom/quest/examples/CMakeLists.txt | 19 + ...est_distributed_distance_query_example.cpp | 598 ++++++++++-------- 2 files changed, 346 insertions(+), 271 deletions(-) diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt index c2d2e9a650..b9be49ab51 100644 --- a/src/axom/quest/examples/CMakeLists.txt +++ b/src/axom/quest/examples/CMakeLists.txt @@ -469,6 +469,25 @@ if(AXOM_ENABLE_MPI AND AXOM_ENABLE_SIDRE AND HDF5_FOUND) PASS_REGULAR_EXPRESSION "Analytic verification for 'circle' found 0 errors") endforeach() + if("seq" IN_LIST AXOM_EXECUTION_POLICIES) + set(_test "quest_distributed_closest_point_memory_instrumentation") + axom_add_test( + NAME ${_test} + COMMAND quest_distributed_distance_query_ex + --mesh-file ${quest_data_dir}/mdmesh.2x1.root + --object-mesh-file ${quest_data_dir}/dcp_object_circle.root + --dist-threshold .3 + --policy seq + --track-memory + --trim-after-query + --sample-memory-ms 1 + --distance-file dcp_closest_point_2d_seq_memory_instrumentation + NUM_MPI_TASKS ${_nranks}) + set_tests_properties(${_test} + PROPERTIES + PASS_REGULAR_EXPRESSION "peak RSS during phase") + endif() + unset(_dim) unset(_mesh) unset(_meshes) diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index b74898f458..71bd2b2e34 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -51,6 +51,323 @@ #include // mallinfo2 / mallinfo / malloc_trim #endif +namespace +{ + +constexpr long long INVALID_BYTES = -1; + +struct ProcRss +{ + long long current {INVALID_BYTES}; + long long peak {INVALID_BYTES}; +}; + +struct ReducedBytes +{ + long long maxValue {INVALID_BYTES}; + long long sumValue {INVALID_BYTES}; + int maxRank {-1}; +}; + +struct MemorySnapshot +{ + ProcRss rss; + long long mallocLive {INVALID_BYTES}; + long long mallocArena {INVALID_BYTES}; + long long umpireCurrent {INVALID_BYTES}; + long long umpireHighWatermark {INVALID_BYTES}; +}; + +/// Read current (VmRSS) and peak (VmHWM) resident set size in bytes. +ProcRss readProcRss() +{ + ProcRss result; +#if defined(__linux__) + std::ifstream status("/proc/self/status"); + std::string line; + while(std::getline(status, line)) + { + long long kb = 0; + if(std::sscanf(line.c_str(), "VmRSS: %lld kB", &kb) == 1) + { + result.current = kb * 1024; + } + else if(std::sscanf(line.c_str(), "VmHWM: %lld kB", &kb) == 1) + { + result.peak = kb * 1024; + } + + if(result.current >= 0 && result.peak >= 0) + { + break; + } + } +#endif + return result; +} + +/// Reset VmHWM so the next RSS-peak read reflects only the following phase. +void resetPeakRss() +{ +#if defined(__linux__) + std::ofstream clear("/proc/self/clear_refs"); + if(clear) + { + clear << "5\n"; + } +#endif +} + +std::string humanBytes(long long bytes) +{ + if(bytes < 0) + { + return "n/a"; + } + + const char* units[] = {"B", "KiB", "MiB", "GiB", "TiB"}; + double value = static_cast(bytes); + int unitIdx = 0; + while(value >= 1024.0 && unitIdx < 4) + { + value /= 1024.0; + ++unitIdx; + } + + return axom::fmt::format("{:.2f} {}", value, units[unitIdx]); +} + +/// Reduce a per-rank byte count to the communicator total and hottest rank. +ReducedBytes reduceBytes(long long localValue, MPI_Comm comm, int rank, int commSize) +{ + const bool hasLocalValue = localValue >= 0; + long long localMax = hasLocalValue ? localValue : INVALID_BYTES; + long long localSum = hasLocalValue ? localValue : 0; + int localCount = hasLocalValue ? 1 : 0; + + ReducedBytes result; + MPI_Allreduce(&localMax, &result.maxValue, 1, MPI_LONG_LONG, MPI_MAX, comm); + + int candidateRank = (hasLocalValue && localValue == result.maxValue) ? rank : commSize; + MPI_Allreduce(&candidateRank, &result.maxRank, 1, MPI_INT, MPI_MIN, comm); + + int validCount = 0; + MPI_Allreduce(&localCount, &validCount, 1, MPI_INT, MPI_SUM, comm); + MPI_Allreduce(&localSum, &result.sumValue, 1, MPI_LONG_LONG, MPI_SUM, comm); + + if(validCount == 0) + { + result = ReducedBytes {}; + } + else if(validCount != commSize) + { + result.sumValue = INVALID_BYTES; + } + + return result; +} + +MemorySnapshot takeMemorySnapshot(int umpireAllocatorId) +{ + MemorySnapshot snapshot; + snapshot.rss = readProcRss(); + +#if defined(__GLIBC__) + #if defined(__GLIBC_PREREQ) && __GLIBC_PREREQ(2, 33) + struct mallinfo2 mi = mallinfo2(); // size_t fields: safe above 2 GiB + snapshot.mallocLive = static_cast(mi.uordblks) + static_cast(mi.hblkhd); + snapshot.mallocArena = static_cast(mi.arena); + #else + struct mallinfo mi = mallinfo(); // NOTE: int fields saturate above ~2 GiB + snapshot.mallocLive = static_cast(mi.uordblks) + static_cast(mi.hblkhd); + snapshot.mallocArena = static_cast(mi.arena); + #endif +#endif + +#if defined(AXOM_USE_UMPIRE) + if(umpireAllocatorId >= 0) + { + auto& rm = umpire::ResourceManager::getInstance(); + umpire::Allocator alloc = rm.getAllocator(umpireAllocatorId); + snapshot.umpireCurrent = static_cast(alloc.getCurrentSize()); + snapshot.umpireHighWatermark = static_cast(alloc.getHighWatermark()); + } +#else + AXOM_UNUSED_VAR(umpireAllocatorId); +#endif + + return snapshot; +} + +bool trimMallocArena() +{ +#if defined(__GLIBC__) + ::malloc_trim(0); + return true; +#else + return false; +#endif +} + +/*! + * \brief Opt-in per-run memory probe for the closest-point query. + * + * Reports RSS, peak RSS, glibc live/arena bytes, and Umpire current/high-water bytes when available. + * Values are reduced across the given communicator as a total and a hottest-rank maximum. + * The optional sampler captures transient RSS spikes during the query phase. + */ +class MemoryProbe +{ +public: + MemoryProbe(bool enabled, int sampleMs, MPI_Comm comm, int umpireAllocatorId = -1) + : m_enabled(enabled) + , m_sampleMs(sampleMs) + , m_comm(comm) + , m_umpireAllocatorId(umpireAllocatorId) + { + MPI_Comm_rank(m_comm, &m_rank); + MPI_Comm_size(m_comm, &m_commSize); + } + + ~MemoryProbe() { stopSamplerThread(); } + + void resetPeak() + { + if(m_enabled) + { + resetPeakRss(); + } + } + + void startSampler() + { + if(!m_enabled || m_sampleMs <= 0 || m_samplerThread.joinable()) + { + return; + } + + m_samplerPeak = INVALID_BYTES; + m_stopSampler.store(false, std::memory_order_relaxed); + m_samplerThread = std::thread([this]() { + while(!m_stopSampler.load(std::memory_order_relaxed)) + { + recordSamplerValue(readProcRss().current); + std::this_thread::sleep_for(std::chrono::milliseconds(m_sampleMs)); + } + }); + } + + void stopSampler(const std::string& label) + { + if(!m_enabled || m_sampleMs <= 0 || !m_samplerThread.joinable()) + { + return; + } + + stopSamplerThread(); + recordSamplerValue(readProcRss().current); + + const ReducedBytes sampled = reduceBytes(m_samplerPeak, m_comm, m_rank, m_commSize); + if(m_rank == 0) + { + SLIC_INFO(axom::fmt::format( + "[mem] {}: peak RSS during phase (sampled @ {} ms): max/rank={} (rank {}), total={}", + label, + m_sampleMs, + humanBytes(sampled.maxValue), + sampled.maxRank, + humanBytes(sampled.sumValue))); + } + } + + void report(const std::string& label) + { + if(!m_enabled) + { + return; + } + + const MemorySnapshot snapshot = takeMemorySnapshot(m_umpireAllocatorId); + const ReducedBytes rss = reduceBytes(snapshot.rss.current, m_comm, m_rank, m_commSize); + const ReducedBytes peakRss = reduceBytes(snapshot.rss.peak, m_comm, m_rank, m_commSize); + const ReducedBytes mallocLive = reduceBytes(snapshot.mallocLive, m_comm, m_rank, m_commSize); + const ReducedBytes mallocArena = reduceBytes(snapshot.mallocArena, m_comm, m_rank, m_commSize); + const ReducedBytes umpireCurrent = + reduceBytes(snapshot.umpireCurrent, m_comm, m_rank, m_commSize); + const ReducedBytes umpireHighWatermark = + reduceBytes(snapshot.umpireHighWatermark, m_comm, m_rank, m_commSize); + + if(m_rank == 0) + { + std::string msg = axom::fmt::format( + "[mem] {} (total over {} ranks | max on one rank)\n" + " RSS : {:>11} | {:>11} (rank {})\n" + " RSS peak : {:>11} | {:>11} (rank {})\n" + " malloc live : {:>11} | {:>11} (rank {})\n" + " malloc arena : {:>11} | {:>11} (rank {})", + label, + m_commSize, + humanBytes(rss.sumValue), + humanBytes(rss.maxValue), + rss.maxRank, + humanBytes(peakRss.sumValue), + humanBytes(peakRss.maxValue), + peakRss.maxRank, + humanBytes(mallocLive.sumValue), + humanBytes(mallocLive.maxValue), + mallocLive.maxRank, + humanBytes(mallocArena.sumValue), + humanBytes(mallocArena.maxValue), + mallocArena.maxRank); +#if defined(AXOM_USE_UMPIRE) + if(umpireHighWatermark.maxValue >= 0) + { + msg += axom::fmt::format( + "\n umpire current: {:>11} | {:>11} (rank {})" + "\n umpire hi-water: {:>10} | {:>11} (rank {})", + humanBytes(umpireCurrent.sumValue), + humanBytes(umpireCurrent.maxValue), + umpireCurrent.maxRank, + humanBytes(umpireHighWatermark.sumValue), + humanBytes(umpireHighWatermark.maxValue), + umpireHighWatermark.maxRank); + } +#endif + SLIC_INFO(msg); + } + } + +private: + void stopSamplerThread() + { + if(m_samplerThread.joinable()) + { + m_stopSampler.store(true, std::memory_order_relaxed); + m_samplerThread.join(); + } + } + + void recordSamplerValue(long long bytes) + { + if(bytes > m_samplerPeak) + { + m_samplerPeak = bytes; + } + } + + bool m_enabled {false}; + int m_sampleMs {0}; + MPI_Comm m_comm {MPI_COMM_NULL}; + int m_rank {-1}; + int m_commSize {-1}; + int m_umpireAllocatorId {-1}; + std::atomic m_stopSampler {false}; + std::thread m_samplerThread; + long long m_samplerPeak {INVALID_BYTES}; +}; + +} // namespace + namespace quest = axom::quest; namespace slic = axom::slic; namespace sidre = axom::sidre; @@ -161,6 +478,7 @@ struct Input "If > 0, run a background thread sampling RSS at this interval (ms) " "during the query and report the peak. Useful for observing in-flight " "send-buffer accumulation. Implies --track-memory.") + ->check(axom::CLI::NonNegativeNumber) ->capture_default_str(); app.add_flag("-v,--verbose,!--no-verbose", m_verboseOutput) @@ -1334,270 +1652,6 @@ void make_coords_interleaved(conduit::Node& coordValues) } } -//----------------------------------------------------------------------------- -// Optional memory instrumentation -// -// Purpose: distinguish a genuine leak (memory still referenced after the call) -// from arena retention (memory freed at the allocator but not returned to the OS). -// RSS shows what the process actually holds; glibc "malloc live" shows what is still allocated; -// "malloc arena" shows free memory kept in the arena. -// * RSS after call ~= baseline -> transient, no problem -// * RSS high, malloc-live ~= baseline -> arena retention (malloc_trim should reclaim it) -// * malloc-live high after call -> genuine leak -// Everything here is Linux/glibc-specific and degrades to "n/a" elsewhere. -// All of it is opt-in (--track-memory) and off by default. -//----------------------------------------------------------------------------- - -/// Read current (VmRSS) and peak (VmHWM) resident set size in bytes; -1 if n/a. -inline void readProcRss(long long& rssBytes, long long& peakRssBytes) -{ - rssBytes = -1; - peakRssBytes = -1; -#if defined(__linux__) - std::ifstream status("/proc/self/status"); - std::string line; - while(std::getline(status, line)) - { - long long kb = 0; - if(std::sscanf(line.c_str(), "VmRSS: %lld kB", &kb) == 1) - { - rssBytes = kb * 1024; - } - else if(std::sscanf(line.c_str(), "VmHWM: %lld kB", &kb) == 1) - { - peakRssBytes = kb * 1024; - } - } -#endif -} - -/// Reset the kernel's peak-RSS high-water mark (VmHWM) to the current RSS, so a -/// later VmHWM read reflects the peak of just the intervening phase. -/// Best-effort: requires Linux clear_refs type 5 (kernel >= 4.0). -inline void resetPeakRss() -{ -#if defined(__linux__) - std::ofstream clear("/proc/self/clear_refs"); - if(clear) - { - clear << "5\n"; - } -#endif -} - -/// Human-readable byte count. -inline std::string humanBytes(long long b) -{ - if(b < 0) - { - return "n/a"; - } - const char* units[] = {"B", "KiB", "MiB", "GiB", "TiB"}; - double v = static_cast(b); - int i = 0; - while(v >= 1024.0 && i < 4) - { - v /= 1024.0; - ++i; - } - return axom::fmt::format("{:.2f} {}", v, units[i]); -} - -/// MPI-reduce a per-rank byte count to its max (with the rank achieving it) and -/// its sum across ranks. A negative local value means "unavailable". -inline void reduceBytes(long long local, long long& maxVal, int& maxRank, long long& sumVal) -{ - struct - { - long val; - int rank; - } in {static_cast(local), my_rank}, out {0, 0}; - MPI_Allreduce(&in, &out, 1, MPI_LONG_INT, MPI_MAXLOC, MPI_COMM_WORLD); - maxVal = out.val; - maxRank = out.rank; - MPI_Allreduce(&local, &sumVal, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); -} - -/*! - * \brief Opt-in per-run memory probe for the closest-point query. - * - * Reports RSS, peak RSS, glibc live/arena bytes (and Umpire high-water when built with Umpire), - * reduced across ranks (total and hottest rank). Also offers a background sampler to capture - * the peak RSS *during* a phase (useful for observing in-flight send-buffer accumulation), - * and a resetPeak() to make VmHWM phase-local. - */ -class MemoryProbe -{ -public: - MemoryProbe(bool enabled, int sampleMs, int umpireAllocatorId = -1) - : m_enabled(enabled) - , m_sampleMs(sampleMs) - , m_umpireAllocatorId(umpireAllocatorId) - { } - - bool enabled() const { return m_enabled; } - - /// Reset VmHWM so the next report()'s peak reflects only the next phase. - void resetPeak() - { - if(m_enabled) - { - resetPeakRss(); - } - } - - /// Start a background thread sampling RSS every sampleMs; no-op if disabled or sampleMs <= 0 - /// Records the peak RSS seen until stopSampler(). - void startSampler() - { - if(!m_enabled || m_sampleMs <= 0) - { - return; - } - m_samplerPeak = 0; - m_stopSampler.store(false, std::memory_order_relaxed); - m_samplerThread = std::thread([this]() { - while(!m_stopSampler.load(std::memory_order_relaxed)) - { - long long rss = -1, peak = -1; - readProcRss(rss, peak); - if(rss > m_samplerPeak) - { - m_samplerPeak = rss; - } - std::this_thread::sleep_for(std::chrono::milliseconds(m_sampleMs)); - } - }); - } - - /// Stop the sampler and report the peak RSS during the sampled phase. - void stopSampler(const std::string& label) - { - if(!m_enabled || m_sampleMs <= 0) - { - return; - } - m_stopSampler.store(true, std::memory_order_relaxed); - if(m_samplerThread.joinable()) - { - m_samplerThread.join(); - } - long long rss = -1, peak = -1; // final read to catch the tail - readProcRss(rss, peak); - if(rss > m_samplerPeak) - { - m_samplerPeak = rss; - } - - long long maxVal, sumVal; - int maxRank; - reduceBytes(m_samplerPeak, maxVal, maxRank, sumVal); - if(my_rank == 0) - { - SLIC_INFO(axom::fmt::format( - "[mem] {}: peak RSS during phase (sampled @ {} ms): max/rank={} (rank {}), total={}", - label, - m_sampleMs, - humanBytes(maxVal), - maxRank, - humanBytes(sumVal))); - } - } - - /// Take a snapshot on every rank, reduce it, and print an aggregate (rank 0). - void report(const std::string& label) - { - if(!m_enabled) - { - return; - } - - long long rss = -1, peakRss = -1; - readProcRss(rss, peakRss); - - long long mallocLive = -1, mallocArena = -1; -#if defined(__GLIBC__) - #if defined(__GLIBC_PREREQ) && __GLIBC_PREREQ(2, 33) - struct mallinfo2 mi = mallinfo2(); // size_t fields: safe above 2 GiB - mallocLive = static_cast(mi.uordblks) + static_cast(mi.hblkhd); - mallocArena = static_cast(mi.arena); - #else - struct mallinfo mi = mallinfo(); // NOTE: int fields saturate above ~2 GiB - mallocLive = static_cast(mi.uordblks) + static_cast(mi.hblkhd); - mallocArena = static_cast(mi.arena); - #endif -#endif - - long long umpireCur = -1, umpireHwm = -1; -#if defined(AXOM_USE_UMPIRE) - if(m_umpireAllocatorId >= 0) - { - auto& rm = umpire::ResourceManager::getInstance(); - umpire::Allocator alloc = rm.getAllocator(m_umpireAllocatorId); - umpireCur = static_cast(alloc.getCurrentSize()); - umpireHwm = static_cast(alloc.getHighWatermark()); - } -#endif - - long long rssMax, rssSum, peakMax, peakSum, liveMax, liveSum, arenaMax, arenaSum; - long long umpCurMax, umpCurSum, umpHwmMax, umpHwmSum; - int rssMaxRank, peakMaxRank, liveMaxRank, arenaMaxRank, umpCurMaxRank, umpHwmMaxRank; - reduceBytes(rss, rssMax, rssMaxRank, rssSum); - reduceBytes(peakRss, peakMax, peakMaxRank, peakSum); - reduceBytes(mallocLive, liveMax, liveMaxRank, liveSum); - reduceBytes(mallocArena, arenaMax, arenaMaxRank, arenaSum); - reduceBytes(umpireCur, umpCurMax, umpCurMaxRank, umpCurSum); - reduceBytes(umpireHwm, umpHwmMax, umpHwmMaxRank, umpHwmSum); - - if(my_rank == 0) - { - std::string msg = axom::fmt::format( - "[mem] {} (total over {} ranks | max on one rank)\n" - " RSS : {:>11} | {:>11} (rank {})\n" - " RSS peak : {:>11} | {:>11} (rank {})\n" - " malloc live : {:>11} | {:>11} (rank {})\n" - " malloc arena : {:>11} | {:>11} (rank {})", - label, - num_ranks, - humanBytes(rssSum), - humanBytes(rssMax), - rssMaxRank, - humanBytes(peakSum), - humanBytes(peakMax), - peakMaxRank, - humanBytes(liveSum), - humanBytes(liveMax), - liveMaxRank, - humanBytes(arenaSum), - humanBytes(arenaMax), - arenaMaxRank); -#if defined(AXOM_USE_UMPIRE) - if(umpHwmMax >= 0) - { - msg += axom::fmt::format( - "\n umpire current: {:>11} | {:>11} (rank {})" - "\n umpire hi-water: {:>10} | {:>11} (rank {})", - humanBytes(umpCurSum), - humanBytes(umpCurMax), - umpCurMaxRank, - humanBytes(umpHwmSum), - humanBytes(umpHwmMax), - umpHwmMaxRank); - } -#endif - SLIC_INFO(msg); - } - } - -private: - bool m_enabled {false}; - int m_sampleMs {0}; - int m_umpireAllocatorId {-1}; - std::atomic m_stopSampler {false}; - std::thread m_samplerThread; - long long m_samplerPeak {0}; -}; - /// Utility function to initialize the logger void initializeLogger() { @@ -1778,7 +1832,7 @@ int main(int argc, char** argv) memUmpireId = umpireAllocator.getId(); #endif const bool trackMem = params.trackMemory || params.trimAfterQuery || params.sampleMemoryMs > 0; - MemoryProbe memProbe(trackMem, params.sampleMemoryMs, memUmpireId); + MemoryProbe memProbe(trackMem, params.sampleMemoryMs, MPI_COMM_WORLD, memUmpireId); memProbe.report("baseline (meshes read, before BVH)"); // Build the spatial index over the object on each rank @@ -1810,12 +1864,14 @@ int main(int argc, char** argv) // (This is a diagnostic; the library-side fix would trim inside computeClosestPoints after releasing transfer buffers.) if(params.trimAfterQuery) { -#if defined(__GLIBC__) - ::malloc_trim(0); - memProbe.report("after malloc_trim(0)"); -#else - SLIC_WARNING("--trim-after-query requested but malloc_trim is glibc-only; skipping."); -#endif + if(trimMallocArena()) + { + memProbe.report("after malloc_trim(0)"); + } + else + { + SLIC_WARNING("--trim-after-query requested but malloc_trim is glibc-only; skipping."); + } } auto getDoubleMinMax = [](double inVal, double& minVal, double& maxVal, double& sumVal) { From 8ad04aa5f711e906510c2e70bc97bdbbd70ab895 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sun, 16 Aug 2026 16:49:17 -0700 Subject: [PATCH 13/16] Quest: Adds allReduce helper to dcp driver --- ...est_distributed_distance_query_example.cpp | 116 +++++++++--------- 1 file changed, 61 insertions(+), 55 deletions(-) diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index 71bd2b2e34..78d5d8fa52 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -41,7 +41,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -54,28 +56,47 @@ namespace { -constexpr long long INVALID_BYTES = -1; +using ByteCount = std::int64_t; + +constexpr ByteCount INVALID_BYTES = -1; +constexpr ByteCount BYTES_PER_KIB = 1024; + +template +T allReduce(T localValue, MPI_Op op, MPI_Comm comm) +{ + T result {}; + MPI_Allreduce(&localValue, &result, 1, axom::mpi_traits::type, op, comm); + return result; +} + +template +void allReduceMinMaxSum(T localValue, T& minValue, T& maxValue, T& sumValue, MPI_Comm comm) +{ + minValue = allReduce(localValue, MPI_MIN, comm); + maxValue = allReduce(localValue, MPI_MAX, comm); + sumValue = allReduce(localValue, MPI_SUM, comm); +} struct ProcRss { - long long current {INVALID_BYTES}; - long long peak {INVALID_BYTES}; + ByteCount current {INVALID_BYTES}; + ByteCount peak {INVALID_BYTES}; }; struct ReducedBytes { - long long maxValue {INVALID_BYTES}; - long long sumValue {INVALID_BYTES}; + ByteCount maxValue {INVALID_BYTES}; + ByteCount sumValue {INVALID_BYTES}; int maxRank {-1}; }; struct MemorySnapshot { ProcRss rss; - long long mallocLive {INVALID_BYTES}; - long long mallocArena {INVALID_BYTES}; - long long umpireCurrent {INVALID_BYTES}; - long long umpireHighWatermark {INVALID_BYTES}; + ByteCount mallocLive {INVALID_BYTES}; + ByteCount mallocArena {INVALID_BYTES}; + ByteCount umpireCurrent {INVALID_BYTES}; + ByteCount umpireHighWatermark {INVALID_BYTES}; }; /// Read current (VmRSS) and peak (VmHWM) resident set size in bytes. @@ -87,14 +108,14 @@ ProcRss readProcRss() std::string line; while(std::getline(status, line)) { - long long kb = 0; - if(std::sscanf(line.c_str(), "VmRSS: %lld kB", &kb) == 1) + ByteCount kb = 0; + if(std::sscanf(line.c_str(), "VmRSS: %" SCNd64 " kB", &kb) == 1) { - result.current = kb * 1024; + result.current = kb * BYTES_PER_KIB; } - else if(std::sscanf(line.c_str(), "VmHWM: %lld kB", &kb) == 1) + else if(std::sscanf(line.c_str(), "VmHWM: %" SCNd64 " kB", &kb) == 1) { - result.peak = kb * 1024; + result.peak = kb * BYTES_PER_KIB; } if(result.current >= 0 && result.peak >= 0) @@ -118,7 +139,7 @@ void resetPeakRss() #endif } -std::string humanBytes(long long bytes) +std::string humanBytes(ByteCount bytes) { if(bytes < 0) { @@ -138,22 +159,21 @@ std::string humanBytes(long long bytes) } /// Reduce a per-rank byte count to the communicator total and hottest rank. -ReducedBytes reduceBytes(long long localValue, MPI_Comm comm, int rank, int commSize) +ReducedBytes reduceBytes(ByteCount localValue, MPI_Comm comm, int rank, int commSize) { const bool hasLocalValue = localValue >= 0; - long long localMax = hasLocalValue ? localValue : INVALID_BYTES; - long long localSum = hasLocalValue ? localValue : 0; + ByteCount localMax = hasLocalValue ? localValue : INVALID_BYTES; + ByteCount localSum = hasLocalValue ? localValue : 0; int localCount = hasLocalValue ? 1 : 0; ReducedBytes result; - MPI_Allreduce(&localMax, &result.maxValue, 1, MPI_LONG_LONG, MPI_MAX, comm); + result.maxValue = allReduce(localMax, MPI_MAX, comm); int candidateRank = (hasLocalValue && localValue == result.maxValue) ? rank : commSize; - MPI_Allreduce(&candidateRank, &result.maxRank, 1, MPI_INT, MPI_MIN, comm); + result.maxRank = allReduce(candidateRank, MPI_MIN, comm); - int validCount = 0; - MPI_Allreduce(&localCount, &validCount, 1, MPI_INT, MPI_SUM, comm); - MPI_Allreduce(&localSum, &result.sumValue, 1, MPI_LONG_LONG, MPI_SUM, comm); + const int validCount = allReduce(localCount, MPI_SUM, comm); + result.sumValue = allReduce(localSum, MPI_SUM, comm); if(validCount == 0) { @@ -175,12 +195,12 @@ MemorySnapshot takeMemorySnapshot(int umpireAllocatorId) #if defined(__GLIBC__) #if defined(__GLIBC_PREREQ) && __GLIBC_PREREQ(2, 33) struct mallinfo2 mi = mallinfo2(); // size_t fields: safe above 2 GiB - snapshot.mallocLive = static_cast(mi.uordblks) + static_cast(mi.hblkhd); - snapshot.mallocArena = static_cast(mi.arena); + snapshot.mallocLive = static_cast(mi.uordblks) + static_cast(mi.hblkhd); + snapshot.mallocArena = static_cast(mi.arena); #else struct mallinfo mi = mallinfo(); // NOTE: int fields saturate above ~2 GiB - snapshot.mallocLive = static_cast(mi.uordblks) + static_cast(mi.hblkhd); - snapshot.mallocArena = static_cast(mi.arena); + snapshot.mallocLive = static_cast(mi.uordblks) + static_cast(mi.hblkhd); + snapshot.mallocArena = static_cast(mi.arena); #endif #endif @@ -189,8 +209,8 @@ MemorySnapshot takeMemorySnapshot(int umpireAllocatorId) { auto& rm = umpire::ResourceManager::getInstance(); umpire::Allocator alloc = rm.getAllocator(umpireAllocatorId); - snapshot.umpireCurrent = static_cast(alloc.getCurrentSize()); - snapshot.umpireHighWatermark = static_cast(alloc.getHighWatermark()); + snapshot.umpireCurrent = static_cast(alloc.getCurrentSize()); + snapshot.umpireHighWatermark = static_cast(alloc.getHighWatermark()); } #else AXOM_UNUSED_VAR(umpireAllocatorId); @@ -347,7 +367,7 @@ class MemoryProbe } } - void recordSamplerValue(long long bytes) + void recordSamplerValue(ByteCount bytes) { if(bytes > m_samplerPeak) { @@ -363,7 +383,7 @@ class MemoryProbe int m_umpireAllocatorId {-1}; std::atomic m_stopSampler {false}; std::thread m_samplerThread; - long long m_samplerPeak {INVALID_BYTES}; + ByteCount m_samplerPeak {INVALID_BYTES}; }; } // namespace @@ -649,7 +669,7 @@ struct BlueprintParticleMesh } int verificationRank = m_hasVerification ? m_rank : m_nranks; - MPI_Allreduce(MPI_IN_PLACE, &verificationRank, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); + verificationRank = allReduce(verificationRank, MPI_MIN, MPI_COMM_WORLD); if(verificationRank < m_nranks) { std::string description = m_rank == verificationRank ? m_description : std::string {}; @@ -693,7 +713,7 @@ struct BlueprintParticleMesh m_dimension = conduit::blueprint::mesh::coordset::dims(coordsetNode); } - MPI_Allreduce(MPI_IN_PLACE, &m_dimension, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD); + m_dimension = allReduce(m_dimension, MPI_MAX, MPI_COMM_WORLD); SLIC_ASSERT(m_dimension > 0); if(domCount > 0) @@ -823,16 +843,10 @@ struct BlueprintParticleMesh numPoints(), domain_count())); - auto getIntMinMax = [](int inVal, int& minVal, int& maxVal, int& sumVal) { - MPI_Allreduce(&inVal, &minVal, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); - MPI_Allreduce(&inVal, &maxVal, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD); - MPI_Allreduce(&inVal, &sumVal, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); - }; - // Output some global mesh size stats { int mn, mx, sum; - getIntMinMax(numPoints(), mn, mx, sum); + allReduceMinMaxSum(numPoints(), mn, mx, sum, MPI_COMM_WORLD); SLIC_INFO(axom::fmt::format("{} has {{min:{}, max:{}, sum:{}, avg:{}}} points", meshLabel, mn, @@ -842,7 +856,7 @@ struct BlueprintParticleMesh } { int mn, mx, sum; - getIntMinMax(domain_count(), mn, mx, sum); + allReduceMinMaxSum(static_cast(domain_count()), mn, mx, sum, MPI_COMM_WORLD); SLIC_INFO(axom::fmt::format("{} has {{min:{}, max:{}, sum:{}, avg:{}}} domains", meshLabel, mn, @@ -1618,10 +1632,8 @@ int verifyAnalyticClosestPoints(BlueprintParticleMesh& queryMesh, } } - int globalErrCount = 0; - MPI_Allreduce(&localErrCount, &globalErrCount, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); - int globalLogCount = 0; - MPI_Allreduce(&localLogCount, &globalLogCount, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + int globalErrCount = allReduce(localErrCount, MPI_SUM, MPI_COMM_WORLD); + int globalLogCount = allReduce(localLogCount, MPI_SUM, MPI_COMM_WORLD); SLIC_INFO( axom::fmt::format("Analytic verification for '{}' found {} errors " "({} detailed messages{}).", @@ -1792,8 +1804,8 @@ int main(int argc, char** argv) // Initialize spatial index for object points, and run query //--------------------------------------------------------------------------- - int globalObjectPointCount = objectMeshWrapper.getParticleMesh().numPoints(); - MPI_Allreduce(MPI_IN_PLACE, &globalObjectPointCount, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + int globalObjectPointCount = + allReduce(objectMeshWrapper.getParticleMesh().numPoints(), MPI_SUM, MPI_COMM_WORLD); auto init_str = banner(axom::fmt::format("Initializing BVH tree over {} object points", globalObjectPointCount)); @@ -1874,19 +1886,13 @@ int main(int argc, char** argv) } } - auto getDoubleMinMax = [](double inVal, double& minVal, double& maxVal, double& sumVal) { - MPI_Allreduce(&inVal, &minVal, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); - MPI_Allreduce(&inVal, &maxVal, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - MPI_Allreduce(&inVal, &sumVal, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - }; - // Output some timing stats { double minInit, maxInit, sumInit; - getDoubleMinMax(initTimer.elapsedTimeInSec(), minInit, maxInit, sumInit); + allReduceMinMaxSum(initTimer.elapsedTimeInSec(), minInit, maxInit, sumInit, MPI_COMM_WORLD); double minQuery, maxQuery, sumQuery; - getDoubleMinMax(queryTimer.elapsedTimeInSec(), minQuery, maxQuery, sumQuery); + allReduceMinMaxSum(queryTimer.elapsedTimeInSec(), minQuery, maxQuery, sumQuery, MPI_COMM_WORLD); SLIC_INFO( axom::fmt::format("Initialization with policy {} took {{avg:{}, min:{}, max:{}}} seconds", From 2b03546833dd62b4d7ae3a27c4b8e2597bdc22e3 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sun, 16 Aug 2026 19:53:57 -0700 Subject: [PATCH 14/16] Minor fixups * In convert_sidre_protocol, run w/ correct protocol * In new point mesh generator script, ensure MPI/conduit are available when we expect them * Fix check for multi-domain runs --- ...uest_distributed_distance_query_example.cpp | 6 ++---- src/tools/convert_sidre_protocol.py | 18 ++++++++++++++---- src/tools/gen-multidom-point-mesh.py | 7 +++++-- 3 files changed, 21 insertions(+), 10 deletions(-) diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index 78d5d8fa52..7d764c099b 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -1184,11 +1184,9 @@ class QueryMeshWrapper { sidre::Group* dstDomains = m_queryMesh.root_group(); bool isMultidomain = conduit::blueprint::mesh::is_multi_domain(node); - if(!isMultidomain) - { - SLIC_ASSERT(!isMultidomain || dstDomains->getNumGroups() == node.number_of_children()); - } const int domainCount = dstDomains->getNumGroups(); + const int srcDomainCount = static_cast(conduit::blueprint::mesh::number_of_domains(node)); + SLIC_ASSERT(domainCount == srcDomainCount); for(int d = 0; d < domainCount; ++d) { sidre::Group& domGroup = *dstDomains->getGroup(d); diff --git a/src/tools/convert_sidre_protocol.py b/src/tools/convert_sidre_protocol.py index 044eae2e8d..2085a68d18 100644 --- a/src/tools/convert_sidre_protocol.py +++ b/src/tools/convert_sidre_protocol.py @@ -83,7 +83,7 @@ def parse_args() -> argparse.Namespace: "--protocol", default="json", choices=VALID_PROTOCOLS, - help="Desired protocol for output datastore", + help="Desired output protocol; valid protocols depend on --input-type", ) parser.add_argument( "-s", @@ -244,6 +244,16 @@ def blueprint_protocol(protocol: str) -> str: f"Use one of: {valid}", ) +def sidre_protocol(protocol: str) -> str: + if protocol in SIDRE_PROTOCOLS: + return protocol + + valid = ", ".join(SIDRE_PROTOCOLS) + raise RuntimeError( + f"Protocol '{protocol}' is not valid for Sidre datastore output. " + f"Use one of: {valid}", ) + + def blueprint_domain_count(mesh) -> int: if mesh.has_path("coordsets") and mesh.has_path("topologies"): return 1 @@ -374,6 +384,7 @@ def convert_sidre_datastore(args: argparse.Namespace, comm_size: int) -> int: if not sidre.AXOM_ENABLE_MPI: raise RuntimeError("sidre.IOManager bindings require an MPI-enabled Axom build") + protocol = sidre_protocol(args.protocol) input_path = Path(args.input) manager = sidre.IOManager() datastore = sidre.DataStore() @@ -420,9 +431,8 @@ def convert_sidre_datastore(args: argparse.Namespace, comm_size: int) -> int: root.createViewString("Note", strip_note("datastore", args.strip)) print( - f"Writing out datastore in '{args.protocol}' protocol to file(s) with base name {args.output}", - ) - manager.write(root, num_files, args.output, args.protocol) + f"Writing out datastore in '{protocol}' protocol to file(s) with base name {args.output}", ) + manager.write(root, num_files, args.output, protocol) return 0 diff --git a/src/tools/gen-multidom-point-mesh.py b/src/tools/gen-multidom-point-mesh.py index 73f0b250de..8423116234 100755 --- a/src/tools/gen-multidom-point-mesh.py +++ b/src/tools/gen-multidom-point-mesh.py @@ -132,8 +132,11 @@ def get_mpi_context(): import conduit.relay.mpi import conduit.relay.mpi.io import conduit.relay.mpi.io.blueprint - except ModuleNotFoundError: - return {'enabled': False, 'comm': None, 'rank': 0, 'size': 1} + except ModuleNotFoundError as exc: + raise RuntimeError('MPI launch detected, but mpi4py and Conduit MPI relay Python modules ' + 'are not available. Run this generator serially, or use Axom\'s ' + 'run_python_with_axom.sh from an MPI-enabled build with Conduit MPI ' + 'Python support.') from exc comm = MPI.COMM_WORLD.py2f() return { From 764538b27231354c38dfc4c36bb869ab9cf308d5 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Tue, 18 Aug 2026 21:06:59 -0700 Subject: [PATCH 15/16] Bugfix for configs that do not have a 64-bit int type --- ...est_distributed_distance_query_example.cpp | 26 ++++++++++++------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index 7d764c099b..9ad7e8cc10 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -12,7 +12,6 @@ // Axom includes #include "axom/config.hpp" #include "axom/core.hpp" -#include "axom/core/NumericLimits.hpp" #include "axom/slic.hpp" #include "axom/primal.hpp" #include "axom/sidre.hpp" @@ -56,9 +55,16 @@ namespace { +#if defined(AXOM_NO_INT64_T) +using ByteCount = std::uint32_t; +constexpr ByteCount INVALID_BYTES = axom::numeric_limits::max(); + #define AXOM_DCP_SCN_BYTE_COUNT SCNu32 +#else using ByteCount = std::int64_t; - constexpr ByteCount INVALID_BYTES = -1; + #define AXOM_DCP_SCN_BYTE_COUNT SCNd64 +#endif + constexpr ByteCount BYTES_PER_KIB = 1024; template @@ -109,16 +115,16 @@ ProcRss readProcRss() while(std::getline(status, line)) { ByteCount kb = 0; - if(std::sscanf(line.c_str(), "VmRSS: %" SCNd64 " kB", &kb) == 1) + if(std::sscanf(line.c_str(), "VmRSS: %" AXOM_DCP_SCN_BYTE_COUNT " kB", &kb) == 1) { result.current = kb * BYTES_PER_KIB; } - else if(std::sscanf(line.c_str(), "VmHWM: %" SCNd64 " kB", &kb) == 1) + else if(std::sscanf(line.c_str(), "VmHWM: %" AXOM_DCP_SCN_BYTE_COUNT " kB", &kb) == 1) { result.peak = kb * BYTES_PER_KIB; } - if(result.current >= 0 && result.peak >= 0) + if(result.current != INVALID_BYTES && result.peak != INVALID_BYTES) { break; } @@ -141,7 +147,7 @@ void resetPeakRss() std::string humanBytes(ByteCount bytes) { - if(bytes < 0) + if(bytes == INVALID_BYTES) { return "n/a"; } @@ -161,8 +167,8 @@ std::string humanBytes(ByteCount bytes) /// Reduce a per-rank byte count to the communicator total and hottest rank. ReducedBytes reduceBytes(ByteCount localValue, MPI_Comm comm, int rank, int commSize) { - const bool hasLocalValue = localValue >= 0; - ByteCount localMax = hasLocalValue ? localValue : INVALID_BYTES; + const bool hasLocalValue = localValue != INVALID_BYTES; + ByteCount localMax = hasLocalValue ? localValue : 0; ByteCount localSum = hasLocalValue ? localValue : 0; int localCount = hasLocalValue ? 1 : 0; @@ -369,7 +375,7 @@ class MemoryProbe void recordSamplerValue(ByteCount bytes) { - if(bytes > m_samplerPeak) + if(bytes != INVALID_BYTES && (m_samplerPeak == INVALID_BYTES || bytes > m_samplerPeak)) { m_samplerPeak = bytes; } @@ -388,6 +394,8 @@ class MemoryProbe } // namespace +#undef AXOM_DCP_SCN_BYTE_COUNT + namespace quest = axom::quest; namespace slic = axom::slic; namespace sidre = axom::sidre; From f97c4769f327e61c31d12aed80adf2d8ee485b2c Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sun, 16 Aug 2026 16:57:42 -0700 Subject: [PATCH 16/16] Updates data submodule to include dcp point meshes --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 8ac544afdc..c446495931 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 8ac544afdc0d75e9cfe0681f9eaa8f2150534dea +Subproject commit c446495931576ffda8017633f683118f791c66f0