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..9347faa18df2 --- /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..292b86a10be6 --- /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..aaec67cfa8e3 --- /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..c04727b9d53c --- /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..d63aad8cd4eb --- /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..6483264102d5 --- /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..76a9b3186aaa --- /dev/null +++ b/example/monte_carlo_pi/src/mccirc/mccirc.cpp @@ -0,0 +1,118 @@ +#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..a7d8c1e1bb59 --- /dev/null +++ b/example/monte_carlo_pi/src/mccirc/mccirc.hpp @@ -0,0 +1,27 @@ +#ifndef _MCCIRC_MCCIRC_HPP_ +#define _MCCIRC_MCCIRC_HPP_ + +#include "Kokkos_Random.hpp" +#include +#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..bb5c7f10bb12 --- /dev/null +++ b/example/monte_carlo_pi/src/pgen.cpp @@ -0,0 +1,87 @@ +//======================================================================================== +// (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..9901c0ba0fdc --- /dev/null +++ b/example/monte_carlo_pi/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 GenerateCircle(parthenon::MeshBlock *pmb, parthenon::ParameterInput *pin); + +#endif // _PGEN_HPP_