diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index 9baba4e..6a54f23 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -27,6 +27,13 @@ jobs: sudo apt-get install -yq cmake cmake --version + # There is a known bug using MPICH ubuntu-24.04 (ubuntu-latest as of 11/10/2025) + # https://bugs.launchpad.net/ubuntu/+source/mpich/+bug/2072338 + - name: Install mpi + run: | + sudo apt-get update -yq + sudo apt-get install -yq openmpi-bin libopenmpi-dev + ## Kokkos - name: Kokkos Checkout repo uses: actions/checkout@v4 diff --git a/CMakeLists.txt b/CMakeLists.txt index c042d51..208ccd7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,6 +9,7 @@ option(MeshFields_USE_Cabana "Build with the Cabana storage backend" OFF) find_package(Kokkos REQUIRED) find_package(Omega_h REQUIRED) +find_package(MPI REQUIRED) # Verify Omega_h was built with Kokkos enabled if(NOT Omega_h_USE_Kokkos) @@ -234,6 +235,18 @@ add_test_property(OmegahSprAdapt PASS_REGULAR_EXPRESSION "afterAdapt nTri: 14439; solution_1 min -4084.78 max 213.322") +test_func(OmegahSprAdapt_p2 + ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 ${MPIEXEC_PREFLAGS} + ./OmegahSprAdapt ${CMAKE_SOURCE_DIR}/meshes/thwaites_basal_effectiveStrain.osh/ + thwaitesAdapted_p2 #output prefix + 0.1 # adapt ratio + 1.0 # min size + 16.0 # max size +) +add_test_property(OmegahSprAdapt_p2 + PASS_REGULAR_EXPRESSION + "afterAdapt nTri: 14468; solution_1 min -4098.9 max 213.322") + #Code Coverage set up ------------------------------------------------------- option(meshfields_ENABLE_COVERAGE_BUILD "Do a coverage build" OFF) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index 100f268..b2b1ebe 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -196,9 +196,12 @@ auto getJacobianDeterminants(FieldElement &fes, * - `post`: Called after the integration process ends. * - `atPoints`: Called for each integration point to perform user-defined * computations. + * - `parallelReduce`: Called after `atPoints` to reduce process-local + * integrations into a global mesh integration. * * To use this class, derive from it and implement the `atPoints` method. - * Optionally, override `pre` and `post` for additional setup or cleanup. + * Optionally, override `parallelReduce`, `pre` and `post` for distributed + * memory support, additional setup, or cleanup. */ class Integrator { public: @@ -227,9 +230,18 @@ class Integrator { */ virtual void atPoints(Kokkos::View p, Kokkos::View w, Kokkos::View dV) = 0; + + /** \brief User callback: distributed memory accumulation. + * + * \details If needed, this function supports the use of + * inter-process communication (e.g., MPI) to reduce + * process-local integrations into a global + * mesh integration. + */ + virtual void parallelReduce(){}; + /** \brief Run the Integrator over the local field elements. - * \param fes FieldElement - * FIXME make the sensible + * \param fes FieldElement * */ template void process(FieldElement &fes) { constexpr auto topo = decltype(FieldElement::elm2dof)::getTopology(); @@ -239,9 +251,8 @@ class Integrator { auto weights = getIntegrationPointWeights(fes, ip); auto dV = getJacobianDeterminants(fes, localCoords, ip.size()); atPoints(localCoords, weights, dV); + parallelReduce(); post(); - // TODO support distributed meshes by running a parallel reduction with user - // functor } protected: diff --git a/src/MeshField_SPR_ErrorEstimator.hpp b/src/MeshField_SPR_ErrorEstimator.hpp index 97eb726..6503e12 100644 --- a/src/MeshField_SPR_ErrorEstimator.hpp +++ b/src/MeshField_SPR_ErrorEstimator.hpp @@ -19,9 +19,12 @@ static double getNaN() { return std::numeric_limits::quiet_NaN(); } /* common base for Scalar Integrator. */ class SInt : public MeshField::Integrator { public: - SInt(int order) : MeshField::Integrator(order), r(0) {} + SInt(int order, Omega_h::CommPtr comm_in) + : MeshField::Integrator(order), r(0), comm(comm_in) {} void reset() { r = 0; } MeshField::Real r; + Omega_h::CommPtr comm; + void parallelReduce() { r = comm->allreduce(r, OMEGA_H_SUM); } }; template class Estimation { @@ -83,8 +86,8 @@ template class SelfProduct : public SInt { public: SelfProduct(EstimationT &estimation_in, OmegahMeshField &omf_in) - : SInt(estimation_in.integration_order), estimation(estimation_in), - omf(omf_in) {} + : SInt(estimation_in.integration_order, estimation_in.mesh.comm()), + estimation(estimation_in), omf(omf_in) {} void atPoints(Kokkos::View p, Kokkos::View w, Kokkos::View dV) { @@ -155,8 +158,9 @@ template class Error : public SInt { public: Error(EstimationT &estimation_in, OmegahMeshField &omf_in) - : SInt(estimation_in.integration_order), estimation(estimation_in), - omf(omf_in), errorNorm("errorNorm", estimation_in.mesh.nelems()) {} + : SInt(estimation_in.integration_order, estimation_in.mesh.comm()), + estimation(estimation_in), omf(omf_in), + errorNorm("errorNorm", estimation_in.mesh.nelems()) {} void atPoints(Kokkos::View p, Kokkos::View w, Kokkos::View dV) { @@ -270,7 +274,8 @@ template std::tuple, MeshField::Real> getSprSizeField(EstimationT &e, OmegahMeshField &omf, FieldElement &coordFe) { Error errorIntegrator(e, omf); - errorIntegrator.process(coordFe); + errorIntegrator.process( + coordFe); // FIXME - needs to sync results across processes computeSizeFactor(e, omf, coordFe, errorIntegrator); getElementSizeField(e, errorIntegrator); return {averageToVertex(e.mesh, e.element_size), errorIntegrator.r}; diff --git a/test/testSprThwaitesAdapt.cpp b/test/testSprThwaitesAdapt.cpp index 4ed0c32..55374ec 100644 --- a/test/testSprThwaitesAdapt.cpp +++ b/test/testSprThwaitesAdapt.cpp @@ -157,6 +157,7 @@ int main(int argc, char **argv) { auto world = lib.world(); Omega_h::Mesh mesh(&lib); Omega_h::binary::read(argv[1], world, &mesh); + const auto prefix = std::string(argv[2]); // adaptRatio = 0.1 is used in scorec/core:cws/sprThwaites test/spr_test.cc const MeshField::Real adaptRatio = std::stof(argv[3]); @@ -166,13 +167,21 @@ int main(int argc, char **argv) { argv[5]); // a value of 8 or 16 roughly maintains the boundary shape // in the downstream parts of the domain // where there is significant coarsening - std::cout << "input mesh: " << argv[1] << " outputMeshPrefix: " << prefix - << " adaptRatio: " << adaptRatio << " max_size: " << max_size - << " min_size: " << min_size << "\n"; + if (!mesh.comm()->rank()) { + std::cout << "input mesh: " << argv[1] << " outputMeshPrefix: " << prefix + << " adaptRatio: " << adaptRatio << " max_size: " << max_size + << " min_size: " << min_size << "\n"; + } const auto outname = prefix + "_adaptRatio_" + std::string(argv[3]) + "_maxSz_" + std::to_string(max_size) + "_minSz_" + std::to_string(min_size); + // equally distribute elements among processes; required when loading a serial + // mesh + mesh.balance(); + // Omega_h::project_by_fit used by spr requires ghosts + mesh.set_parting(Omega_h_Parting::OMEGA_H_GHOSTED); + auto effectiveStrain = getEffectiveStrainRate(mesh); auto recoveredStrain = recoverLinearStrain(mesh, effectiveStrain); mesh.add_tag(VERT, "recoveredStrain", 1, recoveredStrain, false, @@ -203,10 +212,10 @@ int main(int argc, char **argv) { Omega_h::ArrayType::VectorND); { // write vtk - const std::string vtkFileName = "beforeAdapt" + outname + ".vtk"; + const std::string vtkFileName = "beforeAdapt_" + outname + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); const std::string vtkFileName_edges = - "beforeAdapt" + outname + "_edges.vtk"; + "beforeAdapt_" + outname + "_edges.vtk"; Omega_h::vtk::write_parallel(vtkFileName_edges, &mesh, 1); } @@ -223,12 +232,14 @@ int main(int argc, char **argv) { printTriCount(&mesh, "afterAdapt"); const auto sol = mesh.get_array(VERT, "solution_1"); - const auto sol_min = Omega_h::get_min(sol); - const auto sol_max = Omega_h::get_max(sol); - std::cout << "solution_1 min " << sol_min << " max " << sol_max << '\n'; + const auto sol_min = Omega_h::get_min(mesh.comm(), sol); + const auto sol_max = Omega_h::get_max(mesh.comm(), sol); + if (!mesh.comm()->rank()) { + std::cout << "solution_1 min " << sol_min << " max " << sol_max << '\n'; + } { // write vtk and osh for adapted mesh - const std::string outfilename = "afterAdapt" + outname; + const std::string outfilename = "afterAdapt_" + outname; Omega_h::vtk::write_parallel(outfilename + ".vtk", &mesh, 2); Omega_h::binary::write(outfilename + ".osh", &mesh); std::cout << "wrote adapted mesh: " << outfilename + ".osh"