From fddb84e8c86cce88bd44e46c4e5262a91e6c2859 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 29 Aug 2024 13:37:11 -0600 Subject: [PATCH 1/2] new examples --- example/CMakeLists.txt | 2 + example/ellipse/inputs/parthinput.ellipse | 59 +++++++++ example/ellipse/src/CMakeLists.txt | 43 +++++++ example/ellipse/src/driver.hpp | 31 +++++ example/ellipse/src/ellipse/ellipse.cpp | 35 ++++++ example/ellipse/src/ellipse/ellipse.hpp | 27 +++++ example/ellipse/src/indicator/indicator.cpp | 33 +++++ example/ellipse/src/indicator/indicator.hpp | 34 ++++++ example/ellipse/src/main.cpp | 65 ++++++++++ example/ellipse/src/pgen.cpp | 60 ++++++++++ example/ellipse/src/pgen.hpp | 21 ++++ example/monte_carlo_pi/inputs/parthinput.toy | 50 ++++++++ example/monte_carlo_pi/src/CMakeLists.txt | 22 ++++ example/monte_carlo_pi/src/main.cpp | 57 +++++++++ example/monte_carlo_pi/src/mccirc/mccirc.cpp | 120 +++++++++++++++++++ example/monte_carlo_pi/src/mccirc/mccirc.hpp | 27 +++++ example/monte_carlo_pi/src/pgen.cpp | 89 ++++++++++++++ example/monte_carlo_pi/src/pgen.hpp | 22 ++++ 18 files changed, 797 insertions(+) create mode 100644 example/ellipse/inputs/parthinput.ellipse create mode 100644 example/ellipse/src/CMakeLists.txt create mode 100644 example/ellipse/src/driver.hpp create mode 100644 example/ellipse/src/ellipse/ellipse.cpp create mode 100644 example/ellipse/src/ellipse/ellipse.hpp create mode 100644 example/ellipse/src/indicator/indicator.cpp create mode 100644 example/ellipse/src/indicator/indicator.hpp create mode 100644 example/ellipse/src/main.cpp create mode 100644 example/ellipse/src/pgen.cpp create mode 100644 example/ellipse/src/pgen.hpp create mode 100644 example/monte_carlo_pi/inputs/parthinput.toy create mode 100644 example/monte_carlo_pi/src/CMakeLists.txt create mode 100644 example/monte_carlo_pi/src/main.cpp create mode 100644 example/monte_carlo_pi/src/mccirc/mccirc.cpp create mode 100644 example/monte_carlo_pi/src/mccirc/mccirc.hpp create mode 100644 example/monte_carlo_pi/src/pgen.cpp create mode 100644 example/monte_carlo_pi/src/pgen.hpp diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index 4df4f0ae1fed..2be4d678a729 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -23,3 +23,5 @@ add_subdirectory(particle_tracers) add_subdirectory(poisson) add_subdirectory(poisson_gmg) add_subdirectory(sparse_advection) +add_subdirectory(ellipse/src ellipse) +add_subdirectory(monte_carlo_pi/src monte_carlo_pi) diff --git a/example/ellipse/inputs/parthinput.ellipse b/example/ellipse/inputs/parthinput.ellipse new file mode 100644 index 000000000000..5669ed7560c8 --- /dev/null +++ b/example/ellipse/inputs/parthinput.ellipse @@ -0,0 +1,59 @@ +# ======================================================================================== +# (C) (or copyright) 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 +# for the U.S. Department of Energy/National Nuclear Security Administration. All rights +# in the program are reserved by Triad National Security, LLC, and the U.S. Department +# of Energy/National Nuclear Security Administration. The Government is granted for +# itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +# license in this material to reproduce, prepare derivative works, distribute copies to +# the public, perform publicly and display publicly, and to permit others to do so. +# ======================================================================================== + + +problem_id = ellipse + + +refinement = adaptive +numlevel = 3 + +nx1 = 64 +x1min = -2.0 +x1max = 2.0 +ix1_bc = outflow +ox1_bc = outflow + +nx2 = 64 +x2min = -2.0 +x2max = 2.0 +ix2_bc = outflow +ox2_bc = outflow + +nx3 = 1 +x3min = -0.5 +x3max = 0.5 + +# How many meshblocks to use in a premade default kernel. +# A value of <1 means use the whole mesh. +pack_size = 1 + + +nx1 = 8 +nx2 = 8 +nx3 = 1 + + +field = phi # the name of the variable we want to refine on +method = derivative_order_1 # selects the first derivative method +refine_tol = 0.5 # tag for refinement if |(dfield/dx)/field| > refine_tol +derefine_tol = 0.05 # tag for derefinement if |(dfield/dx)/field| < derefine_tol +max_level = 3 # if set, limits refinement level from this criteria to no greater than max_level + + +major_axis = 1.5 +minor_axis = 1.0 + + +file_type = hdf5 +variables = phi diff --git a/example/ellipse/src/CMakeLists.txt b/example/ellipse/src/CMakeLists.txt new file mode 100644 index 000000000000..eb5ff175e1be --- /dev/null +++ b/example/ellipse/src/CMakeLists.txt @@ -0,0 +1,43 @@ +#========================================================================================= +# (C) (or copyright) 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 +# for the U.S. Department of Energy/National Nuclear Security Administration. All rights +# in the program are reserved by Triad National Security, LLC, and the U.S. Department +# of Energy/National Nuclear Security Administration. The Government is granted for +# itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +# license in this material to reproduce, prepare derivative works, distribute copies to +# the public, perform publicly and display publicly, and to permit others to do so. +#========================================================================================= + +# There is no test but we build it if we build examples +if (NOT PARTHENON_DISABLE_EXAMPLES) + set (SRC_LIST + main.cpp + + driver.hpp + + ellipse/ellipse.cpp + ellipse/ellipse.hpp + + indicator/indicator.cpp + indicator/indicator.hpp + + pgen.cpp + pgen.hpp + ) + + add_executable(ellipse ${SRC_LIST}) + target_include_directories(ellipse PRIVATE + $ + $ + ) + target_compile_features(ellipse PRIVATE cxx_std_17) + + if (Kokkos_ENABLE_CUDA) + target_compile_options(ellipse --expt-relaxed-constexpr) + endif() + + target_link_libraries(ellipse PRIVATE parthenon) +endif() diff --git a/example/ellipse/src/driver.hpp b/example/ellipse/src/driver.hpp new file mode 100644 index 000000000000..6ab2c5107e49 --- /dev/null +++ b/example/ellipse/src/driver.hpp @@ -0,0 +1,31 @@ +//======================================================================================== +// (C) (or copyright) 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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#ifndef _DRIVER_HPP_ +#define _DRIVER_HPP_ + +#include + +class ToyDriver : parthenon::Driver { +public: + ToyDriver(ParameterInput *pin, ApplicationInput *app_in, Mesh *pm) + : parthenon::Driver(pin, app_in, pm) { + InitializeOutputs(); + } + parthenon::DriverStatus Execute() override { + pouts->MakeOutputs(pmesh, pinput); + return parthenon::DriverStatus::complete; + } +}; + +#endif // _DRIVER_HPP_ diff --git a/example/ellipse/src/ellipse/ellipse.cpp b/example/ellipse/src/ellipse/ellipse.cpp new file mode 100644 index 000000000000..36f3f4781294 --- /dev/null +++ b/example/ellipse/src/ellipse/ellipse.cpp @@ -0,0 +1,35 @@ +//======================================================================================== +// (C) (or copyright) 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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#include + +#include +using namespace parthenon::package::prelude; + +#include "ellipse.hpp" + +namespace Ellipse { + +std::shared_ptr Initialize(ParameterInput *pin) { + auto pkg = std::make_shared("ellipse"); + + // parse input deck and add params for ellipse shape. + const Real major_axis = pin->GetOrAddReal("ellipse","major_axis", 1.0); + const Real minor_axis = pin->GetOrAddReal("ellipse","minor_axis", 1.0); + pkg->AddParam("major_axis", major_axis); + pkg->AddParam("minor_axis", minor_axis); + + return pkg; +} + +} // namespace Ellipse diff --git a/example/ellipse/src/ellipse/ellipse.hpp b/example/ellipse/src/ellipse/ellipse.hpp new file mode 100644 index 000000000000..d81fc3f04c57 --- /dev/null +++ b/example/ellipse/src/ellipse/ellipse.hpp @@ -0,0 +1,27 @@ +//======================================================================================== +// (C) (or copyright) 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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#ifndef _ELLIPSE_ELLIPSE_HPP_ +#define _ELLIPSE_ELLIPSE_HPP_ + +#include +#include + +namespace Ellipse { +using namespace parthenon::package::prelude; + +std::shared_ptr Initialize(ParameterInput *pin); + +} // namespace Ellipse + +#endif // _ELLIPSE_ELLIPSE_HPP_ diff --git a/example/ellipse/src/indicator/indicator.cpp b/example/ellipse/src/indicator/indicator.cpp new file mode 100644 index 000000000000..0fff0ea183a9 --- /dev/null +++ b/example/ellipse/src/indicator/indicator.cpp @@ -0,0 +1,33 @@ +//======================================================================================== +// (C) (or copyright) 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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#include + +#include +using namespace parthenon::package::prelude; + +#include "indicator.hpp" + +namespace Indicator { + +std::shared_ptr Initialize(ParameterInput *pin) { + auto pkg = std::make_shared("indicator"); + + // register the phi variable + Metadata m({Metadata::OneCopy, Metadata::Cell}); + pkg->AddField(m); + + return pkg; +} + +} // namespace Indicator diff --git a/example/ellipse/src/indicator/indicator.hpp b/example/ellipse/src/indicator/indicator.hpp new file mode 100644 index 000000000000..d47fa9b744b8 --- /dev/null +++ b/example/ellipse/src/indicator/indicator.hpp @@ -0,0 +1,34 @@ +//======================================================================================== +// (C) (or copyright) 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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#ifndef _INDICATOR_INDICATOR_HPP_ +#define _INDICATOR_INDICATOR_HPP_ + +#include +#include + +namespace Indicator { +using namespace parthenon::package::prelude; + +struct phi : public parthenon::variable_names::base_t { + template + KOKKOS_INLINE_FUNCTION phi(Ts &&...args) + : parthenon::variable_names::base_t(std::forward(args)...) {} + static std::string name() { return "phi"; } +}; + +std::shared_ptr Initialize(ParameterInput *pin); + +} // namespace Indicator + +#endif // _INDICATOR_INDICATOR_HPP_ diff --git a/example/ellipse/src/main.cpp b/example/ellipse/src/main.cpp new file mode 100644 index 000000000000..08bab06ef1fc --- /dev/null +++ b/example/ellipse/src/main.cpp @@ -0,0 +1,65 @@ +//======================================================================================== +// (C) (or copyright) 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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#include +#include +using namespace parthenon::driver::prelude; + +#include "driver.hpp" +#include "ellipse/ellipse.hpp" +#include "indicator/indicator.hpp" +#include "pgen.hpp" + +int main(int argc, char *argv[]) { + parthenon::ParthenonManager pman; + + // Set up kokkos and read pin + auto manager_status = pman.ParthenonInitEnv(argc, argv); + if (manager_status == ParthenonStatus::complete) { + pman.ParthenonFinalize(); + return 0; + } + if (manager_status == ParthenonStatus::error) { + pman.ParthenonFinalize(); + return 1; + } + + // Redefine parthenon defaults + pman.app_input->ProcessPackages = [](std::unique_ptr &pin) { + Packages_t packages; + packages.Add(Indicator::Initialize(pin.get())); + packages.Add(Ellipse::Initialize(pin.get())); + return packages; + }; + pman.app_input->ProblemGenerator = SetupEllipse; + + // call ParthenonInit to set up the mesh + // scope so that the mesh object, kokkos views, etc, all get cleaned + // up before kokkos::finalize + pman.ParthenonInitPackagesAndMesh(); + { + + // Initialize the driver + ToyDriver driver(pman.pinput.get(), pman.app_input.get(), pman.pmesh.get()); + + // This line actually runs the simulation + auto driver_status = driver.Execute(); // unneeded here + + // call MPI_Finalize and Kokkos::finalize if necessary + } + pman.ParthenonFinalize(); + + // MPI and Kokkos can no longer be used + + return (0); +} diff --git a/example/ellipse/src/pgen.cpp b/example/ellipse/src/pgen.cpp new file mode 100644 index 000000000000..cccf6ab62a0e --- /dev/null +++ b/example/ellipse/src/pgen.cpp @@ -0,0 +1,60 @@ +//======================================================================================== +// (C) (or copyright) 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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#include +#include +using namespace parthenon::package::prelude; + +#include "indicator/indicator.hpp" + +void SetupEllipse(MeshBlock *pmb, ParameterInput *pin) { + + // get meshblock data object + auto &data = pmb->meshblock_data.Get(); + + // pull out ellipse data + auto pkg = pmb->packages.Get("ellipse"); + const auto major_axis = pkg->Param("major_axis"); + const auto minor_axis = pkg->Param("minor_axis"); + const Real a2 = major_axis*major_axis; + const Real b2 = minor_axis*minor_axis; + + // loop bounds for interior of meshblock + auto cellbounds = pmb->cellbounds; + IndexRange ib = cellbounds.GetBoundsI(IndexDomain::interior); + IndexRange jb = cellbounds.GetBoundsJ(IndexDomain::interior); + IndexRange kb = cellbounds.GetBoundsK(IndexDomain::interior); + + // coordinates object + auto coords = pmb->coords; + + // build a type-based variable pack + static auto desc = parthenon::MakePackDescriptor(data.get()); + auto pack = desc.GetPack(data.get()); + + const int ndim = pmb->pmy_mesh->ndim; + PARTHENON_REQUIRE_THROWS(ndim >= 2, "Calculate area must be at least 2d"); + + const int b = 0; + 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) { + // are we on the ellipse? + Real x = coords.Xc(k, j, i); + Real y = coords.Xc(k, j, i); + Real condition = ((x*x)/(a2 + 1e-20) + (y*y)/(b2 + 1e-20)) <= 1; + // set indicator function appropriately + pack(b, Indicator::phi(), k, j, i) = condition; + }); + return; +} diff --git a/example/ellipse/src/pgen.hpp b/example/ellipse/src/pgen.hpp new file mode 100644 index 000000000000..8d3042f61523 --- /dev/null +++ b/example/ellipse/src/pgen.hpp @@ -0,0 +1,21 @@ +//======================================================================================== +// (C) (or copyright) 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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#ifndef _PGEN_HPP_ +#define _PGEN_HPP_ + +#include + +void SetupEllipse(parthenon::MeshBlock *pmb, parthenon::ParameterInput *pin); + +#endif // _PGEN_HPP_ diff --git a/example/monte_carlo_pi/inputs/parthinput.toy b/example/monte_carlo_pi/inputs/parthinput.toy new file mode 100644 index 000000000000..be5ed0333fd5 --- /dev/null +++ b/example/monte_carlo_pi/inputs/parthinput.toy @@ -0,0 +1,50 @@ +# ======================================================================================== +# (C) (or copyright) 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 +# for the U.S. Department of Energy/National Nuclear Security Administration. All rights +# in the program are reserved by Triad National Security, LLC, and the U.S. Department +# of Energy/National Nuclear Security Administration. The Government is granted for +# itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +# license in this material to reproduce, prepare derivative works, distribute copies to +# the public, perform publicly and display publicly, and to permit others to do so. +# ======================================================================================== + + +problem_id = montecarlo + + +nx1 = 64 +x1min = -1.0 +x1max = 1.0 +ix1_bc = periodic +ox1_bc = periodic + +nx2 = 64 +x2min = -1.0 +x2max = 1.0 +ix2_bc = periodic +ox2_bc = periodic + +nx3 = 1 +x3min = -0.5 +x3max = 0.5 +ix3_bc = periodic +ox3_bc = periodic + +# How many meshblocks to use in a premade default kernel. +# A value of <1 means use the whole mesh. +pack_size = -1 + + +nx1 = 8 +nx2 = 8 +nx3 = 1 + + +radius = 1 + + +num_particles_per_block = 1000 +rng_seed = 1234 \ No newline at end of file diff --git a/example/monte_carlo_pi/src/CMakeLists.txt b/example/monte_carlo_pi/src/CMakeLists.txt new file mode 100644 index 000000000000..f2056e8c0ab9 --- /dev/null +++ b/example/monte_carlo_pi/src/CMakeLists.txt @@ -0,0 +1,22 @@ +if (NOT PARTHENON_DISABLE_EXAMPLES) + set (SRC_LIST + main.cpp + mccirc/mccirc.hpp + mccirc/mccirc.cpp + pgen.cpp + pgen.hpp + ) + + add_executable(mcpi ${SRC_LIST}) + target_include_directories(mcpi PRIVATE + $ + $ + ) + target_compile_features(mcpi PRIVATE cxx_std_17) + + if (Kokkos_ENABLE_CUDA) + target_compile_options(mcpi --expt-relaxed-constexpr) + endif() + + target_link_libraries(mcpi PRIVATE parthenon) +endif() diff --git a/example/monte_carlo_pi/src/main.cpp b/example/monte_carlo_pi/src/main.cpp new file mode 100644 index 000000000000..7458333549e0 --- /dev/null +++ b/example/monte_carlo_pi/src/main.cpp @@ -0,0 +1,57 @@ +//======================================================================================== +// (C) (or copyright) 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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#include +#include +using namespace parthenon::driver::prelude; + +#include "mccirc/mccirc.hpp" +#include "pgen.hpp" + +int main(int argc, char *argv[]) { + parthenon::ParthenonManager pman; + + // Set up kokkos and read pin + auto manager_status = pman.ParthenonInitEnv(argc, argv); + if (manager_status == ParthenonStatus::complete) { + pman.ParthenonFinalize(); + return 0; + } + if (manager_status == ParthenonStatus::error) { + pman.ParthenonFinalize(); + return 1; + } + + // Redefine parthenon defaults + pman.app_input->ProcessPackages = [](std::unique_ptr &pin) { + Packages_t packages; + packages.Add(MCCirc::Initialize(pin.get())); + return packages; + }; + pman.app_input->ProblemGenerator = GenerateCircle; + + // call ParthenonInit to set up the mesh + // scope so that the mesh object, kokkos views, etc, all get cleaned + // up before kokkos::finalize + pman.ParthenonInitPackagesAndMesh(); + { + + MCCirc::ComputeParticleCounts(pman.pmesh.get()); + // call MPI_Finalize and Kokkos::finalize if necessary + } + pman.ParthenonFinalize(); + + // MPI and Kokkos can no longer be used + + return (0); +} diff --git a/example/monte_carlo_pi/src/mccirc/mccirc.cpp b/example/monte_carlo_pi/src/mccirc/mccirc.cpp new file mode 100644 index 000000000000..d0200bbafa2c --- /dev/null +++ b/example/monte_carlo_pi/src/mccirc/mccirc.cpp @@ -0,0 +1,120 @@ +#include +#include +using namespace parthenon::package::prelude; + +#include "mccirc.hpp" + +namespace MCCirc { + +std::shared_ptr Initialize(ParameterInput *pin) { + auto pkg = std::make_shared("MCCirc"); + + const Real radius = pin->GetOrAddReal("circle", "radius", 1.0); + pkg->AddParam("radius", radius); + + const int npart = pin->GetOrAddInteger("MonteCarlo", "num_particles_per_block", 1000); + pkg->AddParam("num_particles", npart); + + // Initialize random number generator pool + int rng_seed = pin->GetOrAddInteger("MonteCarlo", "rng_seed", 1234); + pkg->AddParam("rng_seed", rng_seed); + RNGPool rng_pool(rng_seed); + pkg->AddParam("rng_pool", rng_pool); + + // TO access later + // pkg->Param("name"); + // e.g., pkg->Param("rng_pool"); + + Metadata swarm_metadata({Metadata::Provides, Metadata::None}); + pkg->AddSwarm("samples", swarm_metadata); + + Metadata real_swarmvalue_metadata({Metadata::Real}); + pkg->AddSwarmValue(weight::name(), "samples", real_swarmvalue_metadata); + + Metadata m({Metadata::Cell}); + pkg->AddField(m); + + return pkg; +} + +// in task list should return task status but I'm being lazy. +void ComputeParticleCounts(Mesh *pm) { + // get mesh data + auto md = pm->mesh_data.Get(); + + // Make a SwarmPack via types to get positions + static auto desc_swarm = + parthenon::MakeSwarmPackDescriptor("samples"); + auto pack_swarm = desc_swarm.GetPack(md.get()); + + // pull out circle radius from params + auto pkg = pm->packages.Get("MCCirc"); + auto r = pkg->Param("radius"); + + // build a type-based variable pack + static auto desc = parthenon::MakePackDescriptor(md.get()); + auto pack = desc.GetPack(md.get()); + + IndexRange ib = md->GetBoundsI(IndexDomain::interior); + IndexRange jb = md->GetBoundsJ(IndexDomain::interior); + IndexRange kb = md->GetBoundsK(IndexDomain::interior); + 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, MCCirc::NumParticles(), k, j, i) = 0; + }); + + parthenon::par_for(DEFAULT_LOOP_PATTERN, PARTHENON_AUTO_LABEL, + DevExecSpace(), 0, + pack_swarm.GetMaxFlatIndex(), + // loop over all particles + KOKKOS_LAMBDA(const int idx) { + // block and particle indices + auto [b, n] = pack_swarm.GetBlockParticleIndices(idx); + const auto swarm_d = pack_swarm.GetContext(b); + if (swarm_d.IsActive(n)) { + // computes block-local cell indices of this particle + int i, j, k; + Real x = pack_swarm(b, swarm_position::x(), n); + Real y = pack_swarm(b, swarm_position::y(), n); + Real z = pack_swarm(b, swarm_position::z(), n); + swarm_d.Xtoijk(x, y, z, i, j, k); + + Kokkos::atomic_add(&pack(b, MCCirc::NumParticles(), k, j, i), + pack_swarm(b, MCCirc::weight(), n)); + } + }); + + // local reductions over all meshblocks in meshdata object + Real total_particles; + parthenon::par_reduce(parthenon::LoopPatternMDRange(), + PARTHENON_AUTO_LABEL, DevExecSpace(), + 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, Real &tot) { + tot += pack(b, MCCirc::NumParticles(), k, j, i); + }, total_particles); + Real total_particles_in_circle; + parthenon::par_reduce(parthenon::LoopPatternMDRange(), + PARTHENON_AUTO_LABEL, DevExecSpace(), + 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, Real &tot) { + auto &coords = pack.GetCoordinates(b); + Real x = coords.Xc(k, j, i); + Real y = coords.Xc(k, j, i); + bool in_circle = x*x + y*y < r*r; + tot += in_circle*pack(b, MCCirc::NumParticles(), k, j, i); + }, total_particles_in_circle); + + // just print for simplicity but if we were doing this right, we would call parthenon's reductions + // which take the above data and reduce accross MPI ranks and task lists + printf("particles in circle, particles total, pi = %.14e %.14e %.14e\n", + total_particles_in_circle, total_particles, 4.*total_particles_in_circle/total_particles); + +} + +} // namespace MCCirc diff --git a/example/monte_carlo_pi/src/mccirc/mccirc.hpp b/example/monte_carlo_pi/src/mccirc/mccirc.hpp new file mode 100644 index 000000000000..a73964589f7c --- /dev/null +++ b/example/monte_carlo_pi/src/mccirc/mccirc.hpp @@ -0,0 +1,27 @@ +#ifndef _MCCIRC_MCCIRC_HPP_ +#define _MCCIRC_MCCIRC_HPP_ + +#include +#include "Kokkos_Random.hpp" +#include +#include + +namespace MCCirc { +using namespace parthenon::package::prelude; + +typedef Kokkos::Random_XorShift64_Pool<> RNGPool; + +struct NumParticles : public parthenon::variable_names::base_t { + template + KOKKOS_INLINE_FUNCTION NumParticles(Ts &&...args) + : parthenon::variable_names::base_t(std::forward(args)...) {} + static std::string name() { return "particles_per_cell"; } +}; +SWARM_VARIABLE(Real, mc, weight); + +std::shared_ptr Initialize(ParameterInput *pin); +void ComputeParticleCounts(Mesh *pm); + +} // namespace MCCirc + +#endif // _MCCIRC_MCCIRC_HPP_ diff --git a/example/monte_carlo_pi/src/pgen.cpp b/example/monte_carlo_pi/src/pgen.cpp new file mode 100644 index 000000000000..50bf2f1f2aab --- /dev/null +++ b/example/monte_carlo_pi/src/pgen.cpp @@ -0,0 +1,89 @@ +//======================================================================================== +// (C) (or copyright) 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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#include +#include +using namespace parthenon::package::prelude; + +#include "mccirc/mccirc.hpp" + +void GenerateCircle(parthenon::MeshBlock *pmb, parthenon::ParameterInput *pin) { + auto &data = pmb->meshblock_data.Get(); + + // pull out information/global params from package + auto pkg = pmb->packages.Get("MCCirc"); + auto rng_pool = pkg->Param("rng_pool"); + const int N = pkg->Param("num_particles"); + + // Pull out swarm object + auto swarm = data->GetSwarmData()->Get("samples"); + + // Create an accessor to particles, allocate particles + auto newParticlesContext = swarm->AddEmptyParticles(N); + + // Meshblock geometry + const IndexRange &ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); + const IndexRange &jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); + const IndexRange &kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); + const int &nx_i = pmb->cellbounds.ncellsi(IndexDomain::interior); + const int &nx_j = pmb->cellbounds.ncellsj(IndexDomain::interior); + const int &nx_k = pmb->cellbounds.ncellsk(IndexDomain::interior); + const Real &dx_i = pmb->coords.Dxf<1>(pmb->cellbounds.is(IndexDomain::interior)); + const Real &dx_j = pmb->coords.Dxf<2>(pmb->cellbounds.js(IndexDomain::interior)); + const Real &dx_k = pmb->coords.Dxf<3>(pmb->cellbounds.ks(IndexDomain::interior)); + const Real &minx_i = pmb->coords.Xf<1>(ib.s); + const Real &minx_j = pmb->coords.Xf<2>(jb.s); + const Real &minx_k = pmb->coords.Xf<3>(kb.s); + + auto &x = swarm->Get(swarm_position::x::name()).Get(); + auto &y = swarm->Get(swarm_position::y::name()).Get(); + auto &z = swarm->Get(swarm_position::z::name()).Get(); + auto &weight = swarm->Get(MCCirc::weight::name()).Get(); + + // Make a SwarmPack via types to get positions + static auto desc_swarm = + parthenon::MakeSwarmPackDescriptor("samples"); + auto pack_swarm = desc_swarm.GetPack(data.get()); + auto swarm_d = swarm->GetDeviceContext(); + + // loop over new particles created + parthenon::par_for(DEFAULT_LOOP_PATTERN, PARTHENON_AUTO_LABEL, + DevExecSpace(), 0, + newParticlesContext.GetNewParticlesMaxIndex(), + // new_n ranges from 0 to N_new_particles + KOKKOS_LAMBDA(const int new_n) { + // this is the particle index inside the swarm + const int n = newParticlesContext.GetNewParticleIndex(new_n); + auto rng_gen = rng_pool.get_state(); + + // Normally b would be free-floating and set by pack.GetBlockparticleIndices + // but since we're on a single meshblock for this loop, it's just 0 + // because block index = 0 + const int b = 0; + //auto [b, n] = pack_swarm.GetBlockparticleIndices(idx); + + // randomly sample particle positions + pack_swarm(b, swarm_position::x(), n) = minx_i + nx_i * dx_i * rng_gen.drand(); + pack_swarm(b, swarm_position::y(), n) = minx_j + nx_j * dx_j * rng_gen.drand(); + pack_swarm(b, swarm_position::z(), n) = minx_k + nx_k * dx_k * rng_gen.drand(); + + // set weights to 1 + pack_swarm(b, MCCirc::weight(), n) = 1.0; + + // release random number generator + rng_pool.free_state(rng_gen); + }); +} diff --git a/example/monte_carlo_pi/src/pgen.hpp b/example/monte_carlo_pi/src/pgen.hpp new file mode 100644 index 000000000000..c221648ba67f --- /dev/null +++ b/example/monte_carlo_pi/src/pgen.hpp @@ -0,0 +1,22 @@ +//======================================================================================== +// (C) (or copyright) 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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#ifndef _PGEN_HPP_ +#define _PGEN_HPP_ + +#include + +void GenerateCircle(parthenon::MeshBlock *pmb, parthenon::ParameterInput *pin); + + +#endif // _PGEN_HPP_ From 12e47123329821419306cf7435f7bca216e95fe5 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 29 Aug 2024 13:37:44 -0600 Subject: [PATCH 2/2] formatting --- example/ellipse/src/driver.hpp | 4 +- example/ellipse/src/ellipse/ellipse.cpp | 4 +- example/ellipse/src/indicator/indicator.hpp | 2 +- example/ellipse/src/main.cpp | 14 +-- example/ellipse/src/pgen.cpp | 6 +- example/monte_carlo_pi/src/main.cpp | 6 +- example/monte_carlo_pi/src/mccirc/mccirc.cpp | 102 +++++++++---------- example/monte_carlo_pi/src/mccirc/mccirc.hpp | 8 +- example/monte_carlo_pi/src/pgen.cpp | 50 +++++---- example/monte_carlo_pi/src/pgen.hpp | 1 - 10 files changed, 96 insertions(+), 101 deletions(-) diff --git a/example/ellipse/src/driver.hpp b/example/ellipse/src/driver.hpp index 6ab2c5107e49..9347faa18df2 100644 --- a/example/ellipse/src/driver.hpp +++ b/example/ellipse/src/driver.hpp @@ -17,9 +17,9 @@ #include class ToyDriver : parthenon::Driver { -public: + public: ToyDriver(ParameterInput *pin, ApplicationInput *app_in, Mesh *pm) - : parthenon::Driver(pin, app_in, pm) { + : parthenon::Driver(pin, app_in, pm) { InitializeOutputs(); } parthenon::DriverStatus Execute() override { diff --git a/example/ellipse/src/ellipse/ellipse.cpp b/example/ellipse/src/ellipse/ellipse.cpp index 36f3f4781294..292b86a10be6 100644 --- a/example/ellipse/src/ellipse/ellipse.cpp +++ b/example/ellipse/src/ellipse/ellipse.cpp @@ -24,8 +24,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto pkg = std::make_shared("ellipse"); // parse input deck and add params for ellipse shape. - const Real major_axis = pin->GetOrAddReal("ellipse","major_axis", 1.0); - const Real minor_axis = pin->GetOrAddReal("ellipse","minor_axis", 1.0); + const Real major_axis = pin->GetOrAddReal("ellipse", "major_axis", 1.0); + const Real minor_axis = pin->GetOrAddReal("ellipse", "minor_axis", 1.0); pkg->AddParam("major_axis", major_axis); pkg->AddParam("minor_axis", minor_axis); diff --git a/example/ellipse/src/indicator/indicator.hpp b/example/ellipse/src/indicator/indicator.hpp index d47fa9b744b8..aaec67cfa8e3 100644 --- a/example/ellipse/src/indicator/indicator.hpp +++ b/example/ellipse/src/indicator/indicator.hpp @@ -23,7 +23,7 @@ using namespace parthenon::package::prelude; struct phi : public parthenon::variable_names::base_t { template KOKKOS_INLINE_FUNCTION phi(Ts &&...args) - : parthenon::variable_names::base_t(std::forward(args)...) {} + : parthenon::variable_names::base_t(std::forward(args)...) {} static std::string name() { return "phi"; } }; diff --git a/example/ellipse/src/main.cpp b/example/ellipse/src/main.cpp index 08bab06ef1fc..c04727b9d53c 100644 --- a/example/ellipse/src/main.cpp +++ b/example/ellipse/src/main.cpp @@ -38,7 +38,7 @@ int main(int argc, char *argv[]) { pman.app_input->ProcessPackages = [](std::unique_ptr &pin) { Packages_t packages; packages.Add(Indicator::Initialize(pin.get())); - packages.Add(Ellipse::Initialize(pin.get())); + packages.Add(Ellipse::Initialize(pin.get())); return packages; }; pman.app_input->ProblemGenerator = SetupEllipse; @@ -47,15 +47,15 @@ int main(int argc, char *argv[]) { // scope so that the mesh object, kokkos views, etc, all get cleaned // up before kokkos::finalize pman.ParthenonInitPackagesAndMesh(); - { + { - // Initialize the driver - ToyDriver driver(pman.pinput.get(), pman.app_input.get(), pman.pmesh.get()); + // Initialize the driver + ToyDriver driver(pman.pinput.get(), pman.app_input.get(), pman.pmesh.get()); - // This line actually runs the simulation - auto driver_status = driver.Execute(); // unneeded here + // This line actually runs the simulation + auto driver_status = driver.Execute(); // unneeded here - // call MPI_Finalize and Kokkos::finalize if necessary + // call MPI_Finalize and Kokkos::finalize if necessary } pman.ParthenonFinalize(); diff --git a/example/ellipse/src/pgen.cpp b/example/ellipse/src/pgen.cpp index cccf6ab62a0e..d63aad8cd4eb 100644 --- a/example/ellipse/src/pgen.cpp +++ b/example/ellipse/src/pgen.cpp @@ -26,8 +26,8 @@ void SetupEllipse(MeshBlock *pmb, ParameterInput *pin) { auto pkg = pmb->packages.Get("ellipse"); const auto major_axis = pkg->Param("major_axis"); const auto minor_axis = pkg->Param("minor_axis"); - const Real a2 = major_axis*major_axis; - const Real b2 = minor_axis*minor_axis; + const Real a2 = major_axis * major_axis; + const Real b2 = minor_axis * minor_axis; // loop bounds for interior of meshblock auto cellbounds = pmb->cellbounds; @@ -52,7 +52,7 @@ void SetupEllipse(MeshBlock *pmb, ParameterInput *pin) { // are we on the ellipse? Real x = coords.Xc(k, j, i); Real y = coords.Xc(k, j, i); - Real condition = ((x*x)/(a2 + 1e-20) + (y*y)/(b2 + 1e-20)) <= 1; + Real condition = ((x * x) / (a2 + 1e-20) + (y * y) / (b2 + 1e-20)) <= 1; // set indicator function appropriately pack(b, Indicator::phi(), k, j, i) = condition; }); diff --git a/example/monte_carlo_pi/src/main.cpp b/example/monte_carlo_pi/src/main.cpp index 7458333549e0..6483264102d5 100644 --- a/example/monte_carlo_pi/src/main.cpp +++ b/example/monte_carlo_pi/src/main.cpp @@ -32,7 +32,7 @@ int main(int argc, char *argv[]) { return 1; } - // Redefine parthenon defaults + // Redefine parthenon defaults pman.app_input->ProcessPackages = [](std::unique_ptr &pin) { Packages_t packages; packages.Add(MCCirc::Initialize(pin.get())); @@ -44,8 +44,8 @@ int main(int argc, char *argv[]) { // scope so that the mesh object, kokkos views, etc, all get cleaned // up before kokkos::finalize pman.ParthenonInitPackagesAndMesh(); - { - + { + MCCirc::ComputeParticleCounts(pman.pmesh.get()); // call MPI_Finalize and Kokkos::finalize if necessary } diff --git a/example/monte_carlo_pi/src/mccirc/mccirc.cpp b/example/monte_carlo_pi/src/mccirc/mccirc.cpp index d0200bbafa2c..76a9b3186aaa 100644 --- a/example/monte_carlo_pi/src/mccirc/mccirc.cpp +++ b/example/monte_carlo_pi/src/mccirc/mccirc.cpp @@ -44,10 +44,8 @@ void ComputeParticleCounts(Mesh *pm) { // Make a SwarmPack via types to get positions static auto desc_swarm = - parthenon::MakeSwarmPackDescriptor("samples"); + parthenon::MakeSwarmPackDescriptor("samples"); auto pack_swarm = desc_swarm.GetPack(md.get()); // pull out circle radius from params @@ -62,59 +60,59 @@ void ComputeParticleCounts(Mesh *pm) { IndexRange jb = md->GetBoundsJ(IndexDomain::interior); IndexRange kb = md->GetBoundsK(IndexDomain::interior); 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, MCCirc::NumParticles(), k, j, i) = 0; - }); - - parthenon::par_for(DEFAULT_LOOP_PATTERN, PARTHENON_AUTO_LABEL, - DevExecSpace(), 0, - pack_swarm.GetMaxFlatIndex(), - // loop over all particles - KOKKOS_LAMBDA(const int idx) { - // block and particle indices - auto [b, n] = pack_swarm.GetBlockParticleIndices(idx); - const auto swarm_d = pack_swarm.GetContext(b); - if (swarm_d.IsActive(n)) { - // computes block-local cell indices of this particle - int i, j, k; - Real x = pack_swarm(b, swarm_position::x(), n); - Real y = pack_swarm(b, swarm_position::y(), n); - Real z = pack_swarm(b, swarm_position::z(), n); - swarm_d.Xtoijk(x, y, z, i, j, k); - - Kokkos::atomic_add(&pack(b, MCCirc::NumParticles(), k, j, i), - pack_swarm(b, MCCirc::weight(), n)); - } - }); + 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, MCCirc::NumParticles(), k, j, i) = 0; + }); + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, PARTHENON_AUTO_LABEL, DevExecSpace(), 0, + pack_swarm.GetMaxFlatIndex(), + // loop over all particles + KOKKOS_LAMBDA(const int idx) { + // block and particle indices + auto [b, n] = pack_swarm.GetBlockParticleIndices(idx); + const auto swarm_d = pack_swarm.GetContext(b); + if (swarm_d.IsActive(n)) { + // computes block-local cell indices of this particle + int i, j, k; + Real x = pack_swarm(b, swarm_position::x(), n); + Real y = pack_swarm(b, swarm_position::y(), n); + Real z = pack_swarm(b, swarm_position::z(), n); + swarm_d.Xtoijk(x, y, z, i, j, k); + + Kokkos::atomic_add(&pack(b, MCCirc::NumParticles(), k, j, i), + pack_swarm(b, MCCirc::weight(), n)); + } + }); // local reductions over all meshblocks in meshdata object Real total_particles; - parthenon::par_reduce(parthenon::LoopPatternMDRange(), - PARTHENON_AUTO_LABEL, DevExecSpace(), - 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, Real &tot) { - tot += pack(b, MCCirc::NumParticles(), k, j, i); - }, total_particles); + parthenon::par_reduce( + parthenon::LoopPatternMDRange(), PARTHENON_AUTO_LABEL, DevExecSpace(), 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, Real &tot) { + tot += pack(b, MCCirc::NumParticles(), k, j, i); + }, + total_particles); Real total_particles_in_circle; - parthenon::par_reduce(parthenon::LoopPatternMDRange(), - PARTHENON_AUTO_LABEL, DevExecSpace(), - 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, Real &tot) { - auto &coords = pack.GetCoordinates(b); - Real x = coords.Xc(k, j, i); - Real y = coords.Xc(k, j, i); - bool in_circle = x*x + y*y < r*r; - tot += in_circle*pack(b, MCCirc::NumParticles(), k, j, i); - }, total_particles_in_circle); - - // just print for simplicity but if we were doing this right, we would call parthenon's reductions - // which take the above data and reduce accross MPI ranks and task lists + parthenon::par_reduce( + parthenon::LoopPatternMDRange(), PARTHENON_AUTO_LABEL, DevExecSpace(), 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, Real &tot) { + auto &coords = pack.GetCoordinates(b); + Real x = coords.Xc(k, j, i); + Real y = coords.Xc(k, j, i); + bool in_circle = x * x + y * y < r * r; + tot += in_circle * pack(b, MCCirc::NumParticles(), k, j, i); + }, + total_particles_in_circle); + + // just print for simplicity but if we were doing this right, we would call parthenon's + // reductions which take the above data and reduce accross MPI ranks and task lists printf("particles in circle, particles total, pi = %.14e %.14e %.14e\n", - total_particles_in_circle, total_particles, 4.*total_particles_in_circle/total_particles); - + total_particles_in_circle, total_particles, + 4. * total_particles_in_circle / total_particles); } } // namespace MCCirc diff --git a/example/monte_carlo_pi/src/mccirc/mccirc.hpp b/example/monte_carlo_pi/src/mccirc/mccirc.hpp index a73964589f7c..a7d8c1e1bb59 100644 --- a/example/monte_carlo_pi/src/mccirc/mccirc.hpp +++ b/example/monte_carlo_pi/src/mccirc/mccirc.hpp @@ -1,10 +1,10 @@ #ifndef _MCCIRC_MCCIRC_HPP_ #define _MCCIRC_MCCIRC_HPP_ -#include #include "Kokkos_Random.hpp" -#include #include +#include +#include namespace MCCirc { using namespace parthenon::package::prelude; @@ -12,9 +12,9 @@ using namespace parthenon::package::prelude; typedef Kokkos::Random_XorShift64_Pool<> RNGPool; struct NumParticles : public parthenon::variable_names::base_t { - template + template KOKKOS_INLINE_FUNCTION NumParticles(Ts &&...args) - : parthenon::variable_names::base_t(std::forward(args)...) {} + : parthenon::variable_names::base_t(std::forward(args)...) {} static std::string name() { return "particles_per_cell"; } }; SWARM_VARIABLE(Real, mc, weight); diff --git a/example/monte_carlo_pi/src/pgen.cpp b/example/monte_carlo_pi/src/pgen.cpp index 50bf2f1f2aab..bb5c7f10bb12 100644 --- a/example/monte_carlo_pi/src/pgen.cpp +++ b/example/monte_carlo_pi/src/pgen.cpp @@ -52,38 +52,36 @@ void GenerateCircle(parthenon::MeshBlock *pmb, parthenon::ParameterInput *pin) { // Make a SwarmPack via types to get positions static auto desc_swarm = - parthenon::MakeSwarmPackDescriptor("samples"); + parthenon::MakeSwarmPackDescriptor("samples"); auto pack_swarm = desc_swarm.GetPack(data.get()); auto swarm_d = swarm->GetDeviceContext(); // loop over new particles created - parthenon::par_for(DEFAULT_LOOP_PATTERN, PARTHENON_AUTO_LABEL, - DevExecSpace(), 0, - newParticlesContext.GetNewParticlesMaxIndex(), - // new_n ranges from 0 to N_new_particles - KOKKOS_LAMBDA(const int new_n) { - // this is the particle index inside the swarm - const int n = newParticlesContext.GetNewParticleIndex(new_n); - auto rng_gen = rng_pool.get_state(); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, PARTHENON_AUTO_LABEL, DevExecSpace(), 0, + newParticlesContext.GetNewParticlesMaxIndex(), + // new_n ranges from 0 to N_new_particles + KOKKOS_LAMBDA(const int new_n) { + // this is the particle index inside the swarm + const int n = newParticlesContext.GetNewParticleIndex(new_n); + auto rng_gen = rng_pool.get_state(); - // Normally b would be free-floating and set by pack.GetBlockparticleIndices - // but since we're on a single meshblock for this loop, it's just 0 - // because block index = 0 - const int b = 0; - //auto [b, n] = pack_swarm.GetBlockparticleIndices(idx); + // Normally b would be free-floating and set by pack.GetBlockparticleIndices + // but since we're on a single meshblock for this loop, it's just 0 + // because block index = 0 + const int b = 0; + // auto [b, n] = pack_swarm.GetBlockparticleIndices(idx); - // randomly sample particle positions - pack_swarm(b, swarm_position::x(), n) = minx_i + nx_i * dx_i * rng_gen.drand(); - pack_swarm(b, swarm_position::y(), n) = minx_j + nx_j * dx_j * rng_gen.drand(); - pack_swarm(b, swarm_position::z(), n) = minx_k + nx_k * dx_k * rng_gen.drand(); + // randomly sample particle positions + pack_swarm(b, swarm_position::x(), n) = minx_i + nx_i * dx_i * rng_gen.drand(); + pack_swarm(b, swarm_position::y(), n) = minx_j + nx_j * dx_j * rng_gen.drand(); + pack_swarm(b, swarm_position::z(), n) = minx_k + nx_k * dx_k * rng_gen.drand(); - // set weights to 1 - pack_swarm(b, MCCirc::weight(), n) = 1.0; + // set weights to 1 + pack_swarm(b, MCCirc::weight(), n) = 1.0; - // release random number generator - rng_pool.free_state(rng_gen); - }); + // release random number generator + rng_pool.free_state(rng_gen); + }); } diff --git a/example/monte_carlo_pi/src/pgen.hpp b/example/monte_carlo_pi/src/pgen.hpp index c221648ba67f..9901c0ba0fdc 100644 --- a/example/monte_carlo_pi/src/pgen.hpp +++ b/example/monte_carlo_pi/src/pgen.hpp @@ -18,5 +18,4 @@ void GenerateCircle(parthenon::MeshBlock *pmb, parthenon::ParameterInput *pin); - #endif // _PGEN_HPP_