From 5cd505dd9afa17d3b0f6095d018cd7a5e70b2a8e Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 15 Apr 2024 17:55:01 -0600 Subject: [PATCH 01/22] unified initializer list --- src/mesh/mesh.cpp | 58 ++++++++--------------------------------------- src/mesh/mesh.hpp | 4 +++- 2 files changed, 12 insertions(+), 50 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 5a105dd3aaad..3d74fb9add8a 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -57,11 +57,8 @@ namespace parthenon { -//---------------------------------------------------------------------------------------- -// Mesh constructor, builds mesh at start of calculation using parameters in input file - Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, - int mesh_test) + private_t) : // public members: modified(true), is_restart(false), // aggregate initialization of RegionSize struct: @@ -103,7 +100,13 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // private members: num_mesh_threads_(pin->GetOrAddInteger("parthenon/mesh", "num_threads", 1)), use_uniform_meshgen_fn_{true, true, true, true}, lb_flag_(true), lb_automatic_(), - lb_manual_(), MeshBndryFnctn{nullptr, nullptr, nullptr, nullptr, nullptr, nullptr} { + lb_manual_(), MeshBndryFnctn{nullptr, nullptr, nullptr, nullptr, nullptr, nullptr} {} +//---------------------------------------------------------------------------------------- +// Mesh constructor, builds mesh at start of calculation using parameters in input file + +Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, + int mesh_test) + : Mesh(pin, app_in, packages, private_t()) { std::stringstream msg; RegionSize block_size; BoundaryFlag block_bcs[6]; @@ -434,50 +437,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // Mesh constructor for restarts. Load the restart file Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, Packages_t &packages, int mesh_test) - : // public members: - // aggregate initialization of RegionSize struct: - // (will be overwritten by memcpy from restart file, in this case) - modified(true), is_restart(true), - // aggregate initialization of RegionSize struct: - mesh_size({pin->GetReal("parthenon/mesh", "x1min"), - pin->GetReal("parthenon/mesh", "x2min"), - pin->GetReal("parthenon/mesh", "x3min")}, - {pin->GetReal("parthenon/mesh", "x1max"), - pin->GetReal("parthenon/mesh", "x2max"), - pin->GetReal("parthenon/mesh", "x3max")}, - {pin->GetOrAddReal("parthenon/mesh", "x1rat", 1.0), - pin->GetOrAddReal("parthenon/mesh", "x2rat", 1.0), - pin->GetOrAddReal("parthenon/mesh", "x3rat", 1.0)}, - {pin->GetInteger("parthenon/mesh", "nx1"), - pin->GetInteger("parthenon/mesh", "nx2"), - pin->GetInteger("parthenon/mesh", "nx3")}, - {false, pin->GetInteger("parthenon/mesh", "nx2") == 1, - pin->GetInteger("parthenon/mesh", "nx3") == 1}), - mesh_bcs{ - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix1_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox1_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix2_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox2_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix3_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox3_bc", "reflecting"))}, - ndim((mesh_size.nx(X3DIR) > 1) ? 3 : ((mesh_size.nx(X2DIR) > 1) ? 2 : 1)), - adaptive(pin->GetOrAddString("parthenon/mesh", "refinement", "none") == "adaptive" - ? true - : false), - multilevel( - (adaptive || - pin->GetOrAddString("parthenon/mesh", "refinement", "none") == "static" || - pin->GetOrAddString("parthenon/mesh", "multigrid", "false") == "true") - ? true - : false), - multigrid(pin->GetOrAddString("parthenon/mesh", "multigrid", "false") == "true" - ? true - : false), - nbnew(), nbdel(), step_since_lb(), gflag(), packages(packages), - // private members: - num_mesh_threads_(pin->GetOrAddInteger("parthenon/mesh", "num_threads", 1)), - use_uniform_meshgen_fn_{true, true, true, true}, lb_flag_(true), lb_automatic_(), - lb_manual_(), MeshBndryFnctn{nullptr, nullptr, nullptr, nullptr, nullptr, nullptr} { + : Mesh(pin, app_in, packages, private_t()) { std::stringstream msg; RegionSize block_size; BoundaryFlag block_bcs[6]; diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 0aeebe8faeef..018f7b675337 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -76,9 +76,11 @@ class Mesh { friend class MeshBlock; friend class MeshBlockTree; friend class MeshRefinement; - + + struct private_t{}; public: // 2x function overloads of ctor: normal and restarted simulation + Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, private_t); Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, int test_flag = 0); Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &resfile, From ed404d3bc3f404d23997a3b46dc3624586304e04 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 15 Apr 2024 17:59:02 -0600 Subject: [PATCH 02/22] works... --- src/mesh/mesh.cpp | 38 +++++++++++--------------------------- 1 file changed, 11 insertions(+), 27 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 3d74fb9add8a..c5b8c1b9a700 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -100,7 +100,17 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // private members: num_mesh_threads_(pin->GetOrAddInteger("parthenon/mesh", "num_threads", 1)), use_uniform_meshgen_fn_{true, true, true, true}, lb_flag_(true), lb_automatic_(), - lb_manual_(), MeshBndryFnctn{nullptr, nullptr, nullptr, nullptr, nullptr, nullptr} {} + lb_manual_(), MeshBndryFnctn{nullptr, nullptr, nullptr, nullptr, nullptr, nullptr}, + nslist(Globals::nranks), + nblist(Globals::nranks), + nref(Globals::nranks), + nderef(Globals::nranks), + rdisp(Globals::nranks), + ddisp(Globals::nranks), + bnref(Globals::nranks), + bnderef(Globals::nranks), + brdisp(Globals::nranks), + bddisp(Globals::nranks){} //---------------------------------------------------------------------------------------- // Mesh constructor, builds mesh at start of calculation using parameters in input file @@ -378,19 +388,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, ranklist = std::vector(nbtotal); - nslist = std::vector(Globals::nranks); - nblist = std::vector(Globals::nranks); - if (adaptive) { // allocate arrays for AMR - nref = std::vector(Globals::nranks); - nderef = std::vector(Globals::nranks); - rdisp = std::vector(Globals::nranks); - ddisp = std::vector(Globals::nranks); - bnref = std::vector(Globals::nranks); - bnderef = std::vector(Globals::nranks); - brdisp = std::vector(Globals::nranks); - bddisp = std::vector(Globals::nranks); - } - // initialize cost array with the simplest estimate; all the blocks are equal costlist = std::vector(nbtotal, 1.0); @@ -588,19 +585,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, #endif costlist = std::vector(nbtotal, 1.0); ranklist = std::vector(nbtotal); - nslist = std::vector(Globals::nranks); - nblist = std::vector(Globals::nranks); - - if (adaptive) { // allocate arrays for AMR - nref = std::vector(Globals::nranks); - nderef = std::vector(Globals::nranks); - rdisp = std::vector(Globals::nranks); - ddisp = std::vector(Globals::nranks); - bnref = std::vector(Globals::nranks); - bnderef = std::vector(Globals::nranks); - brdisp = std::vector(Globals::nranks); - bddisp = std::vector(Globals::nranks); - } CalculateLoadBalance(costlist, ranklist, nslist, nblist); From 216f9963efb9d7ec642409af78fa39619280b7ff Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 15 Apr 2024 18:11:56 -0600 Subject: [PATCH 03/22] working --- src/mesh/mesh.cpp | 77 ++++++++++++++++++++++++----------------------- src/mesh/mesh.hpp | 1 + 2 files changed, 40 insertions(+), 38 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index c5b8c1b9a700..64b55ca95831 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -110,7 +110,21 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, bnref(Globals::nranks), bnderef(Globals::nranks), brdisp(Globals::nranks), - bddisp(Globals::nranks){} + bddisp(Globals::nranks) { + for (auto &[dir, label] : std::vector>{ + {X1DIR, "nx1"}, {X2DIR, "nx2"}, {X3DIR, "nx3"}}) { + base_block_size.xrat(dir) = mesh_size.xrat(dir); + base_block_size.symmetry(dir) = mesh_size.symmetry(dir); + if (!base_block_size.symmetry(dir)) { + base_block_size.nx(dir) = + pin->GetOrAddInteger("parthenon/meshblock", label, mesh_size.nx(dir)); + } else { + base_block_size.nx(dir) = mesh_size.nx(dir); + } + nrbx[dir - 1] = mesh_size.nx(dir) / base_block_size.nx(dir); + } +} + //---------------------------------------------------------------------------------------- // Mesh constructor, builds mesh at start of calculation using parameters in input file @@ -119,7 +133,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, : Mesh(pin, app_in, packages, private_t()) { std::stringstream msg; RegionSize block_size; - BoundaryFlag block_bcs[6]; // mesh test if (mesh_test > 0) Globals::nranks = mesh_test; @@ -211,30 +224,18 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, } EnrollBndryFncts_(app_in); - for (auto &[dir, label] : std::vector>{ - {X1DIR, "nx1"}, {X2DIR, "nx2"}, {X3DIR, "nx3"}}) { - block_size.xrat(dir) = mesh_size.xrat(dir); - block_size.symmetry(dir) = mesh_size.symmetry(dir); - if (!block_size.symmetry(dir)) { - block_size.nx(dir) = - pin->GetOrAddInteger("parthenon/meshblock", label, mesh_size.nx(dir)); - } else { - block_size.nx(dir) = mesh_size.nx(dir); - } - nrbx[dir - 1] = mesh_size.nx(dir) / block_size.nx(dir); - } - base_block_size = block_size; + // check consistency of the block and mesh - if (mesh_size.nx(X1DIR) % block_size.nx(X1DIR) != 0 || - mesh_size.nx(X2DIR) % block_size.nx(X2DIR) != 0 || - mesh_size.nx(X3DIR) % block_size.nx(X3DIR) != 0) { + if (mesh_size.nx(X1DIR) % base_block_size.nx(X1DIR) != 0 || + mesh_size.nx(X2DIR) % base_block_size.nx(X2DIR) != 0 || + mesh_size.nx(X3DIR) % base_block_size.nx(X3DIR) != 0) { msg << "### FATAL ERROR in Mesh constructor" << std::endl << "the Mesh must be evenly divisible by the MeshBlock" << std::endl; PARTHENON_FAIL(msg); } - if (block_size.nx(X1DIR) < 4 || (block_size.nx(X2DIR) < 4 && (ndim >= 2)) || - (block_size.nx(X3DIR) < 4 && (ndim >= 3))) { + if (base_block_size.nx(X1DIR) < 4 || (base_block_size.nx(X2DIR) < 4 && (ndim >= 2)) || + (base_block_size.nx(X3DIR) < 4 && (ndim >= 3))) { msg << "### FATAL ERROR in Mesh constructor" << std::endl << "block_size must be larger than or equal to 4 cells." << std::endl; PARTHENON_FAIL(msg); @@ -243,7 +244,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // initialize user-enrollable functions default_pack_size_ = pin->GetOrAddInteger("parthenon/mesh", "pack_size", -1); - forest = forest::Forest::HyperRectangular(mesh_size, block_size, mesh_bcs); + forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); root_level = forest.root_level; current_level = root_level; @@ -266,8 +267,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, InitUserMeshData(this, pin); if (multilevel) { - if (block_size.nx(X1DIR) % 2 == 1 || (block_size.nx(X2DIR) % 2 == 1 && (ndim >= 2)) || - (block_size.nx(X3DIR) % 2 == 1 && (ndim >= 3))) { + if (base_block_size.nx(X1DIR) % 2 == 1 || (base_block_size.nx(X2DIR) % 2 == 1 && (ndim >= 2)) || + (base_block_size.nx(X3DIR) % 2 == 1 && (ndim >= 3))) { msg << "### FATAL ERROR in Mesh constructor" << std::endl << "The size of MeshBlock must be divisible by 2 in order to use SMR or AMR." << std::endl; @@ -418,6 +419,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, block_list.clear(); block_list.resize(nbe - nbs + 1); for (int i = nbs; i <= nbe; i++) { + RegionSize block_size; + BoundaryFlag block_bcs[6]; SetBlockSizeAndBoundaries(loclist[i], block_size, block_bcs); // create a block and add into the link list block_list[i - nbs] = @@ -436,8 +439,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, Packages_t &packages, int mesh_test) : Mesh(pin, app_in, packages, private_t()) { std::stringstream msg; - RegionSize block_size; - BoundaryFlag block_bcs[6]; // mesh test if (mesh_test > 0) Globals::nranks = mesh_test; @@ -460,11 +461,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, nbtotal = mesh_info.nbtotal; root_level = mesh_info.root_level; - const auto bc = mesh_info.bound_cond; - for (int i = 0; i < 6; i++) { - block_bcs[i] = GetBoundaryFlag(bc[i]); - } - // Allow for user overrides to default Parthenon functions if (app_in->InitUserMeshData != nullptr) { InitUserMeshData = app_in->InitUserMeshData; @@ -503,21 +499,20 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, loclist = std::vector(nbtotal); for (auto &dir : {X1DIR, X2DIR, X3DIR}) { - block_size.xrat(dir) = mesh_size.xrat(dir); - block_size.nx(dir) = mesh_info.block_size[dir - 1] - + base_block_size.xrat(dir) = mesh_size.xrat(dir); + base_block_size.nx(dir) = mesh_info.block_size[dir - 1] - (mesh_info.block_size[dir - 1] > 1) * mesh_info.includes_ghost * 2 * mesh_info.n_ghost; - if (block_size.nx(dir) == 1) { - block_size.symmetry(dir) = true; + if (base_block_size.nx(dir) == 1) { + base_block_size.symmetry(dir) = true; mesh_size.symmetry(dir) = true; } else { - block_size.symmetry(dir) = false; + base_block_size.symmetry(dir) = false; mesh_size.symmetry(dir) = false; } // calculate the number of the blocks - nrbx[dir - 1] = mesh_size.nx(dir) / block_size.nx(dir); + nrbx[dir - 1] = mesh_size.nx(dir) / base_block_size.nx(dir); } - base_block_size = block_size; // Load balancing flag and parameters RegisterLoadBalancing_(pin); @@ -553,7 +548,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, } // rebuild the Block Tree - forest = forest::Forest::HyperRectangular(mesh_size, block_size, mesh_bcs); + forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); for (int i = 0; i < nbtotal; i++) { forest.AddMeshBlock(forest.GetForestLocationFromLegacyTreeLocation(loclist[i]), false); @@ -613,6 +608,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, block_list.clear(); block_list.resize(nbe - nbs + 1); for (int i = nbs; i <= nbe; i++) { + RegionSize block_size; + BoundaryFlag block_bcs[6]; for (auto &v : block_bcs) { v = parthenon::BoundaryFlag::undef; } @@ -1210,4 +1207,8 @@ void Mesh::SetupMPIComms() { #endif } +void Mesh::CheckMeshValidity() const { + +} + } // namespace parthenon diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 018f7b675337..816c38bc9291 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -301,6 +301,7 @@ class Mesh { #endif // functions + void CheckMeshValidity() const; void CalculateLoadBalance(std::vector const &costlist, std::vector &ranklist, std::vector &nslist, std::vector &nblist); From 35d3941e944b33b8b5b103e73d1b20ab31926cf4 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 15 Apr 2024 18:16:31 -0600 Subject: [PATCH 04/22] working --- src/mesh/mesh.cpp | 16 +++------------- 1 file changed, 3 insertions(+), 13 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 64b55ca95831..afeb87f5e2f6 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -499,19 +499,9 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, loclist = std::vector(nbtotal); for (auto &dir : {X1DIR, X2DIR, X3DIR}) { - base_block_size.xrat(dir) = mesh_size.xrat(dir); - base_block_size.nx(dir) = mesh_info.block_size[dir - 1] - - (mesh_info.block_size[dir - 1] > 1) * mesh_info.includes_ghost * - 2 * mesh_info.n_ghost; - if (base_block_size.nx(dir) == 1) { - base_block_size.symmetry(dir) = true; - mesh_size.symmetry(dir) = true; - } else { - base_block_size.symmetry(dir) = false; - mesh_size.symmetry(dir) = false; - } - // calculate the number of the blocks - nrbx[dir - 1] = mesh_size.nx(dir) / base_block_size.nx(dir); + PARTHENON_REQUIRE(base_block_size.nx(dir) == mesh_info.block_size[dir - 1] - + (mesh_info.block_size[dir - 1] > 1) * mesh_info.includes_ghost * + 2 * mesh_info.n_ghost, "Block size not consistent on restart."); } // Load balancing flag and parameters From ec1a9e5a2bbd3d884cf1a4f112ee1f7cb55a293a Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 15 Apr 2024 18:24:24 -0600 Subject: [PATCH 05/22] still works --- src/mesh/mesh.cpp | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index afeb87f5e2f6..cc1321745278 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -111,6 +111,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, bnderef(Globals::nranks), brdisp(Globals::nranks), bddisp(Globals::nranks) { + for (auto &[dir, label] : std::vector>{ {X1DIR, "nx1"}, {X2DIR, "nx2"}, {X3DIR, "nx3"}}) { base_block_size.xrat(dir) = mesh_size.xrat(dir); @@ -123,6 +124,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, } nrbx[dir - 1] = mesh_size.nx(dir) / base_block_size.nx(dir); } + + mesh_data.SetMeshPointer(this); } //---------------------------------------------------------------------------------------- @@ -402,8 +405,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, return; } - mesh_data.SetMeshPointer(this); - resolved_packages = ResolvePackages(packages); // Register user defined boundary conditions @@ -451,6 +452,10 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, PARTHENON_FAIL(msg); } + + forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); + root_level = forest.root_level; + // read the restart file // the file is already open and the pointer is set to after @@ -459,7 +464,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, nbnew = mesh_info.nbnew; nbdel = mesh_info.nbdel; nbtotal = mesh_info.nbtotal; - root_level = mesh_info.root_level; // Allow for user overrides to default Parthenon functions if (app_in->InitUserMeshData != nullptr) { @@ -483,20 +487,17 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, EnrollBndryFncts_(app_in); const auto grid_dim = mesh_info.grid_dim; - mesh_size.xmin(X1DIR) = grid_dim[0]; - mesh_size.xmax(X1DIR) = grid_dim[1]; - mesh_size.xrat(X1DIR) = grid_dim[2]; + PARTHENON_REQUIRE(mesh_size.xmin(X1DIR) == grid_dim[0], "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xmax(X1DIR) == grid_dim[1], "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xrat(X1DIR) == grid_dim[2], "Mesh size shouldn't change on restart."); - mesh_size.xmin(X2DIR) = grid_dim[3]; - mesh_size.xmax(X2DIR) = grid_dim[4]; - mesh_size.xrat(X2DIR) = grid_dim[5]; + PARTHENON_REQUIRE(mesh_size.xmin(X2DIR) == grid_dim[3], "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xmax(X2DIR) == grid_dim[4], "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xrat(X2DIR) == grid_dim[5], "Mesh size shouldn't change on restart."); - mesh_size.xmin(X3DIR) = grid_dim[6]; - mesh_size.xmax(X3DIR) = grid_dim[7]; - mesh_size.xrat(X3DIR) = grid_dim[8]; - - // initialize - loclist = std::vector(nbtotal); + PARTHENON_REQUIRE(mesh_size.xmin(X3DIR) == grid_dim[6], "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xmax(X3DIR) == grid_dim[7], "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xrat(X3DIR) == grid_dim[8], "Mesh size shouldn't change on restart."); for (auto &dir : {X1DIR, X2DIR, X3DIR}) { PARTHENON_REQUIRE(base_block_size.nx(dir) == mesh_info.block_size[dir - 1] - @@ -525,6 +526,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, InitUserMeshData(this, pin); // Populate logical locations + loclist = std::vector(nbtotal); auto lx123 = mesh_info.lx123; auto locLevelGidLidCnghostGflag = mesh_info.level_gid_lid_cnghost_gflag; current_level = -1; @@ -538,11 +540,9 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, } // rebuild the Block Tree - forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); - for (int i = 0; i < nbtotal; i++) { + for (int i = 0; i < nbtotal; i++) forest.AddMeshBlock(forest.GetForestLocationFromLegacyTreeLocation(loclist[i]), false); - } loclist = forest.GetMeshBlockListAndResolveGids(); int nnb = loclist.size(); @@ -584,7 +584,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, int nbs = nslist[Globals::my_rank]; int nbe = nbs + nb - 1; - mesh_data.SetMeshPointer(this); + resolved_packages = ResolvePackages(packages); From ac6cbd17bce2a7017971a5076e108716d9fd6745 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 15 Apr 2024 18:36:42 -0600 Subject: [PATCH 06/22] working... --- src/mesh/mesh.cpp | 272 +++++++++++++++++----------------------------- src/mesh/mesh.hpp | 2 + 2 files changed, 102 insertions(+), 172 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index cc1321745278..49cda0db42a7 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -111,7 +111,35 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, bnderef(Globals::nranks), brdisp(Globals::nranks), bddisp(Globals::nranks) { - + // Allow for user overrides to default Parthenon functions + if (app_in->InitUserMeshData != nullptr) { + InitUserMeshData = app_in->InitUserMeshData; + } + if (app_in->MeshProblemGenerator != nullptr) { + ProblemGenerator = app_in->MeshProblemGenerator; + } + if (app_in->MeshPostInitialization != nullptr) { + PostInitialization = app_in->MeshPostInitialization; + } + if (app_in->PreStepMeshUserWorkInLoop != nullptr) { + PreStepUserWorkInLoop = app_in->PreStepMeshUserWorkInLoop; + } + if (app_in->PostStepMeshUserWorkInLoop != nullptr) { + PostStepUserWorkInLoop = app_in->PostStepMeshUserWorkInLoop; + } + if (app_in->UserMeshWorkBeforeOutput != nullptr) { + UserMeshWorkBeforeOutput = app_in->UserMeshWorkBeforeOutput; + } + if (app_in->PreStepDiagnosticsInLoop != nullptr) { + PreStepUserDiagnosticsInLoop = app_in->PreStepDiagnosticsInLoop; + } + if (app_in->PostStepDiagnosticsInLoop != nullptr) { + PostStepUserDiagnosticsInLoop = app_in->PostStepDiagnosticsInLoop; + } + if (app_in->UserWorkAfterLoop != nullptr) { + UserWorkAfterLoop = app_in->UserWorkAfterLoop; + } + for (auto &[dir, label] : std::vector>{ {X1DIR, "nx1"}, {X2DIR, "nx2"}, {X3DIR, "nx3"}}) { base_block_size.xrat(dir) = mesh_size.xrat(dir); @@ -128,6 +156,74 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, mesh_data.SetMeshPointer(this); } +void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, + int mesh_test) { + std::stringstream msg; + nbtotal = loclist.size(); + current_level = -1; + for (const auto &loc : loclist) + if (loc.level() > current_level) current_level = loc.level(); + +#ifdef MPI_PARALLEL + // check if there are sufficient blocks + if (nbtotal < Globals::nranks) { + if (mesh_test == 0) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks (" + << Globals::nranks << ")" << std::endl; + PARTHENON_FAIL(msg); + } else { // test + std::cout << "### Warning in Mesh constructor" << std::endl + << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks (" + << Globals::nranks << ")" << std::endl; + } + } +#endif + + ranklist = std::vector(nbtotal); + + // initialize cost array with the simplest estimate; all the blocks are equal + costlist = std::vector(nbtotal, 1.0); + + CalculateLoadBalance(costlist, ranklist, nslist, nblist); + + // Output some diagnostic information to terminal + + // Output MeshBlock list and quit (mesh test only); do not create meshes + if (mesh_test > 0) { + if (Globals::my_rank == 0) OutputMeshStructure(ndim); + return; + } + + resolved_packages = ResolvePackages(packages); + + // Register user defined boundary conditions + UserBoundaryFunctions = resolved_packages->UserBoundaryFunctions; + + // Setup unique comms for each variable and swarm + SetupMPIComms(); + + // create MeshBlock list for this process + int nbs = nslist[Globals::my_rank]; + int nbe = nbs + nblist[Globals::my_rank] - 1; + // create MeshBlock list for this process + block_list.clear(); + block_list.resize(nbe - nbs + 1); + for (int i = nbs; i <= nbe; i++) { + RegionSize block_size; + BoundaryFlag block_bcs[6]; + SetBlockSizeAndBoundaries(loclist[i], block_size, block_bcs); + // create a block and add into the link list + block_list[i - nbs] = + MeshBlock::Make(i, i - nbs, loclist[i], block_size, block_bcs, this, pin, app_in, + packages, resolved_packages, gflag, costlist[i]); + } + BuildGMGBlockLists(pin, app_in); + SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist); + SetGMGNeighbors(); + ResetLoadBalanceVariables(); +} + //---------------------------------------------------------------------------------------- // Mesh constructor, builds mesh at start of calculation using parameters in input file @@ -176,34 +272,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, PARTHENON_FAIL(msg); } - // Allow for user overrides to default Parthenon functions - if (app_in->InitUserMeshData != nullptr) { - InitUserMeshData = app_in->InitUserMeshData; - } - if (app_in->MeshProblemGenerator != nullptr) { - ProblemGenerator = app_in->MeshProblemGenerator; - } - if (app_in->MeshPostInitialization != nullptr) { - PostInitialization = app_in->MeshPostInitialization; - } - if (app_in->PreStepMeshUserWorkInLoop != nullptr) { - PreStepUserWorkInLoop = app_in->PreStepMeshUserWorkInLoop; - } - if (app_in->PostStepMeshUserWorkInLoop != nullptr) { - PostStepUserWorkInLoop = app_in->PostStepMeshUserWorkInLoop; - } - if (app_in->UserMeshWorkBeforeOutput != nullptr) { - UserMeshWorkBeforeOutput = app_in->UserMeshWorkBeforeOutput; - } - if (app_in->PreStepDiagnosticsInLoop != nullptr) { - PreStepUserDiagnosticsInLoop = app_in->PreStepDiagnosticsInLoop; - } - if (app_in->PostStepDiagnosticsInLoop != nullptr) { - PostStepUserDiagnosticsInLoop = app_in->PostStepDiagnosticsInLoop; - } - if (app_in->UserWorkAfterLoop != nullptr) { - UserWorkAfterLoop = app_in->UserWorkAfterLoop; - } + // check the consistency of the periodic boundaries if (((mesh_bcs[BoundaryFace::inner_x1] == BoundaryFlag::periodic && @@ -369,69 +438,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // initial mesh hierarchy construction is completed here loclist = forest.GetMeshBlockListAndResolveGids(); - nbtotal = loclist.size(); - current_level = -1; - for (const auto &loc : loclist) { - if (loc.level() > current_level) current_level = loc.level(); - } -#ifdef MPI_PARALLEL - // check if there are sufficient blocks - if (nbtotal < Globals::nranks) { - if (mesh_test == 0) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks (" - << Globals::nranks << ")" << std::endl; - PARTHENON_FAIL(msg); - } else { // test - std::cout << "### Warning in Mesh constructor" << std::endl - << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks (" - << Globals::nranks << ")" << std::endl; - } - } -#endif - - ranklist = std::vector(nbtotal); - - // initialize cost array with the simplest estimate; all the blocks are equal - costlist = std::vector(nbtotal, 1.0); - - CalculateLoadBalance(costlist, ranklist, nslist, nblist); - - // Output some diagnostic information to terminal - - // Output MeshBlock list and quit (mesh test only); do not create meshes - if (mesh_test > 0) { - if (Globals::my_rank == 0) OutputMeshStructure(ndim); - return; - } - - resolved_packages = ResolvePackages(packages); - - // Register user defined boundary conditions - UserBoundaryFunctions = resolved_packages->UserBoundaryFunctions; - - // Setup unique comms for each variable and swarm - SetupMPIComms(); - - // create MeshBlock list for this process - int nbs = nslist[Globals::my_rank]; - int nbe = nbs + nblist[Globals::my_rank] - 1; - // create MeshBlock list for this process - block_list.clear(); - block_list.resize(nbe - nbs + 1); - for (int i = nbs; i <= nbe; i++) { - RegionSize block_size; - BoundaryFlag block_bcs[6]; - SetBlockSizeAndBoundaries(loclist[i], block_size, block_bcs); - // create a block and add into the link list - block_list[i - nbs] = - MeshBlock::Make(i, i - nbs, loclist[i], block_size, block_bcs, this, pin, app_in, - packages, resolved_packages, gflag); - } - BuildGMGBlockLists(pin, app_in); - SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist); - SetGMGNeighbors(); - ResetLoadBalanceVariables(); + BuildBlockList(pin, app_in, packages, mesh_test); } //---------------------------------------------------------------------------------------- @@ -465,25 +472,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, nbdel = mesh_info.nbdel; nbtotal = mesh_info.nbtotal; - // Allow for user overrides to default Parthenon functions - if (app_in->InitUserMeshData != nullptr) { - InitUserMeshData = app_in->InitUserMeshData; - } - if (app_in->PreStepMeshUserWorkInLoop != nullptr) { - PreStepUserWorkInLoop = app_in->PreStepMeshUserWorkInLoop; - } - if (app_in->PostStepMeshUserWorkInLoop != nullptr) { - PostStepUserWorkInLoop = app_in->PostStepMeshUserWorkInLoop; - } - if (app_in->PreStepDiagnosticsInLoop != nullptr) { - PreStepUserDiagnosticsInLoop = app_in->PreStepDiagnosticsInLoop; - } - if (app_in->PostStepDiagnosticsInLoop != nullptr) { - PostStepUserDiagnosticsInLoop = app_in->PostStepDiagnosticsInLoop; - } - if (app_in->UserWorkAfterLoop != nullptr) { - UserWorkAfterLoop = app_in->UserWorkAfterLoop; - } EnrollBndryFncts_(app_in); const auto grid_dim = mesh_info.grid_dim; @@ -553,67 +541,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, PARTHENON_FAIL(msg); } -#ifdef MPI_PARALLEL - if (nbtotal < Globals::nranks) { - if (mesh_test == 0) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks (" - << Globals::nranks << ")" << std::endl; - PARTHENON_FAIL(msg); - } else { // test - std::cout << "### Warning in Mesh constructor" << std::endl - << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks (" - << Globals::nranks << ")" << std::endl; - return; - } - } -#endif - costlist = std::vector(nbtotal, 1.0); - ranklist = std::vector(nbtotal); - - CalculateLoadBalance(costlist, ranklist, nslist, nblist); - - // Output MeshBlock list and quit (mesh test only); do not create meshes - if (mesh_test > 0) { - if (Globals::my_rank == 0) OutputMeshStructure(ndim); - return; - } - - // allocate data buffer - int nb = nblist[Globals::my_rank]; - int nbs = nslist[Globals::my_rank]; - int nbe = nbs + nb - 1; - - - - resolved_packages = ResolvePackages(packages); - - // Register user defined boundary conditions - UserBoundaryFunctions = resolved_packages->UserBoundaryFunctions; - - // Setup unique comms for each variable and swarm - SetupMPIComms(); - - // Create MeshBlocks (parallel) - block_list.clear(); - block_list.resize(nbe - nbs + 1); - for (int i = nbs; i <= nbe; i++) { - RegionSize block_size; - BoundaryFlag block_bcs[6]; - for (auto &v : block_bcs) { - v = parthenon::BoundaryFlag::undef; - } - SetBlockSizeAndBoundaries(loclist[i], block_size, block_bcs); - - // create a block and add into the link list - block_list[i - nbs] = - MeshBlock::Make(i, i - nbs, loclist[i], block_size, block_bcs, this, pin, app_in, - packages, resolved_packages, gflag, costlist[i]); - } - BuildGMGBlockLists(pin, app_in); - SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist); - SetGMGNeighbors(); - ResetLoadBalanceVariables(); + BuildBlockList(pin, app_in, packages, mesh_test); } //---------------------------------------------------------------------------------------- diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 816c38bc9291..deaf90956416 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -302,6 +302,8 @@ class Mesh { // functions void CheckMeshValidity() const; + void BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, + int mesh_test); void CalculateLoadBalance(std::vector const &costlist, std::vector &ranklist, std::vector &nslist, std::vector &nblist); From 09a0b2ef16b5ceff19fedf253b6e4d3ae547fd7b Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 15 Apr 2024 18:42:54 -0600 Subject: [PATCH 07/22] factor out mesh validity --- src/mesh/mesh.cpp | 205 ++++++++++++++++++++++------------------------ 1 file changed, 97 insertions(+), 108 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 49cda0db42a7..2e674c4ac98a 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -152,8 +152,17 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, } nrbx[dir - 1] = mesh_size.nx(dir) / base_block_size.nx(dir); } + + // SMR / AMR: + if (adaptive) { + max_level = pin->GetOrAddInteger("parthenon/mesh", "numlevel", 1) + root_level - 1; + } else { + max_level = 63; + } mesh_data.SetMeshPointer(this); + + CheckMeshValidity(); } void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, @@ -231,88 +240,12 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, int mesh_test) : Mesh(pin, app_in, packages, private_t()) { std::stringstream msg; - RegionSize block_size; // mesh test if (mesh_test > 0) Globals::nranks = mesh_test; - // check number of OpenMP threads for mesh - if (num_mesh_threads_ < 1) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "Number of OpenMP threads must be >= 1, but num_threads=" << num_mesh_threads_ - << std::endl; - PARTHENON_FAIL(msg); - } - - for (auto &[dir, label] : std::vector>{ - {X1DIR, "1"}, {X2DIR, "2"}, {X3DIR, "3"}}) { - // check number of grid cells in root level of mesh from input file. - if (mesh_size.nx(dir) < 1) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "In mesh block in input file nx" + label + " must be >= 1, but nx" + label + - "=" - << mesh_size.nx(dir) << std::endl; - PARTHENON_FAIL(msg); - } - // check physical size of mesh (root level) from input file. - if (mesh_size.xmax(dir) <= mesh_size.xmin(dir)) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "Input x" + label + "max must be larger than x" + label + "min: x" + label + - "min=" - << mesh_size.xmin(dir) << " x" + label + "max=" << mesh_size.xmax(dir) - << std::endl; - PARTHENON_FAIL(msg); - } - } - - if (mesh_size.nx(X2DIR) == 1 && mesh_size.nx(X3DIR) > 1) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "In mesh block in input file: nx2=1, nx3=" << mesh_size.nx(X3DIR) - << ", 2D problems in x1-x3 plane not supported" << std::endl; - PARTHENON_FAIL(msg); - } - - - - // check the consistency of the periodic boundaries - if (((mesh_bcs[BoundaryFace::inner_x1] == BoundaryFlag::periodic && - mesh_bcs[BoundaryFace::outer_x1] != BoundaryFlag::periodic) || - (mesh_bcs[BoundaryFace::inner_x1] != BoundaryFlag::periodic && - mesh_bcs[BoundaryFace::outer_x1] == BoundaryFlag::periodic)) || - (mesh_size.nx(X2DIR) > 1 && - ((mesh_bcs[BoundaryFace::inner_x2] == BoundaryFlag::periodic && - mesh_bcs[BoundaryFace::outer_x2] != BoundaryFlag::periodic) || - (mesh_bcs[BoundaryFace::inner_x2] != BoundaryFlag::periodic && - mesh_bcs[BoundaryFace::outer_x2] == BoundaryFlag::periodic))) || - (mesh_size.nx(X3DIR) > 1 && - ((mesh_bcs[BoundaryFace::inner_x3] == BoundaryFlag::periodic && - mesh_bcs[BoundaryFace::outer_x3] != BoundaryFlag::periodic) || - (mesh_bcs[BoundaryFace::inner_x3] != BoundaryFlag::periodic && - mesh_bcs[BoundaryFace::outer_x3] == BoundaryFlag::periodic)))) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "When periodic boundaries are in use, both sides must be periodic." - << std::endl; - PARTHENON_FAIL(msg); - } - EnrollBndryFncts_(app_in); - - // check consistency of the block and mesh - if (mesh_size.nx(X1DIR) % base_block_size.nx(X1DIR) != 0 || - mesh_size.nx(X2DIR) % base_block_size.nx(X2DIR) != 0 || - mesh_size.nx(X3DIR) % base_block_size.nx(X3DIR) != 0) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "the Mesh must be evenly divisible by the MeshBlock" << std::endl; - PARTHENON_FAIL(msg); - } - if (base_block_size.nx(X1DIR) < 4 || (base_block_size.nx(X2DIR) < 4 && (ndim >= 2)) || - (base_block_size.nx(X3DIR) < 4 && (ndim >= 3))) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "block_size must be larger than or equal to 4 cells." << std::endl; - PARTHENON_FAIL(msg); - } - // initialize user-enrollable functions default_pack_size_ = pin->GetOrAddInteger("parthenon/mesh", "pack_size", -1); @@ -339,14 +272,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, InitUserMeshData(this, pin); if (multilevel) { - if (base_block_size.nx(X1DIR) % 2 == 1 || (base_block_size.nx(X2DIR) % 2 == 1 && (ndim >= 2)) || - (base_block_size.nx(X3DIR) % 2 == 1 && (ndim >= 3))) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "The size of MeshBlock must be divisible by 2 in order to use SMR or AMR." - << std::endl; - PARTHENON_FAIL(msg); - } - InputBlock *pib = pin->pfirst_block; while (pib != nullptr) { if (pib->block_name.compare(0, 27, "parthenon/static_refinement") == 0) { @@ -451,15 +376,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, // mesh test if (mesh_test > 0) Globals::nranks = mesh_test; - // check the number of OpenMP threads for mesh - if (num_mesh_threads_ < 1) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "Number of OpenMP threads must be >= 1, but num_threads=" << num_mesh_threads_ - << std::endl; - PARTHENON_FAIL(msg); - } - - forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); root_level = forest.root_level; @@ -496,21 +412,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, // Load balancing flag and parameters RegisterLoadBalancing_(pin); - // SMR / AMR - if (adaptive) { - // read from file or from input? input for now. - // max_level = rr.GetAttr("Info", "MaxLevel"); - max_level = pin->GetOrAddInteger("parthenon/mesh", "numlevel", 1) + root_level - 1; - if (max_level > 63) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "The number of the refinement level must be smaller than " - << 63 - root_level + 1 << "." << std::endl; - PARTHENON_FAIL(msg); - } - } else { - max_level = 63; - } - InitUserMeshData(this, pin); // Populate logical locations @@ -1126,7 +1027,95 @@ void Mesh::SetupMPIComms() { } void Mesh::CheckMeshValidity() const { + std::stringstream msg; + // check number of OpenMP threads for mesh + if (num_mesh_threads_ < 1) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "Number of OpenMP threads must be >= 1, but num_threads=" << num_mesh_threads_ + << std::endl; + PARTHENON_FAIL(msg); + } + for (auto &[dir, label] : std::vector>{ + {X1DIR, "1"}, {X2DIR, "2"}, {X3DIR, "3"}}) { + // check number of grid cells in root level of mesh from input file. + if (mesh_size.nx(dir) < 1) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "In mesh block in input file nx" + label + " must be >= 1, but nx" + label + + "=" + << mesh_size.nx(dir) << std::endl; + PARTHENON_FAIL(msg); + } + // check physical size of mesh (root level) from input file. + if (mesh_size.xmax(dir) <= mesh_size.xmin(dir)) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "Input x" + label + "max must be larger than x" + label + "min: x" + label + + "min=" + << mesh_size.xmin(dir) << " x" + label + "max=" << mesh_size.xmax(dir) + << std::endl; + PARTHENON_FAIL(msg); + } + } + + if (mesh_size.nx(X2DIR) == 1 && mesh_size.nx(X3DIR) > 1) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "In mesh block in input file: nx2=1, nx3=" << mesh_size.nx(X3DIR) + << ", 2D problems in x1-x3 plane not supported" << std::endl; + PARTHENON_FAIL(msg); + } + + // check the consistency of the periodic boundaries + if (((mesh_bcs[BoundaryFace::inner_x1] == BoundaryFlag::periodic && + mesh_bcs[BoundaryFace::outer_x1] != BoundaryFlag::periodic) || + (mesh_bcs[BoundaryFace::inner_x1] != BoundaryFlag::periodic && + mesh_bcs[BoundaryFace::outer_x1] == BoundaryFlag::periodic)) || + (mesh_size.nx(X2DIR) > 1 && + ((mesh_bcs[BoundaryFace::inner_x2] == BoundaryFlag::periodic && + mesh_bcs[BoundaryFace::outer_x2] != BoundaryFlag::periodic) || + (mesh_bcs[BoundaryFace::inner_x2] != BoundaryFlag::periodic && + mesh_bcs[BoundaryFace::outer_x2] == BoundaryFlag::periodic))) || + (mesh_size.nx(X3DIR) > 1 && + ((mesh_bcs[BoundaryFace::inner_x3] == BoundaryFlag::periodic && + mesh_bcs[BoundaryFace::outer_x3] != BoundaryFlag::periodic) || + (mesh_bcs[BoundaryFace::inner_x3] != BoundaryFlag::periodic && + mesh_bcs[BoundaryFace::outer_x3] == BoundaryFlag::periodic)))) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "When periodic boundaries are in use, both sides must be periodic." + << std::endl; + PARTHENON_FAIL(msg); + } + + // check consistency of the block and mesh + if (mesh_size.nx(X1DIR) % base_block_size.nx(X1DIR) != 0 || + mesh_size.nx(X2DIR) % base_block_size.nx(X2DIR) != 0 || + mesh_size.nx(X3DIR) % base_block_size.nx(X3DIR) != 0) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "the Mesh must be evenly divisible by the MeshBlock" << std::endl; + PARTHENON_FAIL(msg); + } + if (base_block_size.nx(X1DIR) < 4 || (base_block_size.nx(X2DIR) < 4 && (ndim >= 2)) || + (base_block_size.nx(X3DIR) < 4 && (ndim >= 3))) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "block_size must be larger than or equal to 4 cells." << std::endl; + PARTHENON_FAIL(msg); + } + + if (max_level > 63) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "The number of the refinement level must be smaller than " + << 63 - root_level + 1 << "." << std::endl; + PARTHENON_FAIL(msg); + } + + if (multilevel) { + if (base_block_size.nx(X1DIR) % 2 == 1 || (base_block_size.nx(X2DIR) % 2 == 1 && (ndim >= 2)) || + (base_block_size.nx(X3DIR) % 2 == 1 && (ndim >= 3))) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "The size of MeshBlock must be divisible by 2 in order to use SMR or AMR." + << std::endl; + PARTHENON_FAIL(msg); + } + } } } // namespace parthenon From a1d92ab87e4a7066150e4f4c9a865a5a58f870ee Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 15 Apr 2024 19:32:54 -0600 Subject: [PATCH 08/22] working --- src/mesh/mesh.cpp | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 2e674c4ac98a..54427a5cc756 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -153,12 +153,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, nrbx[dir - 1] = mesh_size.nx(dir) / base_block_size.nx(dir); } - // SMR / AMR: - if (adaptive) { - max_level = pin->GetOrAddInteger("parthenon/mesh", "numlevel", 1) + root_level - 1; - } else { - max_level = 63; - } mesh_data.SetMeshPointer(this); @@ -412,6 +406,19 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, // Load balancing flag and parameters RegisterLoadBalancing_(pin); + // SMR / AMR: + if (adaptive) { + max_level = pin->GetOrAddInteger("parthenon/mesh", "numlevel", 1) + root_level - 1; + if (max_level > 63) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "The number of the refinement level must be smaller than " + << 63 - root_level + 1 << "." << std::endl; + PARTHENON_FAIL(msg); + } + } else { + max_level = 63; + } + InitUserMeshData(this, pin); // Populate logical locations From eeee5d58b9acc172e5938107331620bbf214b53d Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 15 Apr 2024 19:54:07 -0600 Subject: [PATCH 09/22] working, but gross --- src/mesh/mesh.cpp | 59 ++++++++++++----------------------------------- 1 file changed, 15 insertions(+), 44 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 54427a5cc756..e7b59b6f7632 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -155,6 +155,21 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, mesh_data.SetMeshPointer(this); + + // Load balancing flag and parameters + EnrollBndryFncts_(app_in); + RegisterLoadBalancing_(pin); + + forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); + root_level = forest.root_level; + // SMR / AMR: + if (adaptive) { + max_level = pin->GetOrAddInteger("parthenon/mesh", "numlevel", 1) + root_level - 1; + } else { + max_level = 63; + } + + InitUserMeshData(this, pin); CheckMeshValidity(); } @@ -238,33 +253,10 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // mesh test if (mesh_test > 0) Globals::nranks = mesh_test; - EnrollBndryFncts_(app_in); // initialize user-enrollable functions default_pack_size_ = pin->GetOrAddInteger("parthenon/mesh", "pack_size", -1); - forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); - root_level = forest.root_level; - current_level = root_level; - - // Load balancing flag and parameters - RegisterLoadBalancing_(pin); - - // SMR / AMR: - if (adaptive) { - max_level = pin->GetOrAddInteger("parthenon/mesh", "numlevel", 1) + root_level - 1; - if (max_level > 63) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "The number of the refinement level must be smaller than " - << 63 - root_level + 1 << "." << std::endl; - PARTHENON_FAIL(msg); - } - } else { - max_level = 63; - } - - InitUserMeshData(this, pin); - if (multilevel) { InputBlock *pib = pin->pfirst_block; while (pib != nullptr) { @@ -370,8 +362,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, // mesh test if (mesh_test > 0) Globals::nranks = mesh_test; - forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); - root_level = forest.root_level; // read the restart file // the file is already open and the pointer is set to after @@ -382,8 +372,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, nbdel = mesh_info.nbdel; nbtotal = mesh_info.nbtotal; - EnrollBndryFncts_(app_in); - const auto grid_dim = mesh_info.grid_dim; PARTHENON_REQUIRE(mesh_size.xmin(X1DIR) == grid_dim[0], "Mesh size shouldn't change on restart."); PARTHENON_REQUIRE(mesh_size.xmax(X1DIR) == grid_dim[1], "Mesh size shouldn't change on restart."); @@ -403,23 +391,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, 2 * mesh_info.n_ghost, "Block size not consistent on restart."); } - // Load balancing flag and parameters - RegisterLoadBalancing_(pin); - - // SMR / AMR: - if (adaptive) { - max_level = pin->GetOrAddInteger("parthenon/mesh", "numlevel", 1) + root_level - 1; - if (max_level > 63) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "The number of the refinement level must be smaller than " - << 63 - root_level + 1 << "." << std::endl; - PARTHENON_FAIL(msg); - } - } else { - max_level = 63; - } - - InitUserMeshData(this, pin); // Populate logical locations loclist = std::vector(nbtotal); From 459f2e0478c7643b67e05f9e7268fb971f4068d0 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 15 Apr 2024 20:11:30 -0600 Subject: [PATCH 10/22] working... --- src/mesh/forest/forest.cpp | 4 + src/mesh/mesh.cpp | 186 +++++++++++++++++++------------------ src/mesh/mesh.hpp | 1 + 3 files changed, 100 insertions(+), 91 deletions(-) diff --git a/src/mesh/forest/forest.cpp b/src/mesh/forest/forest.cpp index e321304e073c..a548b731e55c 100644 --- a/src/mesh/forest/forest.cpp +++ b/src/mesh/forest/forest.cpp @@ -83,6 +83,10 @@ Forest Forest::HyperRectangular(RegionSize mesh_size, RegionSize block_size, max_common_power2_divisor = std::min(max_common_power2_divisor, MaximumPowerOf2Divisor(nblock[dir - 1])); } + + // Every root block is a tree if the next line is set to one + max_common_power2_divisor = 1; + int max_ntree = 0; for (auto dir : {X1DIR, X2DIR, X3DIR}) { if (mesh_size.symmetry(dir)) { diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index e7b59b6f7632..047b2876c02d 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -160,8 +160,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, EnrollBndryFncts_(app_in); RegisterLoadBalancing_(pin); - forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); - root_level = forest.root_level; + + root_level = 0; // SMR / AMR: if (adaptive) { max_level = pin->GetOrAddInteger("parthenon/mesh", "numlevel", 1) + root_level - 1; @@ -257,97 +257,10 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // initialize user-enrollable functions default_pack_size_ = pin->GetOrAddInteger("parthenon/mesh", "pack_size", -1); - if (multilevel) { - InputBlock *pib = pin->pfirst_block; - while (pib != nullptr) { - if (pib->block_name.compare(0, 27, "parthenon/static_refinement") == 0) { - RegionSize ref_size; - ref_size.xmin(X1DIR) = pin->GetReal(pib->block_name, "x1min"); - ref_size.xmax(X1DIR) = pin->GetReal(pib->block_name, "x1max"); - if (ndim >= 2) { - ref_size.xmin(X2DIR) = pin->GetReal(pib->block_name, "x2min"); - ref_size.xmax(X2DIR) = pin->GetReal(pib->block_name, "x2max"); - } else { - ref_size.xmin(X2DIR) = mesh_size.xmin(X2DIR); - ref_size.xmax(X2DIR) = mesh_size.xmax(X2DIR); - } - if (ndim == 3) { - ref_size.xmin(X3DIR) = pin->GetReal(pib->block_name, "x3min"); - ref_size.xmax(X3DIR) = pin->GetReal(pib->block_name, "x3max"); - } else { - ref_size.xmin(X3DIR) = mesh_size.xmin(X3DIR); - ref_size.xmax(X3DIR) = mesh_size.xmax(X3DIR); - } - int ref_lev = pin->GetInteger(pib->block_name, "level"); - int lrlev = ref_lev + root_level; - // range check - if (ref_lev < 1) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "Refinement level must be larger than 0 (root level = 0)" << std::endl; - PARTHENON_FAIL(msg); - } - if (lrlev > max_level) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "Refinement level exceeds the maximum level (specify " - << "'maxlevel' parameter in input block if adaptive)." - << std::endl; - - PARTHENON_FAIL(msg); - } - if (ref_size.xmin(X1DIR) > ref_size.xmax(X1DIR) || - ref_size.xmin(X2DIR) > ref_size.xmax(X2DIR) || - ref_size.xmin(X3DIR) > ref_size.xmax(X3DIR)) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "Invalid refinement region is specified." << std::endl; - PARTHENON_FAIL(msg); - } - if (ref_size.xmin(X1DIR) < mesh_size.xmin(X1DIR) || - ref_size.xmax(X1DIR) > mesh_size.xmax(X1DIR) || - ref_size.xmin(X2DIR) < mesh_size.xmin(X2DIR) || - ref_size.xmax(X2DIR) > mesh_size.xmax(X2DIR) || - ref_size.xmin(X3DIR) < mesh_size.xmin(X3DIR) || - ref_size.xmax(X3DIR) > mesh_size.xmax(X3DIR)) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "Refinement region must be smaller than the whole mesh." << std::endl; - PARTHENON_FAIL(msg); - } - std::int64_t l_region_min[3]{0, 0, 0}; - std::int64_t l_region_max[3]{1, 1, 1}; - for (auto dir : {X1DIR, X2DIR, X3DIR}) { - if (!mesh_size.symmetry(dir)) { - l_region_min[dir - 1] = - GetLLFromMeshCoordinate(dir, lrlev, ref_size.xmin(dir)); - l_region_max[dir - 1] = - GetLLFromMeshCoordinate(dir, lrlev, ref_size.xmax(dir)); - l_region_min[dir - 1] = - std::max(l_region_min[dir - 1], static_cast(0)); - l_region_max[dir - 1] = - std::min(l_region_max[dir - 1], - static_cast(nrbx[dir - 1] * (1LL << ref_lev) - 1)); - auto current_loc = - LogicalLocation(lrlev, l_region_max[0], l_region_max[1], l_region_max[2]); - // Remove last block if it just it's boundary overlaps with the region - if (GetMeshCoordinate(dir, BlockLocation::Left, current_loc) == - ref_size.xmax(dir)) - l_region_max[dir - 1]--; - if (l_region_min[dir - 1] % 2 == 1) l_region_min[dir - 1]--; - if (l_region_max[dir - 1] % 2 == 0) l_region_max[dir - 1]++; - } - } - for (std::int64_t k = l_region_min[2]; k < l_region_max[2]; k += 2) { - for (std::int64_t j = l_region_min[1]; j < l_region_max[1]; j += 2) { - for (std::int64_t i = l_region_min[0]; i < l_region_max[0]; i += 2) { - LogicalLocation nloc(lrlev, i, j, k); - forest.AddMeshBlock(forest.GetForestLocationFromLegacyTreeLocation(nloc)); - } - } - } - } - pib = pib->pnext; - } - } + if (multilevel) DoStaticRefinement(pin); // initial mesh hierarchy construction is completed here + forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); loclist = forest.GetMeshBlockListAndResolveGids(); BuildBlockList(pin, app_in, packages, mesh_test); } @@ -407,6 +320,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, } // rebuild the Block Tree + forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); for (int i = 0; i < nbtotal; i++) forest.AddMeshBlock(forest.GetForestLocationFromLegacyTreeLocation(loclist[i]), false); @@ -1096,4 +1010,94 @@ void Mesh::CheckMeshValidity() const { } } +void Mesh::DoStaticRefinement(ParameterInput *pin) { + std::stringstream msg; + InputBlock *pib = pin->pfirst_block; + while (pib != nullptr) { + if (pib->block_name.compare(0, 27, "parthenon/static_refinement") == 0) { + RegionSize ref_size; + ref_size.xmin(X1DIR) = pin->GetReal(pib->block_name, "x1min"); + ref_size.xmax(X1DIR) = pin->GetReal(pib->block_name, "x1max"); + if (ndim >= 2) { + ref_size.xmin(X2DIR) = pin->GetReal(pib->block_name, "x2min"); + ref_size.xmax(X2DIR) = pin->GetReal(pib->block_name, "x2max"); + } else { + ref_size.xmin(X2DIR) = mesh_size.xmin(X2DIR); + ref_size.xmax(X2DIR) = mesh_size.xmax(X2DIR); + } + if (ndim == 3) { + ref_size.xmin(X3DIR) = pin->GetReal(pib->block_name, "x3min"); + ref_size.xmax(X3DIR) = pin->GetReal(pib->block_name, "x3max"); + } else { + ref_size.xmin(X3DIR) = mesh_size.xmin(X3DIR); + ref_size.xmax(X3DIR) = mesh_size.xmax(X3DIR); + } + int ref_lev = pin->GetInteger(pib->block_name, "level"); + int lrlev = ref_lev + root_level; + // range check + if (ref_lev < 1) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "Refinement level must be larger than 0 (root level = 0)" << std::endl; + PARTHENON_FAIL(msg); + } + if (lrlev > max_level) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "Refinement level exceeds the maximum level (specify " + << "'maxlevel' parameter in input block if adaptive)." + << std::endl; + + PARTHENON_FAIL(msg); + } + if (ref_size.xmin(X1DIR) > ref_size.xmax(X1DIR) || + ref_size.xmin(X2DIR) > ref_size.xmax(X2DIR) || + ref_size.xmin(X3DIR) > ref_size.xmax(X3DIR)) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "Invalid refinement region is specified." << std::endl; + PARTHENON_FAIL(msg); + } + if (ref_size.xmin(X1DIR) < mesh_size.xmin(X1DIR) || + ref_size.xmax(X1DIR) > mesh_size.xmax(X1DIR) || + ref_size.xmin(X2DIR) < mesh_size.xmin(X2DIR) || + ref_size.xmax(X2DIR) > mesh_size.xmax(X2DIR) || + ref_size.xmin(X3DIR) < mesh_size.xmin(X3DIR) || + ref_size.xmax(X3DIR) > mesh_size.xmax(X3DIR)) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "Refinement region must be smaller than the whole mesh." << std::endl; + PARTHENON_FAIL(msg); + } + std::int64_t l_region_min[3]{0, 0, 0}; + std::int64_t l_region_max[3]{1, 1, 1}; + for (auto dir : {X1DIR, X2DIR, X3DIR}) { + if (!mesh_size.symmetry(dir)) { + l_region_min[dir - 1] = + GetLLFromMeshCoordinate(dir, lrlev, ref_size.xmin(dir)); + l_region_max[dir - 1] = + GetLLFromMeshCoordinate(dir, lrlev, ref_size.xmax(dir)); + l_region_min[dir - 1] = + std::max(l_region_min[dir - 1], static_cast(0)); + l_region_max[dir - 1] = + std::min(l_region_max[dir - 1], + static_cast(nrbx[dir - 1] * (1LL << ref_lev) - 1)); + auto current_loc = + LogicalLocation(lrlev, l_region_max[0], l_region_max[1], l_region_max[2]); + // Remove last block if it just it's boundary overlaps with the region + if (GetMeshCoordinate(dir, BlockLocation::Left, current_loc) == + ref_size.xmax(dir)) + l_region_max[dir - 1]--; + if (l_region_min[dir - 1] % 2 == 1) l_region_min[dir - 1]--; + if (l_region_max[dir - 1] % 2 == 0) l_region_max[dir - 1]++; + } + } + for (std::int64_t k = l_region_min[2]; k < l_region_max[2]; k += 2) { + for (std::int64_t j = l_region_min[1]; j < l_region_max[1]; j += 2) { + for (std::int64_t i = l_region_min[0]; i < l_region_max[0]; i += 2) { + LogicalLocation nloc(lrlev, i, j, k); + forest.AddMeshBlock(forest.GetForestLocationFromLegacyTreeLocation(nloc)); + } + } + } + } + pib = pib->pnext; + } +} } // namespace parthenon diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index deaf90956416..ca903ee0c6a8 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -304,6 +304,7 @@ class Mesh { void CheckMeshValidity() const; void BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, int mesh_test); + void DoStaticRefinement(ParameterInput *pin); void CalculateLoadBalance(std::vector const &costlist, std::vector &ranklist, std::vector &nslist, std::vector &nblist); From cd046ac0d0e79c7ab3a6d9824c347bff2d1273d9 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 15 Apr 2024 20:12:20 -0600 Subject: [PATCH 11/22] format and lint --- src/mesh/forest/forest.cpp | 2 +- src/mesh/mesh.cpp | 81 +++++++++++++++++++------------------- src/mesh/mesh.hpp | 7 ++-- 3 files changed, 45 insertions(+), 45 deletions(-) diff --git a/src/mesh/forest/forest.cpp b/src/mesh/forest/forest.cpp index a548b731e55c..632c09db842e 100644 --- a/src/mesh/forest/forest.cpp +++ b/src/mesh/forest/forest.cpp @@ -85,7 +85,7 @@ Forest Forest::HyperRectangular(RegionSize mesh_size, RegionSize block_size, } // Every root block is a tree if the next line is set to one - max_common_power2_divisor = 1; + max_common_power2_divisor = 1; int max_ntree = 0; for (auto dir : {X1DIR, X2DIR, X3DIR}) { diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 047b2876c02d..827a38f17c06 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -101,16 +101,10 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, num_mesh_threads_(pin->GetOrAddInteger("parthenon/mesh", "num_threads", 1)), use_uniform_meshgen_fn_{true, true, true, true}, lb_flag_(true), lb_automatic_(), lb_manual_(), MeshBndryFnctn{nullptr, nullptr, nullptr, nullptr, nullptr, nullptr}, - nslist(Globals::nranks), - nblist(Globals::nranks), - nref(Globals::nranks), - nderef(Globals::nranks), - rdisp(Globals::nranks), - ddisp(Globals::nranks), - bnref(Globals::nranks), - bnderef(Globals::nranks), - brdisp(Globals::nranks), - bddisp(Globals::nranks) { + nslist(Globals::nranks), nblist(Globals::nranks), nref(Globals::nranks), + nderef(Globals::nranks), rdisp(Globals::nranks), ddisp(Globals::nranks), + bnref(Globals::nranks), bnderef(Globals::nranks), brdisp(Globals::nranks), + bddisp(Globals::nranks) { // Allow for user overrides to default Parthenon functions if (app_in->InitUserMeshData != nullptr) { InitUserMeshData = app_in->InitUserMeshData; @@ -152,15 +146,13 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, } nrbx[dir - 1] = mesh_size.nx(dir) / base_block_size.nx(dir); } - mesh_data.SetMeshPointer(this); - + // Load balancing flag and parameters EnrollBndryFncts_(app_in); RegisterLoadBalancing_(pin); - - + root_level = 0; // SMR / AMR: if (adaptive) { @@ -174,8 +166,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, CheckMeshValidity(); } -void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, - int mesh_test) { +void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, + Packages_t &packages, int mesh_test) { std::stringstream msg; nbtotal = loclist.size(); current_level = -1; @@ -229,7 +221,7 @@ void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, Package block_list.resize(nbe - nbs + 1); for (int i = nbs; i <= nbe; i++) { RegionSize block_size; - BoundaryFlag block_bcs[6]; + BoundaryFlag block_bcs[6]; SetBlockSizeAndBoundaries(loclist[i], block_size, block_bcs); // create a block and add into the link list block_list[i - nbs] = @@ -253,7 +245,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // mesh test if (mesh_test > 0) Globals::nranks = mesh_test; - // initialize user-enrollable functions default_pack_size_ = pin->GetOrAddInteger("parthenon/mesh", "pack_size", -1); @@ -269,13 +260,12 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // Mesh constructor for restarts. Load the restart file Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, Packages_t &packages, int mesh_test) - : Mesh(pin, app_in, packages, private_t()) { + : Mesh(pin, app_in, packages, private_t()) { std::stringstream msg; // mesh test if (mesh_test > 0) Globals::nranks = mesh_test; - // read the restart file // the file is already open and the pointer is set to after @@ -286,25 +276,35 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, nbtotal = mesh_info.nbtotal; const auto grid_dim = mesh_info.grid_dim; - PARTHENON_REQUIRE(mesh_size.xmin(X1DIR) == grid_dim[0], "Mesh size shouldn't change on restart."); - PARTHENON_REQUIRE(mesh_size.xmax(X1DIR) == grid_dim[1], "Mesh size shouldn't change on restart."); - PARTHENON_REQUIRE(mesh_size.xrat(X1DIR) == grid_dim[2], "Mesh size shouldn't change on restart."); - - PARTHENON_REQUIRE(mesh_size.xmin(X2DIR) == grid_dim[3], "Mesh size shouldn't change on restart."); - PARTHENON_REQUIRE(mesh_size.xmax(X2DIR) == grid_dim[4], "Mesh size shouldn't change on restart."); - PARTHENON_REQUIRE(mesh_size.xrat(X2DIR) == grid_dim[5], "Mesh size shouldn't change on restart."); - - PARTHENON_REQUIRE(mesh_size.xmin(X3DIR) == grid_dim[6], "Mesh size shouldn't change on restart."); - PARTHENON_REQUIRE(mesh_size.xmax(X3DIR) == grid_dim[7], "Mesh size shouldn't change on restart."); - PARTHENON_REQUIRE(mesh_size.xrat(X3DIR) == grid_dim[8], "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xmin(X1DIR) == grid_dim[0], + "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xmax(X1DIR) == grid_dim[1], + "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xrat(X1DIR) == grid_dim[2], + "Mesh size shouldn't change on restart."); + + PARTHENON_REQUIRE(mesh_size.xmin(X2DIR) == grid_dim[3], + "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xmax(X2DIR) == grid_dim[4], + "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xrat(X2DIR) == grid_dim[5], + "Mesh size shouldn't change on restart."); + + PARTHENON_REQUIRE(mesh_size.xmin(X3DIR) == grid_dim[6], + "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xmax(X3DIR) == grid_dim[7], + "Mesh size shouldn't change on restart."); + PARTHENON_REQUIRE(mesh_size.xrat(X3DIR) == grid_dim[8], + "Mesh size shouldn't change on restart."); for (auto &dir : {X1DIR, X2DIR, X3DIR}) { PARTHENON_REQUIRE(base_block_size.nx(dir) == mesh_info.block_size[dir - 1] - - (mesh_info.block_size[dir - 1] > 1) * mesh_info.includes_ghost * - 2 * mesh_info.n_ghost, "Block size not consistent on restart."); + (mesh_info.block_size[dir - 1] > 1) * + mesh_info.includes_ghost * 2 * + mesh_info.n_ghost, + "Block size not consistent on restart."); } - // Populate logical locations loclist = std::vector(nbtotal); auto lx123 = mesh_info.lx123; @@ -918,7 +918,7 @@ void Mesh::SetupMPIComms() { #endif } -void Mesh::CheckMeshValidity() const { +void Mesh::CheckMeshValidity() const { std::stringstream msg; // check number of OpenMP threads for mesh if (num_mesh_threads_ < 1) { @@ -955,7 +955,7 @@ void Mesh::CheckMeshValidity() const { << ", 2D problems in x1-x3 plane not supported" << std::endl; PARTHENON_FAIL(msg); } - + // check the consistency of the periodic boundaries if (((mesh_bcs[BoundaryFace::inner_x1] == BoundaryFlag::periodic && mesh_bcs[BoundaryFace::outer_x1] != BoundaryFlag::periodic) || @@ -1000,7 +1000,8 @@ void Mesh::CheckMeshValidity() const { } if (multilevel) { - if (base_block_size.nx(X1DIR) % 2 == 1 || (base_block_size.nx(X2DIR) % 2 == 1 && (ndim >= 2)) || + if (base_block_size.nx(X1DIR) % 2 == 1 || + (base_block_size.nx(X2DIR) % 2 == 1 && (ndim >= 2)) || (base_block_size.nx(X3DIR) % 2 == 1 && (ndim >= 3))) { msg << "### FATAL ERROR in Mesh constructor" << std::endl << "The size of MeshBlock must be divisible by 2 in order to use SMR or AMR." @@ -1069,10 +1070,8 @@ void Mesh::DoStaticRefinement(ParameterInput *pin) { std::int64_t l_region_max[3]{1, 1, 1}; for (auto dir : {X1DIR, X2DIR, X3DIR}) { if (!mesh_size.symmetry(dir)) { - l_region_min[dir - 1] = - GetLLFromMeshCoordinate(dir, lrlev, ref_size.xmin(dir)); - l_region_max[dir - 1] = - GetLLFromMeshCoordinate(dir, lrlev, ref_size.xmax(dir)); + l_region_min[dir - 1] = GetLLFromMeshCoordinate(dir, lrlev, ref_size.xmin(dir)); + l_region_max[dir - 1] = GetLLFromMeshCoordinate(dir, lrlev, ref_size.xmax(dir)); l_region_min[dir - 1] = std::max(l_region_min[dir - 1], static_cast(0)); l_region_max[dir - 1] = diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index ca903ee0c6a8..d8da4a2f38ba 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -76,8 +76,9 @@ class Mesh { friend class MeshBlock; friend class MeshBlockTree; friend class MeshRefinement; - - struct private_t{}; + + struct private_t {}; + public: // 2x function overloads of ctor: normal and restarted simulation Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, private_t); @@ -303,7 +304,7 @@ class Mesh { // functions void CheckMeshValidity() const; void BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, - int mesh_test); + int mesh_test); void DoStaticRefinement(ParameterInput *pin); void CalculateLoadBalance(std::vector const &costlist, std::vector &ranklist, std::vector &nslist, From 4b7790e5c61ee4aac8a88d33b1603744553f64d0 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 17 Apr 2024 13:11:58 -0600 Subject: [PATCH 12/22] split up constructors more rationally --- src/mesh/mesh.cpp | 225 ++++++++++++++++++++++++---------------------- src/mesh/mesh.hpp | 11 +-- 2 files changed, 122 insertions(+), 114 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 827a38f17c06..dab1ef83cce8 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -56,34 +56,10 @@ #include "utils/partition_stl_containers.hpp" namespace parthenon { - Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, - private_t) + base_constructor_selector_t) : // public members: modified(true), is_restart(false), - // aggregate initialization of RegionSize struct: - mesh_size({pin->GetReal("parthenon/mesh", "x1min"), - pin->GetReal("parthenon/mesh", "x2min"), - pin->GetReal("parthenon/mesh", "x3min")}, - {pin->GetReal("parthenon/mesh", "x1max"), - pin->GetReal("parthenon/mesh", "x2max"), - pin->GetReal("parthenon/mesh", "x3max")}, - {pin->GetOrAddReal("parthenon/mesh", "x1rat", 1.0), - pin->GetOrAddReal("parthenon/mesh", "x2rat", 1.0), - pin->GetOrAddReal("parthenon/mesh", "x3rat", 1.0)}, - {pin->GetInteger("parthenon/mesh", "nx1"), - pin->GetInteger("parthenon/mesh", "nx2"), - pin->GetInteger("parthenon/mesh", "nx3")}, - {false, pin->GetInteger("parthenon/mesh", "nx2") == 1, - pin->GetInteger("parthenon/mesh", "nx3") == 1}), - mesh_bcs{ - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix1_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox1_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix2_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox2_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix3_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox3_bc", "reflecting"))}, - ndim((mesh_size.nx(X3DIR) > 1) ? 3 : ((mesh_size.nx(X2DIR) > 1) ? 2 : 1)), adaptive(pin->GetOrAddString("parthenon/mesh", "refinement", "none") == "adaptive" ? true : false), @@ -96,7 +72,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, multigrid(pin->GetOrAddString("parthenon/mesh", "multigrid", "false") == "true" ? true : false), - nbnew(), nbdel(), step_since_lb(), gflag(), packages(packages), + nbnew(), nbdel(), step_since_lb(), gflag(), packages(packages), + resolved_packages(ResolvePackages(packages)), // private members: num_mesh_threads_(pin->GetOrAddInteger("parthenon/mesh", "num_threads", 1)), use_uniform_meshgen_fn_{true, true, true, true}, lb_flag_(true), lb_automatic_(), @@ -133,6 +110,45 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, if (app_in->UserWorkAfterLoop != nullptr) { UserWorkAfterLoop = app_in->UserWorkAfterLoop; } + + // Default root level, may be overwritten by another constructor + root_level = 0; + // SMR / AMR: + if (adaptive) { + max_level = pin->GetOrAddInteger("parthenon/mesh", "numlevel", 1) + root_level - 1; + } else { + max_level = 63; + } + + // Setup unique comms for each variable and swarm + SetupMPIComms(); +} + +Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, + hyper_rectangular_constructor_selector_t) + : Mesh(pin, app_in, packages, base_constructor_selector_t()) { + mesh_size = RegionSize({pin->GetReal("parthenon/mesh", "x1min"), + pin->GetReal("parthenon/mesh", "x2min"), + pin->GetReal("parthenon/mesh", "x3min")}, + {pin->GetReal("parthenon/mesh", "x1max"), + pin->GetReal("parthenon/mesh", "x2max"), + pin->GetReal("parthenon/mesh", "x3max")}, + {pin->GetOrAddReal("parthenon/mesh", "x1rat", 1.0), + pin->GetOrAddReal("parthenon/mesh", "x2rat", 1.0), + pin->GetOrAddReal("parthenon/mesh", "x3rat", 1.0)}, + {pin->GetInteger("parthenon/mesh", "nx1"), + pin->GetInteger("parthenon/mesh", "nx2"), + pin->GetInteger("parthenon/mesh", "nx3")}, + {false, pin->GetInteger("parthenon/mesh", "nx2") == 1, + pin->GetInteger("parthenon/mesh", "nx3") == 1}); + mesh_bcs = { + GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix1_bc", "reflecting")), + GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox1_bc", "reflecting")), + GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix2_bc", "reflecting")), + GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox2_bc", "reflecting")), + GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix3_bc", "reflecting")), + GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox3_bc", "reflecting"))}; + ndim = (mesh_size.nx(X3DIR) > 1) ? 3 : ((mesh_size.nx(X2DIR) > 1) ? 2 : 1); for (auto &[dir, label] : std::vector>{ {X1DIR, "nx1"}, {X2DIR, "nx2"}, {X3DIR, "nx3"}}) { @@ -152,94 +168,29 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // Load balancing flag and parameters EnrollBndryFncts_(app_in); RegisterLoadBalancing_(pin); - - root_level = 0; + + forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); + root_level = forest.root_level; // SMR / AMR: if (adaptive) { max_level = pin->GetOrAddInteger("parthenon/mesh", "numlevel", 1) + root_level - 1; } else { max_level = 63; } - - InitUserMeshData(this, pin); - - CheckMeshValidity(); -} - -void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, - Packages_t &packages, int mesh_test) { - std::stringstream msg; - nbtotal = loclist.size(); - current_level = -1; - for (const auto &loc : loclist) - if (loc.level() > current_level) current_level = loc.level(); - -#ifdef MPI_PARALLEL - // check if there are sufficient blocks - if (nbtotal < Globals::nranks) { - if (mesh_test == 0) { - msg << "### FATAL ERROR in Mesh constructor" << std::endl - << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks (" - << Globals::nranks << ")" << std::endl; - PARTHENON_FAIL(msg); - } else { // test - std::cout << "### Warning in Mesh constructor" << std::endl - << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks (" - << Globals::nranks << ")" << std::endl; - } - } -#endif - - ranklist = std::vector(nbtotal); - - // initialize cost array with the simplest estimate; all the blocks are equal - costlist = std::vector(nbtotal, 1.0); - - CalculateLoadBalance(costlist, ranklist, nslist, nblist); - - // Output some diagnostic information to terminal - - // Output MeshBlock list and quit (mesh test only); do not create meshes - if (mesh_test > 0) { - if (Globals::my_rank == 0) OutputMeshStructure(ndim); - return; - } - - resolved_packages = ResolvePackages(packages); - + // Register user defined boundary conditions UserBoundaryFunctions = resolved_packages->UserBoundaryFunctions; + + InitUserMeshData(this, pin); - // Setup unique comms for each variable and swarm - SetupMPIComms(); - - // create MeshBlock list for this process - int nbs = nslist[Globals::my_rank]; - int nbe = nbs + nblist[Globals::my_rank] - 1; - // create MeshBlock list for this process - block_list.clear(); - block_list.resize(nbe - nbs + 1); - for (int i = nbs; i <= nbe; i++) { - RegionSize block_size; - BoundaryFlag block_bcs[6]; - SetBlockSizeAndBoundaries(loclist[i], block_size, block_bcs); - // create a block and add into the link list - block_list[i - nbs] = - MeshBlock::Make(i, i - nbs, loclist[i], block_size, block_bcs, this, pin, app_in, - packages, resolved_packages, gflag, costlist[i]); - } - BuildGMGBlockLists(pin, app_in); - SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist); - SetGMGNeighbors(); - ResetLoadBalanceVariables(); + CheckMeshValidity(); } //---------------------------------------------------------------------------------------- // Mesh constructor, builds mesh at start of calculation using parameters in input file - Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, int mesh_test) - : Mesh(pin, app_in, packages, private_t()) { + : Mesh(pin, app_in, packages, hyper_rectangular_constructor_selector_t()) { std::stringstream msg; // mesh test @@ -251,8 +202,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, if (multilevel) DoStaticRefinement(pin); // initial mesh hierarchy construction is completed here - forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); - loclist = forest.GetMeshBlockListAndResolveGids(); BuildBlockList(pin, app_in, packages, mesh_test); } @@ -260,7 +209,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // Mesh constructor for restarts. Load the restart file Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, Packages_t &packages, int mesh_test) - : Mesh(pin, app_in, packages, private_t()) { + : Mesh(pin, app_in, packages, hyper_rectangular_constructor_selector_t()) { + is_restart = true; std::stringstream msg; // mesh test @@ -313,20 +263,14 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, for (int i = 0; i < nbtotal; i++) { loclist[i] = LogicalLocation(locLevelGidLidCnghostGflag[5 * i], lx123[3 * i], lx123[3 * i + 1], lx123[3 * i + 2]); - - if (loclist[i].level() > current_level) { - current_level = loclist[i].level(); - } } // rebuild the Block Tree - forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); for (int i = 0; i < nbtotal; i++) forest.AddMeshBlock(forest.GetForestLocationFromLegacyTreeLocation(loclist[i]), false); - loclist = forest.GetMeshBlockListAndResolveGids(); - int nnb = loclist.size(); + int nnb = forest.CountMeshBlock(); if (nnb != nbtotal) { msg << "### FATAL ERROR in Mesh constructor" << std::endl << "Tree reconstruction failed. The total numbers of the blocks do not match. (" @@ -337,6 +281,69 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, BuildBlockList(pin, app_in, packages, mesh_test); } +void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, + Packages_t &packages, int mesh_test) { + // LFR: This routine should work for general block lists + std::stringstream msg; + + loclist = forest.GetMeshBlockListAndResolveGids(); + nbtotal = loclist.size(); + current_level = -1; + for (const auto &loc : loclist) + if (loc.level() > current_level) current_level = loc.level(); + +#ifdef MPI_PARALLEL + // check if there are sufficient blocks + if (nbtotal < Globals::nranks) { + if (mesh_test == 0) { + msg << "### FATAL ERROR in Mesh constructor" << std::endl + << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks (" + << Globals::nranks << ")" << std::endl; + PARTHENON_FAIL(msg); + } else { // test + std::cout << "### Warning in Mesh constructor" << std::endl + << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks (" + << Globals::nranks << ")" << std::endl; + } + } +#endif + + ranklist = std::vector(nbtotal); + + // initialize cost array with the simplest estimate; all the blocks are equal + costlist = std::vector(nbtotal, 1.0); + + CalculateLoadBalance(costlist, ranklist, nslist, nblist); + + // Output MeshBlock list and quit (mesh test only); do not create meshes + if (mesh_test > 0) { + if (Globals::my_rank == 0) OutputMeshStructure(ndim); + return; + } + + // create MeshBlock list for this process + int nbs = nslist[Globals::my_rank]; + int nbe = nbs + nblist[Globals::my_rank] - 1; + // create MeshBlock list for this process + block_list.clear(); + block_list.resize(nbe - nbs + 1); + for (int i = nbs; i <= nbe; i++) { + RegionSize block_size; + BoundaryFlag block_bcs[6]; + SetBlockSizeAndBoundaries(loclist[i], block_size, block_bcs); + // create a block and add into the link list + block_list[i - nbs] = + MeshBlock::Make(i, i - nbs, loclist[i], block_size, block_bcs, this, pin, app_in, + packages, resolved_packages, gflag, costlist[i]); + } + BuildGMGBlockLists(pin, app_in); + SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist); + SetGMGNeighbors(); + ResetLoadBalanceVariables(); +} + + + //---------------------------------------------------------------------------------------- // destructor diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index d8da4a2f38ba..554b0d6ff9ed 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -74,14 +74,15 @@ class Mesh { friend class RestartOutput; friend class HistoryOutput; friend class MeshBlock; - friend class MeshBlockTree; friend class MeshRefinement; - struct private_t {}; + struct base_constructor_selector_t{}; + Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, base_constructor_selector_t); + struct hyper_rectangular_constructor_selector_t{}; + Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, hyper_rectangular_constructor_selector_t); public: // 2x function overloads of ctor: normal and restarted simulation - Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, private_t); Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, int test_flag = 0); Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &resfile, @@ -106,11 +107,11 @@ class Mesh { // data bool modified; - const bool is_restart; + bool is_restart; RegionSize mesh_size; RegionSize base_block_size; std::array mesh_bcs; - const int ndim; // number of dimensions + int ndim; // number of dimensions const bool adaptive, multilevel, multigrid; int nbtotal, nbnew, nbdel; std::uint64_t mbcnt; From 0c6f3d610a07ad946a014815ad27157b30245664 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 17 Apr 2024 13:25:44 -0600 Subject: [PATCH 13/22] format --- src/mesh/mesh.cpp | 40 ++++++++++++++++++---------------------- src/mesh/mesh.hpp | 10 ++++++---- 2 files changed, 24 insertions(+), 26 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 40aef6e061ff..9319a0aa1a1e 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -72,7 +72,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, multigrid(pin->GetOrAddString("parthenon/mesh", "multigrid", "false") == "true" ? true : false), - nbnew(), nbdel(), step_since_lb(), gflag(), packages(packages), + nbnew(), nbdel(), step_since_lb(), gflag(), packages(packages), resolved_packages(ResolvePackages(packages)), // private members: num_mesh_threads_(pin->GetOrAddInteger("parthenon/mesh", "num_threads", 1)), @@ -110,7 +110,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, if (app_in->UserWorkAfterLoop != nullptr) { UserWorkAfterLoop = app_in->UserWorkAfterLoop; } - + // Default root level, may be overwritten by another constructor root_level = 0; // SMR / AMR: @@ -127,20 +127,18 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, hyper_rectangular_constructor_selector_t) : Mesh(pin, app_in, packages, base_constructor_selector_t()) { - mesh_size = RegionSize({pin->GetReal("parthenon/mesh", "x1min"), - pin->GetReal("parthenon/mesh", "x2min"), - pin->GetReal("parthenon/mesh", "x3min")}, - {pin->GetReal("parthenon/mesh", "x1max"), - pin->GetReal("parthenon/mesh", "x2max"), - pin->GetReal("parthenon/mesh", "x3max")}, - {pin->GetOrAddReal("parthenon/mesh", "x1rat", 1.0), - pin->GetOrAddReal("parthenon/mesh", "x2rat", 1.0), - pin->GetOrAddReal("parthenon/mesh", "x3rat", 1.0)}, - {pin->GetInteger("parthenon/mesh", "nx1"), - pin->GetInteger("parthenon/mesh", "nx2"), - pin->GetInteger("parthenon/mesh", "nx3")}, - {false, pin->GetInteger("parthenon/mesh", "nx2") == 1, - pin->GetInteger("parthenon/mesh", "nx3") == 1}); + mesh_size = RegionSize( + {pin->GetReal("parthenon/mesh", "x1min"), pin->GetReal("parthenon/mesh", "x2min"), + pin->GetReal("parthenon/mesh", "x3min")}, + {pin->GetReal("parthenon/mesh", "x1max"), pin->GetReal("parthenon/mesh", "x2max"), + pin->GetReal("parthenon/mesh", "x3max")}, + {pin->GetOrAddReal("parthenon/mesh", "x1rat", 1.0), + pin->GetOrAddReal("parthenon/mesh", "x2rat", 1.0), + pin->GetOrAddReal("parthenon/mesh", "x3rat", 1.0)}, + {pin->GetInteger("parthenon/mesh", "nx1"), pin->GetInteger("parthenon/mesh", "nx2"), + pin->GetInteger("parthenon/mesh", "nx3")}, + {false, pin->GetInteger("parthenon/mesh", "nx2") == 1, + pin->GetInteger("parthenon/mesh", "nx3") == 1}); mesh_bcs = { GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix1_bc", "reflecting")), GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox1_bc", "reflecting")), @@ -168,7 +166,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // Load balancing flag and parameters EnrollBndryFncts_(app_in); RegisterLoadBalancing_(pin); - + forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); root_level = forest.root_level; // SMR / AMR: @@ -177,10 +175,10 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, } else { max_level = 63; } - + // Register user defined boundary conditions UserBoundaryFunctions = resolved_packages->UserBoundaryFunctions; - + InitUserMeshData(this, pin); CheckMeshValidity(); @@ -285,7 +283,7 @@ void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, int mesh_test) { // LFR: This routine should work for general block lists std::stringstream msg; - + loclist = forest.GetMeshBlockListAndResolveGids(); nbtotal = loclist.size(); current_level = -1; @@ -342,8 +340,6 @@ void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, ResetLoadBalanceVariables(); } - - //---------------------------------------------------------------------------------------- // destructor diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 7922f9a6bb0f..cc48b6b859b7 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -76,10 +76,12 @@ class Mesh { friend class MeshBlock; friend class MeshRefinement; - struct base_constructor_selector_t{}; - Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, base_constructor_selector_t); - struct hyper_rectangular_constructor_selector_t{}; - Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, hyper_rectangular_constructor_selector_t); + struct base_constructor_selector_t {}; + Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, + base_constructor_selector_t); + struct hyper_rectangular_constructor_selector_t {}; + Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, + hyper_rectangular_constructor_selector_t); public: // 2x function overloads of ctor: normal and restarted simulation From 3f5dcf4f52444a1d51462a4b378561905eca1ee6 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 17 Apr 2024 14:04:24 -0600 Subject: [PATCH 14/22] Remove nrbx and specific coordinate location stuff --- src/mesh/mesh.cpp | 27 ++++++++++++++++++++++++--- src/mesh/mesh.hpp | 21 --------------------- 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 9319a0aa1a1e..1541d1a7b5ab 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -158,7 +158,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, } else { base_block_size.nx(dir) = mesh_size.nx(dir); } - nrbx[dir - 1] = mesh_size.nx(dir) / base_block_size.nx(dir); } mesh_data.SetMeshPointer(this); @@ -364,8 +363,7 @@ void Mesh::OutputMeshStructure(const int ndim, // Write overall Mesh structure to stdout and file std::cout << std::endl; - std::cout << "Root grid = " << nrbx[0] << " x " << nrbx[1] << " x " << nrbx[2] - << " MeshBlocks" << std::endl; + std::cout << "Number of Trees = " << forest.CountTrees() << std::endl; std::cout << "Total number of MeshBlocks = " << nbtotal << std::endl; std::cout << "Number of physical refinement levels = " << (current_level - root_level) << std::endl; @@ -1016,6 +1014,29 @@ void Mesh::CheckMeshValidity() const { void Mesh::DoStaticRefinement(ParameterInput *pin) { std::stringstream msg; + + // TODO(LFR): Static refinement currently only works for hyper-rectangular meshes + std::array nrbx; + for (auto dir : {X1DIR, X2DIR, X3DIR}) + nrbx[dir - 1] = mesh_size.nx(dir) / base_block_size.nx(dir); + + auto GetMeshCoordinate = [this, nrbx](CoordinateDirection dir, BlockLocation bloc, + const LogicalLocation &loc) -> Real { + auto xll = loc.LLCoord(dir, bloc); + auto root_fac = static_cast(1 << this->root_level) / static_cast(nrbx[dir - 1]); + xll *= root_fac; + return this->mesh_size.xmin(dir) * (1.0 - xll) + this->mesh_size.xmax(dir) * xll; + }; + + auto GetLLFromMeshCoordinate = [this, nrbx](CoordinateDirection dir, int level, + Real xmesh) -> std::int64_t { + auto root_fac = static_cast(1 << this->root_level) / static_cast(nrbx[dir - 1]); + auto xLL = (xmesh - this->mesh_size.xmin(dir)) / + (this->mesh_size.xmax(dir) - this->mesh_size.xmin(dir)) / root_fac; + return static_cast((1 << std::max(level, 0)) * xLL); + }; + + InputBlock *pib = pin->pfirst_block; while (pib != nullptr) { if (pib->block_name.compare(0, 27, "parthenon/static_refinement") == 0) { diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index cc48b6b859b7..27652f02c468 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -270,9 +270,6 @@ class Mesh { std::vector loclist; forest::Forest forest; - // number of MeshBlocks in the x1, x2, x3 directions of the root grid: - // (unlike LogicalLocation.lxi, nrbxi don't grow w/ AMR # of levels, so keep 32-bit int) - std::array nrbx; // flags are false if using non-uniform or user meshgen function bool use_uniform_meshgen_fn_[4]; @@ -330,24 +327,6 @@ class Mesh { void CommunicateBoundaries(std::string md_name = "base"); void PreCommFillDerived(); void FillDerived(); - - // Transform from logical location coordinates to uniform mesh coordinates accounting - // for root grid - Real GetMeshCoordinate(CoordinateDirection dir, BlockLocation bloc, - const LogicalLocation &loc) const { - auto xll = loc.LLCoord(dir, bloc); - auto root_fac = static_cast(1 << root_level) / static_cast(nrbx[dir - 1]); - xll *= root_fac; - return mesh_size.xmin(dir) * (1.0 - xll) + mesh_size.xmax(dir) * xll; - } - - std::int64_t GetLLFromMeshCoordinate(CoordinateDirection dir, int level, - Real xmesh) const { - auto root_fac = static_cast(1 << root_level) / static_cast(nrbx[dir - 1]); - auto xLL = (xmesh - mesh_size.xmin(dir)) / - (mesh_size.xmax(dir) - mesh_size.xmin(dir)) / root_fac; - return static_cast((1 << std::max(level, 0)) * xLL); - } }; } // namespace parthenon From 9c4721778f84d1df202a393466550e573f9a75f8 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 17 Apr 2024 14:05:00 -0600 Subject: [PATCH 15/22] format --- src/mesh/mesh.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 1541d1a7b5ab..c86d931c5a42 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -1014,29 +1014,30 @@ void Mesh::CheckMeshValidity() const { void Mesh::DoStaticRefinement(ParameterInput *pin) { std::stringstream msg; - + // TODO(LFR): Static refinement currently only works for hyper-rectangular meshes - std::array nrbx; + std::array nrbx; for (auto dir : {X1DIR, X2DIR, X3DIR}) nrbx[dir - 1] = mesh_size.nx(dir) / base_block_size.nx(dir); auto GetMeshCoordinate = [this, nrbx](CoordinateDirection dir, BlockLocation bloc, - const LogicalLocation &loc) -> Real { + const LogicalLocation &loc) -> Real { auto xll = loc.LLCoord(dir, bloc); - auto root_fac = static_cast(1 << this->root_level) / static_cast(nrbx[dir - 1]); + auto root_fac = + static_cast(1 << this->root_level) / static_cast(nrbx[dir - 1]); xll *= root_fac; return this->mesh_size.xmin(dir) * (1.0 - xll) + this->mesh_size.xmax(dir) * xll; }; - + auto GetLLFromMeshCoordinate = [this, nrbx](CoordinateDirection dir, int level, - Real xmesh) -> std::int64_t { - auto root_fac = static_cast(1 << this->root_level) / static_cast(nrbx[dir - 1]); + Real xmesh) -> std::int64_t { + auto root_fac = + static_cast(1 << this->root_level) / static_cast(nrbx[dir - 1]); auto xLL = (xmesh - this->mesh_size.xmin(dir)) / (this->mesh_size.xmax(dir) - this->mesh_size.xmin(dir)) / root_fac; return static_cast((1 << std::max(level, 0)) * xLL); }; - InputBlock *pib = pin->pfirst_block; while (pib != nullptr) { if (pib->block_name.compare(0, 27, "parthenon/static_refinement") == 0) { From 3d958da0ca78c8e1fd01ffce427f2252b68c163a Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 17 Apr 2024 14:14:25 -0600 Subject: [PATCH 16/22] changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6bcd1a67830a..68501b9cbeef 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,7 @@ - [[PR 1031]](https://github.com/parthenon-hpc-lab/parthenon/pull/1031) Fix bug in non-cell centered AMR ### Infrastructure (changes irrelevant to downstream codes) +- [[PR 1055]](https://github.com/parthenon-hpc-lab/parthenon/pull/1055) Refactor mesh constructors - [[PR 1035]](https://github.com/parthenon-hpc-lab/parthenon/pull/1035) Fix multigrid infrastructure to work with forest - [[PR 1048]](https://github.com/parthenon-hpc-lab/parthenon/pull/1048) Tiny fixes to custom coords logic - [[PR 1028]](https://github.com/parthenon-hpc-lab/parthenon/pull/1028) Internal reorganization of LogicalLocation files From e3dc06548f09e4b84af641624c3cb41aae2f0a6e Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 17 Apr 2024 14:15:09 -0600 Subject: [PATCH 17/22] use old root grid --- src/mesh/forest/forest.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/mesh/forest/forest.cpp b/src/mesh/forest/forest.cpp index 632c09db842e..e321304e073c 100644 --- a/src/mesh/forest/forest.cpp +++ b/src/mesh/forest/forest.cpp @@ -83,10 +83,6 @@ Forest Forest::HyperRectangular(RegionSize mesh_size, RegionSize block_size, max_common_power2_divisor = std::min(max_common_power2_divisor, MaximumPowerOf2Divisor(nblock[dir - 1])); } - - // Every root block is a tree if the next line is set to one - max_common_power2_divisor = 1; - int max_ntree = 0; for (auto dir : {X1DIR, X2DIR, X3DIR}) { if (mesh_size.symmetry(dir)) { From 9e2330ed315f0f7d40a8a12b9babba16275e3165 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 17 Apr 2024 18:54:29 -0600 Subject: [PATCH 18/22] move stuff to where it belongs --- src/mesh/mesh.cpp | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index c86d931c5a42..e99e97cfb8b9 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -74,6 +74,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, : false), nbnew(), nbdel(), step_since_lb(), gflag(), packages(packages), resolved_packages(ResolvePackages(packages)), + default_pack_size_(pin->GetOrAddInteger("parthenon/mesh", "pack_size", -1)), // private members: num_mesh_threads_(pin->GetOrAddInteger("parthenon/mesh", "num_threads", 1)), use_uniform_meshgen_fn_{true, true, true, true}, lb_flag_(true), lb_automatic_(), @@ -120,8 +121,13 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, max_level = 63; } - // Setup unique comms for each variable and swarm SetupMPIComms(); + + RegisterLoadBalancing_(pin); + + mesh_data.SetMeshPointer(this); + + InitUserMeshData(this, pin); } Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, @@ -160,11 +166,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, } } - mesh_data.SetMeshPointer(this); - // Load balancing flag and parameters EnrollBndryFncts_(app_in); - RegisterLoadBalancing_(pin); forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); root_level = forest.root_level; @@ -178,8 +181,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // Register user defined boundary conditions UserBoundaryFunctions = resolved_packages->UserBoundaryFunctions; - InitUserMeshData(this, pin); - CheckMeshValidity(); } @@ -188,14 +189,9 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, int mesh_test) : Mesh(pin, app_in, packages, hyper_rectangular_constructor_selector_t()) { - std::stringstream msg; - // mesh test if (mesh_test > 0) Globals::nranks = mesh_test; - // initialize user-enrollable functions - default_pack_size_ = pin->GetOrAddInteger("parthenon/mesh", "pack_size", -1); - if (multilevel) DoStaticRefinement(pin); // initial mesh hierarchy construction is completed here From 284b54e0ca11b55e75bd2915bebbbc2bca54441d Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 13 May 2024 21:51:16 -0600 Subject: [PATCH 19/22] fix issue --- src/mesh/mesh.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index ff32a45a1230..a50b38c5e330 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -130,7 +130,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, mesh_data.SetMeshPointer(this); - InitUserMeshData(this, pin); + if (InitUserMeshData) InitUserMeshData(this, pin); } Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, @@ -184,7 +184,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // Register user defined boundary conditions UserBoundaryFunctions = resolved_packages->UserBoundaryFunctions; UserSwarmBoundaryFunctions = resolved_packages->UserSwarmBoundaryFunctions; - + CheckMeshValidity(); } From d1e96ce4957c80fa54b0f3f88050c5b44edb7966 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Tue, 14 May 2024 12:45:25 -0600 Subject: [PATCH 20/22] correctly pass dealloc count --- src/mesh/mesh.cpp | 17 +++++++++++------ src/mesh/mesh.hpp | 3 ++- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index a00b95c72018..a55dcd7f6a24 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -254,12 +254,14 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, // Populate logical locations loclist = std::vector(nbtotal); + std::unordered_map dealloc_count; auto lx123 = mesh_info.lx123; auto locLevelGidLidCnghostGflag = mesh_info.level_gid_lid_cnghost_gflag; current_level = -1; for (int i = 0; i < nbtotal; i++) { - loclist[i] = LogicalLocation(locLevelGidLidCnghostGflag[5 * i], lx123[3 * i], + loclist[i] = LogicalLocation(locLevelGidLidCnghostGflag[6 * i], lx123[3 * i], lx123[3 * i + 1], lx123[3 * i + 2]); + dealloc_count[loclist[i]] = locLevelGidLidCnghostGflag[6 * i + 5]; } // rebuild the Block Tree @@ -275,11 +277,12 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, PARTHENON_FAIL(msg); } - BuildBlockList(pin, app_in, packages, mesh_test); + BuildBlockList(pin, app_in, packages, mesh_test, dealloc_count); } void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, - Packages_t &packages, int mesh_test) { + Packages_t &packages, int mesh_test, + const std::unordered_map &dealloc_count) { // LFR: This routine should work for general block lists std::stringstream msg; @@ -334,7 +337,7 @@ void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, packages, resolved_packages, gflag, costlist[i]); if (block_list[i - nbs]->pmr) block_list[i - nbs]->pmr->DerefinementCount() = - locLevelGidLidCnghostGflag[6 * i + 5]; + dealloc_count.count(loclist[i]) ? dealloc_count.at(loclist[i]) : 0; } BuildGMGBlockLists(pin, app_in); SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist); @@ -1138,8 +1141,10 @@ void Mesh::DoStaticRefinement(ParameterInput *pin) { std::int64_t l_region_max[3]{1, 1, 1}; for (auto dir : {X1DIR, X2DIR, X3DIR}) { if (!mesh_size.symmetry(dir)) { - l_region_min[dir - 1] = GetLegacyLLFromMeshCoordinate(dir, lrlev, ref_size.xmin(dir)); - l_region_max[dir - 1] = GetLegacyLLFromMeshCoordinate(dir, lrlev, ref_size.xmax(dir)); + l_region_min[dir - 1] = + GetLegacyLLFromMeshCoordinate(dir, lrlev, ref_size.xmin(dir)); + l_region_max[dir - 1] = + GetLegacyLLFromMeshCoordinate(dir, lrlev, ref_size.xmax(dir)); l_region_min[dir - 1] = std::max(l_region_min[dir - 1], static_cast(0)); l_region_max[dir - 1] = diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 0422308a1d61..0e1a25c0b76d 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -295,7 +295,8 @@ class Mesh { // functions void CheckMeshValidity() const; void BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, - int mesh_test); + int mesh_test, + const std::unordered_map &dealloc_count = {}); void DoStaticRefinement(ParameterInput *pin); void CalculateLoadBalance(std::vector const &costlist, std::vector &ranklist, std::vector &nslist, From 0b85e7faae395ebf3790d961f89646a9ce94bfcc Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 16 May 2024 15:20:36 -0600 Subject: [PATCH 21/22] format --- src/mesh/mesh.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index fb640e394fb6..b54b09a7cb2c 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -259,8 +259,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, auto locLevelGidLidCnghostGflag = mesh_info.level_gid_lid_cnghost_gflag; current_level = -1; for (int i = 0; i < nbtotal; i++) { - loclist[i] = LogicalLocation(locLevelGidLidCnghostGflag[NumIDsAndFlags * i], lx123[3 * i], - lx123[3 * i + 1], lx123[3 * i + 2]); + loclist[i] = LogicalLocation(locLevelGidLidCnghostGflag[NumIDsAndFlags * i], + lx123[3 * i], lx123[3 * i + 1], lx123[3 * i + 2]); dealloc_count[loclist[i]] = locLevelGidLidCnghostGflag[NumIDsAndFlags * i + 5]; } From d8cc662ed10f2bc7fb389ea8d2b1f5a11b5e0653 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 29 May 2024 12:08:27 -0600 Subject: [PATCH 22/22] Use forest locations for mapping dealloc_count --- src/mesh/mesh.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index cd309175f20e..913ac28a8b13 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -257,17 +257,16 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, std::unordered_map dealloc_count; auto lx123 = mesh_info.lx123; auto locLevelGidLidCnghostGflag = mesh_info.level_gid_lid_cnghost_gflag; - current_level = -1; for (int i = 0; i < nbtotal; i++) { - loclist[i] = LogicalLocation(locLevelGidLidCnghostGflag[NumIDsAndFlags * i], - lx123[3 * i], lx123[3 * i + 1], lx123[3 * i + 2]); + loclist[i] = forest.GetForestLocationFromLegacyTreeLocation( + LogicalLocation(locLevelGidLidCnghostGflag[NumIDsAndFlags * i], lx123[3 * i], + lx123[3 * i + 1], lx123[3 * i + 2])); dealloc_count[loclist[i]] = mesh_info.derefinement_count[i]; } // rebuild the Block Tree for (int i = 0; i < nbtotal; i++) - forest.AddMeshBlock(forest.GetForestLocationFromLegacyTreeLocation(loclist[i]), - false); + forest.AddMeshBlock(loclist[i], false); int nnb = forest.CountMeshBlock(); if (nnb != nbtotal) {