From 93f34c199a5e2fe95da776035799514c563b95d3 Mon Sep 17 00:00:00 2001 From: Hugh Carson Date: Wed, 26 Jul 2023 10:04:59 -0400 Subject: [PATCH] Change to allow for partitioning of Nonconformal meshes. - The partitioning is done using the MeshPartitioner where possible, but for nonconformal meshes will now use the standard ParMesh constructor. - Change to use Mesh& where a parameter out value is not allowed to be invalidated - Move #if 0 to behind constexpr false to ensure can continue to compile --- palace/utils/configfile.cpp | 51 ++++ palace/utils/configfile.hpp | 32 ++ palace/utils/geodata.cpp | 493 +++++++++++++++++-------------- scripts/schema/config/model.json | 19 ++ 4 files changed, 374 insertions(+), 221 deletions(-) diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index 34803c400..d417d248e 100644 --- a/palace/utils/configfile.cpp +++ b/palace/utils/configfile.cpp @@ -397,10 +397,61 @@ void RefinementData::SetUp(json &model) } } + auto adapt = refinement->find("Adaptation"); + if (adapt != refinement->end()) + { + MFEM_ABORT("Placeholder: Not currently supported!"); + + // Load Values + adaptation.tolerance = adapt->value("Tol", adaptation.tolerance); + adaptation.max_its = adapt->value("MaxIts", adaptation.max_its); + adaptation.min_its = adapt->value("MinIts", adaptation.min_its); + adaptation.update_fraction = adapt->value("UpdateFraction", adaptation.update_fraction); + adaptation.max_nc_levels = adapt->value("MaxNCLevels", adaptation.max_nc_levels); + adaptation.dof_limit = adapt->value("DOFLimit", adaptation.dof_limit); + adaptation.coarsening_fraction = + adapt->value("CoarseningFraction", adaptation.coarsening_fraction); + adaptation.save_step = adapt->value("SaveStep", adaptation.save_step); + adaptation.nonconformal = adapt->value("Nonconformal", adaptation.nonconformal); + adaptation.maximum_imbalance = + adapt->value("MaximumImbalance", adaptation.maximum_imbalance); + + // Perform Checks + MFEM_VERIFY(adaptation.tolerance > 0, "\"Tol\" must be strictly positive"); + MFEM_VERIFY(adaptation.max_its >= 0, "\"MaxIts\" must be non-negative"); + MFEM_VERIFY(adaptation.min_its >= 0, "\"MinIts\" must be non-negative"); + MFEM_VERIFY(adaptation.min_its <= adaptation.max_its, + "\"MinIts\" must be smaller than \"MaxIts\": " << adaptation.min_its << "," + << adaptation.max_its); + MFEM_VERIFY(adaptation.update_fraction > 0 && adaptation.update_fraction < 1, + "\"UpdateFraction\" must be in (0,1)"); + MFEM_VERIFY(adaptation.coarsening_fraction >= 0 && adaptation.coarsening_fraction < 1, + "\"CoarseningFraction\" must be in [0, 1)"); + MFEM_VERIFY(adaptation.max_nc_levels >= 0, "\"MaxNCLevels\" must non-negative"); + MFEM_VERIFY(adaptation.dof_limit >= 0, "\"DOFLimit\" must be non-negative"); + MFEM_VERIFY(adaptation.save_step >= 0, "\"SaveStep\" must be non-negative"); + MFEM_VERIFY(adaptation.maximum_imbalance >= 1, + "\"MaximumImbalance\" must be greater than or equal to 1"); + + // Cleanup + const auto fields = { + "Tol", "MaxIts", "MinIts", "UpdateFraction", "CoarseningFraction", + "MaxNCLevels", "DOFLimit", "SaveStep", "Nonconformal", "MaximumImbalance"}; + for (const auto &f : fields) + { + adapt->erase(f); + } + + MFEM_VERIFY(adapt->empty(), + "Found an unsupported configuration file keyword under \"Adaptation\"!\n" + << adapt->dump(2)); + } + // Cleanup refinement->erase("UniformLevels"); refinement->erase("Boxes"); refinement->erase("Spheres"); + refinement->erase("Adaptation"); MFEM_VERIFY(refinement->empty(), "Found an unsupported configuration file keyword under \"Refinement\"!\n" << refinement->dump(2)); diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp index 17c55ae17..fe85d1b5a 100644 --- a/palace/utils/configfile.hpp +++ b/palace/utils/configfile.hpp @@ -128,11 +128,43 @@ struct SphereRefinementData std::vector center = {}; }; +// Stores data specifying the adaptive mesh refinement algorithm. +struct AdaptiveRefinementData +{ + // Non-dimensional tolerance used to specify convergence of the AMR. + double tolerance = 1e-2; + // Maximum number of iterations to perform during the AMR. + int max_its = 0; + // Minimum number of iterations to perform during the AMR. + int min_its = 0; + // Dörfler update fraction. The set of marked elements is the minimum set + // that contains update_fraction of the total error. + double update_fraction = 0.4; + // Whether or not to perform coarsening during the AMR. + double coarsening_fraction = 0.0; + // Maximum difference in non-conformal refinements between two adjacent + // elements. Default = 0 implies there is no constraint on local non-conformity. + int max_nc_levels = 0; + // If a refinement results in a greater number of DOFs than this value, no + // future refinement will be allowed unless coarsening is allowed to occur. + int dof_limit = 0; + // Frequency with which to store the post processing results for a given + // adaptation, e.g. save_step = 3 means save every third adaptation. + int save_step = 0; + // Whether or not to perform nonconformal adaptation. + bool nonconformal = true; + // Maximum allowable ratio of number of elements across processors before + // rebalancing is performed. + double maximum_imbalance = 1.5; +}; + struct RefinementData { public: // Parallel uniform mesh refinement levels. int uniform_ref_levels = 0; + // Adaptive refinement configuration data. + AdaptiveRefinementData adaptation; private: // Refinement data for mesh regions. diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index ce51a140e..e6d6d6b88 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -27,35 +27,52 @@ namespace // Floating point precision for mesh IO. This precision is important, make sure nothing is // lost! -const auto MSH_FLT_PRECISION = std::numeric_limits::max_digits10; +constexpr auto MSH_FLT_PRECISION = std::numeric_limits::max_digits10; // Load the serial mesh from disk. -std::unique_ptr LoadMesh(const std::string &, bool); +std::unique_ptr LoadMesh(const std::string &path, bool remove_curvature); // Optionally reorder mesh elements based on MFEM's internal reordeing tools for improved // cache usage. -void ReorderMesh(mfem::Mesh &); +void ReorderMesh(mfem::Mesh &mesh); // Generate element-based mesh partitioning, using either a provided file or METIS. -std::unique_ptr GetMeshPartitioning(mfem::Mesh &, int, const std::string &); +std::unique_ptr GetMeshPartitioning(mfem::Mesh &mesh, int size, + const std::string &partition = ""); // Cleanup the provided serial mesh by removing unnecessary domain and elements, adding // boundary elements for material interfaces and exterior boundaries, and adding boundary // elements for subdomain interfaces. -std::map> CheckMesh(std::unique_ptr &, - const std::unique_ptr &, const IoData &, - bool, bool, bool); +std::map> CheckMesh(mfem::Mesh &orig_mesh, + const std::unique_ptr &partitioning, + const IoData &iodata, bool clean_elem, + bool add_bdr, bool add_subdomain); -// Given a serial mesh on the root processor and element partitioning, create a parallel -// mesh oer the given communicator. -std::unique_ptr DistributeMesh(MPI_Comm, std::unique_ptr &, - std::unique_ptr &, - const std::string &); +// Given a serial mesh on the root processor, and element partitioning, create a parallel +// mesh over the given communicator. +std::unique_ptr +DistributeMesh(MPI_Comm comm, std::unique_ptr &mesh, + const std::unique_ptr &partitioning = nullptr, + const std::string &output_dir = ""); // Get list of domain and boundary attribute markers used in configuration file for mesh // cleaning. -void GetUsedAttributeMarkers(const IoData &, int, int, mfem::Array &, - mfem::Array &); +void GetUsedAttributeMarkers(const IoData &iodata, int n_mat, int n_bdr, + mfem::Array &mat_marker, mfem::Array &bdr_marker); + +// Simplified helper for describing the element types in a mesh +struct ElementTypeInfo +{ + bool has_simplices; + bool has_tensors; + bool has_wedges; + bool has_pyramids; +}; +ElementTypeInfo CheckElements(mfem::Mesh &mesh) +{ + auto meshgen = mesh.MeshGenerator(); + return {bool(meshgen & 1), bool(meshgen & 2), bool(meshgen & 4), bool(meshgen & 8)}; +} } // namespace @@ -65,43 +82,73 @@ namespace mesh std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, bool reorder, bool clean, bool add_bdr, bool unassembled) { - // On root, read the serial mesh (converting format if necessary), and do all necessary - // serial preprocessing. When finished, distribute the mesh to all processes. Count disk - // I/O time separately for the mesh read from file. + // If possible on root, read the serial mesh (converting format if necessary), and do all + // necessary serial preprocessing. When finished, distribute the mesh to all processes. + // Count disk I/O time separately for the mesh read from file. + + // If not adapting, or performing conformal adaptation, can use the mesh partitioner. std::unique_ptr smesh; + std::unique_ptr partitioning; + const auto use_amr = iodata.model.refinement.adaptation.max_its > 0; + const bool use_mesh_partitioner = + !use_amr || !iodata.model.refinement.adaptation.nonconformal; { BlockTimer bt(Timer::IO); - if (Mpi::Root(comm)) + + if (Mpi::Root(comm) || !use_mesh_partitioner) { + // If using the mesh partitioner, only the root node needs to load the mesh. // Optionally reorder elements (and vertices) based on spatial location after loading // the serial mesh. smesh = LoadMesh(iodata.model.mesh, iodata.model.remove_curvature); if (reorder) { - ReorderMesh(*smesh); + // Optionally reorder elements (and vertices) based on spatial location after loading + // the serial mesh. + smesh = LoadMesh(iodata.model.mesh, iodata.model.remove_curvature); + if (reorder) + { + ReorderMesh(*smesh); + } } + Mpi::Barrier(comm); } - Mpi::Barrier(comm); } - std::unique_ptr partitioning; - if (Mpi::Root(comm)) + if (Mpi::Root(comm) || !use_mesh_partitioner) { - // Generate the parallel mesh partitioning on the root process. + const auto element_types = CheckElements(*smesh); + + MFEM_VERIFY(!use_amr || !element_types.has_tensors || + iodata.model.refinement.adaptation.nonconformal, + "If there are tensor elements, AMR must be nonconformal"); + MFEM_VERIFY(!use_amr || !element_types.has_pyramids || + iodata.model.refinement.adaptation.nonconformal, + "If there are pyramid elements, AMR must be nonconformal"); + MFEM_VERIFY(!use_amr || !element_types.has_wedges || + iodata.model.refinement.adaptation.nonconformal, + "If there are wedge elements, AMR must be nonconformal"); + partitioning = GetMeshPartitioning(*smesh, Mpi::Size(comm), iodata.model.partition); // Clean up unused domain elements from the mesh, add new boundary elements for material // interfaces if not present, and optionally (when running unassembled) add subdomain // interface boundary elements. - std::map> attr_map = - CheckMesh(smesh, partitioning, iodata, clean, add_bdr, unassembled); + static_cast(CheckMesh(*smesh, partitioning, iodata, clean, add_bdr, unassembled)); } - // Construct the parallel mesh data structure by distributing the serial mesh from the - // root process. The serial mesh and partitioning are deleted inside. - std::unique_ptr mesh = - DistributeMesh(comm, smesh, partitioning, iodata.problem.output); -#if 0 + std::unique_ptr pmesh; + if (iodata.model.refinement.adaptation.nonconformal && use_amr) + { + smesh->EnsureNCMesh(true); + pmesh = std::make_unique(comm, *smesh, partitioning.get()); + } + else + { + pmesh = DistributeMesh(comm, smesh, partitioning, iodata.problem.output); + } + + if constexpr (false) { std::string tmp = iodata.problem.output; if (tmp.back() != '/') @@ -113,28 +160,30 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, boo { std::filesystem::create_directories(tmp); } - int width = 1 + static_cast(std::log10(Mpi::Size(comm)-1)); - std::unique_ptr gsmesh = LoadMesh(iodata.model.mesh); + int width = 1 + static_cast(std::log10(Mpi::Size(comm) - 1)); + std::unique_ptr gsmesh = + LoadMesh(iodata.model.mesh, iodata.model.remove_curvature); std::unique_ptr gpartitioning = GetMeshPartitioning(*gsmesh, Mpi::Size(comm)); mfem::ParMesh gpmesh(comm, *gsmesh, gpartitioning.get(), 0); { - std::string pfile = mfem::MakeParFilename(tmp + "part.", Mpi::Rank(comm), ".mesh", width); + std::string pfile = + mfem::MakeParFilename(tmp + "part.", Mpi::Rank(comm), ".mesh", width); std::ofstream fo(pfile); // mfem::ofgzstream fo(pfile, true); // Use zlib compression if available fo.precision(MSH_FLT_PRECISION); gpmesh.ParPrint(fo); } { - std::string pfile = mfem::MakeParFilename(tmp + "final.", Mpi::Rank(comm), ".mesh", width); + std::string pfile = + mfem::MakeParFilename(tmp + "final.", Mpi::Rank(comm), ".mesh", width); std::ofstream fo(pfile); // mfem::ofgzstream fo(pfile, true); // Use zlib compression if available fo.precision(MSH_FLT_PRECISION); - mesh->ParPrint(fo); + pmesh->ParPrint(fo); } } -#endif - return mesh; + return pmesh; } void RefineMesh(const IoData &iodata, std::vector> &mesh) @@ -183,7 +232,7 @@ void RefineMesh(const IoData &iodata, std::vector { if (mesh.capacity() > 1) { - mesh.push_back(std::make_unique(*mesh.back())); + mesh.emplace_back(std::make_unique(*mesh.back())); } mesh.back()->UniformRefinement(); } @@ -302,7 +351,7 @@ void RefineMesh(const IoData &iodata, std::vector // (adds hanging nodes). if (mesh.capacity() > 1) { - mesh.push_back(std::make_unique(*mesh.back())); + mesh.emplace_back(std::make_unique(*mesh.back())); } mesh.back()->GeneralRefinement(refs, -1); region_ref_level++; @@ -1168,24 +1217,25 @@ void ReorderMesh(mfem::Mesh &mesh) { mfem::Array ordering; -#if 0 - // Gecko reordering. - Mpi::Print(\n); - mfem::Array tentative; - int outer = 3, inner = 3, window = 4, period = 2; - double best_cost = mfem::infinity(); - for (int i = 0; i < outer; i++) - { - int seed = i+1; - double cost = mesh.GetGeckoElementOrdering(tentative, inner, window, eriod, seed, true); - if (cost < best_cost) + if constexpr (false) + { + // Gecko reordering. + mfem::Array tentative; + int outer = 3, inner = 3, window = 4, period = 2; + double best_cost = mfem::infinity(); + for (int i = 0; i < outer; i++) { - ordering = tentative; - best_cost = cost; + int seed = i + 1; + double cost = + mesh.GetGeckoElementOrdering(tentative, inner, window, period, seed, true); + if (cost < best_cost) + { + ordering = tentative; + best_cost = cost; + } } + Mpi::Print("Final cost: {:e}\n", best_cost); } - Mpi::Print("Final cost: {:e}\n", best_cost); -#endif // (Faster) Hilbert reordering. mesh.GetHilbertElementOrdering(ordering); @@ -1239,7 +1289,7 @@ std::unique_ptr GetMeshPartitioning(mfem::Mesh &mesh, int size, return partitioning; } -std::map> CheckMesh(std::unique_ptr &orig_mesh, +std::map> CheckMesh(mfem::Mesh &orig_mesh, const std::unique_ptr &partitioning, const IoData &iodata, bool clean_elem, bool add_bdr, bool add_subdomain) @@ -1250,22 +1300,24 @@ std::map> CheckMesh(std::unique_ptr &orig_me // interfaces if these elements do not yet exist. // - If desired, create a new mesh which has removed all domain elements which do not have // an associated material property specified in the input file. - MFEM_VERIFY(orig_mesh->Dimension() == 3 && !orig_mesh->Nonconforming(), + MFEM_VERIFY(orig_mesh.Dimension() == 3 && !orig_mesh.Nonconforming(), "Nonconforming or 2D meshes have not been tested yet!"); + MFEM_VERIFY(dynamic_cast(&orig_mesh) == nullptr, + "This function does not work for ParMesh"); mfem::Array mat_marker, bdr_marker; GetUsedAttributeMarkers( - iodata, orig_mesh->attributes.Size() ? orig_mesh->attributes.Max() : 0, - orig_mesh->bdr_attributes.Size() ? orig_mesh->bdr_attributes.Max() : 0, mat_marker, + iodata, orig_mesh.attributes.Size() ? orig_mesh.attributes.Max() : 0, + orig_mesh.bdr_attributes.Size() ? orig_mesh.bdr_attributes.Max() : 0, mat_marker, bdr_marker); bool warn = false; - for (int be = 0; be < orig_mesh->GetNBE(); be++) + for (int be = 0; be < orig_mesh.GetNBE(); be++) { - int attr = orig_mesh->GetBdrAttribute(be); + int attr = orig_mesh.GetBdrAttribute(be); if (!bdr_marker[attr - 1]) { int f, o, e1, e2; - orig_mesh->GetBdrElementFace(be, &f, &o); - orig_mesh->GetFaceElements(f, &e1, &e2); + orig_mesh.GetBdrElementFace(be, &f, &o); + orig_mesh.GetFaceElements(f, &e1, &e2); if (e1 < 0 || e2 < 0) // Internal boundary elements are allowed to have no BC { warn = true; @@ -1288,33 +1340,33 @@ std::map> CheckMesh(std::unique_ptr &orig_me } // Count deleted or added domain and boundary elements. - int new_ne = orig_mesh->GetNE(); - int new_nbdr = orig_mesh->GetNBE(); + int new_ne = orig_mesh.GetNE(); + int new_nbdr = orig_mesh.GetNBE(); mfem::Array elem_delete, bdr_delete; mfem::Array orig_bdr_faces, add_bdr_faces; - elem_delete.SetSize(orig_mesh->GetNE(), false); - bdr_delete.SetSize(orig_mesh->GetNBE(), false); - orig_bdr_faces.SetSize(orig_mesh->GetNumFaces(), -1); - for (int be = 0; be < orig_mesh->GetNBE(); be++) + elem_delete.SetSize(orig_mesh.GetNE(), false); + bdr_delete.SetSize(orig_mesh.GetNBE(), false); + orig_bdr_faces.SetSize(orig_mesh.GetNumFaces(), -1); + for (int be = 0; be < orig_mesh.GetNBE(); be++) { int f, o; - orig_mesh->GetBdrElementFace(be, &f, &o); + orig_mesh.GetBdrElementFace(be, &f, &o); MFEM_VERIFY(orig_bdr_faces[f] < 0, "Mesh should not define boundary elements multiple times!"); orig_bdr_faces[f] = be; } if (add_bdr || add_subdomain) { - add_bdr_faces.SetSize(orig_mesh->GetNumFaces(), -1); + add_bdr_faces.SetSize(orig_mesh.GetNumFaces(), -1); } if (clean_elem) { // Delete domain and boundary elements which have no associated material or BC attribute // from the mesh. - for (int e = 0; e < orig_mesh->GetNE(); e++) + for (int e = 0; e < orig_mesh.GetNE(); e++) { - int attr = orig_mesh->GetAttribute(e); + int attr = orig_mesh.GetAttribute(e); if (!mat_marker[attr - 1]) { elem_delete[e] = true; @@ -1324,13 +1376,13 @@ std::map> CheckMesh(std::unique_ptr &orig_me // Make sure to remove any boundary elements which are no longer attached to elements in // the domain. - for (int f = 0; f < orig_mesh->GetNumFaces(); f++) + for (int f = 0; f < orig_mesh.GetNumFaces(); f++) { const int &be = orig_bdr_faces[f]; if (be >= 0) { int e1, e2; - orig_mesh->GetFaceElements(f, &e1, &e2); + orig_mesh.GetFaceElements(f, &e1, &e2); if ((e1 < 0 || elem_delete[e1]) && (e2 < 0 || elem_delete[e2])) { // Mpi::Print("Deleting an unattached boundary element!\n"); @@ -1339,15 +1391,15 @@ std::map> CheckMesh(std::unique_ptr &orig_me } } } - if (new_ne < orig_mesh->GetNE()) + if (new_ne < orig_mesh.GetNE()) { Mpi::Print("Removed {:d} unmarked domain elements from the mesh\n", - orig_mesh->GetNE() - new_ne); + orig_mesh.GetNE() - new_ne); } - if (new_nbdr < orig_mesh->GetNBE()) + if (new_nbdr < orig_mesh.GetNBE()) { Mpi::Print("Removed {:d} unattached boundary elements from the mesh\n", - orig_mesh->GetNBE() - new_nbdr); + orig_mesh.GetNBE() - new_nbdr); } } int new_ne_step1 = new_ne; @@ -1357,16 +1409,17 @@ std::map> CheckMesh(std::unique_ptr &orig_me { // Add new boundary elements at material interfaces or on the exterior boundary of the // simulation domain, if there is not already a boundary element present. - MFEM_VERIFY(!orig_mesh->Nonconforming(), "Adding material interface boundary elements " - "is not supported for nonconforming meshes!"); + MFEM_VERIFY(!orig_mesh.Nonconforming(), "Adding material interface boundary elements " + "is not supported for nonconforming meshes!"); int add_bdr_ext = 0, add_bdr_int = 0; - for (int f = 0; f < orig_mesh->GetNumFaces(); f++) + for (int f = 0; f < orig_mesh.GetNumFaces(); f++) { const int &be = orig_bdr_faces[f]; if (be < 0 && add_bdr_faces[f] < 0) { int e1, e2; - orig_mesh->GetFaceElements(f, &e1, &e2); + orig_mesh.GetFaceElements(f, &e1, &e2); + bool no_e1 = (e1 < 0 || elem_delete[e1]); bool no_e2 = (e2 < 0 || elem_delete[e2]); if ((no_e1 || no_e2) && !(no_e1 && no_e2)) @@ -1375,7 +1428,7 @@ std::map> CheckMesh(std::unique_ptr &orig_me add_bdr_faces[f] = 1; add_bdr_ext++; } - else if (orig_mesh->GetAttribute(e1) != orig_mesh->GetAttribute(e2)) + else if (orig_mesh.GetAttribute(e1) != orig_mesh.GetAttribute(e2)) { // Add new boundary element at material interface between two domains. // Mpi::Print("Adding material interface boundary element!\n"); @@ -1401,19 +1454,20 @@ std::map> CheckMesh(std::unique_ptr &orig_me if (add_subdomain) { - // Add new boundary elements at interfaces between elements beloning to different + // Add new boundary elements at interfaces between elements belonging to different // subdomains. This uses similar code to mfem::Mesh::PrintWithPartitioning. MFEM_VERIFY(partitioning, "Cannot add subdomain interface boundary elements without " "supplied mesh partitioning!"); - MFEM_VERIFY(!orig_mesh->Nonconforming(), "Adding subdomain interface boundary elements " - "is not supported for nonconforming meshes!"); - for (int f = 0; f < orig_mesh->GetNumFaces(); f++) + MFEM_VERIFY(!orig_mesh.Nonconforming(), "Adding subdomain interface boundary elements " + "is not supported for nonconforming meshes!"); + for (int f = 0; f < orig_mesh.GetNumFaces(); f++) { const int &be = orig_bdr_faces[f]; if (be < 0 && add_bdr_faces[f] < 0) { int e1, e2; - orig_mesh->GetFaceElements(f, &e1, &e2); + + orig_mesh.GetFaceElements(f, &e1, &e2); bool no_e1 = (e1 < 0 || elem_delete[e1]); bool no_e2 = (e2 < 0 || elem_delete[e2]); if (!no_e1 && !no_e2 && partitioning[e1] != partitioning[e2]) @@ -1440,62 +1494,54 @@ std::map> CheckMesh(std::unique_ptr &orig_me // Create the new mesh. if (new_ne == new_ne_step1 && new_ne_step1 == new_ne_step2 && - new_ne_step2 == orig_mesh->GetNE() && new_nbdr == new_nbdr_step1 && - new_nbdr_step1 == new_nbdr_step2 && new_nbdr_step2 == orig_mesh->GetNBE()) + new_ne_step2 == orig_mesh.GetNE() && new_nbdr == new_nbdr_step1 && + new_nbdr_step1 == new_nbdr_step2 && new_nbdr_step2 == orig_mesh.GetNBE()) { return new_attr_map; } - std::unique_ptr new_mesh = - std::make_unique(orig_mesh->Dimension(), orig_mesh->GetNV(), new_ne, - new_nbdr, orig_mesh->SpaceDimension()); + + auto new_mesh = mfem::Mesh(orig_mesh.Dimension(), orig_mesh.GetNV(), new_ne, new_nbdr, + orig_mesh.SpaceDimension()); // Copy vertices and non-deleted domain and boundary elements. - for (int v = 0; v < orig_mesh->GetNV(); v++) + for (int v = 0; v < orig_mesh.GetNV(); v++) { - new_mesh->AddVertex(orig_mesh->GetVertex(v)); + new_mesh.AddVertex(orig_mesh.GetVertex(v)); } - for (int e = 0; e < orig_mesh->GetNE(); e++) + for (int e = 0; e < orig_mesh.GetNE(); e++) { if (!elem_delete[e]) { - mfem::Element *ne = orig_mesh->GetElement(e)->Duplicate(new_mesh.get()); - new_mesh->AddElement(ne); + mfem::Element *ne = orig_mesh.GetElement(e)->Duplicate(&new_mesh); + new_mesh.AddElement(ne); } } - for (int be = 0; be < orig_mesh->GetNBE(); be++) + for (int be = 0; be < orig_mesh.GetNBE(); be++) { if (!bdr_delete[be]) { - mfem::Element *ne = orig_mesh->GetBdrElement(be)->Duplicate(new_mesh.get()); - new_mesh->AddBdrElement(ne); + mfem::Element *ne = orig_mesh.GetBdrElement(be)->Duplicate(&new_mesh); + new_mesh.AddBdrElement(ne); } } // Add new boundary elements. if (add_bdr || add_subdomain) { + auto FlipVertices = [](mfem::Element *e) { mfem::Array v; e->GetVertices(v); - int start = 0, end = v.Size() - 1; - while (start < end) - { - int t = v[start]; - v[start] = v[end]; - v[end] = t; - start++; - end--; - } + std::reverse(v.begin(), v.end()); e->SetVertices(v.HostRead()); }; // 1-based, some boundary attributes may be empty since they were removed from the // original mesh, but to keep indices the same as config file we don't compact the // list. - int max_bdr_attr = - orig_mesh->bdr_attributes.Size() ? orig_mesh->bdr_attributes.Max() : 0; - for (int f = 0; f < orig_mesh->GetNumFaces(); f++) + int max_bdr_attr = orig_mesh.bdr_attributes.Size() ? orig_mesh.bdr_attributes.Max() : 0; + for (int f = 0; f < orig_mesh.GetNumFaces(); f++) { if (add_bdr_faces[f] > 0) { @@ -1505,23 +1551,23 @@ std::map> CheckMesh(std::unique_ptr &orig_me // inverse so that the attributes of e1 and e2 can be easily referenced using the // new attribute. Since attributes are in 1-based indexing, a, b > 0. int e1, e2, a = 0, b = 0; - orig_mesh->GetFaceElements(f, &e1, &e2); + orig_mesh.GetFaceElements(f, &e1, &e2); bool no_e1 = (e1 < 0 || elem_delete[e1]); bool no_e2 = (e2 < 0 || elem_delete[e2]); if (!no_e1 && !no_e2) { - a = std::max(orig_mesh->GetAttribute(e1), orig_mesh->GetAttribute(e2)); - b = (a == orig_mesh->GetAttribute(e1)) ? orig_mesh->GetAttribute(e2) - : orig_mesh->GetAttribute(e1); + a = std::max(orig_mesh.GetAttribute(e1), orig_mesh.GetAttribute(e2)); + b = (a == orig_mesh.GetAttribute(e1)) ? orig_mesh.GetAttribute(e2) + : orig_mesh.GetAttribute(e1); } else if (!no_e1) { - a = orig_mesh->GetAttribute(e1); + a = orig_mesh.GetAttribute(e1); b = 0; } else if (!no_e2) { - a = orig_mesh->GetAttribute(e2); + a = orig_mesh.GetAttribute(e2); b = 0; } MFEM_VERIFY(a + b > 0, "Invalid new boundary element attribute!"); @@ -1533,16 +1579,16 @@ std::map> CheckMesh(std::unique_ptr &orig_me } // Add the boundary elements with the new boundary attribute. - mfem::Element *ne = orig_mesh->GetFace(f)->Duplicate(new_mesh.get()); + mfem::Element *ne = orig_mesh.GetFace(f)->Duplicate(&new_mesh); ne->SetAttribute(new_attr); - new_mesh->AddBdrElement(ne); + new_mesh.AddBdrElement(ne); if (add_bdr_faces[f] > 1) { // Flip order of vertices to reverse normal direction of second added element. - ne = orig_mesh->GetFace(f)->Duplicate(new_mesh.get()); + ne = orig_mesh.GetFace(f)->Duplicate(&new_mesh); FlipVertices(ne); ne->SetAttribute(new_attr); - new_mesh->AddBdrElement(ne); + new_mesh.AddBdrElement(ne); // Mpi::Print("Adding two BE with attr {:d} from elements {:d} and {:d}\n", // new_attr, a, b); } @@ -1556,21 +1602,21 @@ std::map> CheckMesh(std::unique_ptr &orig_me // reference. After we have copied the high-order nodes information, topological changes // in Mesh::Finalize are OK (with refine = true). constexpr bool generate_bdr = false, refine = true, fix_orientation = true; - new_mesh->FinalizeTopology(generate_bdr); - new_mesh->RemoveUnusedVertices(); - if (orig_mesh->GetNodes()) + new_mesh.FinalizeTopology(generate_bdr); + new_mesh.RemoveUnusedVertices(); + if (orig_mesh.GetNodes()) { - const mfem::GridFunction *nodes = orig_mesh->GetNodes(); + const mfem::GridFunction *nodes = orig_mesh.GetNodes(); const mfem::FiniteElementSpace *fespace = nodes->FESpace(); mfem::Ordering::Type ordering = fespace->GetOrdering(); int order = fespace->GetMaxElementOrder(); - int sdim = orig_mesh->SpaceDimension(); + int sdim = orig_mesh.SpaceDimension(); bool discont = dynamic_cast(fespace->FEColl()) != nullptr; - new_mesh->SetCurvature(order, discont, sdim, ordering); - mfem::GridFunction *new_nodes = new_mesh->GetNodes(); + new_mesh.SetCurvature(order, discont, sdim, ordering); + mfem::GridFunction *new_nodes = new_mesh.GetNodes(); const mfem::FiniteElementSpace *new_fespace = new_nodes->FESpace(); // The element loop works because we know the mapping from old_mesh to new_mesh element @@ -1578,7 +1624,7 @@ std::map> CheckMesh(std::unique_ptr &orig_me mfem::Array vdofs, new_vdofs; mfem::Vector loc_vec; int te = 0; - for (int e = 0; e < orig_mesh->GetNE(); e++) + for (int e = 0; e < orig_mesh.GetNE(); e++) { if (!elem_delete[e]) { @@ -1590,14 +1636,14 @@ std::map> CheckMesh(std::unique_ptr &orig_me } } } - new_mesh->Finalize(refine, fix_orientation); + new_mesh.Finalize(refine, fix_orientation); orig_mesh = std::move(new_mesh); return new_attr_map; } std::unique_ptr DistributeMesh(MPI_Comm comm, std::unique_ptr &smesh, - std::unique_ptr &partitioning, + const std::unique_ptr &partitioning, const std::string &output_dir) { // Take a serial mesh and partitioning on the root process and construct the global @@ -1605,107 +1651,112 @@ std::unique_ptr DistributeMesh(MPI_Comm comm, // pass arguments to ensure no topological changes (this isn't required since the serial // mesh was marked for refinement). constexpr bool generate_edges = true, refine = true, fix_orientation = true; -#if 0 - // Write each processor's component to file. - std::string tmp = output_dir; - if (tmp.back() != '/') + if constexpr (false) { - tmp += '/'; - } - tmp += "tmp/"; - int width = 1 + static_cast(std::log10(Mpi::Size(comm) - 1)); - if (Mpi::Root(comm)) - { - if (!std::filesystem::exists(tmp)) + // Write each processor's component to file. + std::string tmp = output_dir; + if (tmp.back() != '/') { - std::filesystem::create_directories(tmp); + tmp += '/'; } - mfem::MeshPartitioner partitioner(*smesh, Mpi::Size(comm), partitioning.get()); - for (int i = 0; i < Mpi::Size(comm); i++) + tmp += "tmp/"; + int width = 1 + static_cast(std::log10(Mpi::Size(comm) - 1)); + if (Mpi::Root(comm)) { - mfem::MeshPart part; - partitioner.ExtractPart(i, part); - std::string pfile = mfem::MakeParFilename(tmp + "part.", i, ".mesh", width); - std::ofstream fo(pfile); - // mfem::ofgzstream fo(pfile, true); // Use zlib compression if available - // fo << std::fixed; - fo << std::scientific; - fo.precision(MSH_FLT_PRECISION); - part.Print(fo); + if (!std::filesystem::exists(tmp)) + { + std::filesystem::create_directories(tmp); + } + mfem::MeshPartitioner partitioner(*smesh, Mpi::Size(comm), partitioning.get()); + for (int i = 0; i < Mpi::Size(comm); i++) + { + mfem::MeshPart part; + partitioner.ExtractPart(i, part); + std::string pfile = mfem::MakeParFilename(tmp + "part.", i, ".mesh", width); + std::ofstream fo(pfile); + // mfem::ofgzstream fo(pfile, true); // Use zlib compression if available + // fo << std::fixed; + fo << std::scientific; + fo.precision(MSH_FLT_PRECISION); + part.Print(fo); + } } - } - // Each process loads its own partitioned mesh file and constructs the parallel mesh. - std::string pfile = mfem::MakeParFilename(tmp + "part.", Mpi::Rank(comm), ".mesh", width); - int exists = 0; - while (!exists) // Wait for root to finish writing all files - { - exists = std::filesystem::exists(pfile); - Mpi::GlobalMax(1, &exists, comm); - } - std::ifstream fi(pfile); - // mfem::ifgzstream fi(pfile); - if (!fi.good()) - { - MFEM_ABORT("Unable to open partitioned mesh file \"" << pfile << "\"!"); - } - auto pmesh = - std::make_unique(comm, fi, generate_edges, refine, fix_orientation); - Mpi::Barrier(comm); - if (Mpi::Root(comm)) - { - std::filesystem::remove_all(tmp); // Remove the temporary directory - } - return pmesh; -#endif - // Send each processor's component as a byte string. - std::unique_ptr pmesh; - if (Mpi::Root(comm)) - { - mfem::MeshPartitioner partitioner(*smesh, Mpi::Size(comm), partitioning.get()); - std::vector send_requests(Mpi::Size(comm) - 1, MPI_REQUEST_NULL); - std::vector so; - so.reserve(Mpi::Size(comm)); - for (int i = 0; i < Mpi::Size(comm); i++) + // Each process loads its own partitioned mesh file and constructs the parallel mesh. + std::string pfile = + mfem::MakeParFilename(tmp + "part.", Mpi::Rank(comm), ".mesh", width); + int exists = 0; + while (!exists) // Wait for root to finish writing all files { - mfem::MeshPart part; - partitioner.ExtractPart(i, part); - std::ostringstream fo(std::stringstream::out); - // fo << std::fixed; - fo << std::scientific; - fo.precision(MSH_FLT_PRECISION); - part.Print(fo); - so.push_back(fo.str()); - // so.push_back((i > 0) ? zlib::CompressString(fo.str()) : fo.str()); - if (i > 0) - { - int slen = static_cast(so[i].length()); - MFEM_VERIFY(so[i].length() == (std::size_t)slen, - "Overflow error distributing parallel mesh!"); - MPI_Isend(so[i].data(), slen, MPI_CHAR, i, i, comm, &send_requests[i - 1]); - } + exists = std::filesystem::exists(pfile); + Mpi::GlobalMax(1, &exists, comm); + } + std::ifstream fi(pfile); + // mfem::ifgzstream fi(pfile); + if (!fi.good()) + { + MFEM_ABORT("Unable to open partitioned mesh file \"" << pfile << "\"!"); } - std::istringstream fi(so[0]); // This is never compressed - pmesh = + auto pmesh = std::make_unique(comm, fi, generate_edges, refine, fix_orientation); - MPI_Waitall(static_cast(send_requests.size()), send_requests.data(), - MPI_STATUSES_IGNORE); + Mpi::Barrier(comm); + if (Mpi::Root(comm)) + { + std::filesystem::remove_all(tmp); // Remove the temporary directory + } + return pmesh; } else { - MPI_Status status; - int rlen; - std::string si; - MPI_Probe(0, Mpi::Rank(comm), comm, &status); - MPI_Get_count(&status, MPI_CHAR, &rlen); - si.resize(rlen); - MPI_Recv(si.data(), rlen, MPI_CHAR, 0, Mpi::Rank(comm), comm, MPI_STATUS_IGNORE); - std::istringstream fi(si); - // std::istringstream fi(zlib::DecompressString(si)); - pmesh = - std::make_unique(comm, fi, generate_edges, refine, fix_orientation); + // Send each processor's component as a byte string. + std::unique_ptr pmesh; + if (Mpi::Root(comm)) + { + mfem::MeshPartitioner partitioner(*smesh, Mpi::Size(comm), partitioning.get()); + std::vector send_requests(Mpi::Size(comm) - 1, MPI_REQUEST_NULL); + std::vector so; + so.reserve(Mpi::Size(comm)); + for (int i = 0; i < Mpi::Size(comm); i++) + { + mfem::MeshPart part; + partitioner.ExtractPart(i, part); + std::ostringstream fo(std::stringstream::out); + // fo << std::fixed; + fo << std::scientific; + fo.precision(MSH_FLT_PRECISION); + part.Print(fo); + so.push_back(fo.str()); + // so.push_back((i > 0) ? zlib::CompressString(fo.str()) : fo.str()); + if (i > 0) + { + int slen = static_cast(so[i].length()); + MFEM_VERIFY(so[i].length() == (std::size_t)slen, + "Overflow error distributing parallel mesh!"); + MPI_Isend(so[i].data(), slen, MPI_CHAR, i, i, comm, &send_requests[i - 1]); + } + } + std::istringstream fi(so[0]); // This is never compressed + pmesh = std::make_unique(comm, fi, generate_edges, refine, + fix_orientation); + MPI_Waitall(static_cast(send_requests.size()), send_requests.data(), + MPI_STATUSES_IGNORE); + } + else + { + MPI_Status status; + int rlen; + std::string si; + MPI_Probe(0, Mpi::Rank(comm), comm, &status); + MPI_Get_count(&status, MPI_CHAR, &rlen); + si.resize(rlen); + MPI_Recv(si.data(), rlen, MPI_CHAR, 0, Mpi::Rank(comm), comm, MPI_STATUS_IGNORE); + std::istringstream fi(si); + // std::istringstream fi(zlib::DecompressString(si)); + pmesh = std::make_unique(comm, fi, generate_edges, refine, + fix_orientation); + } + return pmesh; } - return pmesh; } void GetUsedAttributeMarkers(const IoData &iodata, int n_mat, int n_bdr, diff --git a/scripts/schema/config/model.json b/scripts/schema/config/model.json index dd7c0055a..1ff8204a1 100644 --- a/scripts/schema/config/model.json +++ b/scripts/schema/config/model.json @@ -20,6 +20,25 @@ "properties": { "UniformLevels": { "type": "integer", "minimum": 0 }, + "Adaptation": + { + "type": "object", + "additionalProperties": false, + "required": ["MaxIts"], + "properties": + { + "Tol": {"type": "number", "exclusiveMinimum": 0.0}, + "MaxIts": {"type": "integer", "inclusiveMinimum": 0}, + "MinIts": {"type": "integer", "inclusiveMinimum": 0}, + "UpdateFraction": {"type": "number", "exclusiveMinimum": 0.0, "exclusiveMaximum": 1.0}, + "CoarseningFraction": {"type": "number", "exclusiveMaximum": 1.0, "inclusiveMinimum": 0.0}, + "DOFLimit": {"type": "number", "inclusiveMinimum": 0}, + "SaveStep": {"type": "integer", "inclusiveMinimum": 0}, + "MaxNCLevels": {"type": "integer", "inclusiveMinimum": 0}, + "Nonconformal": {"type": "boolean"}, + "MaximumImbalance": {"type": "number", "inclusiveMinimum": 1.0} + } + }, "Boxes": { "type": "array",