diff --git a/CHANGELOG.md b/CHANGELOG.md index 6a00ac911ba5..8336392611a8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[PR 1103]](https://github.com/parthenon-hpc-lab/parthenon/pull/1103) Add sparsity to vector wave equation test - [[PR 1185]](https://github.com/parthenon-hpc-lab/parthenon/pull/1185) Bugfix to particle defragmentation - [[PR 1184]](https://github.com/parthenon-hpc-lab/parthenon/pull/1184) Fix swarm block neighbor indexing in 1D, 2D - [[PR 1183]](https://github.com/parthenon-hpc-lab/parthenon/pull/1183) Fix particle leapfrog example initialization data @@ -12,6 +13,7 @@ ### Changed (changing behavior/API/variables/...) - [[PR 1191]](https://github.com/parthenon-hpc-lab/parthenon/pull/1191) Update Kokkos version to 4.4.1 +- [[PR 1209]](https://github.com/parthenon-hpc-lab/parthenon/pull/1209) Ordered history output - [[PR 1206]](https://github.com/parthenon-hpc-lab/parthenon/pull/1206) Leapfrog fix - [[PR 1203]](https://github.com/parthenon-hpc-lab/parthenon/pull/1203) Pin Ubuntu CI image - [[PR 1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag @@ -20,6 +22,7 @@ - [[PR 1172]](https://github.com/parthenon-hpc-lab/parthenon/pull/1172) Make parthenon manager robust against external MPI init and finalize calls ### Fixed (not changing behavior/API/variables/...) +- [[PR 1188]](https://github.com/parthenon-hpc-lab/parthenon/pull/1188) Fix hdf5 output issue for metadata none variables, update test. - [[PR 1170]](https://github.com/parthenon-hpc-lab/parthenon/pull/1170) Fixed incorrect initialization of array by a const not constexpr - [[PR 1189]](https://github.com/parthenon-hpc-lab/parthenon/pull/1189) Address CUDA MPI/ICP issue with Kokkos <=4.4.1 - [[PR 1178]](https://github.com/parthenon-hpc-lab/parthenon/pull/1178) Fix issue with mesh pointer when using relative residual tolerance in BiCGSTAB solver. diff --git a/CMakeLists.txt b/CMakeLists.txt index 61fd4fd3cfff..590f277ee8d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -75,9 +75,9 @@ include(cmake/Format.cmake) include(cmake/Lint.cmake) # regression test reference data -set(REGRESSION_GOLD_STANDARD_VER 24 CACHE STRING "Version of gold standard to download and use") +set(REGRESSION_GOLD_STANDARD_VER 25 CACHE STRING "Version of gold standard to download and use") set(REGRESSION_GOLD_STANDARD_HASH - "SHA512=e220df92a335131131e42ddb52dc221a6dbd6bb56361483b4af0292620eeb82ffb21ef3b95fd9a7c5cc158fb754da0bf1a1015bec98b5bbad05f4bceb1ee99bc" + "SHA512=314dc8312366d81ba33d1fde25812e9a7697b2f529de29e22662df0d458f1c4bc5b5bb4e649888170f66ffec0df1be20a9cf401944531c1c1ad835e26eaad28f" CACHE STRING "Hash of default gold standard file to download") option(REGRESSION_GOLD_STANDARD_SYNC "Automatically sync gold standard files." ON) diff --git a/doc/sphinx/src/outputs.rst b/doc/sphinx/src/outputs.rst index 7196e4777bb3..fa44e27da298 100644 --- a/doc/sphinx/src/outputs.rst +++ b/doc/sphinx/src/outputs.rst @@ -194,7 +194,9 @@ block might look like This will produce a text file (``.hst``) output file every 1 units of simulation time. The content of the file is determined by the functions -enrolled by specific packages, see :ref:`state history output`. +enrolled by specific packages, see :ref:`state history output`. Per-package history +outputs will always be in alphabetical order by package name, which may not match +the order in which packages were added to a simulation. Histograms ---------- diff --git a/example/advection/advection_package.cpp b/example/advection/advection_package.cpp index f799a2909262..64e4ea30712d 100644 --- a/example/advection/advection_package.cpp +++ b/example/advection/advection_package.cpp @@ -196,6 +196,20 @@ std::shared_ptr Initialize(ParameterInput *pin) { m = Metadata({Metadata::Cell, Metadata::OneCopy}, std::vector({1})); pkg->AddField("my_derived_var", m); + // Create a Metadata::None variable for IO testing purposes. + // Only load if test_metadata_none is specified in the Advection block + auto test_metadata_none = + pin->GetOrAddBoolean("Advection", "test_metadata_none", false); + pkg->AddParam("test_metadata_none", test_metadata_none); + if (test_metadata_none) { + const int nx1 = pin->GetOrAddInteger("parthenon/meshblock", "nx1", 1); + const int nx2 = pin->GetOrAddInteger("parthenon/meshblock", "nx2", 1); + const int nx3 = pin->GetOrAddInteger("parthenon/meshblock", "nx3", 1); + std::vector test_shape = {nx1 + 1, nx2 + 1, nx3 + 1, 3}; + m = Metadata({Metadata::OneCopy, Metadata::None}, test_shape); + pkg->AddField("metadata_none_var", m); + } + // List (vector) of HistoryOutputVar that will all be enrolled as output variables parthenon::HstVar_list hst_vars = {}; // Now we add a couple of callback functions @@ -281,6 +295,7 @@ AmrTag CheckRefinement(MeshBlockData *rc) { void PreFill(MeshBlockData *rc) { auto pmb = rc->GetBlockPointer(); auto pkg = pmb->packages.Get("advection_package"); + const bool test_metadata_none = pkg->Param("test_metadata_none"); bool fill_derived = pkg->Param("fill_derived"); if (fill_derived) { @@ -302,6 +317,24 @@ void PreFill(MeshBlockData *rc) { v(out + n, k, j, i) = 1.0 - v(in + n, k, j, i); }); } + + // Fill the metadata::None var with index gymnastics. + if (test_metadata_none) { + const int nx1 = pmb->cellbounds.ncellsi(IndexDomain::interior); + const int nx2 = pmb->cellbounds.ncellsj(IndexDomain::interior); + const int nx3 = pmb->cellbounds.ncellsk(IndexDomain::interior); + + // packing in principle unnecessary/convoluted here and just done for demonstration + std::vector vars({"metadata_none_var"}); + PackIndexMap imap; + const auto &v = rc->PackVariables(vars, imap); + + pmb->par_for( + PARTHENON_AUTO_LABEL, 0, 2, 0, nx3, 0, nx2, 0, nx1, + KOKKOS_LAMBDA(const int n, const int k, const int j, const int i) { + v(n, k, j, i) = n + k + j + i; + }); + } } // this is the package registered function to fill derived diff --git a/example/fine_advection/advection_driver.cpp b/example/fine_advection/advection_driver.cpp index c02b4cd555ff..ca55c04f8278 100644 --- a/example/fine_advection/advection_driver.cpp +++ b/example/fine_advection/advection_driver.cpp @@ -21,6 +21,7 @@ #include "amr_criteria/refinement_package.hpp" #include "bvals/comms/bvals_in_one.hpp" #include "interface/metadata.hpp" +#include "interface/state_descriptor.hpp" #include "interface/update.hpp" #include "mesh/meshblock_pack.hpp" #include "parthenon/driver.hpp" @@ -72,6 +73,12 @@ TaskCollection AdvectionDriver::MakeTaskCollection(BlockList_t &blocks, const in TaskCollection tc; TaskID none(0); + std::shared_ptr pkg = + pmesh->packages.Get("advection_package"); + const auto do_regular_advection = pkg->Param("do_regular_advection"); + const auto do_fine_advection = pkg->Param("do_fine_advection"); + const auto do_CT_advection = pkg->Param("do_CT_advection"); + // Build MeshBlockData containers that will be included in MeshData containers. It is // gross that this has to be done by hand. const auto &stage_name = integrator->stage_name; @@ -118,35 +125,49 @@ TaskCollection AdvectionDriver::MakeTaskCollection(BlockList_t &blocks, const in auto flx = none; auto flx_fine = none; for (auto face : faces) { - flx = flx | tl.AddTask(none, advection_package::CalculateFluxes, desc, - face, parthenon::CellLevel::same, mc0.get()); - flx_fine = flx_fine | - tl.AddTask(none, advection_package::CalculateFluxes, - desc_fine, face, parthenon::CellLevel::fine, mc0.get()); + if (do_regular_advection) { + flx = flx | tl.AddTask(none, advection_package::CalculateFluxes, + desc, face, parthenon::CellLevel::same, mc0.get()); + } + if (do_fine_advection) { + flx_fine = flx_fine | + tl.AddTask(none, advection_package::CalculateFluxes, + desc_fine, face, parthenon::CellLevel::fine, mc0.get()); + } } auto vf_dep = none; - for (auto edge : std::vector{TE::E1, TE::E2, TE::E3}) { - vf_dep = tl.AddTask(vf_dep, advection_package::CalculateVectorFluxes, edge, - parthenon::CellLevel::same, 1.0, mc0.get()); - vf_dep = tl.AddTask(vf_dep, advection_package::CalculateVectorFluxes, edge, - parthenon::CellLevel::same, -1.0, mc0.get()); + if (do_CT_advection) { + for (auto edge : std::vector{TE::E1, TE::E2, TE::E3}) { + vf_dep = tl.AddTask(vf_dep, advection_package::CalculateVectorFluxes, edge, + parthenon::CellLevel::same, 1.0, mc0.get()); + vf_dep = tl.AddTask(vf_dep, advection_package::CalculateVectorFluxes, edge, + parthenon::CellLevel::same, -1.0, mc0.get()); + } } auto set_flx = parthenon::AddFluxCorrectionTasks( start_flxcor | flx | flx_fine | vf_dep, tl, mc0, pmesh->multilevel); - auto update = - AddUpdateTasks(set_flx, tl, parthenon::CellLevel::same, TT::Cell, beta, dt, desc, - mbase.get(), mc0.get(), mdudt.get(), mc1.get()); + auto update = set_flx; + if (do_regular_advection) { + update = AddUpdateTasks(set_flx, tl, parthenon::CellLevel::same, TT::Cell, beta, dt, + desc, mbase.get(), mc0.get(), mdudt.get(), mc1.get()); + } - auto update_fine = - AddUpdateTasks(set_flx, tl, parthenon::CellLevel::fine, TT::Cell, beta, dt, - desc_fine, mbase.get(), mc0.get(), mdudt.get(), mc1.get()); + auto update_fine = set_flx; + if (do_fine_advection) { + update_fine = + AddUpdateTasks(set_flx, tl, parthenon::CellLevel::fine, TT::Cell, beta, dt, + desc_fine, mbase.get(), mc0.get(), mdudt.get(), mc1.get()); + } - auto update_vec = - AddUpdateTasks(set_flx, tl, parthenon::CellLevel::same, TT::Face, beta, dt, - desc_vec, mbase.get(), mc0.get(), mdudt.get(), mc1.get()); + auto update_vec = set_flx; + if (do_CT_advection) { + update_vec = + AddUpdateTasks(set_flx, tl, parthenon::CellLevel::same, TT::Face, beta, dt, + desc_vec, mbase.get(), mc0.get(), mdudt.get(), mc1.get()); + } auto boundaries = parthenon::AddBoundaryExchangeTasks( update | update_vec | update_fine | start_send, tl, mc1, pmesh->multilevel); @@ -155,7 +176,8 @@ TaskCollection AdvectionDriver::MakeTaskCollection(BlockList_t &blocks, const in tl.AddTask(boundaries, parthenon::Update::FillDerived>, mc1.get()); if (stage == integrator->nstages) { - auto new_dt = tl.AddTask(fill_derived, EstimateTimestep>, mc1.get()); + auto dealloc = tl.AddTask(fill_derived, SparseDealloc, mc1.get()); + auto new_dt = tl.AddTask(dealloc, EstimateTimestep>, mc1.get()); if (pmesh->adaptive) { auto tag_refine = tl.AddTask(new_dt, parthenon::Refinement::Tag>, mc1.get()); diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp index 8ca6c6458b20..428d91b9245f 100644 --- a/example/fine_advection/advection_package.cpp +++ b/example/fine_advection/advection_package.cpp @@ -65,34 +65,58 @@ std::shared_ptr Initialize(ParameterInput *pin) { } pkg->AddParam<>("profile", profile_str); - pkg->AddField( - Metadata({Metadata::Cell, Metadata::Fine, Metadata::Independent, - Metadata::WithFluxes, Metadata::FillGhost})); - - pkg->AddField(Metadata({Metadata::Cell, Metadata::Independent, - Metadata::WithFluxes, Metadata::FillGhost})); - pkg->AddField( - Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy})); - - Metadata m( - {Metadata::Face, Metadata::Independent, Metadata::WithFluxes, Metadata::FillGhost}); - m.RegisterRefinementOps(); - pkg->AddField(m); - pkg->AddField(m); - pkg->AddField(Metadata( - {Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, std::vector{4})); - pkg->AddField(Metadata( - {Metadata::Face, Metadata::Derived, Metadata::OneCopy}, std::vector{2})); - pkg->AddField(Metadata( - {Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, std::vector{3})); - pkg->AddField(Metadata( - {Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, std::vector{3})); - pkg->AddField( - Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy})); - pkg->AddField( - Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy})); + bool do_regular_advection = + pin->GetOrAddBoolean("Advection", "do_regular_advection", true); + pkg->AddParam<>("do_regular_advection", do_regular_advection); + if (do_regular_advection) { + int shape_size = pin->GetOrAddInteger("Advection", "shape_size", 1); + int sparse_size = pin->GetOrAddInteger("Advection", "sparse_size", 1); + pkg->AddParam<>("sparse_size", sparse_size); + Real alloc_threshold = pin->GetOrAddReal("Advection", "alloc_threshold", 1.e-6); + Real dealloc_threshold = pin->GetOrAddReal("Advection", "dealloc_threshold", 5.e-7); + Metadata m({Metadata::Cell, Metadata::Independent, Metadata::WithFluxes, + Metadata::FillGhost, Metadata::Sparse}, + std::vector{shape_size}); + m.SetSparseThresholds(alloc_threshold, dealloc_threshold, 0.0); + std::vector sparse_idxs(sparse_size); + std::iota(sparse_idxs.begin(), sparse_idxs.end(), 0); + pkg->AddSparsePool(m, sparse_idxs); + } + + bool do_fine_advection = pin->GetOrAddBoolean("Advection", "do_fine_advection", true); + pkg->AddParam<>("do_fine_advection", do_fine_advection); + if (do_fine_advection) { + pkg->AddField( + Metadata({Metadata::Cell, Metadata::Fine, Metadata::Independent, + Metadata::WithFluxes, Metadata::FillGhost})); + + pkg->AddField( + Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy})); + } + + bool do_CT_advection = pin->GetOrAddBoolean("Advection", "do_CT_advection", true); + pkg->AddParam<>("do_CT_advection", do_CT_advection); + if (do_CT_advection) { + auto m = Metadata({Metadata::Face, Metadata::Independent, Metadata::WithFluxes, + Metadata::FillGhost}); + m.RegisterRefinementOps(); + pkg->AddField(m); + pkg->AddField(m); + pkg->AddField(Metadata( + {Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, std::vector{4})); + pkg->AddField(Metadata( + {Metadata::Face, Metadata::Derived, Metadata::OneCopy}, std::vector{2})); + pkg->AddField(Metadata( + {Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, std::vector{3})); + pkg->AddField(Metadata( + {Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, std::vector{3})); + pkg->AddField( + Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy})); + pkg->AddField( + Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy})); + } pkg->CheckRefinementBlock = CheckRefinement; pkg->EstimateTimestepMesh = EstimateTimestep; @@ -101,35 +125,42 @@ std::shared_ptr Initialize(ParameterInput *pin) { } AmrTag CheckRefinement(MeshBlockData *rc) { - // refine on advected, for example. could also be a derived quantity - static auto desc = parthenon::MakePackDescriptor(rc); - auto pack = desc.GetPack(rc); + std::shared_ptr pkg = + rc->GetMeshPointer()->packages.Get("advection_package"); + auto do_regular_advection = pkg->Param("do_regular_advection"); + if (do_regular_advection) { + // refine on advected, for example. could also be a derived quantity + static auto desc = parthenon::MakePackDescriptor(rc); + auto pack = desc.GetPack(rc); - auto pmb = rc->GetBlockPointer(); - IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); - IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); - IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); + auto pmb = rc->GetBlockPointer(); + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); - typename Kokkos::MinMax::value_type minmax; - parthenon::par_reduce( - parthenon::loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), - pack.GetLowerBoundHost(0), pack.GetUpperBoundHost(0), kb.s, kb.e, jb.s, jb.e, ib.s, - ib.e, - KOKKOS_LAMBDA(const int n, const int k, const int j, const int i, - typename Kokkos::MinMax::value_type &lminmax) { - lminmax.min_val = - (pack(n, k, j, i) < lminmax.min_val ? pack(n, k, j, i) : lminmax.min_val); - lminmax.max_val = - (pack(n, k, j, i) > lminmax.max_val ? pack(n, k, j, i) : lminmax.max_val); - }, - Kokkos::MinMax(minmax)); + typename Kokkos::MinMax::value_type minmax; + parthenon::par_reduce( + parthenon::loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), 0, + pack.GetNBlocks() - 1, // Runs from [0, 0] since pack built from MeshBlockData + pack.GetLowerBoundHost(0), pack.GetUpperBoundHost(0), kb.s, kb.e, jb.s, jb.e, + ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int n, const int k, const int j, const int i, + typename Kokkos::MinMax::value_type &lminmax) { + lminmax.min_val = (pack(b, n, k, j, i) < lminmax.min_val ? pack(b, n, k, j, i) + : lminmax.min_val); + lminmax.max_val = (pack(b, n, k, j, i) > lminmax.max_val ? pack(b, n, k, j, i) + : lminmax.max_val); + }, + Kokkos::MinMax(minmax)); - auto pkg = pmb->packages.Get("advection_package"); - const auto &refine_tol = pkg->Param("refine_tol"); - const auto &derefine_tol = pkg->Param("derefine_tol"); + auto pkg = pmb->packages.Get("advection_package"); + const auto &refine_tol = pkg->Param("refine_tol"); + const auto &derefine_tol = pkg->Param("derefine_tol"); - if (minmax.max_val > refine_tol && minmax.min_val < derefine_tol) return AmrTag::refine; - if (minmax.max_val < derefine_tol) return AmrTag::derefine; + if (minmax.max_val > refine_tol && minmax.min_val < derefine_tol) + return AmrTag::refine; + if (minmax.max_val < derefine_tol) return AmrTag::derefine; + } return AmrTag::same; } @@ -175,74 +206,84 @@ TaskStatus FillDerived(MeshData *md) { md); auto pack = desc.GetPack(md); + std::shared_ptr pkg = + md->GetMeshPointer()->packages.Get("advection_package"); + IndexRange ib = md->GetBoundsI(IndexDomain::interior); IndexRange jb = md->GetBoundsJ(IndexDomain::interior); IndexRange kb = md->GetBoundsK(IndexDomain::interior); const int ndim = md->GetMeshPointer()->ndim; const int nghost = parthenon::Globals::nghost; - parthenon::par_for( - PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { - const int kf = (ndim > 2) ? (k - nghost) * 2 + nghost : k; - const int jf = (ndim > 1) ? (j - nghost) * 2 + nghost : j; - const int fi = (ndim > 0) ? (i - nghost) * 2 + nghost : i; - pack(b, Conserved::phi_fine_restricted(), k, j, i) = 0.0; - Real ntot = 0.0; - for (int ioff = 0; ioff <= (ndim > 0); ++ioff) - for (int joff = 0; joff <= (ndim > 1); ++joff) - for (int koff = 0; koff <= (ndim > 2); ++koff) { - ntot += 1.0; - pack(b, Conserved::phi_fine_restricted(), k, j, i) += - pack(b, Conserved::phi_fine(), kf + koff, jf + joff, fi + ioff); - } - pack(b, Conserved::phi_fine_restricted(), k, j, i) /= ntot; - }); - - using TE = parthenon::TopologicalElement; - parthenon::par_for( - PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { - pack(b, Conserved::C_cc(0), k, j, i) = - 0.5 * (pack(b, TE::F1, Conserved::C(), k, j, i) + - pack(b, TE::F1, Conserved::C(), k, j, i + (ndim > 0))); - pack(b, Conserved::C_cc(1), k, j, i) = - 0.5 * (pack(b, TE::F2, Conserved::C(), k, j, i) + - pack(b, TE::F2, Conserved::C(), k, j + (ndim > 1), i)); - pack(b, Conserved::C_cc(2), k, j, i) = - 0.5 * (pack(b, TE::F3, Conserved::C(), k, j, i) + - pack(b, TE::F3, Conserved::C(), k + (ndim > 2), j, i)); - auto &coords = pack.GetCoordinates(b); - pack(b, Conserved::divC(), k, j, i) = - (pack(b, TE::F1, Conserved::C(), k, j, i + (ndim > 0)) - - pack(b, TE::F1, Conserved::C(), k, j, i)) / - coords.Dxc(k, j, i) + - (pack(b, TE::F2, Conserved::C(), k, j + (ndim > 1), i) - - pack(b, TE::F2, Conserved::C(), k, j, i)) / - coords.Dxc(k, j, i) + - (pack(b, TE::F3, Conserved::C(), k + (ndim > 2), j, i) - - pack(b, TE::F3, Conserved::C(), k, j, i)) / - coords.Dxc(k, j, i); - - pack(b, Conserved::D_cc(0), k, j, i) = - 0.5 * (pack(b, TE::F1, Conserved::D(), k, j, i) + - pack(b, TE::F1, Conserved::D(), k, j, i + (ndim > 0))); - pack(b, Conserved::D_cc(1), k, j, i) = - 0.5 * (pack(b, TE::F2, Conserved::D(), k, j, i) + - pack(b, TE::F2, Conserved::D(), k, j + (ndim > 1), i)); - pack(b, Conserved::D_cc(2), k, j, i) = - 0.5 * (pack(b, TE::F3, Conserved::D(), k, j, i) + - pack(b, TE::F3, Conserved::D(), k + (ndim > 2), j, i)); - pack(b, Conserved::divD(), k, j, i) = - (pack(b, TE::F1, Conserved::D(), k, j, i + (ndim > 0)) - - pack(b, TE::F1, Conserved::D(), k, j, i)) / - coords.Dxc(k, j, i) + - (pack(b, TE::F2, Conserved::D(), k, j + (ndim > 1), i) - - pack(b, TE::F2, Conserved::D(), k, j, i)) / - coords.Dxc(k, j, i) + - (pack(b, TE::F3, Conserved::D(), k + (ndim > 2), j, i) - - pack(b, TE::F3, Conserved::D(), k, j, i)) / - coords.Dxc(k, j, i); - }); + + auto do_fine_advection = pkg->Param("do_fine_advection"); + if (do_fine_advection) { + parthenon::par_for( + PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, + ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { + const int kf = (ndim > 2) ? (k - nghost) * 2 + nghost : k; + const int jf = (ndim > 1) ? (j - nghost) * 2 + nghost : j; + const int fi = (ndim > 0) ? (i - nghost) * 2 + nghost : i; + pack(b, Conserved::phi_fine_restricted(), k, j, i) = 0.0; + Real ntot = 0.0; + for (int koff = 0; koff <= (ndim > 2); ++koff) + for (int joff = 0; joff <= (ndim > 1); ++joff) + for (int ioff = 0; ioff <= (ndim > 0); ++ioff) { + ntot += 1.0; + pack(b, Conserved::phi_fine_restricted(), k, j, i) += + pack(b, Conserved::phi_fine(), kf + koff, jf + joff, fi + ioff); + } + pack(b, Conserved::phi_fine_restricted(), k, j, i) /= ntot; + }); + } + + auto do_CT_advection = pkg->Param("do_CT_advection"); + if (do_CT_advection) { + using TE = parthenon::TopologicalElement; + parthenon::par_for( + PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, + ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { + pack(b, Conserved::C_cc(0), k, j, i) = + 0.5 * (pack(b, TE::F1, Conserved::C(), k, j, i) + + pack(b, TE::F1, Conserved::C(), k, j, i + (ndim > 0))); + pack(b, Conserved::C_cc(1), k, j, i) = + 0.5 * (pack(b, TE::F2, Conserved::C(), k, j, i) + + pack(b, TE::F2, Conserved::C(), k, j + (ndim > 1), i)); + pack(b, Conserved::C_cc(2), k, j, i) = + 0.5 * (pack(b, TE::F3, Conserved::C(), k, j, i) + + pack(b, TE::F3, Conserved::C(), k + (ndim > 2), j, i)); + auto &coords = pack.GetCoordinates(b); + pack(b, Conserved::divC(), k, j, i) = + (pack(b, TE::F1, Conserved::C(), k, j, i + (ndim > 0)) - + pack(b, TE::F1, Conserved::C(), k, j, i)) / + coords.Dxc(k, j, i) + + (pack(b, TE::F2, Conserved::C(), k, j + (ndim > 1), i) - + pack(b, TE::F2, Conserved::C(), k, j, i)) / + coords.Dxc(k, j, i) + + (pack(b, TE::F3, Conserved::C(), k + (ndim > 2), j, i) - + pack(b, TE::F3, Conserved::C(), k, j, i)) / + coords.Dxc(k, j, i); + + pack(b, Conserved::D_cc(0), k, j, i) = + 0.5 * (pack(b, TE::F1, Conserved::D(), k, j, i) + + pack(b, TE::F1, Conserved::D(), k, j, i + (ndim > 0))); + pack(b, Conserved::D_cc(1), k, j, i) = + 0.5 * (pack(b, TE::F2, Conserved::D(), k, j, i) + + pack(b, TE::F2, Conserved::D(), k, j + (ndim > 1), i)); + pack(b, Conserved::D_cc(2), k, j, i) = + 0.5 * (pack(b, TE::F3, Conserved::D(), k, j, i) + + pack(b, TE::F3, Conserved::D(), k + (ndim > 2), j, i)); + pack(b, Conserved::divD(), k, j, i) = + (pack(b, TE::F1, Conserved::D(), k, j, i + (ndim > 0)) - + pack(b, TE::F1, Conserved::D(), k, j, i)) / + coords.Dxc(k, j, i) + + (pack(b, TE::F2, Conserved::D(), k, j + (ndim > 1), i) - + pack(b, TE::F2, Conserved::D(), k, j, i)) / + coords.Dxc(k, j, i) + + (pack(b, TE::F3, Conserved::D(), k + (ndim > 2), j, i) - + pack(b, TE::F3, Conserved::D(), k, j, i)) / + coords.Dxc(k, j, i); + }); + } return TaskStatus::complete; } } // namespace advection_package diff --git a/example/fine_advection/advection_package.hpp b/example/fine_advection/advection_package.hpp index 87e0c0b2edab..12b6b9b65aad 100644 --- a/example/fine_advection/advection_package.hpp +++ b/example/fine_advection/advection_package.hpp @@ -19,6 +19,7 @@ #include #include +#include #include #define VARIABLE(ns, varname) \ @@ -79,13 +80,20 @@ TaskStatus CalculateFluxes(pack_desc_t &desc, parthenon::TopologicalElement FACE IndexRange ib = md->GetBoundsI(cl, IndexDomain::interior, FACE); IndexRange jb = md->GetBoundsJ(cl, IndexDomain::interior, FACE); IndexRange kb = md->GetBoundsK(cl, IndexDomain::interior, FACE); - parthenon::par_for( - PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, pack.GetLowerBoundHost(0), - pack.GetUpperBoundHost(0), // Warning: only works for dense variables - kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int b, const int l, const int k, const int j, const int i) { - // Calculate the flux using upwind donor cell reconstruction - pack.flux(b, FACE, l, k, j, i) = v * pack(b, l, k + koff, j + joff, i + ioff); + constexpr int unused_scratch_size = 0; + constexpr int unused_scratch_level = 1; + parthenon::par_for_outer( + PARTHENON_AUTO_LABEL, unused_scratch_size, unused_scratch_level, 0, + pack.GetNBlocks() - 1, kb.s, kb.e, + KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) { + parthenon::Indexer2D idxer({jb.s, jb.e}, {ib.s, ib.e}); + for (int l = pack.GetLowerBound(b); l <= pack.GetUpperBound(b); ++l) { + parthenon::par_for_inner(member, 0, idxer.size() - 1, [&](const int idx) { + const auto [j, i] = idxer(idx); + // Calculate the flux using upwind donor cell reconstruction + pack.flux(b, FACE, l, k, j, i) = v * pack(b, l, k + koff, j + joff, i + ioff); + }); + } }); return TaskStatus::complete; } @@ -116,23 +124,33 @@ TaskStatus CalculateVectorFluxes(parthenon::TopologicalElement edge, int koff = edge == TE::E3 ? ndim > 2 : 0; int joff = edge == TE::E2 ? ndim > 1 : 0; int ioff = edge == TE::E1 ? ndim > 0 : 0; - parthenon::par_for( - PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, kb.s - (ndim > 2), - kb.e + (ndim > 2), jb.s - (ndim > 1), jb.e + (ndim > 1), ib.s - 1, ib.e + 1, - KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { - // Piecewise linear in the orthogonal directions, could do something better here - pack(b, TE::CC, recon(0), k, j, i) = - 0.5 * (pack(b, fe, flux_var(), k, j, i) + - pack(b, fe, flux_var(), k + koff, j + joff, i + ioff)); - pack(b, TE::CC, recon(1), k, j, i) = - 0.5 * (pack(b, fe, flux_var(), k, j, i) + - pack(b, fe, flux_var(), k + koff, j + joff, i + ioff)); - pack(b, TE::CC, recon(2), k, j, i) = - 0.5 * (pack(b, fe, flux_var(), k, j, i) + - pack(b, fe, flux_var(), k + koff, j + joff, i + ioff)); - pack(b, TE::CC, recon(3), k, j, i) = - 0.5 * (pack(b, fe, flux_var(), k, j, i) + - pack(b, fe, flux_var(), k + koff, j + joff, i + ioff)); + constexpr int scratch_size = 0; + constexpr int scratch_level = 1; + parthenon::par_for_outer( + PARTHENON_AUTO_LABEL, scratch_size, scratch_level, 0, pack.GetNBlocks() - 1, + kb.s - (ndim > 2), kb.e + (ndim > 2), + KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) { + parthenon::Indexer2D idxer({jb.s - (ndim > 1), jb.e + (ndim > 1)}, + {ib.s - (ndim > 0), ib.e + (ndim > 0)}); + if (pack.GetLowerBound(b, flux_var()) <= pack.GetUpperBound(b, flux_var())) { + parthenon::par_for_inner(member, 0, idxer.size() - 1, [&](const int idx) { + const auto [j, i] = idxer(idx); + // Piecewise linear in the orthogonal directions, could do something better + // here + pack(b, TE::CC, recon(0), k, j, i) = + 0.5 * (pack(b, fe, flux_var(), k, j, i) + + pack(b, fe, flux_var(), k + koff, j + joff, i + ioff)); + pack(b, TE::CC, recon(1), k, j, i) = + 0.5 * (pack(b, fe, flux_var(), k, j, i) + + pack(b, fe, flux_var(), k + koff, j + joff, i + ioff)); + pack(b, TE::CC, recon(2), k, j, i) = + 0.5 * (pack(b, fe, flux_var(), k, j, i) + + pack(b, fe, flux_var(), k + koff, j + joff, i + ioff)); + pack(b, TE::CC, recon(3), k, j, i) = + 0.5 * (pack(b, fe, flux_var(), k, j, i) + + pack(b, fe, flux_var(), k + koff, j + joff, i + ioff)); + }); + } }); // 2. Calculate the quartant averaged flux @@ -143,20 +161,26 @@ TaskStatus CalculateVectorFluxes(parthenon::TopologicalElement edge, ib = md->GetBoundsI(cl, IndexDomain::interior, edge); jb = md->GetBoundsJ(cl, IndexDomain::interior, edge); kb = md->GetBoundsK(cl, IndexDomain::interior, edge); - parthenon::par_for( - PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { - pack.flux(b, edge, var(), k, j, i) = 0.0; - for (int dk = -koff; dk < 1; ++dk) - for (int dj = -joff; dj < 1; ++dj) - for (int di = -ioff; di < 1; ++di) { - // TODO(LFR): Pick out the correct component of the reconstructed flux, - // current version is not an issue for piecewise constant flux though. - pack.flux(b, edge, var(), k, j, i) += - pack(b, TE::CC, recon(0), k + dk, j + dj, i + di); - } - pack.flux(b, edge, var(), k, j, i) /= npoints; - pack.flux(b, edge, var(), k, j, i) *= fac; + parthenon::par_for_outer( + PARTHENON_AUTO_LABEL, scratch_size, scratch_level, 0, pack.GetNBlocks() - 1, kb.s, + kb.e, KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) { + parthenon::Indexer2D idxer({jb.s, jb.e}, {ib.s, ib.e}); + if (pack.GetLowerBound(b, var()) <= pack.GetUpperBound(b, var())) { + parthenon::par_for_inner(member, 0, idxer.size() - 1, [&](const int idx) { + const auto [j, i] = idxer(idx); + pack.flux(b, edge, var(), k, j, i) = 0.0; + for (int dk = -koff; dk < 1; ++dk) + for (int dj = -joff; dj < 1; ++dj) + for (int di = -ioff; di < 1; ++di) { + // TODO(LFR): Pick out the correct component of the reconstructed flux, + // current version is not an issue for piecewise constant flux though. + pack.flux(b, edge, var(), k, j, i) += + pack(b, TE::CC, recon(0), k + dk, j + dj, i + di); + } + pack.flux(b, edge, var(), k, j, i) /= npoints; + pack.flux(b, edge, var(), k, j, i) *= fac; + }); + } }); // 3. Reconstruct the transverse components of the advected field at the edge @@ -167,13 +191,21 @@ TaskStatus CalculateVectorFluxes(parthenon::TopologicalElement edge, ib = md->GetBoundsI(cl, IndexDomain::interior, f); jb = md->GetBoundsJ(cl, IndexDomain::interior, f); kb = md->GetBoundsK(cl, IndexDomain::interior, f); - parthenon::par_for( - PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, kb.s - (ndim > 2), - kb.e + (ndim > 2), jb.s - (ndim > 1), jb.e + (ndim > 1), ib.s - 1, ib.e + 1, - KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { - // Piecewise linear in the orthogonal directions, could do something better here - pack(b, f, recon_f(0), k, j, i) = pack(b, f, var(), k, j, i); - pack(b, f, recon_f(1), k, j, i) = pack(b, f, var(), k, j, i); + parthenon::par_for_outer( + PARTHENON_AUTO_LABEL, scratch_size, scratch_level, 0, pack.GetNBlocks() - 1, + kb.s - (ndim > 2), kb.e + (ndim > 2), + KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) { + parthenon::Indexer2D idxer({jb.s - (ndim > 1), jb.e + (ndim > 1)}, + {ib.s - (ndim > 0), ib.e + (ndim > 0)}); + if (pack.GetLowerBound(b, var()) <= pack.GetUpperBound(b, var())) { + parthenon::par_for_inner(member, 0, idxer.size() - 1, [&](const int idx) { + const auto [j, i] = idxer(idx); + // Piecewise linear in the orthogonal directions, could do something better + // here + pack(b, f, recon_f(0), k, j, i) = pack(b, f, var(), k, j, i); + pack(b, f, recon_f(1), k, j, i) = pack(b, f, var(), k, j, i); + }); + } }); } @@ -189,15 +221,21 @@ TaskStatus CalculateVectorFluxes(parthenon::TopologicalElement edge, d1[static_cast(f1) % 3] = 0; d2[static_cast(edge) % 3] = 0; d2[static_cast(f2) % 3] = 0; - parthenon::par_for( - PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { - pack.flux(b, edge, var(), k, j, i) += - 0.5 * (pack(b, f1, recon_f(0), k, j, i) - - pack(b, f1, recon_f(1), k - d1[2], j - d1[1], i - d1[0])); - pack.flux(b, edge, var(), k, j, i) -= - 0.5 * (pack(b, f2, recon_f(0), k, j, i) - - pack(b, f2, recon_f(1), k - d2[2], j - d2[1], i - d2[0])); + parthenon::par_for_outer( + PARTHENON_AUTO_LABEL, scratch_size, scratch_level, 0, pack.GetNBlocks() - 1, kb.s, + kb.e, KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) { + parthenon::Indexer2D idxer({jb.s, jb.e}, {ib.s, ib.e}); + if (pack.GetLowerBound(b, var()) <= pack.GetUpperBound(b, var())) { + parthenon::par_for_inner(member, 0, idxer.size() - 1, [&](const int idx) { + const auto [j, i] = idxer(idx); + pack.flux(b, edge, var(), k, j, i) += + 0.5 * (pack(b, f1, recon_f(0), k, j, i) - + pack(b, f1, recon_f(1), k - d1[2], j - d1[1], i - d1[0])); + pack.flux(b, edge, var(), k, j, i) -= + 0.5 * (pack(b, f2, recon_f(0), k, j, i) - + pack(b, f2, recon_f(1), k - d2[2], j - d2[1], i - d2[0])); + }); + } }); // Add in the diffusive component diff --git a/example/fine_advection/parthenon_app_inputs.cpp b/example/fine_advection/parthenon_app_inputs.cpp index ab970175e9cf..8025970f35ca 100644 --- a/example/fine_advection/parthenon_app_inputs.cpp +++ b/example/fine_advection/parthenon_app_inputs.cpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2020-2021. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC @@ -50,6 +50,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto pkg = pmb->packages.Get("advection_package"); const auto &profile = pkg->Param("profile"); + auto do_regular_advection = pkg->Param("do_regular_advection"); + auto do_fine_advection = pkg->Param("do_fine_advection"); + auto do_CT_advection = pkg->Param("do_CT_advection"); auto cellbounds = pmb->cellbounds; IndexRange ib = cellbounds.GetBoundsI(IndexDomain::interior); @@ -59,6 +62,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto coords = pmb->coords; using namespace advection_package::Conserved; + if (do_regular_advection) { + const int sparse_size = pkg->Param("sparse_size"); + for (int s = 0; s < sparse_size; ++s) + pmb->AllocSparseID(phi::name(), s); + } static auto desc = parthenon::MakePackDescriptor(data.get()); auto pack = desc.GetPack(data.get()); @@ -68,78 +76,105 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { if (profile == "hard_sphere") profile_type = 2; if (profile == "block") profile_type = 3; Real amp = 1.0; - const int b = 0; const int ndim = pmb->pmy_mesh->ndim; const int nghost = parthenon::Globals::nghost; - parthenon::par_for( - PARTHENON_AUTO_LABEL, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int k, const int j, const int i) { - const int kf = ndim > 2 ? (k - nghost) * 2 + nghost : k; - const int jf = ndim > 1 ? (j - nghost) * 2 + nghost : j; - const int fi = ndim > 0 ? (i - nghost) * 2 + nghost : i; - if (profile_type == 0) { - PARTHENON_FAIL("Profile type wave not implemented."); - } else if (profile_type == 1) { - Real rsq = coords.Xc<1>(i) * coords.Xc<1>(i) + - coords.Xc<2>(j) * coords.Xc<2>(j) + - coords.Xc<3>(k) * coords.Xc<3>(k); - pack(b, phi(), k, j, i) = 1. + amp * exp(-100.0 * rsq); - for (int ioff = 0; ioff <= (ndim > 0); ++ioff) - for (int joff = 0; joff <= (ndim > 1); ++joff) - for (int koff = 0; koff <= (ndim > 2); ++koff) { - pack(b, phi_fine(), kf + koff, jf + joff, fi + ioff) = - 1. + amp * exp(-100.0 * rsq); - } - } else if (profile_type == 2) { - Real rsq = coords.Xc<1>(i) * coords.Xc<1>(i) + - coords.Xc<2>(j) * coords.Xc<2>(j) + - coords.Xc<3>(k) * coords.Xc<3>(k); - pack(b, phi(), k, j, i) = (rsq < 0.15 * 0.15 ? 1.0 : 0.0); - for (int ioff = 0; ioff <= (ndim > 0); ++ioff) - for (int joff = 0; joff <= (ndim > 1); ++joff) - for (int koff = 0; koff <= (ndim > 2); ++koff) { - pack(b, phi_fine(), kf + koff, jf + joff, fi + ioff) = - (rsq < 0.15 * 0.15 ? 1.0 : 0.0); - } - } else { - pack(b, phi(), k, j, i) = 0.0; - for (int ioff = 0; ioff <= (ndim > 0); ++ioff) - for (int joff = 0; joff <= (ndim > 1); ++joff) - for (int koff = 0; koff <= (ndim > 2); ++koff) { - pack(b, phi_fine(), kf + koff, jf + joff, fi + ioff) = 0.0; - } - } - }); - - ib = cellbounds.GetBoundsI(IndexDomain::interior, TE::F2); - jb = cellbounds.GetBoundsJ(IndexDomain::interior, TE::F2); - kb = cellbounds.GetBoundsK(IndexDomain::interior, TE::F2); - const Real x0 = 0.2; - parthenon::par_for( - PARTHENON_AUTO_LABEL, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int k, const int j, const int i) { - auto &coords = pack.GetCoordinates(b); - Real xlo = coords.X(k, j, i); - Real xhi = coords.X(k, j, i + 1); - Real Alo = std::abs(xlo) < x0 ? sign(xlo, xhi) * (x0 - std::abs(xlo)) : 0.0; - Real Ahi = std::abs(xhi) < x0 ? sign(xhi, xlo) * (x0 - std::abs(xhi)) : 0.0; - pack(b, TE::F2, C(), k, j, i) = (Alo - Ahi) / (xhi - xlo); - }); - - ib = cellbounds.GetBoundsI(IndexDomain::interior, TE::F3); - jb = cellbounds.GetBoundsJ(IndexDomain::interior, TE::F3); - kb = cellbounds.GetBoundsK(IndexDomain::interior, TE::F3); - parthenon::par_for( - PARTHENON_AUTO_LABEL, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int k, const int j, const int i) { - auto &coords = pack.GetCoordinates(b); - Real xlo = coords.X(k, j, i); - Real xhi = coords.X(k, j, i + 1); - Real Alo = std::abs(xlo) < x0 ? sign(xlo, xhi) * (x0 - std::abs(xlo)) : 0.0; - Real Ahi = std::abs(xhi) < x0 ? sign(xhi, xlo) * (x0 - std::abs(xhi)) : 0.0; - pack(b, TE::F3, D(), k, j, i) = (Alo - Ahi) / (xhi - xlo); - }); + + if (do_regular_advection) { + parthenon::par_for( + PARTHENON_AUTO_LABEL, pack.GetLowerBoundHost(b, phi()), + pack.GetUpperBoundHost(b, phi()), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int l, const int k, const int j, const int i) { + const int kf = ndim > 2 ? (k - nghost) * 2 + nghost : k; + const int jf = ndim > 1 ? (j - nghost) * 2 + nghost : j; + const int fi = ndim > 0 ? (i - nghost) * 2 + nghost : i; + if (profile_type == 0) { + PARTHENON_FAIL("Profile type wave not implemented."); + } else if (profile_type == 1) { + Real rsq = coords.Xc<1>(i) * coords.Xc<1>(i) + + coords.Xc<2>(j) * coords.Xc<2>(j) + + coords.Xc<3>(k) * coords.Xc<3>(k); + pack(b, phi(l), k, j, i) = 1. + amp * exp(-100.0 * rsq); + } else if (profile_type == 2) { + Real rsq = coords.Xc<1>(i) * coords.Xc<1>(i) + + coords.Xc<2>(j) * coords.Xc<2>(j) + + coords.Xc<3>(k) * coords.Xc<3>(k); + pack(b, phi(l), k, j, i) = (rsq < 0.15 * 0.15 ? 1.0 : 0.0); + } else { + pack(b, phi(l), k, j, i) = 0.0; + } + }); + } + + if (do_fine_advection) { + parthenon::par_for( + PARTHENON_AUTO_LABEL, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + const int kf = ndim > 2 ? (k - nghost) * 2 + nghost : k; + const int jf = ndim > 1 ? (j - nghost) * 2 + nghost : j; + const int fi = ndim > 0 ? (i - nghost) * 2 + nghost : i; + if (profile_type == 0) { + PARTHENON_FAIL("Profile type wave not implemented."); + } else if (profile_type == 1) { + Real rsq = coords.Xc<1>(i) * coords.Xc<1>(i) + + coords.Xc<2>(j) * coords.Xc<2>(j) + + coords.Xc<3>(k) * coords.Xc<3>(k); + for (int koff = 0; koff <= (ndim > 2); ++koff) + for (int joff = 0; joff <= (ndim > 1); ++joff) + for (int ioff = 0; ioff <= (ndim > 0); ++ioff) { + pack(b, phi_fine(), kf + koff, jf + joff, fi + ioff) = + 1. + amp * exp(-100.0 * rsq); + } + } else if (profile_type == 2) { + Real rsq = coords.Xc<1>(i) * coords.Xc<1>(i) + + coords.Xc<2>(j) * coords.Xc<2>(j) + + coords.Xc<3>(k) * coords.Xc<3>(k); + for (int koff = 0; koff <= (ndim > 2); ++koff) + for (int joff = 0; joff <= (ndim > 1); ++joff) + for (int ioff = 0; ioff <= (ndim > 0); ++ioff) { + pack(b, phi_fine(), kf + koff, jf + joff, fi + ioff) = + (rsq < 0.15 * 0.15 ? 1.0 : 0.0); + } + } else { + for (int koff = 0; koff <= (ndim > 2); ++koff) + for (int joff = 0; joff <= (ndim > 1); ++joff) + for (int ioff = 0; ioff <= (ndim > 0); ++ioff) { + pack(b, phi_fine(), kf + koff, jf + joff, fi + ioff) = 0.0; + } + } + }); + } + + if (do_CT_advection) { + ib = cellbounds.GetBoundsI(IndexDomain::interior, TE::F2); + jb = cellbounds.GetBoundsJ(IndexDomain::interior, TE::F2); + kb = cellbounds.GetBoundsK(IndexDomain::interior, TE::F2); + const Real x0 = 0.2; + parthenon::par_for( + PARTHENON_AUTO_LABEL, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + auto &coords = pack.GetCoordinates(b); + Real xlo = coords.X(k, j, i); + Real xhi = coords.X(k, j, i + 1); + Real Alo = std::abs(xlo) < x0 ? sign(xlo, xhi) * (x0 - std::abs(xlo)) : 0.0; + Real Ahi = std::abs(xhi) < x0 ? sign(xhi, xlo) * (x0 - std::abs(xhi)) : 0.0; + pack(b, TE::F2, C(), k, j, i) = (Alo - Ahi) / (xhi - xlo); + }); + + ib = cellbounds.GetBoundsI(IndexDomain::interior, TE::F3); + jb = cellbounds.GetBoundsJ(IndexDomain::interior, TE::F3); + kb = cellbounds.GetBoundsK(IndexDomain::interior, TE::F3); + parthenon::par_for( + PARTHENON_AUTO_LABEL, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + auto &coords = pack.GetCoordinates(b); + Real xlo = coords.X(k, j, i); + Real xhi = coords.X(k, j, i + 1); + Real Alo = std::abs(xlo) < x0 ? sign(xlo, xhi) * (x0 - std::abs(xlo)) : 0.0; + Real Ahi = std::abs(xhi) < x0 ? sign(xhi, xlo) * (x0 - std::abs(xhi)) : 0.0; + pack(b, TE::F3, D(), k, j, i) = (Alo - Ahi) / (xhi - xlo); + }); + } } Packages_t ProcessPackages(std::unique_ptr &pin) { diff --git a/example/fine_advection/parthinput.advection b/example/fine_advection/parthinput.advection index 31ee255c5053..b3bf7ca10d0c 100644 --- a/example/fine_advection/parthinput.advection +++ b/example/fine_advection/parthinput.advection @@ -57,6 +57,13 @@ profile = hard_sphere refine_tol = 0.3 # control the package specific refinement tagging function derefine_tol = 0.03 +# Include advection equation for regular resolution CC fields +do_regular_advection = true +# Include advection equation for fine resolution CC fields +do_fine_advection = true +# Include vector wave equation using CT +do_CT_advection = true + file_type = rst dt = 0.05 diff --git a/example/fine_advection/stokes.hpp b/example/fine_advection/stokes.hpp index 27b529def8a5..fdc3ec408402 100644 --- a/example/fine_advection/stokes.hpp +++ b/example/fine_advection/stokes.hpp @@ -19,6 +19,7 @@ #include #include +#include namespace advection_example { using namespace parthenon::driver::prelude; @@ -36,21 +37,27 @@ TaskStatus WeightedSumDataElement(parthenon::CellLevel cl, IndexRange jb = in1->GetBoundsJ(cl, IndexDomain::entire, te); IndexRange kb = in1->GetBoundsK(cl, IndexDomain::entire, te); - PARTHENON_REQUIRE(pack1.GetLowerBoundHost(0) == pack2.GetLowerBoundHost(0), - "Packs are different size."); - PARTHENON_REQUIRE(pack1.GetLowerBoundHost(0) == pack_out.GetLowerBoundHost(0), - "Packs are different size."); - PARTHENON_REQUIRE(pack1.GetUpperBoundHost(0) == pack2.GetUpperBoundHost(0), - "Packs are different size."); - PARTHENON_REQUIRE(pack1.GetUpperBoundHost(0) == pack_out.GetUpperBoundHost(0), - "Packs are different size."); - parthenon::par_for( - PARTHENON_AUTO_LABEL, 0, pack1.GetNBlocks() - 1, pack1.GetLowerBoundHost(0), - pack1.GetUpperBoundHost(0), // This is safe for dense vars - kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int b, const int l, const int k, const int j, const int i) { - pack_out(b, te, l, k, j, i) = - w1 * pack1(b, te, l, k, j, i) + w2 * pack2(b, te, l, k, j, i); + constexpr int scratch_size = 0; + constexpr int scratch_level = 1; + parthenon::par_for_outer( + PARTHENON_AUTO_LABEL, scratch_size, scratch_level, 0, pack1.GetNBlocks() - 1, kb.s, + kb.e, KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) { + parthenon::Indexer2D idxer({jb.s, jb.e}, {ib.s, ib.e}); + PARTHENON_REQUIRE(pack1.GetLowerBound(b) == pack2.GetLowerBound(b), + "Packs are different size."); + PARTHENON_REQUIRE(pack1.GetLowerBound(b) == pack_out.GetLowerBound(b), + "Packs are different size."); + PARTHENON_REQUIRE(pack1.GetUpperBound(b) == pack2.GetUpperBound(b), + "Packs are different size."); + PARTHENON_REQUIRE(pack1.GetUpperBound(b) == pack_out.GetUpperBound(b), + "Packs are different size."); + for (int l = pack1.GetLowerBound(b); l <= pack1.GetUpperBound(b); ++l) { + parthenon::par_for_inner(member, 0, idxer.size() - 1, [&](const int idx) { + const auto [j, i] = idxer(idx); + pack_out(b, te, l, k, j, i) = + w1 * pack1(b, te, l, k, j, i) + w2 * pack2(b, te, l, k, j, i); + }); + } }); return TaskStatus::complete; } @@ -73,12 +80,18 @@ void StokesZero(parthenon::CellLevel cl, parthenon::TopologicalElement TeVar, IndexRange jb = out->GetBoundsJ(cl, IndexDomain::interior, TeVar); IndexRange kb = out->GetBoundsK(cl, IndexDomain::interior, TeVar); - parthenon::par_for( - PARTHENON_AUTO_LABEL, 0, pack_out.GetNBlocks() - 1, pack_out.GetLowerBoundHost(0), - pack_out.GetUpperBoundHost(0), // This is safe for dense vars only - kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int b, const int l, const int k, const int j, const int i) { - pack_out(b, TeVar, l, k, j, i) = 0.0; + constexpr int scratch_size = 0; + constexpr int scratch_level = 1; + parthenon::par_for_outer( + PARTHENON_AUTO_LABEL, scratch_size, scratch_level, 0, pack_out.GetNBlocks() - 1, + kb.s, kb.e, KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) { + parthenon::Indexer2D idxer({jb.s, jb.e}, {ib.s, ib.e}); + for (int l = pack_out.GetLowerBound(b); l <= pack_out.GetUpperBound(b); ++l) { + parthenon::par_for_inner(member, 0, idxer.size() - 1, [&](const int idx) { + const auto [j, i] = idxer(idx); + pack_out(b, TeVar, l, k, j, i) = 0.0; + }); + } }); } @@ -103,22 +116,29 @@ void StokesComponent(Real fac, parthenon::CellLevel cl, koff = ndim > 2 ? koff : 0; joff = ndim > 1 ? joff : 0; - PARTHENON_REQUIRE(pack_in.GetLowerBoundHost(0) == pack_out.GetLowerBoundHost(0), - "Packs are different size."); - PARTHENON_REQUIRE(pack_in.GetUpperBoundHost(0) == pack_out.GetUpperBoundHost(0), - "Packs are different size."); - parthenon::par_for( - PARTHENON_AUTO_LABEL, 0, pack_in.GetNBlocks() - 1, pack_in.GetLowerBoundHost(0), - pack_in.GetUpperBoundHost(0), // This is safe for dense vars only - kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int b, const int l, const int k, const int j, const int i) { + constexpr int scratch_size = 0; + constexpr int scratch_level = 1; + parthenon::par_for_outer( + PARTHENON_AUTO_LABEL, scratch_size, scratch_level, 0, pack_out.GetNBlocks() - 1, + kb.s, kb.e, KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) { auto &coords = pack_in.GetCoordinates(b); - pack_out(b, TeVar, l, k, j, i) += - fac * - (coords.Volume(cl, TeFlux, k, j, i) * pack_in.flux(b, TeFlux, l, k, j, i) - - coords.Volume(cl, TeFlux, k + koff, j + joff, i + ioff) * - pack_in.flux(b, TeFlux, l, k + koff, j + joff, i + ioff)) / - coords.Volume(cl, TeVar, k, j, i); + parthenon::Indexer2D idxer({jb.s, jb.e}, {ib.s, ib.e}); + PARTHENON_REQUIRE(pack_in.GetLowerBound(b) == pack_out.GetLowerBound(b), + "Packs are different size."); + PARTHENON_REQUIRE(pack_in.GetUpperBound(b) == pack_out.GetUpperBound(b), + "Packs are different size."); + for (int l = pack_out.GetLowerBound(b); l <= pack_out.GetUpperBound(b); ++l) { + parthenon::par_for_inner(member, 0, idxer.size() - 1, [&](const int idx) { + const auto [j, i] = idxer(idx); + pack_out(b, TeVar, l, k, j, i) += + fac * + (coords.Volume(cl, TeFlux, k, j, i) * + pack_in.flux(b, TeFlux, l, k, j, i) - + coords.Volume(cl, TeFlux, k + koff, j + joff, i + ioff) * + pack_in.flux(b, TeFlux, l, k + koff, j + joff, i + ioff)) / + coords.Volume(cl, TeVar, k, j, i); + }); + } }); } diff --git a/src/outputs/history.cpp b/src/outputs/history.cpp index 2e1310062d48..074bea4feb0c 100644 --- a/src/outputs/history.cpp +++ b/src/outputs/history.cpp @@ -18,6 +18,7 @@ // \brief writes history output data, volume-averaged quantities that are output // frequently in time to trace their history. +#include #include #include #include @@ -93,9 +94,16 @@ void HistoryOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, md_base->Initialize(pm->block_list, pm); } - // Loop over all packages of the application - for (const auto &pkg : packages) { - const auto ¶ms = pkg.second->AllParams(); + // Loop over all packages of the application in alphabetical order to ensure consistency + // of ordering of data in columns. + std::vector keys; + for (const auto &pair : packages) { + keys.push_back(pair.first); + } + std::sort(keys.begin(), keys.end()); + for (const auto &key : keys) { + const auto &pkg = packages[key]; + const auto ¶ms = pkg->AllParams(); // Check if the package has enrolled scalar history functions which are stored in the // Params under the `hist_param_key` name. diff --git a/src/outputs/output_utils.hpp b/src/outputs/output_utils.hpp index 5c094ef688d7..e9fcee04c986 100644 --- a/src/outputs/output_utils.hpp +++ b/src/outputs/output_utils.hpp @@ -312,11 +312,11 @@ void PackOrUnpackVar(const VarInfo &info, bool do_ghosts, idx_t &idx, Function_t auto [kb, jb, ib] = info.GetPaddedBoundsKJI(domain); if (info.where == MetadataFlag({Metadata::None})) { kb.s = 0; - kb.e = shape[4]; + kb.e = std::max(0, shape[4] - 1); jb.s = 0; - jb.e = shape[5]; + jb.e = std::max(0, shape[5] - 1); ib.s = 0; - ib.e = shape[6]; + ib.e = std::max(0, shape[6] - 1); } for (int topo = 0; topo < shape[0]; ++topo) { for (int t = 0; t < shape[1]; ++t) { diff --git a/tst/regression/test_suites/output_hdf5/output_hdf5.py b/tst/regression/test_suites/output_hdf5/output_hdf5.py index cbef0aba2aa7..fd68c8e17851 100644 --- a/tst/regression/test_suites/output_hdf5/output_hdf5.py +++ b/tst/regression/test_suites/output_hdf5/output_hdf5.py @@ -190,6 +190,18 @@ def Analyse(self, parameters): ) analyze_status = False + # check contents of matadata_none_var + for dim in [2, 3]: + data = phdf.phdf(f"advection_{dim}d.out0.final.phdf") + profile = data.Get("metadata_none_var", flatten=False) + for index in np.ndindex(profile.shape): + ib, j, k, l, m = index + expected_value = l + k + j + m + actual_value = profile[ib, j, k, l, m] + if expected_value != actual_value: + print("metadata_none_var is incorrect") + analyze_status = False + # Checking Parthenon histograms versus numpy ones for dim in [2, 3]: # 1D histogram with binning of a variable with bins defined by a var diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index 3061d2bff961..118a82b58e56 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -54,6 +54,7 @@ vx = 1.0 vy = 1.0 vz = 1.0 profile = hard_sphere +test_metadata_none = true refine_tol = 0.3 # control the package specific refinement tagging function derefine_tol = 0.03 @@ -64,7 +65,8 @@ file_type = hdf5 dt = 1.0 variables = advected, one_minus_advected, & # comments are ok one_minus_advected_sq, & # on every (& characters are ok in comments) - one_minus_sqrt_one_minus_advected_sq # line + one_minus_sqrt_one_minus_advected_sq, & # line + metadata_none_var file_type = hst