diff --git a/data b/data index 8ac544afdc..c446495931 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 8ac544afdc0d75e9cfe0681f9eaa8f2150534dea +Subproject commit c446495931576ffda8017633f683118f791c66f0 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 +``` 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/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt index 5261e39afd..b9be49ab51 100644 --- a/src/axom/quest/examples/CMakeLists.txt +++ b/src/axom/quest/examples/CMakeLists.txt @@ -381,54 +381,50 @@ if(AXOM_ENABLE_MPI AND AXOM_ENABLE_SIDRE AND HDF5_FOUND) if(AXOM_ENABLE_TESTS AND AXOM_DATA_DIR) set(_nranks 3) - # Run the distributed closest point example on N ranks for each enabled policy - # Non-zero empty-rank probability tests domain underloading case + # Run the distributed closest point example on N ranks for each enabled policy. + # 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") - # 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}) + 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(_shape circle) + set(_object_mesh ${quest_data_dir}/dcp_object_circle.root) + elseif(_ndim EQUAL 3) + set(_dim 3) + set(_shape sphere) + set(_object_mesh ${quest_data_dir}/dcp_object_sphere.root) 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) - - 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 + --object-mesh-file ${_object_mesh} --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 + PASS_REGULAR_EXPRESSION "Analytic verification for '${_shape}' found 0 errors") - if(_pol STREQUAL "seq" AND (_mesh STREQUAL "mdmesh.2x1" OR _mesh STREQUAL "mdmesh.2x2x1")) + if(_pol STREQUAL "seq" AND + (_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") @@ -436,26 +432,72 @@ if(AXOM_ENABLE_MPI AND AXOM_ENABLE_SIDRE AND HDF5_FOUND) 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 + --object-mesh-file ${_object_mesh} --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 + PASS_REGULAR_EXPRESSION "Analytic verification for '${_shape}' found 0 errors") endif() endforeach() endforeach() - unset(optional_dependency) + 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() + + 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) + unset(_ndim) + unset(_num_threads) + unset(_object_mesh) + unset(_pol) unset(_nranks) + unset(_shape) + unset(_sizes) unset(_test) unset(_static_test) endif() 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..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" @@ -36,19 +35,372 @@ #include "mpi.h" // C/C++ includes +#include #include -#include +#include #include #include +#include +#include +#include +#include +#include +#include +#include +#include +#if defined(__GLIBC__) + #include // mallinfo2 / mallinfo / malloc_trim +#endif + +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 +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 +{ + ByteCount current {INVALID_BYTES}; + ByteCount peak {INVALID_BYTES}; +}; + +struct ReducedBytes +{ + ByteCount maxValue {INVALID_BYTES}; + ByteCount sumValue {INVALID_BYTES}; + int maxRank {-1}; +}; + +struct MemorySnapshot +{ + ProcRss rss; + 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. +ProcRss readProcRss() +{ + ProcRss result; +#if defined(__linux__) + std::ifstream status("/proc/self/status"); + std::string line; + while(std::getline(status, line)) + { + ByteCount kb = 0; + 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: %" AXOM_DCP_SCN_BYTE_COUNT " kB", &kb) == 1) + { + result.peak = kb * BYTES_PER_KIB; + } + + if(result.current != INVALID_BYTES && result.peak != INVALID_BYTES) + { + 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(ByteCount bytes) +{ + if(bytes == INVALID_BYTES) + { + 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(ByteCount localValue, MPI_Comm comm, int rank, int commSize) +{ + const bool hasLocalValue = localValue != INVALID_BYTES; + ByteCount localMax = hasLocalValue ? localValue : 0; + ByteCount localSum = hasLocalValue ? localValue : 0; + int localCount = hasLocalValue ? 1 : 0; + + ReducedBytes result; + result.maxValue = allReduce(localMax, MPI_MAX, comm); + + int candidateRank = (hasLocalValue && localValue == result.maxValue) ? rank : commSize; + result.maxRank = allReduce(candidateRank, MPI_MIN, comm); + + const int validCount = allReduce(localCount, MPI_SUM, comm); + result.sumValue = allReduce(localSum, 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(ByteCount bytes) + { + if(bytes != INVALID_BYTES && (m_samplerPeak == INVALID_BYTES || 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; + ByteCount m_samplerPeak {INVALID_BYTES}; +}; + +} // namespace + +#undef AXOM_DCP_SCN_BYTE_COUNT 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; @@ -59,21 +411,29 @@ 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 { public: std::string meshFile; 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}; - - // Latitudinal direction of object mesh - std::vector latRange {-90.0, 90.0}; - int latPointCount {20}; + // 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 RuntimePolicy policy {RuntimePolicy::seq}; @@ -81,12 +441,6 @@ struct Input bool dynamicDistanceFiltering {true}; - bool checkResults {false}; - - bool randomSpacing {true}; - - std::vector objDomainCountRange {1, 1}; - private: bool m_verboseOutput {false}; @@ -127,40 +481,36 @@ 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. " + "Generate this mesh with src/tools/gen-multidom-point-mesh.py.") + ->check(axom::CLI::ExistingFile) + ->required(); - app.add_flag("-v,--verbose,!--no-verbose", m_verboseOutput) - ->description("Enable/disable verbose output") + 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_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") + 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(); - object_options->add_option("-n,--long-point-count", longPointCount) - ->description("Number of points around the longitudinal direction") + 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.") + ->check(axom::CLI::NonNegativeNumber) ->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)") + app.add_flag("-v,--verbose,!--no-verbose", m_verboseOutput) + ->description("Enable/disable verbose output") ->capture_default_str(); app.add_option("-d,--dist-threshold", distThreshold) @@ -179,10 +529,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 @@ -250,7 +596,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; } @@ -286,6 +633,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. @@ -300,21 +650,49 @@ 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) + m_hasVerification = false; + m_verification.reset(); + m_description.clear(); + for(conduit::index_t d = 0; d < domCount; ++d) { - m_coordsAreStrided = - mdMesh[0].fetch_existing("topologies/mesh/elements/dims").has_child("strides"); - if(m_coordsAreStrided) + const conduit::Node& domain = mdMesh.child(d); + if(m_description.empty() && domain.has_path("state/description")) { - SLIC_WARNING( - axom::fmt::format("Mesh '{}' is strided. Stride support is under development.", - meshFilename)); + 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; + verificationRank = allReduce(verificationRank, 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()) @@ -324,13 +702,26 @@ struct BlueprintParticleMesh } 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].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)); + } + m_coordsetName = mdMesh[0].fetch_existing(topologyPath + "/coordset").as_string(); const conduit::Node coordsetNode = mdMesh[0].fetch_existing("coordsets").fetch_existing(m_coordsetName); 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) @@ -460,16 +851,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, @@ -479,7 +864,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, @@ -527,6 +912,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() @@ -554,13 +953,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 @@ -574,7 +973,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(); @@ -595,11 +994,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(); } @@ -634,8 +1033,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 @@ -699,8 +1098,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; @@ -716,9 +1118,14 @@ struct BlueprintParticleMesh class ObjectMeshWrapper { public: - ObjectMeshWrapper(sidre::Group* group) : m_objectMesh(group, "mesh", "coords") + //!@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; } @@ -728,20 +1135,12 @@ 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); - } + 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; - bool m_verbose {false}; }; class QueryMeshWrapper @@ -793,17 +1192,16 @@ 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); 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")) @@ -849,10 +1247,17 @@ 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(); + 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) + { + continue; + } + + conduit::float64_array dst = dstView->getArray(); + const conduit::float64_array src = srcComponent.value(); for(int i = 0; i < nPts; ++i) { dst[i] = src[i]; @@ -862,308 +1267,387 @@ class QueryMeshWrapper } } - /** - * 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) +private: + BlueprintParticleMesh m_queryMesh; +}; + +//--------------------------------------------------------------------------- +// Transform closest points to distances and directions +//--------------------------------------------------------------------------- +template +void computeDistancesAndDirections(BlueprintParticleMesh& queryMesh, + const std::string& cpCoordsField, + const std::string& cpIndexField, + const std::string& distanceField, + const std::string& directionField) +{ + SLIC_ASSERT(queryMesh.dimension() == DIM); + + using primal::squared_distance; + using PointType = primal::Point; + using IndexSet = slam::PositionSet<>; + + PointType nowhere(axom::numeric_limits::signaling_NaN()); + const double nodist = axom::numeric_limits::signaling_NaN(); + + queryMesh.registerNodalScalarField(distanceField); + queryMesh.registerNodalVectorField(directionField); + for(axom::IndexType di = 0; di < queryMesh.domain_count(); ++di) { - using PointType = axom::primal::Point; + auto cpCoords = queryMesh.getNodalVectorField(cpCoordsField, di); - m_queryMesh.registerNodalScalarField("error_flag"); + auto cpIndices = queryMesh.getNodalScalarField(cpIndexField, di); - int sumErrCount = 0; - int sumWarningCount = 0; - for(axom::IndexType dIdx = 0; dIdx < m_queryMesh.domain_count(); ++dIdx) + axom::Array qPts = queryMesh.getVertexPositions(di); + axom::ArrayView distances = queryMesh.getNodalScalarField("distance", di); + axom::ArrayView directions = queryMesh.getNodalVectorField("direction", di); + axom::IndexType ptCount = queryMesh.numPoints(di); + for(auto ptIdx : IndexSet(ptCount)) { - auto queryPts = m_queryMesh.getPoints(dIdx); - - axom::ArrayView cpCoords = - m_queryMesh.getNodalVectorField("cp_coords", dIdx); - SLIC_INFO(axom::fmt::format("Closest points ({}):", cpCoords.size())); + const bool has_cp = cpIndices[ptIdx] >= 0; + const PointType& cp = has_cp ? cpCoords[ptIdx] : nowhere; + const PointType& qPt = has_cp ? PointType(&qPts[ptIdx][0]) : nowhere; + distances[ptIdx] = has_cp ? sqrt(squared_distance(qPt, cp)) : nodist; + directions[ptIdx] = PointType(has_cp ? (cp - qPt).array() : nowhere.array()); + } + } +} - axom::ArrayView cpIndices = - m_queryMesh.getNodalScalarField("cp_index", dIdx); +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; + } - axom::ArrayView errorFlag = - m_queryMesh.getNodalScalarField("error_flag", dIdx); + 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; +} - SLIC_ASSERT(queryPts.size() == cpCoords.size()); - SLIC_ASSERT(queryPts.size() == cpIndices.size()); +double conduitDouble(const conduit::Node& node, const std::string& path, double defaultValue) +{ + return node.has_path(path) ? node.fetch_existing(path).to_double() : defaultValue; +} - if(params.isVerbose()) - { - SLIC_INFO(axom::fmt::format("Closest points ({}):", cpCoords.size())); - } +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; +} - /* - 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; +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; +} - 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) - { - 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)); - } +template +struct AnalyticTorus +{ + using PointType = primal::Point; - 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)); - } + PointType center; + double majorRadius {0.0}; + double minorRadius {0.0}; - double dist = sqrt(primal::squared_distance(qPt, cpCoord)); - if(!axom::utilities::isNearlyEqual(dist, analyticalDist, allowableSlack)) - { - 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)); - } - } - } - errorFlag[i] = errf; - sumErrCount += errf; - } + 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; } + } +}; - SLIC_INFO( - axom::fmt::format("Local partition has {} errors, {} warnings in closest distance results.", - sumErrCount, - sumWarningCount)); +template +using AnalyticPrimitive = + std::variant, primal::Plane, AnalyticTorus>; - return sumErrCount; - } +template +bool hasPrimitive(const AnalyticPrimitive& primitive) +{ + return !std::holds_alternative(primitive); +} -private: - BlueprintParticleMesh m_queryMesh; -}; +template +bool supportsDistanceEnvelope(const AnalyticPrimitive& primitive) +{ + return hasPrimitive(primitive) && !std::holds_alternative>(primitive); +} -/** - * 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) +template +double signedDistance(const std::monostate&, const primal::Point&) { - using axom::utilities::random_real; + return axom::numeric_limits::max(); +} - int rank = particleMesh.getRank(); - int nranks = particleMesh.getNumRanks(); +template +double signedDistance(const primal::Sphere& sphere, const primal::Point& pt) +{ + return sphere.computeSignedDistance(pt); +} - // rank scan to sum longPointCount and determine local range of longitudinal angles. - axom::Array sums(nranks, nranks); +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 { - axom::Array indivDomainCounts(nranks, nranks); - indivDomainCounts.fill(-1); - MPI_Allgather(&localDomainCount, 1, MPI_INT, indivDomainCounts.data(), 1, MPI_INT, MPI_COMM_WORLD); + if(dimension != DIM) + { + return std::monostate {}; + } - SLIC_DEBUG_IF( - params.isVerbose(), - axom::fmt::format("After all gather: [{}]", axom::fmt::join(indivDomainCounts, ","))); + const auto c = pointFromVector(center); + if((shapeName == "circle" || shapeName == "sphere") && radius > 0.0) + { + return primal::Sphere(c, radius); + } - sums[0] = indivDomainCounts[0]; - for(int i = 1; i < nranks; ++i) + if(shapeName == "plane") { - sums[i] = sums[i - 1] + indivDomainCounts[i]; + const auto n = vectorFromVector(normal); + return n.is_zero() ? AnalyticPrimitive {std::monostate {}} + : AnalyticPrimitive {primal::Plane(n, c)}; } - // If no rank has any domains, force last one to have 1 domain. - if(sums[nranks - 1] == 0) + + if constexpr(DIM == 2) { - sums[nranks - 1] = 1; - if(rank == nranks - 1) + if(shapeName == "annulus" && outerRadius > innerRadius && innerRadius > 0.0) { - localDomainCount = 1; + 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}; } } - } - - 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; + return std::monostate {}; } - 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) + static AnalyticVerification fromNode(const conduit::Node& node, const std::string& description) { - 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; + AnalyticVerification result; + result.description = description; + if(!node.has_path("shape")) + { + return result; + } - for(int li = 0; li < latPointCount; ++li) + 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")) { - 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; - } + result.center = conduitVector(node.fetch_existing("center")); } - particleMesh.setPoints(di, pts); + 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; } +}; - axom::slic::flushStreams(); - SLIC_ASSERT(particleMesh.isValid()); -} - -//--------------------------------------------------------------------------- -// Transform closest points to distances and directions -//--------------------------------------------------------------------------- template -void computeDistancesAndDirections(BlueprintParticleMesh& queryMesh, - const std::string& cpCoordsField, - const std::string& cpIndexField, - const std::string& distanceField, - const std::string& directionField) +int verifyAnalyticClosestPoints(BlueprintParticleMesh& queryMesh, + const AnalyticVerification& verification, + double distThreshold) { SLIC_ASSERT(queryMesh.dimension() == DIM); - using primal::squared_distance; using PointType = primal::Point; using IndexSet = slam::PositionSet<>; - PointType nowhere(axom::numeric_limits::signaling_NaN()); - const double nodist = axom::numeric_limits::signaling_NaN(); - - queryMesh.registerNodalScalarField(distanceField); - queryMesh.registerNodalVectorField(directionField); - for(axom::IndexType di = 0; di < queryMesh.domain_count(); ++di) + const AnalyticPrimitive primitive = verification.makePrimitive(); + if(!hasPrimitive(primitive)) { - auto cpCoords = queryMesh.getNodalVectorField(cpCoordsField, di); + SLIC_WARNING( + axom::fmt::format("Skipping unsupported analytic verification '{}'", verification.shapeName)); + return 0; + } + const bool checkDistanceEnvelope = supportsDistanceEnvelope(primitive); - auto cpIndices = queryMesh.getNodalScalarField(cpIndexField, di); + 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); - axom::ArrayView distances = queryMesh.getNodalScalarField("distance", di); - axom::ArrayView directions = queryMesh.getNodalVectorField("direction", di); - axom::IndexType ptCount = queryMesh.numPoints(di); - for(auto ptIdx : IndexSet(ptCount)) + + for(auto ptIdx : IndexSet(queryMesh.numPoints(di))) { - const bool has_cp = cpIndices[ptIdx] >= 0; - const PointType& cp = has_cp ? cpCoords[ptIdx] : nowhere; - const PointType& qPt = has_cp ? PointType(&qPts[ptIdx][0]) : nowhere; - distances[ptIdx] = has_cp ? sqrt(squared_distance(qPt, cp)) : nodist; - directions[ptIdx] = PointType(has_cp ? (cp - qPt).array() : nowhere.array()); + 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 = 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{}).", + verification.shapeName, + globalErrCount, + globalLogCount, + globalLogCount > MAX_LOCAL_LOGS * num_ranks ? " shown partly" : "")); + return globalErrCount; } void make_coords_contiguous(conduit::Node& coordValues) @@ -1253,15 +1737,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. @@ -1303,31 +1778,26 @@ int main(int argc, char** argv) slic::flushStreams(); const size_t spatialDim = queryMeshWrapper.getParticleMesh().dimension(); - SLIC_ASSERT(params.circleCenter.size() == spatialDim); //--------------------------------------------------------------------------- - // Generate object mesh + // Object (second) mesh //--------------------------------------------------------------------------- - ObjectMeshWrapper objectMeshWrapper(dataStore.getRoot()->createGroup("object_mesh", true)); - objectMeshWrapper.setVerbosity(params.isVerbose()); + ObjectMeshWrapper objectMeshWrapper(dataStore.getRoot()->createGroup("object_mesh", true), + params.objectMeshFile); + AnalyticVerification analyticVerification; + const bool hasAnalyticVerification = objectMeshWrapper.hasVerification(); + if(hasAnalyticVerification) { - 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); + 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()) @@ -1336,15 +1806,15 @@ int main(int argc, char** argv) } slic::flushStreams(); - objectMeshWrapper.saveMesh(params.objectFile); - slic::flushStreams(); - //--------------------------------------------------------------------------- // Initialize spatial index for object points, and run query //--------------------------------------------------------------------------- + int globalObjectPointCount = + allReduce(objectMeshWrapper.getParticleMesh().numPoints(), 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); @@ -1360,25 +1830,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); @@ -1393,6 +1844,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, MPI_COMM_WORLD, memUmpireId); + memProbe.report("baseline (meshes read, before BVH)"); + // Build the spatial index over the object on each rank SLIC_INFO(init_str); slic::flushStreams(); @@ -1400,28 +1860,45 @@ 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"); - 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); - }; + 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(trimMallocArena()) + { + memProbe.report("after malloc_trim(0)"); + } + else + { + SLIC_WARNING("--trim-after-query requested but malloc_trim is glibc-only; skipping."); + } + } // 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", @@ -1438,25 +1915,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(), @@ -1474,22 +1932,32 @@ 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(); 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 globalVerificationErrors == 0 ? 0 : 1; } diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt index e8782b9992..cbd6a462f3 100644 --- a/src/tools/CMakeLists.txt +++ b/src/tools/CMakeLists.txt @@ -263,6 +263,24 @@ 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 + --strip 4 + NUM_MPI_TASKS 2 + ) + + set_tests_properties(${_testname} PROPERTIES + PASS_REGULAR_EXPRESSION "Truncated [0-9]+ numeric Blueprint array") else() axom_add_python_test( NAME ${_testname} @@ -272,7 +290,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..2085a68d18 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 @@ -16,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 @@ -31,7 +36,7 @@ import numpy as np import axom.sidre as sidre -VALID_PROTOCOLS = ( +SIDRE_PROTOCOLS = ( "json", "sidre_hdf5", "sidre_conduit_json", @@ -41,27 +46,44 @@ "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", "--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", @@ -79,6 +101,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,26 +232,159 @@ 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 sidre_protocol(protocol: str) -> str: + if protocol in SIDRE_PROTOCOLS: + return protocol - comm_size = MPI.COMM_WORLD.Get_size() + 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 + 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: + 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 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}", ) + 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()}") + + 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}", ) + conduit.relay.io.blueprint.save_mesh(mesh, args.output, protocol) + + 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") + + protocol = sidre_protocol(args.protocol) input_path = Path(args.input) manager = sidre.IOManager() datastore = sidre.DataStore() @@ -254,24 +428,29 @@ def main() -> 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}", - ) - manager.write(root, num_files, args.output, args.protocol) - - if initialized_mpi and not MPI.Is_finalized(): - MPI.Finalize() + 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 +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: + 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()) diff --git a/src/tools/gen-multidom-point-mesh.py b/src/tools/gen-multidom-point-mesh.py new file mode 100755 index 0000000000..8423116234 --- /dev/null +++ b/src/tools/gen-multidom-point-mesh.py @@ -0,0 +1,1123 @@ +#!/usr/bin/env python3 + +# gen-multidom-point-mesh.py +# Write a mesh following the Conduit mesh blueprint for distributed closest point (DCP) testing. +# +# The generated Blueprint hierarchy is: +# (single-domain output) +# / (multidomain output) +# |-- state +# | `-- 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 +# | |-- 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 == +# | |-- 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 (i-fastest node ordering for grids) +# | |-- 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) +# +# 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, +# 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 os +import shlex +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' +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 + import conduit.relay.mpi + import conduit.relay.mpi.io + import conduit.relay.mpi.io.blueprint + 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 { + '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 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('--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') + + +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 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, + default=(-90.0, 90.0), + help='Latitude range in degrees for latitude sampling') + sphere.add_argument('--random-spacing', + action='store_true', + 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') + + 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') + 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: + 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, 'annulus': 2} + if opts.shape in fixed_dims: + return fixed_dims[opts.shape] + + dim = opts.dimension + inferred = [] + 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)) + + 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: + 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) + if len(values) != dim: + raise RuntimeError(f'{name} must have one value or {dim} values') + 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: + 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') + + 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)) + + 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_id, + 'domain_index': domain_index, + 'points': points[begin:end], + 'global_id_start': begin, + 'grid_counts': None, + }) + 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, 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(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 + + return points + + +def generate_torus(point_count, center, major_radius, minor_radius): + 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_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))) + + +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(order='F') for coord in coords]) + + +def generate_grid_domains(grid_counts, lower, upper, domain_counts): + '''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 + 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]) + + 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, normal, stddev, domain_counts): + '''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) + + 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, opts.random_spacing, + opts.seed) + elif opts.shape == 'sphere': + 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': + 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 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 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 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 + + 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]}' + dom[values_path].set(np.ascontiguousarray(points[:, d], dtype=np.float64)) + + topo = dom[f'topologies/{opts.topology_name}'] + topo['coordset'] = opts.coordset_name + + 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(chunk['global_id_start'], + chunk['global_id_start'] + point_count, + dtype=np.int64)) + + +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, metadata) + + +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(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']) + + if np.any(upper <= lower): + 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 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') + + metadata = { + 'description': shape_description(opts, dim, lower, upper, center, normal, stddev), + 'command_line': sanitized_command_line(), + '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], metadata) + else: + for chunk in local_chunks: + create_domain(mesh, opts, chunk, metadata) + + 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(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}') + + return 0 + + +if __name__ == '__main__': + sys.exit(main())