From 7670dd6518ed3522002d2b943cca0136c5000162 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Tue, 17 Mar 2026 23:57:41 -0700 Subject: [PATCH 01/55] templated coordinates type --- src/grid/geometry.hpp | 119 ++++++++++++++++++++++++++++++++++++++++++ src/grid/subpack.hpp | 4 ++ 2 files changed, 123 insertions(+) create mode 100644 src/grid/geometry.hpp diff --git a/src/grid/geometry.hpp b/src/grid/geometry.hpp new file mode 100644 index 0000000..e00f9cb --- /dev/null +++ b/src/grid/geometry.hpp @@ -0,0 +1,119 @@ +#ifndef GRID_GEOMETRY_HPP_ +#define GRID_GEOMETRY_HPP_ + +#include +#include + +#include "Kokkos_Macros.hpp" + +#include "coordinates/coordinates.hpp" + +#include "dispatcher/options.hpp" +#include "grid/grid_types.hpp" +#include "grid/subpack.hpp" + +namespace kamayan { +POLYMORPHIC_PARM(Geometry, cartesian, cylindrical); + +namespace grid { + +// This is a (incomplete) wrapper around parthenon's uniform_coordinate (cartesian) +// to calculate coordinate related items in a different geometry. +template +struct Coordinates { + // distance between cell centers + template + KOKKOS_FORCEINLINE_FUNCTION Real Dx() const { + if constexpr (geom == Geometry::cylindrical && ax == Axis::KAXIS) { + // I guess we're assuming this only is a useful geometry in NDIM < 3 + return 2.0 * std::numbers::pi; + } + return coords_.Dx(); + } + + KOKKOS_FORCEINLINE_FUNCTION Real Dx(const Axis &ax) const { + return AxisOverload([this]() { return Dx(); }, ax); + } + + // position at cell centroids + template + KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int &idx) const { + constexpr auto dir = AxisToInt(ax); + if constexpr (ax == Axis::IAXIS && geom == Geometry::cylindrical) { + const Real r0 = coords_.Xf(idx); + const Real r0sq = r0 * r0; + const Real r1 = coords_.Xf(idx + 1); + const Real r1sq = r1 * r1; + return (2.0 / 3.0) * (r1sq * r1 - r0sq * r0) / (r1sq - r0sq); + } + return coords_.Xc(idx); + } + + KOKKOS_FORCEINLINE_FUNCTION Real Xc(const Axis &ax, const int &idx) const { + return AxisOverload([&, this]() { return Xc(idx); }, ax); + } + + // position at face centers + template + KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int &idx) const { + return coords_.Xf(idx); + } + + KOKKOS_FORCEINLINE_FUNCTION Real Xf(const Axis &ax, const int &idx) const { + return AxisOverload([&, this]() { return Xf(idx); }, ax); + } + + template + KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int &k, const int &j, + const int &i) const { + constexpr auto dir = AxisToInt(ax); + if constexpr (geom == Geometry::cylindrical) { + if constexpr (ax == Axis::IAXIS) { + // 2 * pi * r_i * dz + return Xf(i) * Dx() * Dx(); + } else if constexpr (ax == Axis::JAXIS) { + // pi * (r_i+1/2**2 - r_i-1/2**2) + const Real rp = Xf(i + 1); + const Real rm = Xf(i); + return 0.5 * (rp * rp - rm * rm) * Dx(); + } + } + return coords_.FaceArea(k, j, i); + } + + KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const Axis &ax, const int &k, const int &j, + const int &i) const { + return AxisOverload([&, this]() { return FaceArea(k, j, i); }, ax); + } + + KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int &k, const int &j, const int &i) { + if constexpr (geom == Geometry::cylindrical) { + // dphi / 2 * (r_p^2 - r_m^2) * dz + const Real rp = Xf(i + 1); + const Real rm = Xf(i); + return 0.5 * Dx() * (rp * rp - rm * rm); + } + return coords_.CellVolume(k, j, i); + } + + private: + const parthenon::Coordinates_t coords_; + + template + KOKKOS_FORCEINLINE_FUNCTION Real AxisOverload(Func function, Axis ax, Args &&...args) { + if (ax == Axis::IAXIS) { + return function.template operator()(std::forward(args)...); + } else if (ax == Axis::JAXIS) { + return function.template operator()(std::forward(args)...); + } else if (ax == Axis::KAXIS) { + return function.template operator()(std::forward(args)...); + } + + PARTHENON_FAIL("Axis should only be one of [IJK]AXIS"); + return function.template operator()(std::forward(args)...); + } +}; + +} // namespace grid +} // namespace kamayan +#endif // GRID_GEOMETRY_HPP_ diff --git a/src/grid/subpack.hpp b/src/grid/subpack.hpp index b62074d..6230b06 100644 --- a/src/grid/subpack.hpp +++ b/src/grid/subpack.hpp @@ -9,6 +9,10 @@ namespace kamayan { enum class Axis { KAXIS = 0, JAXIS = 1, IAXIS = 2 }; +constexpr int AxisToInt(Axis ax) { + return ax == Axis::IAXIS ? 1 : ax == Axis::JAXIS ? 2 : 3; +} + template requires(TemplateSpecialization) struct SubPack_impl { From 02773068dcbc732ad4f26d7dc8ecf2c9e647ed7a Mon Sep 17 00:00:00 2001 From: adam reyes Date: Wed, 18 Mar 2026 21:28:52 -0700 Subject: [PATCH 02/55] cli debug help --- python/kamayan/cli/app.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/python/kamayan/cli/app.py b/python/kamayan/cli/app.py index 57d97d7..14f5bd2 100644 --- a/python/kamayan/cli/app.py +++ b/python/kamayan/cli/app.py @@ -57,11 +57,16 @@ def run( input_file: Optional[str] = None, dry_run: bool = False, info: bool = False, + debug: bool = False, **kwargs, ): parthenon_args: typer.Context = kwargs.pop("ctx") km = self.func(*args, **kwargs) + if debug: + import os + typer.echo(f"PID: {os.getpid()}") + breakpoint() if input_file: km.input_file = input_file @@ -139,6 +144,18 @@ def run( ) func_params += [info_param] + debug_param = inspect.Parameter( + "debug", + inspect.Parameter.KEYWORD_ONLY, + default=typer.Option( + False, + "--debug", + help="Pause execution and wait for debugger to attach", + ), + annotation=bool, + ) + func_params += [debug_param] + ctx_param = inspect.Parameter( "ctx", inspect.Parameter.KEYWORD_ONLY, annotation=typer.Context ) From 25d8ea8caaf73e5ee4205b973f2b4fd81d20a402 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Wed, 18 Mar 2026 21:29:23 -0700 Subject: [PATCH 03/55] seems to work for cartesian --- src/grid/geometry.hpp | 68 ++++++++++++++++++++++++++++++++++------ src/grid/grid.cpp | 58 +++++++++++++++++++++++----------- src/grid/grid_update.hpp | 25 ++++++++------- src/grid/subpack.hpp | 5 +++ src/kamayan/config.hpp | 2 ++ 5 files changed, 118 insertions(+), 40 deletions(-) diff --git a/src/grid/geometry.hpp b/src/grid/geometry.hpp index e00f9cb..98774d7 100644 --- a/src/grid/geometry.hpp +++ b/src/grid/geometry.hpp @@ -11,16 +11,26 @@ #include "dispatcher/options.hpp" #include "grid/grid_types.hpp" #include "grid/subpack.hpp" +#include "pack/sparse_pack.hpp" namespace kamayan { POLYMORPHIC_PARM(Geometry, cartesian, cylindrical); namespace grid { +using GeometryOptions = OptList; + // This is a (incomplete) wrapper around parthenon's uniform_coordinate (cartesian) // to calculate coordinate related items in a different geometry. template struct Coordinates { + KOKKOS_INLINE_FUNCTION Coordinates(parthenon::Coordinates_t coords) : coords_(coords) {} + + template + KOKKOS_INLINE_FUNCTION Coordinates(const parthenon::SparsePack &pack, + const int block) + : coords_(pack.GetCoordinates(block)) {} + // distance between cell centers template KOKKOS_FORCEINLINE_FUNCTION Real Dx() const { @@ -37,7 +47,7 @@ struct Coordinates { // position at cell centroids template - KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int &idx) const { + KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int idx) const { constexpr auto dir = AxisToInt(ax); if constexpr (ax == Axis::IAXIS && geom == Geometry::cylindrical) { const Real r0 = coords_.Xf(idx); @@ -49,23 +59,22 @@ struct Coordinates { return coords_.Xc(idx); } - KOKKOS_FORCEINLINE_FUNCTION Real Xc(const Axis &ax, const int &idx) const { + KOKKOS_FORCEINLINE_FUNCTION Real Xc(const Axis &ax, const int idx) const { return AxisOverload([&, this]() { return Xc(idx); }, ax); } // position at face centers template - KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int &idx) const { + KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int idx) const { return coords_.Xf(idx); } - KOKKOS_FORCEINLINE_FUNCTION Real Xf(const Axis &ax, const int &idx) const { + KOKKOS_FORCEINLINE_FUNCTION Real Xf(const Axis &ax, const int idx) const { return AxisOverload([&, this]() { return Xf(idx); }, ax); } template - KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int &k, const int &j, - const int &i) const { + KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int k, const int j, const int i) const { constexpr auto dir = AxisToInt(ax); if constexpr (geom == Geometry::cylindrical) { if constexpr (ax == Axis::IAXIS) { @@ -81,12 +90,28 @@ struct Coordinates { return coords_.FaceArea(k, j, i); } - KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const Axis &ax, const int &k, const int &j, - const int &i) const { + KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const Axis &ax, const int k, const int j, + const int i) const { return AxisOverload([&, this]() { return FaceArea(k, j, i); }, ax); } - KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int &k, const int &j, const int &i) { + template + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, + const int i) const { + if constexpr (geom == Geometry::cylindrical && ax == Axis::KAXIS) { + // phi aligned edges -- dphi * r + return std::abs(Xf(i)) * Dx(); + } + return coords_.EdgeLength(k, j, i); + } + + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const Axis &ax, const int k, const int j, + const int i) const { + return AxisOverload([&, this]() { return EdgeLength(k, j, i); }, ax); + } + + KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, + const int i) const { if constexpr (geom == Geometry::cylindrical) { // dphi / 2 * (r_p^2 - r_m^2) * dz const Real rp = Xf(i + 1); @@ -96,6 +121,31 @@ struct Coordinates { return coords_.CellVolume(k, j, i); } + // Generalized volumes for a topological element + KOKKOS_FORCEINLINE_FUNCTION + Real Volume(TopologicalElement el, const int k, const int j, const int i) const { + using TE = TopologicalElement; + if (el == TE::CC) + return CellVolume(k, j, i); + else if (el == TE::F1) + return FaceArea(k, j, i); + else if (el == TE::F2) + return FaceArea(k, j, i); + else if (el == TE::F3) + return FaceArea(k, j, i); + else if (el == TE::E1) + return EdgeLength(k, j, i); + else if (el == TE::E2) + return EdgeLength(k, j, i); + else if (el == TE::E3) + return EdgeLength(k, j, i); + else if (el == TE::NN) + return 1.0; + PARTHENON_FAIL("If you reach this point, someone has added a new value to the the " + "TopologicalElement enum."); + return 0.0; + } + private: const parthenon::Coordinates_t coords_; diff --git a/src/grid/grid.cpp b/src/grid/grid.cpp index f6d921e..a632841 100644 --- a/src/grid/grid.cpp +++ b/src/grid/grid.cpp @@ -5,11 +5,14 @@ #include #include "driver/kamayan_driver_types.hpp" +#include "grid/geometry.hpp" #include "grid/grid_refinement.hpp" #include "grid/grid_types.hpp" #include "grid/grid_update.hpp" #include "grid/scratch_variables.hpp" #include "kamayan/runtime_parameters.hpp" +#include "physics/hydro/hydro_types.hpp" +#include "utils/instrument.hpp" namespace kamayan::grid { @@ -21,6 +24,12 @@ std::shared_ptr ProcessUnit() { } void SetupParams(KamayanUnit *unit) { + // Geometry + auto &geometry = unit->AddData("geometry"); + geometry.AddParm( + "geometry", "cartesian", "grid geometry", + {{"cartesian", Geometry::cartesian}, {"cylindrical", Geometry::cylindrical}}); + // most of what we're doing here is wrapping the parthenon mesh related // input parameters as runtime parameters for the docs! // @@ -133,27 +142,38 @@ void InitializeData(KamayanUnit *unit) { } } -TaskStatus FluxesToDuDt(MeshData *md, MeshData *dudt) { - static const int ndim = md->GetNDim(); - using TE = TopologicalElement; - switch (ndim) { - case 1: - FluxDivergence(md, dudt); - break; - case 2: - FluxDivergence(md, dudt); - FluxStokes(md, dudt); - FluxStokes(md, dudt); - break; - case 3: - FluxDivergence(md, dudt); - FluxStokes(md, dudt); - FluxStokes(md, dudt); - FluxStokes(md, dudt); - break; +struct FluxesToDuDt_impl { + using options = OptTypeList; + using value = TaskStatus; + + template + value dispatch(MeshData *md, MeshData *dudt) { + static const int ndim = md->GetNDim(); + using TE = TopologicalElement; + switch (ndim) { + case 1: + FluxDivergence(md, dudt); + break; + case 2: + FluxDivergence(md, dudt); + FluxStokes(md, dudt); + FluxStokes(md, dudt); + break; + case 3: + FluxDivergence(md, dudt); + FluxStokes(md, dudt); + FluxStokes(md, dudt); + FluxStokes(md, dudt); + break; + } + + return TaskStatus::complete; } +}; - return TaskStatus::complete; +TaskStatus FluxesToDuDt(MeshData *md, MeshData *dudt) { + auto cfg = GetConfig(md); + return Dispatcher(PARTHENON_AUTO_LABEL, cfg.get()).execute(md, dudt); } template diff --git a/src/grid/grid_update.hpp b/src/grid/grid_update.hpp index 6fe9a2d..f40fb77 100644 --- a/src/grid/grid_update.hpp +++ b/src/grid/grid_update.hpp @@ -3,13 +3,14 @@ #include #include "grid.hpp" +#include "grid/geometry.hpp" #include "grid_types.hpp" #include "kamayan/fields.hpp" #include "kamayan_utils/parallel.hpp" namespace kamayan::grid { -template +template void FluxDivergence(MeshData *md, MeshData *dudt_data) { static auto desc_cc = GetPackDescriptor(md, {Metadata::Cell, Metadata::WithFluxes}, {PDOpt::WithFluxes}); @@ -26,9 +27,9 @@ void FluxDivergence(MeshData *md, MeshData *dudt_data) { PARTHENON_AUTO_LABEL, 0, nblocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int km, const int jm, const int im) { using TE = TopologicalElement; - const auto coords = u0.GetCoordinates(b); - const Real dxi[sizeof...(faces)]{ - 1.0 / coords.Dxc(faces) % 3 + 1>()...}; + const auto coords = Coordinates(u0, b); + const Real dxi[sizeof...(faces)]{1.0 / + coords.template Dx()...}; // we have to check the variable bounds for each block in case // that we are using sparse fields for (int var = u0.GetLowerBound(b); var <= u0.GetUpperBound(b); var++) { @@ -70,7 +71,7 @@ KOKKOS_INLINE_FUNCTION constexpr int AxisFromFaceEdge(const TopologicalElement & return axis; } -template +template void FluxStokes(MeshData *md, MeshData *dudt_data) { static_assert(Face >= TopologicalElement::F1 && Face <= TopologicalElement::F3, "Face must be F1-F3"); @@ -94,7 +95,7 @@ void FluxStokes(MeshData *md, MeshData *dudt_data) { PARTHENON_AUTO_LABEL, 0, nblocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int km, const int jm, const int im) { using TE = TopologicalElement; - const auto coords = u0.GetCoordinates(b); + const auto coords = Coordinates(u0, b); // we have to check the variable bounds for each block in case // that we are using sparse fields for (int var = u0.GetLowerBound(b); var <= u0.GetUpperBound(b); var++) { @@ -110,18 +111,18 @@ void FluxStokes(MeshData *md, MeshData *dudt_data) { // d_t v_i ~ eps_{ijk}d_j E_k // j -- axis // k -- edge + // cylindrical r, z, phi has odd permutations constexpr Real sign = - ((axis + 1) % 3 == static_cast(edges) % 3) ? 1.0 : -1.0; + (geom == Geometry::cylindrical ? -1.0 : 1.0) * + (((axis + 1) % 3 == static_cast(edges) % 3) ? 1.0 : -1.0); dudt(b, Face, var, km, jm, im) += - sign * (coords.Volume(parthenon::CellLevel::same, edges, km, jm, im) * + sign * (coords.Volume(edges, km, jm, im) * u0.flux(b, edges, var, km, jm, im) - - coords.Volume(parthenon::CellLevel::same, edges, ijk[2], - ijk[1], ijk[0]) * + coords.Volume(edges, ijk[2], ijk[1], ijk[0]) * u0.flux(b, edges, var, ijk[2], ijk[1], ijk[0])); }(), ...); - dudt(b, Face, var, km, jm, im) *= - 1. / coords.Volume(parthenon::CellLevel::same, Face, km, jm, im); + dudt(b, Face, var, km, jm, im) *= 1. / coords.Volume(Face, km, jm, im); } }); } diff --git a/src/grid/subpack.hpp b/src/grid/subpack.hpp index 6230b06..c776d73 100644 --- a/src/grid/subpack.hpp +++ b/src/grid/subpack.hpp @@ -13,6 +13,11 @@ constexpr int AxisToInt(Axis ax) { return ax == Axis::IAXIS ? 1 : ax == Axis::JAXIS ? 2 : 3; } +constexpr Axis AxisFromTE(const TopologicalElement el) { + auto dir = static_cast(el) % 3; + return dir == 0 ? Axis::IAXIS : dir == 1 ? Axis::JAXIS : Axis::KAXIS; +} + template requires(TemplateSpecialization) struct SubPack_impl { diff --git a/src/kamayan/config.hpp b/src/kamayan/config.hpp index 0ac4022..f93425f 100644 --- a/src/kamayan/config.hpp +++ b/src/kamayan/config.hpp @@ -34,6 +34,8 @@ class Config { return _params.template Get(OptInfo::key()); } + void List() { _params.list(); } + private: using Mutability = parthenon::Params::Mutability; parthenon::Params _params; From af5979ef3c2698fb3ab8accc4b4af14221e5d409 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Thu, 19 Mar 2026 18:24:58 -0700 Subject: [PATCH 04/55] coordinates variant --- src/grid/geometry.hpp | 140 +++++++++++++++++++++++++++++++++-- src/grid/grid_refinement.cpp | 2 + src/problems/sedov.cpp | 8 +- 3 files changed, 140 insertions(+), 10 deletions(-) diff --git a/src/grid/geometry.hpp b/src/grid/geometry.hpp index 98774d7..a2789fc 100644 --- a/src/grid/geometry.hpp +++ b/src/grid/geometry.hpp @@ -4,6 +4,8 @@ #include #include +#include + #include "Kokkos_Macros.hpp" #include "coordinates/coordinates.hpp" @@ -12,6 +14,7 @@ #include "grid/grid_types.hpp" #include "grid/subpack.hpp" #include "pack/sparse_pack.hpp" +#include "ports-of-call/variant.hpp" namespace kamayan { POLYMORPHIC_PARM(Geometry, cartesian, cylindrical); @@ -78,10 +81,8 @@ struct Coordinates { constexpr auto dir = AxisToInt(ax); if constexpr (geom == Geometry::cylindrical) { if constexpr (ax == Axis::IAXIS) { - // 2 * pi * r_i * dz return Xf(i) * Dx() * Dx(); } else if constexpr (ax == Axis::JAXIS) { - // pi * (r_i+1/2**2 - r_i-1/2**2) const Real rp = Xf(i + 1); const Real rm = Xf(i); return 0.5 * (rp * rp - rm * rm) * Dx(); @@ -99,7 +100,6 @@ struct Coordinates { KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, const int i) const { if constexpr (geom == Geometry::cylindrical && ax == Axis::KAXIS) { - // phi aligned edges -- dphi * r return std::abs(Xf(i)) * Dx(); } return coords_.EdgeLength(k, j, i); @@ -113,7 +113,6 @@ struct Coordinates { KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, const int i) const { if constexpr (geom == Geometry::cylindrical) { - // dphi / 2 * (r_p^2 - r_m^2) * dz const Real rp = Xf(i + 1); const Real rm = Xf(i); return 0.5 * Dx() * (rp * rp - rm * rm); @@ -121,7 +120,6 @@ struct Coordinates { return coords_.CellVolume(k, j, i); } - // Generalized volumes for a topological element KOKKOS_FORCEINLINE_FUNCTION Real Volume(TopologicalElement el, const int k, const int j, const int i) const { using TE = TopologicalElement; @@ -147,10 +145,11 @@ struct Coordinates { } private: - const parthenon::Coordinates_t coords_; + parthenon::Coordinates_t coords_; template - KOKKOS_FORCEINLINE_FUNCTION Real AxisOverload(Func function, Axis ax, Args &&...args) { + KOKKOS_FORCEINLINE_FUNCTION Real AxisOverload(Func function, Axis ax, + Args &&...args) const { if (ax == Axis::IAXIS) { return function.template operator()(std::forward(args)...); } else if (ax == Axis::JAXIS) { @@ -164,6 +163,133 @@ struct Coordinates { } }; +namespace impl { +template +struct CoordinatesVariant {}; + +template +struct CoordinatesVariant> { + using type = PortsOfCall::variant, Coordinates...>; + + static KOKKOS_INLINE_FUNCTION type Get(const Geometry geometry, + const parthenon::Coordinates_t coords) { + type coordinates = Coordinates(coords); + ( + [&]() { + if (geoms == geometry) coordinates = Coordinates(coords); + }(), + ...); + + return coordinates; + } +}; + +} // namespace impl + +using CoordinatesManager = impl::CoordinatesVariant; +using CoordinatesVariant = CoordinatesManager::type; + +namespace impl { +KOKKOS_INLINE_FUNCTION CoordinatesManager::type GetCoordinates(MeshBlock *mb) { + const Geometry geometry = GetConfig(mb)->Get(); + return CoordinatesManager::Get(geometry, mb->coords); +} + +KOKKOS_INLINE_FUNCTION CoordinatesManager::type +GetCoordinates(const Geometry geometry, const parthenon::Coordinates_t &coords) { + return CoordinatesManager::Get(geometry, coords); +} +} // namespace impl + +struct GenericCoordinate { + KOKKOS_INLINE_FUNCTION GenericCoordinate(const CoordinatesVariant coords) + : coords_(coords) {} + + KOKKOS_INLINE_FUNCTION GenericCoordinate(MeshBlock *mb) + : coords_(impl::GetCoordinates(mb)) {} + KOKKOS_INLINE_FUNCTION GenericCoordinate(const Geometry geometry, + const parthenon::Coordinates_t &coords) + : coords_(impl::GetCoordinates(geometry, coords)) {} + + template + KOKKOS_FORCEINLINE_FUNCTION Real Dx() const { + return PortsOfCall::visit([](const auto &coords) { return coords.template Dx(); }, + coords_); + } + + KOKKOS_FORCEINLINE_FUNCTION Real Dx(const Axis &ax) const { + return PortsOfCall::visit([&ax](const auto &coords) { return coords.Dx(ax); }, + coords_); + } + + template + KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int idx) const { + return PortsOfCall::visit( + [idx](const auto &coords) { return coords.template Xc(idx); }, coords_); + } + + KOKKOS_FORCEINLINE_FUNCTION Real Xc(const Axis &ax, const int idx) const { + return PortsOfCall::visit( + [&ax, idx](const auto &coords) { return coords.Xc(ax, idx); }, coords_); + } + + template + KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int idx) const { + return PortsOfCall::visit( + [idx](const auto &coords) { return coords.template Xf(idx); }, coords_); + } + + KOKKOS_FORCEINLINE_FUNCTION Real Xf(const Axis &ax, const int idx) const { + return PortsOfCall::visit( + [&ax, idx](const auto &coords) { return coords.Xf(ax, idx); }, coords_); + } + + template + KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int k, const int j, const int i) const { + return PortsOfCall::visit( + [k, j, i](const auto &coords) { return coords.template FaceArea(k, j, i); }, + coords_); + } + + KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const Axis &ax, const int k, const int j, + const int i) const { + return PortsOfCall::visit( + [&ax, k, j, i](const auto &coords) { return coords.FaceArea(ax, k, j, i); }, + coords_); + } + + template + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, + const int i) const { + return PortsOfCall::visit( + [k, j, i](const auto &coords) { return coords.template EdgeLength(k, j, i); }, + coords_); + } + + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const Axis &ax, const int k, const int j, + const int i) const { + return PortsOfCall::visit( + [&ax, k, j, i](const auto &coords) { return coords.EdgeLength(ax, k, j, i); }, + coords_); + } + + KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, + const int i) const { + return PortsOfCall::visit( + [k, j, i](const auto &coords) { return coords.CellVolume(k, j, i); }, coords_); + } + + KOKKOS_FORCEINLINE_FUNCTION Real Volume(TopologicalElement el, const int k, const int j, + const int i) const { + return PortsOfCall::visit( + [el, k, j, i](const auto &coords) { return coords.Volume(el, k, j, i); }, + coords_); + } + + private: + CoordinatesVariant coords_; +}; + } // namespace grid } // namespace kamayan #endif // GRID_GEOMETRY_HPP_ diff --git a/src/grid/grid_refinement.cpp b/src/grid/grid_refinement.cpp index b1bcd28..38e246a 100644 --- a/src/grid/grid_refinement.cpp +++ b/src/grid/grid_refinement.cpp @@ -62,6 +62,8 @@ void AMRLoehner::operator()(MeshData *md, PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, kb.s - k3d, kb.e + k3d, jb.s - k2d, jb.e + k2d, ib.s - 1, ib.e + 1, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { + // TODO(acreyes): this could be updated to be templated on geometry, but + // would only matter for something more than 2D rz cylindrical to be different const auto coords = pack.GetCoordinates(b); using TE = TopologicalElement; // --8<-- [start:FirstDer] diff --git a/src/problems/sedov.cpp b/src/problems/sedov.cpp index 482466a..6fa7a16 100644 --- a/src/problems/sedov.cpp +++ b/src/problems/sedov.cpp @@ -6,8 +6,10 @@ #include "Kokkos_Macros.hpp" #include "driver/kamayan_driver_types.hpp" +#include "grid/geometry.hpp" #include "grid/grid.hpp" #include "grid/grid_types.hpp" +#include "grid/subpack.hpp" #include "kamayan/fields.hpp" #include "kamayan/kamayan.hpp" #include "kamayan/unit.hpp" @@ -109,7 +111,7 @@ void ProblemGenerator(MeshBlock *mb) { auto jb = cellbounds.GetBoundsJ(IndexDomain::interior); auto kb = cellbounds.GetBoundsK(IndexDomain::interior); - auto coords = mb->coords; + auto coords = grid::GenericCoordinate(mb); auto nspecies = mb->packages.Get("material")->Param("nspecies"); if (nspecies > 1) { @@ -123,8 +125,8 @@ void ProblemGenerator(MeshBlock *mb) { 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 Real r2 = - coords.Xc<1>(i) * coords.Xc<1>(i) + coords.Xc<2>(j) * coords.Xc<2>(j); + const Real r2 = coords.Xc(i) * coords.Xc(i) + + coords.Xc(j) * coords.Xc(j); const auto r = Kokkos::sqrt(r2); auto state = sedov_data.State(r); type_for(SedovData::variables(), [&](const Vars &) { From 1cfd82105fb46c86be34e1be7b38781975f1d816 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Thu, 19 Mar 2026 18:32:41 -0700 Subject: [PATCH 05/55] generic coords in python --- src/grid/pybind/grid_bindings.cpp | 88 +++++++++++++++++++------------ src/problems/shock_tube.py | 6 +-- 2 files changed, 56 insertions(+), 38 deletions(-) diff --git a/src/grid/pybind/grid_bindings.cpp b/src/grid/pybind/grid_bindings.cpp index 16475f9..d01a53f 100644 --- a/src/grid/pybind/grid_bindings.cpp +++ b/src/grid/pybind/grid_bindings.cpp @@ -12,7 +12,7 @@ #include "kamayan/pybind/kamayan_bindings.hpp" #include "kamayan/pybind/kamayan_nanobind.h" -#include "coordinates/coordinates.hpp" +#include "grid/geometry.hpp" #include "grid/grid_types.hpp" #include "kamayan/config.hpp" #include "kamayan_utils/parallel.hpp" @@ -29,6 +29,7 @@ struct SparsePack_py { auto desc = parthenon::MakePackDescriptor(pkg, vars, {}, {PDOpt::WithFluxes}); pack = desc.GetPack(mb->meshblock_data.Get().get()); map = desc.GetMap(); + config_ = GetConfig(mb); } SparsePack_py(MeshData *md, const std::vector vars, @@ -37,6 +38,7 @@ struct SparsePack_py { md->GetMeshPointer()->resolved_packages.get(), vars, {}, options); pack = desc.GetPack(md); map = desc.GetMap(); + config_ = GetConfig(md); } parthenon::ParArray3D GetParArray3D(const int &block, const std::string &var, @@ -46,13 +48,16 @@ struct SparsePack_py { return pack(block, te, idx + comp); } - parthenon::Coordinates_t GetCoordinates(const int b) const { - return pack.GetCoordinates(b); + grid::GenericCoordinate GetCoordinates(const int b) const { + auto cfg = config_.lock(); + const Geometry geometry = cfg->Get(); + return grid::GenericCoordinate(geometry, pack.GetCoordinates(b)); } private: parthenon::SparsePack<> pack; parthenon::SparsePackIdxMap map; + std::weak_ptr config_; }; void sparse_pack_py(nanobind::module_ &m) { @@ -139,62 +144,75 @@ void grid_module(nanobind::module_ &m) { te.value("E3", TopologicalElement::E3); te.value("NN", TopologicalElement::NN); - using Coordinates_t = parthenon::Coordinates_t; - nanobind::class_ coords(m, "Coordinates_t"); - coords.def("Dx", [](Coordinates_t &self, const int dir) { return self.Dx(dir); }); - coords.def("Dx1", &Coordinates_t::Dx<1>); - coords.def("Dx2", &Coordinates_t::Dx<2>); - coords.def("Dx3", &Coordinates_t::Dx<3>); + nanobind::class_ gen_coords(m, "GenericCoordinate"); + + gen_coords.def( + "Dx", [](grid::GenericCoordinate &self, const Axis &ax) { return self.Dx(ax); }); + gen_coords.def("Dx1", + [](grid::GenericCoordinate &self) { return self.Dx(); }); + gen_coords.def("Dx2", + [](grid::GenericCoordinate &self) { return self.Dx(); }); + gen_coords.def("Dx3", + [](grid::GenericCoordinate &self) { return self.Dx(); }); Vectorize( - coords, "Xc", - KOKKOS_LAMBDA(const Coordinates_t &self, const int &idx, const int &dir) { - return (dir == 1) * self.Xc<1>(idx) + (dir == 2) * self.Xc<2>(idx) + - (dir == 3) * self.Xc<3>(idx); + gen_coords, "Xc", + KOKKOS_LAMBDA(const grid::GenericCoordinate &self, const int &idx, const int &dir) { + return (dir == 1) * self.Xc(idx) + + (dir == 2) * self.Xc(idx) + + (dir == 3) * self.Xc(idx); }, TypeList()); Vectorize( - coords, "Xf", - KOKKOS_LAMBDA(const Coordinates_t &self, const int &idx, const int &dir) { - return (dir == 1) * self.Xf<1>(idx) + (dir == 2) * self.Xc<2>(idx) + - (dir == 3) * self.Xf<3>(idx); + gen_coords, "Xf", + KOKKOS_LAMBDA(const grid::GenericCoordinate &self, const int &idx, const int &dir) { + return (dir == 1) * self.Xf(idx) + + (dir == 2) * self.Xf(idx) + + (dir == 3) * self.Xf(idx); }, TypeList()); Vectorize( - coords, "Xc1", - KOKKOS_LAMBDA(const Coordinates_t &self, const int &idx) { - return self.Xc<1>(idx); + gen_coords, "Xc1", + KOKKOS_LAMBDA(const grid::GenericCoordinate &self, const int &idx) { + return self.Xc(idx); }, TypeList()); Vectorize( - coords, "Xc2", - KOKKOS_LAMBDA(const Coordinates_t &self, const int &idx) { - return self.Xc<2>(idx); + gen_coords, "Xc2", + KOKKOS_LAMBDA(const grid::GenericCoordinate &self, const int &idx) { + return self.Xc(idx); }, TypeList()); Vectorize( - coords, "Xc3", - KOKKOS_LAMBDA(const Coordinates_t &self, const int &idx) { - return self.Xc<3>(idx); + gen_coords, "Xc3", + KOKKOS_LAMBDA(const grid::GenericCoordinate &self, const int &idx) { + return self.Xc(idx); }, TypeList()); Vectorize( - coords, "Xf1", - KOKKOS_LAMBDA(const Coordinates_t &self, const int &idx) { - return self.Xf<1>(idx); + gen_coords, "Xf1", + KOKKOS_LAMBDA(const grid::GenericCoordinate &self, const int &idx) { + return self.Xf(idx); }, TypeList()); Vectorize( - coords, "Xf2", - KOKKOS_LAMBDA(const Coordinates_t &self, const int &idx) { - return self.Xf<2>(idx); + gen_coords, "Xf2", + KOKKOS_LAMBDA(const grid::GenericCoordinate &self, const int &idx) { + return self.Xf(idx); }, TypeList()); Vectorize( - coords, "Xf3", - KOKKOS_LAMBDA(const Coordinates_t &self, const int &idx) { - return self.Xf<3>(idx); + gen_coords, "Xf3", + KOKKOS_LAMBDA(const grid::GenericCoordinate &self, const int &idx) { + return self.Xf(idx); }, TypeList()); + + gen_coords.def("FaceArea", [](grid::GenericCoordinate &self, const Axis &ax, int k, + int j, int i) { return self.FaceArea(ax, k, j, i); }); + gen_coords.def("EdgeLength", [](grid::GenericCoordinate &self, const Axis &ax, int k, + int j, int i) { return self.EdgeLength(ax, k, j, i); }); + gen_coords.def("CellVolume", &grid::GenericCoordinate::CellVolume); + gen_coords.def("Volume", &grid::GenericCoordinate::Volume); } } // namespace kamayan diff --git a/src/problems/shock_tube.py b/src/problems/shock_tube.py index 75f5250..b194c28 100644 --- a/src/problems/shock_tube.py +++ b/src/problems/shock_tube.py @@ -9,7 +9,7 @@ from kamayan.code_units import Grid, driver, eos, physics from kamayan.code_units.Grid import AdaptiveGrid import kamayan.pyKamayan as pyKamayan -from kamayan.pyKamayan.Grid import Coordinates_t, TopologicalElement as te +from kamayan.pyKamayan.Grid import GenericCoordinate, TopologicalElement as te from kamayan.code_units.Hydro import Hydro import kamayan.kamayan_manager as kman from kamayan.kamayan_manager import KamayanManager @@ -132,7 +132,7 @@ def _get_state(s: Literal["L", "R"]) -> State: def _get_coords( - indices: np.ndarray, coords: Coordinates_t, el: te = te.CC + indices: np.ndarray, coords: GenericCoordinate, el: te = te.CC ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: xx = np.array(coords.Xc1(indices[2, :, :, :])) yy = np.array(coords.Xc2(indices[1, :, :, :])) @@ -163,7 +163,7 @@ def _x_perp( def _calc_face_mag( face: np.ndarray, el: te, - coords: Coordinates_t, + coords: GenericCoordinate, n_perp: np.ndarray, n_par1: np.ndarray, vL: State, From 40720c36b02a70795a2726cb6c91f20709ef53bd Mon Sep 17 00:00:00 2001 From: adam reyes Date: Thu, 19 Mar 2026 23:32:47 -0700 Subject: [PATCH 06/55] cylindrical source terms --- src/grid/geometry.hpp | 81 +++++++++++++++++++++++++++++- src/grid/grid_types.hpp | 10 ++++ src/grid/subpack.hpp | 44 ++++++++-------- src/physics/hydro/hydro.cpp | 59 ++++++++++++++++++++-- src/physics/hydro/hydro.hpp | 3 +- src/physics/hydro/primconsflux.cpp | 25 +++++++-- src/physics/hydro/primconsflux.hpp | 22 ++++++++ 7 files changed, 211 insertions(+), 33 deletions(-) diff --git a/src/grid/geometry.hpp b/src/grid/geometry.hpp index a2789fc..aa3874b 100644 --- a/src/grid/geometry.hpp +++ b/src/grid/geometry.hpp @@ -12,7 +12,6 @@ #include "dispatcher/options.hpp" #include "grid/grid_types.hpp" -#include "grid/subpack.hpp" #include "pack/sparse_pack.hpp" #include "ports-of-call/variant.hpp" @@ -163,6 +162,86 @@ struct Coordinates { } }; +template +concept CoordinateSystem = requires(T coords) { + // Templated methods (compile-time axis dispatch) + { coords.template Dx() } -> std::same_as; + { coords.template Xc(1) } -> std::same_as; + { coords.template Xf(1) } -> std::same_as; + { coords.template FaceArea(0, 0, 0) } -> std::same_as; + { coords.template EdgeLength(0, 0, 0) } -> std::same_as; + // Runtime methods (runtime axis dispatch) + { coords.Dx(Axis::IAXIS) } -> std::same_as; + { coords.Xc(Axis::IAXIS, 1) } -> std::same_as; + { coords.Xf(Axis::IAXIS, 1) } -> std::same_as; + { coords.FaceArea(Axis::IAXIS, 0, 0, 0) } -> std::same_as; + { coords.EdgeLength(Axis::IAXIS, 0, 0, 0) } -> std::same_as; + { coords.CellVolume(0, 0, 0) } -> std::same_as; + { coords.Volume(TopologicalElement::CC, 0, 0, 0) } -> std::same_as; +}; + +template +struct CoordinateIndexer { + KOKKOS_INLINE_FUNCTION CoordinateIndexer(const T &coords, const int k, const int j, + const int i) + : coords_(coords), kji_({k, j, i}) {} + + template + KOKKOS_FORCEINLINE_FUNCTION Real Dx() const { + return coords_.template Dx(); + } + + KOKKOS_FORCEINLINE_FUNCTION Real Dx(const Axis &ax) const { return coords_.Dx(ax); } + + template + KOKKOS_FORCEINLINE_FUNCTION Real Xc() const { + return coords_.template Xc(kji_[static_cast(ax)]); + } + + KOKKOS_FORCEINLINE_FUNCTION Real Xc(const Axis &ax) const { + return coords_.Xc(ax, kji_[static_cast(ax)]); + } + + template + KOKKOS_FORCEINLINE_FUNCTION Real Xf() const { + return coords_.template Xf(kji_[static_cast(ax)]); + } + + KOKKOS_FORCEINLINE_FUNCTION Real Xf(const Axis &ax) const { + return coords_.Xf(ax, kji_[static_cast(ax)]); + } + + template + KOKKOS_FORCEINLINE_FUNCTION Real FaceArea() const { + return coords_.template FaceArea(kji_[2], kji_[1], kji_[0]); + } + + KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const Axis &ax) const { + return coords_.FaceArea(ax, kji_[2], kji_[1], kji_[0]); + } + + template + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength() const { + return coords_.template EdgeLength(kji_[2], kji_[1], kji_[0]); + } + + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const Axis &ax) const { + return coords_.EdgeLength(ax, kji_[2], kji_[1], kji_[0]); + } + + KOKKOS_FORCEINLINE_FUNCTION Real CellVolume() const { + return coords_.CellVolume(kji_[2], kji_[1], kji_[0]); + } + + KOKKOS_FORCEINLINE_FUNCTION Real Volume(TopologicalElement el) const { + return coords_.Volume(el, kji_[2], kji_[1], kji_[0]); + } + + private: + const T coords_; + const Kokkos::Array kji_; +}; + namespace impl { template struct CoordinatesVariant {}; diff --git a/src/grid/grid_types.hpp b/src/grid/grid_types.hpp index 0d44dab..e2c0c1b 100644 --- a/src/grid/grid_types.hpp +++ b/src/grid/grid_types.hpp @@ -7,6 +7,7 @@ namespace kamayan { using Real = parthenon::Real; using TopologicalElement = parthenon::TopologicalElement; using TopologicalType = parthenon::TopologicalType; +enum class Axis { KAXIS = 0, JAXIS = 1, IAXIS = 2 }; KOKKOS_INLINE_FUNCTION constexpr TopologicalElement IncrementTE(const TopologicalElement &out_te, const TopologicalElement &in_te, @@ -16,6 +17,15 @@ IncrementTE(const TopologicalElement &out_te, const TopologicalElement &in_te, return static_cast(out); } +constexpr int AxisToInt(Axis ax) { + return ax == Axis::IAXIS ? 1 : ax == Axis::JAXIS ? 2 : 3; +} + +constexpr Axis AxisFromTE(const TopologicalElement el) { + auto dir = static_cast(el) % 3; + return dir == 0 ? Axis::IAXIS : dir == 1 ? Axis::JAXIS : Axis::KAXIS; +} + template concept EdgeElement = (edge >= TopologicalElement::E1 && edge <= TopologicalElement::E3); diff --git a/src/grid/subpack.hpp b/src/grid/subpack.hpp index c776d73..e2bef42 100644 --- a/src/grid/subpack.hpp +++ b/src/grid/subpack.hpp @@ -3,26 +3,17 @@ #include +#include "grid/geometry.hpp" #include "grid/grid_types.hpp" #include "kamayan_utils/type_abstractions.hpp" namespace kamayan { -enum class Axis { KAXIS = 0, JAXIS = 1, IAXIS = 2 }; - -constexpr int AxisToInt(Axis ax) { - return ax == Axis::IAXIS ? 1 : ax == Axis::JAXIS ? 2 : 3; -} - -constexpr Axis AxisFromTE(const TopologicalElement el) { - auto dir = static_cast(el) % 3; - return dir == 0 ? Axis::IAXIS : dir == 1 ? Axis::JAXIS : Axis::KAXIS; -} template requires(TemplateSpecialization) struct SubPack_impl { - KOKKOS_INLINE_FUNCTION SubPack_impl(PackType &pack, const int &b, const int &k, - const int &j, const int &i) + KOKKOS_INLINE_FUNCTION SubPack_impl(PackType &pack, const int b, const int k, + const int j, const int i) : pack_(pack), b_(b), k_(k), j_(j), i_(i) {} template @@ -45,6 +36,13 @@ struct SubPack_impl { return pack_.GetSize(b_, var); } + template + KOKKOS_INLINE_FUNCTION grid::CoordinateIndexer> + GetCoordinates() const { + return grid::CoordinateIndexer>(pack_.GetCoordinates(b_), k_, + j_, i_); + } + private: PackType &pack_; const int b_, k_, j_, i_; @@ -53,8 +51,8 @@ struct SubPack_impl { template requires(TemplateSpecialization) struct StencilSubPack_impl { - KOKKOS_INLINE_FUNCTION StencilSubPack_impl(PackType &pack, const int &b, const int &k, - const int &j, const int &i) + KOKKOS_INLINE_FUNCTION StencilSubPack_impl(PackType &pack, const int b, const int k, + const int j, const int i) : pack_(pack), b_(b), kji_({k, j, i}) {} template @@ -100,9 +98,9 @@ struct StencilSubPack_impl { template requires(TemplateSpecialization) struct VarStencilSubPack_impl { - KOKKOS_INLINE_FUNCTION VarStencilSubPack_impl(PackType &pack, const int &b, - const Var_t &var, const int &k, - const int &j, const int &i) + KOKKOS_INLINE_FUNCTION VarStencilSubPack_impl(PackType &pack, const int b, + const Var_t &var, const int k, + const int j, const int i) : pack_(pack), b_(b), var_(var), kji_({k, j, i}) {} template @@ -140,20 +138,20 @@ struct VarStencilSubPack_impl { }; template -KOKKOS_INLINE_FUNCTION auto SubPack(PackType &pack, const int &b, const Var_t &var, - const int &k, const int &j, const int &i) { +KOKKOS_INLINE_FUNCTION auto SubPack(PackType &pack, const int b, const Var_t &var, + const int k, const int j, const int i) { return VarStencilSubPack_impl(pack, b, var, k, j, i); } template -KOKKOS_INLINE_FUNCTION auto SubPack(PackType &pack, const int &b, const int &k, - const int &j, const int &i) { +KOKKOS_INLINE_FUNCTION auto SubPack(PackType &pack, const int b, const int k, const int j, + const int i) { return StencilSubPack_impl(pack, b, k, j, i); } template -KOKKOS_INLINE_FUNCTION auto SubPack(PackType &pack, const int &b, const int &k, - const int &j, const int &i) { +KOKKOS_INLINE_FUNCTION auto SubPack(PackType &pack, const int b, const int k, const int j, + const int i) { return SubPack_impl(pack, b, k, j, i); } diff --git a/src/physics/hydro/hydro.cpp b/src/physics/hydro/hydro.cpp index daf11b9..bebf2bc 100644 --- a/src/physics/hydro/hydro.cpp +++ b/src/physics/hydro/hydro.cpp @@ -9,8 +9,11 @@ #include "dispatcher/options.hpp" #include "driver/kamayan_driver_types.hpp" +#include "grid/geometry.hpp" #include "grid/grid.hpp" +#include "grid/subpack.hpp" #include "hydro_types.hpp" +#include "kamayan/config.hpp" #include "kamayan/fields.hpp" #include "kamayan/unit_data.hpp" #include "kamayan_utils/parallel.hpp" @@ -30,6 +33,7 @@ std::shared_ptr ProcessUnit() { hydro->PreparePrimitive.Register(PreparePrimitive, /*after=*/{}, /*before=*/{"eos"}); hydro->PostMeshInitialization.Register(PostMeshInitialization, {"eos"}); hydro->AddFluxTasks.Register(AddFluxTasks); + hydro->AddTasksOneStep.Register(AddTasksOneStep); // --8<-- [end:register] return hydro; } @@ -140,20 +144,21 @@ struct FillDerived_impl { auto jb = md->GetBoundsJ(IndexDomain::interior); auto kb = md->GetBoundsK(IndexDomain::interior); const auto ndim = md->GetNDim(); + const auto geometry = GetConfig(md)->Get(); parthenon::par_for( PARTHENON_AUTO_LABEL, 0, nblocks - 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) { capture(ndim); - const auto coords = pack.GetCoordinates(b); + const auto coords = grid::GenericCoordinate(geometry, pack.GetCoordinates(b)); // also need to average the face-fields if doing constrained transport if constexpr (hydro_traits::MHD == Mhd::ct) { using te = TopologicalElement; if (ndim > 1) { - pack(b, DIVB(), k, j, i) = 1. / coords.template Dxc<1>() * + pack(b, DIVB(), k, j, i) = 1. / coords.template Dx() * (pack(b, te::F1, MAG(), k, j, i + 1) - pack(b, te::F1, MAG(), k, j, i)) + - 1. / coords.template Dxc<2>() * + 1. / coords.template Dx() * (pack(b, te::F2, MAG(), k, j + 1, i) - pack(b, te::F2, MAG(), k, j, i)); } @@ -169,4 +174,52 @@ TaskStatus FillDerived(MeshData *md) { return Dispatcher(PARTHENON_AUTO_LABEL, cfg.get()).execute(md); } +struct AddSourceTerms { + using options = OptTypeList; + using value = TaskStatus; + + template + requires(geom == Geometry::cartesian) + value dispatch(MeshData *md, MeshData *dudt) { + return TaskStatus::complete; + } + template + value dispatch(MeshData *md, MeshData *dudt) { + using primitives = hydro_traits::Primitive; + using conserved = hydro_traits::Conserved; + + auto pack_dudt = grid::GetPack(conserved(), dudt); + auto pack_prim = grid::GetPack(primitives(), md); + + const int nblocks = pack_prim.GetNBlocks(); + auto ib = md->GetBoundsI(IndexDomain::interior); + auto jb = md->GetBoundsJ(IndexDomain::interior); + auto kb = md->GetBoundsK(IndexDomain::interior); + + par_for( + PARTHENON_AUTO_LABEL, 0, nblocks - 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) { + auto prim = SubPack(pack_prim, b, k, j, i); + auto du = SubPack(pack_dudt, b, k, j, i); + AddGeometricSource( + prim.template GetCoordinates(), prim, du); + }); + return TaskStatus::complete; + } +}; + +TaskID AddTasksOneStep(TaskID prev, TaskList &tl, MeshData *md, MeshData *dudt) { + auto sources = tl.AddTask( + prev, "hydro::AddSourceTerms", + [](MeshData *md, MeshData *dudt) { + auto cfg = GetConfig(md); + return Dispatcher(PARTHENON_AUTO_LABEL, cfg.get()) + .execute(md, dudt); + }, + md, dudt); + + auto final = sources; + return final; +} + } // namespace kamayan::hydro diff --git a/src/physics/hydro/hydro.hpp b/src/physics/hydro/hydro.hpp index bcb3d03..97131d6 100644 --- a/src/physics/hydro/hydro.hpp +++ b/src/physics/hydro/hydro.hpp @@ -5,10 +5,8 @@ #include "driver/kamayan_driver_types.hpp" #include "grid/grid_types.hpp" -#include "kamayan/config.hpp" #include "kamayan/runtime_parameters.hpp" #include "kamayan/unit.hpp" -#include "kamayan/unit_data.hpp" namespace kamayan::hydro { std::shared_ptr ProcessUnit(); @@ -19,6 +17,7 @@ void InitializeData(KamayanUnit *unit); TaskID AddFluxTasks(TaskID prev, TaskList &tl, MeshData *md); TaskID AddTasksOneStep(TaskID prev, TaskList &tl, MeshData *md, MeshData *dudt); Real EstimateTimeStepMesh(MeshData *md); +TaskStatus PrepareConserved(MeshData *md); TaskStatus FillDerived(MeshData *md); diff --git a/src/physics/hydro/primconsflux.cpp b/src/physics/hydro/primconsflux.cpp index 7f195dc..bdb387d 100644 --- a/src/physics/hydro/primconsflux.cpp +++ b/src/physics/hydro/primconsflux.cpp @@ -1,6 +1,7 @@ #include "physics/hydro/primconsflux.hpp" #include "dispatcher/options.hpp" #include "driver/kamayan_driver_types.hpp" +#include "grid/geometry.hpp" #include "grid/grid.hpp" #include "grid/grid_types.hpp" #include "grid/indexer.hpp" @@ -13,10 +14,10 @@ namespace kamayan::hydro { // --8<-- [start:impl] struct PrepareConserved_impl { - using options = OptTypeList; + using options = OptTypeList; using value = TaskStatus; - template + template requires(NonTypeTemplateSpecialization) value dispatch(MeshData *md) { auto pack = grid::GetPack(typename hydro_traits::ConsPrim(), md); @@ -48,6 +49,14 @@ struct PrepareConserved_impl { auto U = SubPack(pack, b, k, j, i); Prim2Cons(U, U); // --8<-- [end:make-idx] + if constexpr (geom == Geometry::cylindrical) { + auto coords = grid::Coordinates(pack, b); + // conserve angular momentum + U(MOMENTUM(2)) *= coords.template Xc(i); + if constexpr (hydro_traits::MHD != Mhd::off) { + U(MAGC(2)) *= 1.0 / coords.template Xc(i); + } + } }); return TaskStatus::complete; } @@ -55,10 +64,10 @@ struct PrepareConserved_impl { // --8<-- [end:impl] struct PreparePrimitive_impl { - using options = OptTypeList; + using options = OptTypeList; using value = TaskStatus; - template + template requires(NonTypeTemplateSpecialization) value dispatch(MeshData *md) { auto pack = grid::GetPack(typename hydro_traits::ConsPrim(), md); @@ -87,6 +96,14 @@ struct PreparePrimitive_impl { } auto U = SubPack(pack, b, k, j, i); Cons2Prim(U, U); + if constexpr (geom == Geometry::cylindrical) { + auto coords = grid::Coordinates(pack, b); + // conserve angular momentum + U(VELOCITY(2)) *= 1.0 / coords.template Xc(i); + if constexpr (hydro_traits::MHD != Mhd::off) { + U(MAGC(2)) *= coords.template Xc(i); + } + } }); return TaskStatus::complete; } diff --git a/src/physics/hydro/primconsflux.hpp b/src/physics/hydro/primconsflux.hpp index 3012962..a38e8c9 100644 --- a/src/physics/hydro/primconsflux.hpp +++ b/src/physics/hydro/primconsflux.hpp @@ -2,8 +2,11 @@ #define PHYSICS_HYDRO_PRIMCONSFLUX_HPP_ #include +#include "Kokkos_Macros.hpp" #include "driver/kamayan_driver_types.hpp" +#include "grid/geometry.hpp" #include "grid/grid_types.hpp" +#include "grid/subpack.hpp" #include "kamayan/fields.hpp" #include "physics/hydro/hydro_types.hpp" #include "physics/physics_types.hpp" @@ -134,5 +137,24 @@ KOKKOS_INLINE_FUNCTION void Prim2Flux(const Prim &V, Flux &F) { } } +template +concept CoordinatePoint = requires(T coords) { + { coords.template Xc() } -> std::same_as; +}; + +template +requires(geom == Geometry::cylindrical) +KOKKOS_FORCEINLINE_FUNCTION void AddGeometricSource(const Coords &coords, const Prim &V, + DuDt &dudt) { + dudt(MOMENTUM(0)) = (V(DENS()) * V(VELOCITY(2)) * V(VELOCITY(2)) + V(PRES())) / + coords.template Xc(); + if constexpr (mhd != Mhd::off) { + const Real pmag = 0.5 * (V(MAGC(0)) * V(MAGC(0)) + V(MAGC(1)) * V(MAGC(1)) + + V(MAGC(2)) * V(MAGC(2))); + dudt(MOMENTUM(0)) += + (pmag - V(MAGC(2)) * V(MAGC(2))) / coords.template Xc(); + } +} + } // namespace kamayan::hydro #endif // PHYSICS_HYDRO_PRIMCONSFLUX_HPP_ From c9c3669bf4d374032292de95cdf60e2176b0b2b7 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 22 Mar 2026 16:09:20 -0700 Subject: [PATCH 07/55] add boundary conditions --- python/kamayan/code_units/Grid.py | 1 + src/CMakeLists.txt | 1 + src/grid/boundary_conditions.cpp | 83 +++++++++++++++++++++++++++++++ src/grid/grid.hpp | 2 + src/kamayan/fields.hpp | 14 ++++-- src/problems/sedov.py | 5 +- 6 files changed, 100 insertions(+), 6 deletions(-) create mode 100644 src/grid/boundary_conditions.cpp diff --git a/python/kamayan/code_units/Grid.py b/python/kamayan/code_units/Grid.py index 23f2503..875b07f 100644 --- a/python/kamayan/code_units/Grid.py +++ b/python/kamayan/code_units/Grid.py @@ -12,6 +12,7 @@ _resolution = None | int _dx = None | float refinement_strategy = Literal["none", "static", "adaptive"] +GEOMETRY = Literal["cartesian", "cylindrical"] def _get_N(N: _resolution, dx: _dx, L: float) -> int: diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5d5386a..1076482 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,5 +1,6 @@ set(_sources driver/kamayan_driver.cpp + grid/boundary_conditions.cpp grid/grid.cpp grid/grid_refinement.cpp kamayan/callback_dag.cpp diff --git a/src/grid/boundary_conditions.cpp b/src/grid/boundary_conditions.cpp new file mode 100644 index 0000000..3d1b293 --- /dev/null +++ b/src/grid/boundary_conditions.cpp @@ -0,0 +1,83 @@ +#include "bvals/boundary_conditions.hpp" +#include +#include +#include +#include + +#include + +#include "basic_types.hpp" +#include "grid/grid.hpp" +#include "grid/grid_types.hpp" +#include "kamayan/fields.hpp" +#include "pack/pack_utils.hpp" + +namespace kamayan::grid { +using TE = TopologicalElement; +// use this as an interface to write boundary conditions per variable topological element +// types that can be registered through parthenon's appliciation input interface at the +// bottom of this file +using BoundaryConditionInterface = + std::function &, bool)>; + +void Axisymmetric(const TE el, std::shared_ptr &mbd, bool coarse) { + const auto ttFlag = [el] { + const auto tt = GetTopologicalType(el); + switch (tt) { + case (TopologicalType::Cell): + return Metadata::Cell; + case (TopologicalType::Face): + return Metadata::Face; + case (TopologicalType::Edge): + return Metadata::Edge; + case (TopologicalType::Node): + return Metadata::Node; + default: + PARTHENON_FAIL("Unknown topological type") + } + }(); + + std::vector flags{Metadata::FillGhost, ttFlag}; + std::set opts = coarse ? std::set{PDOpt::Coarse} : std::set{}; + const auto desc = + MakePackDescriptor(mbd.get(), flags, opts); + + auto q = desc.GetPack(mbd.get()); + const int b = 0; + const int lstart = q.GetLowerBoundHost(b); + const int lend = q.GetUpperBoundHost(b); + // early return if we don't get anything + if (lend < lstart) return; + auto nb = parthenon::IndexRange{lstart, lend}; + + MeshBlock *pmb = mbd->GetBlockPointer(); + const auto &bounds = (coarse ? pmb->c_cellbounds : pmb->cellbounds); + + const auto &range = bounds.GetBoundsI(IndexDomain::interior, el); + const int ref = range.s; + + constexpr IndexDomain domain = IndexDomain::inner_x1; + // used for reflections + const int offset = 2 * ref - 1; + + pmb->par_for_bndry( + PARTHENON_AUTO_LABEL, nb, domain, el, coarse, false, + KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + const auto vec_component = q(b, el, l).vector_component; + const Real factor = (vec_component == 1 || vec_component == 3) ? -1.0 : 1.0; + q(b, el, l, k, j, i) = factor * q(b, el, l, k, j, offset - i); + }); +} + +parthenon::BValFunc GenericBC(BoundaryConditionInterface bndry_func) { + return [=](std::shared_ptr &mbd, bool coarse) { + for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) + bndry_func(el, mbd, coarse); + }; +} + +void RegisterBoundaryConditions(parthenon::ApplicationInput *app) { + app->RegisterBoundaryCondition(parthenon::BoundaryFace::inner_x1, "axisymmetric", + GenericBC(Axisymmetric)); +} +} // namespace kamayan::grid diff --git a/src/grid/grid.hpp b/src/grid/grid.hpp index 3f4f96b..e1f1a88 100644 --- a/src/grid/grid.hpp +++ b/src/grid/grid.hpp @@ -21,6 +21,8 @@ std::shared_ptr ProcessUnit(); void SetupParams(KamayanUnit *unit); void InitializeData(KamayanUnit *unit); +void RegisterBoundaryConditions(parthenon::ApplicationInput *app); + template requires(std::is_same_v || std::is_same_v) auto GetPackDescriptor(Container *md, std::vector m = {}, diff --git a/src/kamayan/fields.hpp b/src/kamayan/fields.hpp index 0981731..a8f2b8f 100644 --- a/src/kamayan/fields.hpp +++ b/src/kamayan/fields.hpp @@ -18,7 +18,10 @@ namespace kamayan { using MetadataFlag = parthenon::MetadataFlag; using Metadata = parthenon::Metadata; -template +enum class VariableRank { scalar, vector }; + +template struct VariableBase : public parthenon::variable_names::base_t { template KOKKOS_INLINE_FUNCTION VariableBase(Ts &&...args) @@ -34,12 +37,14 @@ struct VariableBase : public parthenon::variable_names::base_t static constexpr std::size_t n_comps = (1 * ... * NCOMP); static constexpr std::size_t n_dims = sizeof...(NCOMP); static constexpr auto vname = var_name; + static constexpr VariableRank variable_rank = vr; }; template concept DenseVar = requires { { T::n_comps } -> std::same_as; // number of components { T::Shape() } -> std::same_as>; // shape of the array + { T::variable_rank } -> std::same_as; requires std::same_as().idx), const int>; requires NonTypeTemplateSpecialization; }; @@ -92,6 +97,7 @@ requires(DenseVar) void AddField(parthenon::StateDescriptor *pkg, std::vector m, std::vector shape = T::Shape()) { // can also add refinement ops here depending on the metadata + if (T::variable_rank == VariableRank::vector) m.push_back(Metadata::Vector); pkg->AddField(Metadata(m, shape)); } @@ -152,19 +158,19 @@ constexpr int count_components(T) { // conserved variables using DENS = VariableBase<"dens">; // will register with shape=std::vector{3} -using MOMENTUM = VariableBase<"momentum", 3>; +using MOMENTUM = VariableBase<"momentum", VariableRank::vector, 3>; using ENER = VariableBase<"ener">; using MAG = VariableBase<"mag">; // --8<-- [end:cons] // primitives & Eos should be FillGhost? -using MAGC = VariableBase<"magc", 3>; +using MAGC = VariableBase<"magc", VariableRank::vector, 3>; using EINT = VariableBase<"eint">; using PRES = VariableBase<"pres">; using BMOD = VariableBase<"bulk modulus">; using TEMP = VariableBase<"temp">; -using VELOCITY = VariableBase<"velocity", 3>; +using VELOCITY = VariableBase<"velocity", VariableRank::vector, 3>; // 3T using TELE = VariableBase<"tele">; diff --git a/src/problems/sedov.py b/src/problems/sedov.py index 28b4411..4c14473 100644 --- a/src/problems/sedov.py +++ b/src/problems/sedov.py @@ -16,7 +16,7 @@ from kamayan.cli import kamayan_app from kamayan.code_units import Grid as gr, eos as eos, driver -from kamayan.code_units.Grid import AdaptiveGrid +from kamayan.code_units.Grid import GEOMETRY, AdaptiveGrid from kamayan.code_units.Hydro import Hydro @@ -118,7 +118,7 @@ def initialize(unit: pyKamayan.KamayanUnit): # --8<-- [start:py_sedov] @kamayan_app(description="Sedov blast wave simulation") def sedov( - nxb: int = typer.Option(32, help="cells across block dimension"), + geometry: GEOMETRY = typer.Option(default="cartesian", help="""Geometry setup."""), ) -> KamayanManager: """Build the KamayanManager for Sedov.""" units = kman.process_units( @@ -126,6 +126,7 @@ def sedov( ) km = KamayanManager("sedov", units) + nxb = 32 nblocks = int(128 / nxb) # number of blocks to get 128 zones at coarsest resolution km.grid = AdaptiveGrid( xbnd1=(-0.5, 0.5), # xmin/max From 72a60c987be598bfe8618f00c3876eec845c51a2 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 22 Mar 2026 16:46:09 -0700 Subject: [PATCH 08/55] fix it --- src/kamayan/fields.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/kamayan/fields.hpp b/src/kamayan/fields.hpp index a8f2b8f..42e77e2 100644 --- a/src/kamayan/fields.hpp +++ b/src/kamayan/fields.hpp @@ -44,16 +44,17 @@ template concept DenseVar = requires { { T::n_comps } -> std::same_as; // number of components { T::Shape() } -> std::same_as>; // shape of the array - { T::variable_rank } -> std::same_as; + { T::variable_rank } -> std::same_as; requires std::same_as().idx), const int>; requires NonTypeTemplateSpecialization; }; template -struct SparseBase : public VariableBase { +struct SparseBase : public VariableBase { template KOKKOS_INLINE_FUNCTION SparseBase(Ts &&...args) - : VariableBase(std::forward(args)...) {} + : VariableBase( + std::forward(args)...) {} // tensor indexer for sparse variable at sparse index n, with tensor // index (ds...) From e0dade1a0fbb18e0c7a130afedf6fce8c0d4be16 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 22 Mar 2026 21:32:10 -0700 Subject: [PATCH 09/55] cylindrical looks right --- python/kamayan/code_units/Grid.py | 10 +++++++++- src/grid/geometry.hpp | 11 +++++++++++ src/grid/grid.cpp | 6 +++--- src/grid/grid_update.hpp | 13 +++++++------ src/kamayan/kamayan.cpp | 3 +++ src/physics/hydro/primconsflux.hpp | 5 +++-- src/problems/sedov.py | 17 ++++++++++++++--- 7 files changed, 50 insertions(+), 15 deletions(-) diff --git a/python/kamayan/code_units/Grid.py b/python/kamayan/code_units/Grid.py index 875b07f..574cbc0 100644 --- a/python/kamayan/code_units/Grid.py +++ b/python/kamayan/code_units/Grid.py @@ -57,7 +57,7 @@ def _set_mesh_block(params: KamayanParams, nxb1: int, nxb2: int, nxb3: int) -> N } -_boundary_conditions = Literal["outflow", "periodic", "user", "reflect"] +_boundary_conditions = Literal["outflow", "periodic", "user", "reflect", "axisymmetric"] @dataclass @@ -126,6 +126,7 @@ class KamayanGrid(Node): xbnd2: tuple[float, float] | None = None xbnd3: tuple[float, float] | None = None numlevel: int = 1 + geometry: GEOMETRY = "cartesian" boundary_conditions = auto_property_node(BoundaryConditions, "boundary_conditions") @@ -145,6 +146,7 @@ def set_params(self, params) -> None: ) _set_mesh_block(params, self.nxb1, self.nxb2, self.nxb3) + params["geometry"] = {"geometry": self.geometry} class UniformGrid(KamayanGrid): @@ -161,6 +163,7 @@ def __init__( dx1: _dx = None, dx2: _dx = None, dx3: _dx = None, + **kwargs, ) -> None: """Set parameters for a uniform grid. @@ -177,6 +180,7 @@ def __init__( dx1: zone size along dimensions of the domain. dx2: zone size along dimensions of the domain. dx3: zone size along dimensions of the domain. + kwargs: kwarguments to init KamayanGrid. Raises: ValueError: When a dimension is over-constrained @@ -225,6 +229,7 @@ def __init__( xbnd1=xbnd1, xbnd2=xbnd2, xbnd3=xbnd3, + **kwargs, ) @@ -334,6 +339,7 @@ def __init__( dx2: _dx = None, dx3: _dx = None, num_levels: int = 1, + **kwargs, ) -> None: """Set parameters for a uniform grid. @@ -354,6 +360,7 @@ def __init__( dx2: zone size along dimensions of the domain at highest resoltution. dx3: zone size along dimensions of the domain at highest resoltution. num_levels: number of total amr refinement levels + kwargs: kwarguments to init KamayanGrid. Raises: ValueError: When a dimension is over-constrained @@ -395,6 +402,7 @@ def __init__( xbnd1=xbnd1, xbnd2=xbnd2, xbnd3=xbnd3, + **kwargs, ) self.refinement_fields = RefinementCollection() diff --git a/src/grid/geometry.hpp b/src/grid/geometry.hpp index aa3874b..faf25b2 100644 --- a/src/grid/geometry.hpp +++ b/src/grid/geometry.hpp @@ -61,6 +61,12 @@ struct Coordinates { return coords_.Xc(idx); } + // position at cell centers + template + KOKKOS_FORCEINLINE_FUNCTION Real Xi(const int idx) const { + return coords_.Xc(idx); + } + KOKKOS_FORCEINLINE_FUNCTION Real Xc(const Axis &ax, const int idx) const { return AxisOverload([&, this]() { return Xc(idx); }, ax); } @@ -198,6 +204,11 @@ struct CoordinateIndexer { return coords_.template Xc(kji_[static_cast(ax)]); } + template + KOKKOS_FORCEINLINE_FUNCTION Real Xi() const { + return coords_.template Xi(kji_[static_cast(ax)]); + } + KOKKOS_FORCEINLINE_FUNCTION Real Xc(const Axis &ax) const { return coords_.Xc(ax, kji_[static_cast(ax)]); } diff --git a/src/grid/grid.cpp b/src/grid/grid.cpp index a632841..d1b193c 100644 --- a/src/grid/grid.cpp +++ b/src/grid/grid.cpp @@ -56,9 +56,9 @@ void SetupParams(KamayanUnit *unit) { parthenon_mesh.AddParm("x2max", 1.0, "Maximum x2 value of domain."); parthenon_mesh.AddParm("x3max", 1.0, "Maximum x3 value of domain."); - parthenon_mesh.AddParm("ix1_bc", "outflow", - "Inner boundary condition along x1.", - {"periodic", "outflow", "reflect", "user"}); + parthenon_mesh.AddParm( + "ix1_bc", "outflow", "Inner boundary condition along x1.", + {"periodic", "outflow", "reflect", "user", "axisymmetric"}); parthenon_mesh.AddParm("ix2_bc", "outflow", "Inner boundary condition along x2.", {"periodic", "outflow", "reflect", "user"}); diff --git a/src/grid/grid_update.hpp b/src/grid/grid_update.hpp index f40fb77..0df820f 100644 --- a/src/grid/grid_update.hpp +++ b/src/grid/grid_update.hpp @@ -23,17 +23,16 @@ void FluxDivergence(MeshData *md, MeshData *dudt_data) { auto ib = md->GetBoundsI(IndexDomain::interior); auto jb = md->GetBoundsJ(IndexDomain::interior); auto kb = md->GetBoundsK(IndexDomain::interior); + // could also be done as par_for_outer(b,k,j) par_for_inner(var,i) par_for( PARTHENON_AUTO_LABEL, 0, nblocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int km, const int jm, const int im) { using TE = TopologicalElement; const auto coords = Coordinates(u0, b); - const Real dxi[sizeof...(faces)]{1.0 / - coords.template Dx()...}; // we have to check the variable bounds for each block in case // that we are using sparse fields for (int var = u0.GetLowerBound(b); var <= u0.GetUpperBound(b); var++) { - dudt(b, var, km, jm, im) = 0.; + Real du = 0.; ( [&]() { constexpr int dir = static_cast(faces) % 3; @@ -41,11 +40,13 @@ void FluxDivergence(MeshData *md, MeshData *dudt_data) { const int jp = jm + (dir == 1); const int ip = im + (dir == 0); - dudt(b, var, km, jm, im) -= - dxi[dir] * (u0.flux(b, faces, var, kp, jp, ip) - - u0.flux(b, faces, var, km, jm, im)); + du -= (coords.template FaceArea(kp, jp, ip) * + u0.flux(b, faces, var, kp, jp, ip) - + coords.template FaceArea(km, jm, im) * + u0.flux(b, faces, var, km, jm, im)); }(), ...); + dudt(b, var, km, jm, im) = du / coords.CellVolume(km, jm, im); } }); } diff --git a/src/kamayan/kamayan.cpp b/src/kamayan/kamayan.cpp index 98916ba..9c8eb96 100644 --- a/src/kamayan/kamayan.cpp +++ b/src/kamayan/kamayan.cpp @@ -6,6 +6,7 @@ #include "driver/kamayan_driver.hpp" #include "driver/kamayan_driver_types.hpp" +#include "grid/grid.hpp" #include "kamayan/unit.hpp" #include "parthenon_manager.hpp" #include "physics/material_properties/eos/eos.hpp" @@ -39,6 +40,8 @@ KamayanDriver InitPackages(std::shared_ptr pman, auto runtime_parameters = std::make_shared(pin); + grid::RegisterBoundaryConditions(pman->app_input.get()); + auto config = std::make_shared(); for (auto &kamayan_unit : *units) { kamayan_unit.second->SetUnits(units); diff --git a/src/physics/hydro/primconsflux.hpp b/src/physics/hydro/primconsflux.hpp index a38e8c9..6a94574 100644 --- a/src/physics/hydro/primconsflux.hpp +++ b/src/physics/hydro/primconsflux.hpp @@ -146,8 +146,9 @@ template (); + // ideally we should use the face values and do a better job integrating... + dudt(MOMENTUM(0)) += (V(DENS()) * V(VELOCITY(2)) * V(VELOCITY(2)) + V(PRES())) / + coords.template Xi(); if constexpr (mhd != Mhd::off) { const Real pmag = 0.5 * (V(MAGC(0)) * V(MAGC(0)) + V(MAGC(1)) * V(MAGC(1)) + V(MAGC(2)) * V(MAGC(2))); diff --git a/src/problems/sedov.py b/src/problems/sedov.py index 4c14473..fb5d38b 100644 --- a/src/problems/sedov.py +++ b/src/problems/sedov.py @@ -2,6 +2,7 @@ """Setup for the Sedov blast wave.""" from dataclasses import dataclass +from typing import Literal import numpy as np from numpy.typing import NDArray @@ -115,10 +116,16 @@ def initialize(unit: pyKamayan.KamayanUnit): # --8<-- [end:py_set_param] +DIMENSION_ = Literal[1, 2, 3] + + # --8<-- [start:py_sedov] @kamayan_app(description="Sedov blast wave simulation") def sedov( geometry: GEOMETRY = typer.Option(default="cartesian", help="""Geometry setup."""), + dimension: DIMENSION_ = typer.Option( + default=2, help="""Number of spatial dimensions""" + ), ) -> KamayanManager: """Build the KamayanManager for Sedov.""" units = kman.process_units( @@ -128,17 +135,21 @@ def sedov( nxb = 32 nblocks = int(128 / nxb) # number of blocks to get 128 zones at coarsest resolution + xbnd1 = (-0.5, 0.5) if geometry == "cartesian" else (0.0, 1.0) km.grid = AdaptiveGrid( - xbnd1=(-0.5, 0.5), # xmin/max + xbnd1=xbnd1, # xmin/max xbnd2=(-0.5, 0.5), # ymin/max nxb1=nxb, # zones per block along x - nxb2=nxb, + nxb2=nxb if dimension > 1 else 1, num_levels=3, # 3 levels of refinement nblocks1=nblocks, # number of root blocks in each direction - nblocks2=nblocks, + nblocks2=nblocks if dimension > 1 else 1, + geometry=geometry, ) km.grid.refinement_fields.add("pres") km.grid.boundary_conditions = gr.outflow_box() + if geometry == "cylindrical": + km.grid.boundary_conditions.ix1 = "axisymmetric" km.driver = driver.Driver(integrator="rk2", tlim=0.05) km.outputs.add("restarts", "rst", dt=0.01) From bcd061fb83a43bac58f3cab1ea8164fb7ffd5b6e Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 22 Mar 2026 22:56:56 -0700 Subject: [PATCH 10/55] another bug --- src/grid/geometry.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/grid/geometry.hpp b/src/grid/geometry.hpp index faf25b2..18c89f9 100644 --- a/src/grid/geometry.hpp +++ b/src/grid/geometry.hpp @@ -120,7 +120,7 @@ struct Coordinates { if constexpr (geom == Geometry::cylindrical) { const Real rp = Xf(i + 1); const Real rm = Xf(i); - return 0.5 * Dx() * (rp * rp - rm * rm); + return 0.5 * Dx() * (rp * rp - rm * rm) * Dx(); } return coords_.CellVolume(k, j, i); } From 16eff5edc72f1c4709793d012b2ed5229c80eda7 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 22 Mar 2026 23:18:18 -0700 Subject: [PATCH 11/55] format --- python/kamayan/cli/app.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/kamayan/cli/app.py b/python/kamayan/cli/app.py index 14f5bd2..fb454e3 100644 --- a/python/kamayan/cli/app.py +++ b/python/kamayan/cli/app.py @@ -65,6 +65,7 @@ def run( km = self.func(*args, **kwargs) if debug: import os + typer.echo(f"PID: {os.getpid()}") breakpoint() if input_file: From af92c33e877077aaa57ec262495e748ec10124ef Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 22 Mar 2026 23:21:30 -0700 Subject: [PATCH 12/55] good job codex --- src/grid/geometry.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/grid/geometry.hpp b/src/grid/geometry.hpp index 18c89f9..10c3ea3 100644 --- a/src/grid/geometry.hpp +++ b/src/grid/geometry.hpp @@ -266,7 +266,7 @@ struct CoordinatesVariant> { type coordinates = Coordinates(coords); ( [&]() { - if (geoms == geometry) coordinates = Coordinates(coords); + if (geoms == geometry) coordinates = Coordinates(coords); }(), ...); From 268503964c042d6f191012e167942834352af64b Mon Sep 17 00:00:00 2001 From: adam reyes Date: Tue, 24 Mar 2026 21:21:58 -0700 Subject: [PATCH 13/55] add refinement operations for geometry --- src/grid/refinement_operations.hpp | 367 +++++++++++++++++++++++++++++ 1 file changed, 367 insertions(+) create mode 100644 src/grid/refinement_operations.hpp diff --git a/src/grid/refinement_operations.hpp b/src/grid/refinement_operations.hpp new file mode 100644 index 0000000..c48a0d0 --- /dev/null +++ b/src/grid/refinement_operations.hpp @@ -0,0 +1,367 @@ +#ifndef GRID_REFINEMENT_OPERATIONS_HPP_ +#define GRID_REFINEMENT_OPERATIONS_HPP_ +#include + +#include + +#include "basic_types.hpp" +#include "coordinates/coordinates.hpp" +#include "grid/geometry.hpp" +#include "grid/grid_types.hpp" +#include "interface/variable_state.hpp" +namespace kamayan::grid { + +namespace util { +// compute distances from cell center to the nearest center in the + or - +// coordinate direction. Do so for both coarse and fine grids. +template +KOKKOS_FORCEINLINE_FUNCTION void +GetGridSpacings(const parthenon::Coordinates_t &coords, + const parthenon::Coordinates_t &coarse_coords, + const parthenon::IndexRange &cib, const parthenon::IndexRange &ib, int i, + int fi, Real *dxm, Real *dxp, Real *dxfm, Real *dxfp) { + // here "f" signifies the fine grid, not face locations. + const Real xm = coarse_coords.X(i - 1); + const Real xc = coarse_coords.X(i); + const Real xp = coarse_coords.X(i + 1); + *dxm = xc - xm; + *dxp = xp - xc; + const Real fxm = coords.X(fi); + const Real fxp = coords.X(fi + 1); + *dxfm = xc - fxm; + *dxfp = fxp - xc; +} + +KOKKOS_FORCEINLINE_FUNCTION +Real GradMinMod(const Real fc, const Real fm, const Real fp, const Real dxm, + const Real dxp, Real &gxm, Real &gxp) { + gxm = (fc - fm) / dxm; + gxp = (fp - fc) / dxp; + return 0.5 * (SIGN(gxm) + SIGN(gxp)) * std::min(std::abs(gxm), std::abs(gxp)); +} + +} // namespace util + +template +struct RestrictAverage { + static constexpr bool OperationRequired(TopologicalElement fel, + TopologicalElement cel) { + return fel == cel; + } + + template + KOKKOS_FORCEINLINE_FUNCTION static void + Do(const int l, const int m, const int n, const int ck, const int cj, const int ci, + const parthenon::IndexRange &ckb, const parthenon::IndexRange &cjb, + const parthenon::IndexRange &cib, const parthenon::IndexRange &kb, + const parthenon::IndexRange &jb, const parthenon::IndexRange &ib, + const parthenon::Coordinates_t &pcoords, + const parthenon::Coordinates_t &pcoarse_coords, + const parthenon::ParArrayND *pcoarse, + const parthenon::ParArrayND *pfine) { + using TE = TopologicalElement; + + auto coords = Coordinates(pcoords); + auto coarse_coords = Coordinates(pcoarse_coords); + + constexpr bool INCLUDE_X1 = + (DIM > 0) && (el == TE::CC || el == TE::F2 || el == TE::F3 || el == TE::E1); + constexpr bool INCLUDE_X2 = + (DIM > 1) && (el == TE::CC || el == TE::F3 || el == TE::F1 || el == TE::E2); + constexpr bool INCLUDE_X3 = + (DIM > 2) && (el == TE::CC || el == TE::F1 || el == TE::F2 || el == TE::E3); + constexpr int element_idx = static_cast(el) % 3; + + auto &coarse = *pcoarse; + auto &fine = *pfine; + + const int i = (DIM > 0) ? (ci - cib.s) * 2 + ib.s : ib.s; + const int j = (DIM > 1) ? (cj - cjb.s) * 2 + jb.s : jb.s; + const int k = (DIM > 2) ? (ck - ckb.s) * 2 + kb.s : kb.s; + + // JMM: If dimensionality is wrong, accesses are out of bounds. Only + // access cells if dimensionality is correct. + Real vol[2][2][2], terms[2][2][2]; // memset not available on all accelerators + for (int ok = 0; ok < 2; ++ok) { + for (int oj = 0; oj < 2; ++oj) { + for (int oi = 0; oi < 2; ++oi) { + vol[ok][oj][oi] = terms[ok][oj][oi] = 0; + } + } + } + + for (int ok = 0; ok < 1 + INCLUDE_X3; ++ok) { + for (int oj = 0; oj < 1 + INCLUDE_X2; ++oj) { + for (int oi = 0; oi < 1 + INCLUDE_X1; ++oi) { + vol[ok][oj][oi] = coords.template Volume(k + ok, j + oj, i + oi); + terms[ok][oj][oi] = + vol[ok][oj][oi] * fine(element_idx, l, m, n, k + ok, j + oj, i + oi); + } + } + } + // KGF: add the off-centered quantities first to preserve FP + // symmetry + const Real tvol = ((vol[0][0][0] + vol[0][1][0]) + (vol[0][0][1] + vol[0][1][1])) + + ((vol[1][0][0] + vol[1][1][0]) + (vol[1][0][1] + vol[1][1][1])); + coarse(element_idx, l, m, n, ck, cj, ci) = + tvol > 0.0 + ? (((terms[0][0][0] + terms[0][1][0]) + (terms[0][0][1] + terms[0][1][1])) + + ((terms[1][0][0] + terms[1][1][0]) + (terms[1][0][1] + terms[1][1][1]))) / + tvol + : 0.0; + } +}; + +template +struct ProlongateSharedGeneral { + static constexpr bool OperationRequired(TopologicalElement fel, + TopologicalElement cel) { + return fel == cel; + } + + template + KOKKOS_FORCEINLINE_FUNCTION static void + Do(const int l, const int m, const int n, const int k, const int j, const int i, + const parthenon::IndexRange &ckb, const parthenon::IndexRange &cjb, + const parthenon::IndexRange &cib, const parthenon::IndexRange &kb, + const parthenon::IndexRange &jb, const parthenon::IndexRange &ib, + const parthenon::Coordinates_t &pcoords, + const parthenon::Coordinates_t &pcoarse_coords, + const parthenon::ParArrayND *pcoarse, + const parthenon::ParArrayND *pfine) { + using util::GradMinMod; + using TE = TopologicalElement; + + auto coords = Coordinates(pcoords); + auto coarse_coords = Coordinates(pcoarse_coords); + + auto &coarse = *pcoarse; + auto &fine = *pfine; + + constexpr int element_idx = static_cast(el) % 3; + + const int fi = (DIM > 0) ? (i - cib.s) * 2 + ib.s : ib.s; + const int fj = (DIM > 1) ? (j - cjb.s) * 2 + jb.s : jb.s; + const int fk = (DIM > 2) ? (k - ckb.s) * 2 + kb.s : kb.s; + + constexpr bool INCLUDE_X1 = + (DIM > 0) && (el == TE::CC || el == TE::F2 || el == TE::F3 || el == TE::E1); + constexpr bool INCLUDE_X2 = + (DIM > 1) && (el == TE::CC || el == TE::F3 || el == TE::F1 || el == TE::E2); + constexpr bool INCLUDE_X3 = + (DIM > 2) && (el == TE::CC || el == TE::F1 || el == TE::F2 || el == TE::E3); + + const Real fc = coarse(element_idx, l, m, n, k, j, i); + + Real dx1fm = 0; + [[maybe_unused]] Real dx1fp = 0; + [[maybe_unused]] Real gx1m = 0, gx1p = 0; + if constexpr (INCLUDE_X1) { + Real dx1m, dx1p; + GetGridSpacings<1, el>(pcoords, pcoarse_coords, cib, ib, i, fi, &dx1m, &dx1p, + &dx1fm, &dx1fp); + + Real gx1c = + GradMinMod(fc, coarse(element_idx, l, m, n, k, j, i - 1), + coarse(element_idx, l, m, n, k, j, i + 1), dx1m, dx1p, gx1m, gx1p); + if constexpr (use_minmod_slope) { + gx1m = gx1c; + gx1p = gx1c; + } + } + + Real dx2fm = 0; + [[maybe_unused]] Real dx2fp = 0; + [[maybe_unused]] Real gx2m = 0, gx2p = 0; + if constexpr (INCLUDE_X2) { + Real dx2m, dx2p; + GetGridSpacings<2, el>(pcoords, pcoarse_coords, cjb, jb, j, fj, &dx2m, &dx2p, + &dx2fm, &dx2fp); + Real gx2c = + GradMinMod(fc, coarse(element_idx, l, m, n, k, j - 1, i), + coarse(element_idx, l, m, n, k, j + 1, i), dx2m, dx2p, gx2m, gx2p); + if constexpr (use_minmod_slope) { + gx2m = gx2c; + gx2p = gx2c; + } + } + + Real dx3fm = 0; + [[maybe_unused]] Real dx3fp = 0; + [[maybe_unused]] Real gx3m = 0, gx3p = 0; + if constexpr (INCLUDE_X3) { + Real dx3m, dx3p; + GetGridSpacings<3, el>(pcoords, pcoarse_coords, ckb, kb, k, fk, &dx3m, &dx3p, + &dx3fm, &dx3fp); + Real gx3c = + GradMinMod(fc, coarse(element_idx, l, m, n, k - 1, j, i), + coarse(element_idx, l, m, n, k + 1, j, i), dx3m, dx3p, gx3m, gx3p); + if constexpr (use_minmod_slope) { + gx3m = gx3c; + gx3p = gx3c; + } + } + + if constexpr (piecewise_constant) { + gx1m = 0.0; + gx1p = 0.0; + gx2m = 0.0; + gx2p = 0.0; + gx3m = 0.0; + gx3p = 0.0; + } + + // KGF: add the off-centered quantities first to preserve FP symmetry + // JMM: Extraneous quantities are zero + fine(element_idx, l, m, n, fk, fj, fi) = + fc - (gx1m * dx1fm + gx2m * dx2fm + gx3m * dx3fm); + if constexpr (INCLUDE_X1) + fine(element_idx, l, m, n, fk, fj, fi + 1) = + fc + (gx1p * dx1fp - gx2m * dx2fm - gx3m * dx3fm); + if constexpr (INCLUDE_X2) + fine(element_idx, l, m, n, fk, fj + 1, fi) = + fc - (gx1m * dx1fm - gx2p * dx2fp + gx3m * dx3fm); + if constexpr (INCLUDE_X2 && INCLUDE_X1) + fine(element_idx, l, m, n, fk, fj + 1, fi + 1) = + fc + (gx1p * dx1fp + gx2p * dx2fp - gx3m * dx3fm); + if constexpr (INCLUDE_X3) + fine(element_idx, l, m, n, fk + 1, fj, fi) = + fc - (gx1m * dx1fm + gx2m * dx2fm - gx3p * dx3fp); + if constexpr (INCLUDE_X3 && INCLUDE_X1) + fine(element_idx, l, m, n, fk + 1, fj, fi + 1) = + fc + (gx1p * dx1fp - gx2m * dx2fm + gx3p * dx3fp); + if constexpr (INCLUDE_X3 && INCLUDE_X2) + fine(element_idx, l, m, n, fk + 1, fj + 1, fi) = + fc - (gx1m * dx1fm - gx2p * dx2fp - gx3p * dx3fp); + if constexpr (INCLUDE_X3 && INCLUDE_X2 && INCLUDE_X1) + fine(element_idx, l, m, n, fk + 1, fj + 1, fi + 1) = + fc + (gx1p * dx1fp + gx2p * dx2fp + gx3p * dx3fp); + } +}; + +template +using ProlongateSharedMinMod = ProlongateSharedGeneral; +template +using ProlongateSharedLinear = ProlongateSharedGeneral; +template +using ProlongatePiecewiseConstant = ProlongateSharedGeneral; + +// Implements divergence-free prolongation to internal faces using the method +// described in Toth & Roe (2002). Any prolongation method for faces shared +// between the coarse and the fine grid can be used alongside this internal +// prolongation operation. Obviously, this prolongation operation is only +// defined for face fields. +// +// We modify it here for cylindrical coordinates where the divergence constraint +// is satisfied for r*B +template +struct ProlongateInternalTothAndRoe { + static constexpr bool OperationRequired(TopologicalElement fel, + TopologicalElement cel) { + using TE = TopologicalElement; + return (cel == TE::CC) && (GetTopologicalType(fel) == TopologicalType::Face); + } + // Here, fel is the topological element on which the field is defined and + // cel is the topological element on which we are filling the internal values + // of the field. So, for instance, we could fill the fine cell values of an + // x-face field within the volume of a coarse cell. This is assumes that the + // values of the fine cells on the elements corresponding with the coarse cell + // have been filled. + template + KOKKOS_FORCEINLINE_FUNCTION static void + Do(const int l, const int m, const int n, const int k, const int j, const int i, + const parthenon::IndexRange &ckb, const parthenon::IndexRange &cjb, + const parthenon::IndexRange &cib, const parthenon::IndexRange &kb, + const parthenon::IndexRange &jb, const parthenon::IndexRange &ib, + const parthenon::Coordinates_t &pcoords, + const parthenon::Coordinates_t &pcoarse_coords, + const parthenon::ParArrayND *, + const parthenon::ParArrayND *pfine) { + using TE = TopologicalElement; + using util::GradMinMod; + + auto coords = Coordinates(pcoords); + auto coarse_coords = Coordinates(pcoarse_coords); + + if constexpr (!IsSubmanifold(fel, cel)) { + return; + } else { + const int fi = (DIM > 0) ? (i - cib.s) * 2 + ib.s : ib.s; + const int fj = (DIM > 1) ? (j - cjb.s) * 2 + jb.s : jb.s; + const int fk = (DIM > 2) ? (k - ckb.s) * 2 + kb.s : kb.s; + + // Here, we write the update for the x-component of the B-field and recover the + // other components by cyclic permutation + constexpr int element_idx = static_cast(fel) % 3; + auto get_fine_permuted_kji = [&](int ok, int oj, int oi) -> std::array { + // Guard against offsetting in symmetry dimensions + constexpr int g3 = (DIM > 2); + constexpr int g2 = (DIM > 1); + if constexpr (fel == TE::F1) { + return {fk + ok * g3, fj + oj * g2, fi + oi}; + } else if constexpr (fel == TE::F2) { + return {fk + oj * g3, fj + oi * g2, fi + ok}; + } else { + return {fk + oi * g3, fj + ok * g2, fi + oj}; + } + }; + auto get_fine_permuted = [&](int eidx, int ok, int oj, int oi) -> Real & { + eidx = (element_idx + eidx) % 3; + const auto [kk, jj, ii] = get_fine_permuted_kji(ok, oj, oi); + return (*pfine)(eidx, l, m, n, kk, jj, ii); + }; + auto get_geometric_factor = [&](int v, int u, int t) { + if constexpr (geom == Geometry::cylindrical) { + return coords.template Xf(fi + t); + } else { + return 1.0; + } + }; + + using iarr2 = std::array; + auto sg = [](int offset) -> Real { return offset == 0 ? -1.0 : 1.0; }; + Real Uxx{0.0}; + Real Vxyz{0.0}; + Real Wxyz{0.0}; + for (const int v : iarr2{0, 1}) { + // Note step size of 2 for the direction normal to the eidx2/eidx3 + for (const int u : iarr2{0, 2}) { + for (const int t : iarr2{0, 1}) { + // we should multiply by the radial coordinate in cylindrical + // coordinates for r-faces + const auto fine2 = + get_fine_permuted(1, v, u, t) * get_geometric_factor(v, u, t); + const auto fine3 = + get_fine_permuted(2, u, v, t) * get_geometric_factor(v, u, t); + Uxx += sg(t) * sg(u) * (fine2 + fine3); + Vxyz += sg(t) * sg(u) * sg(v) * fine2; + Wxyz += sg(t) * sg(u) * sg(v) * fine3; + } + } + } + Uxx *= 0.125; + const int dir1 = element_idx + 1; + const int dir2 = (element_idx + 1) % 3 + 1; + const int dir3 = (element_idx + 2) % 3 + 1; + const auto dx2 = std::pow(coarse_coords.Dxc(dir1, k, j, i), 2); + const auto dy2 = std::pow(coarse_coords.Dxc(dir2, k, j, i), 2); + const auto dz2 = std::pow(coarse_coords.Dxc(dir3, k, j, i), 2); + Vxyz *= 0.125 * dz2 / (dx2 + dz2); + Wxyz *= 0.125 * dy2 / (dx2 + dy2); + + for (int ok : iarr2{0, 1}) { + for (int oj : iarr2{0, 1}) { + get_fine_permuted(0, ok, oj, 1) = + 0.5 * (get_fine_permuted(0, ok, oj, 0) + get_fine_permuted(0, ok, oj, 2)) + + Uxx + sg(ok) * Vxyz + sg(oj) * Wxyz; + get_fine_permuted(0, ok, oj, 1) *= 1. / get_geometric_factor(ok, oj, 1); + } + } + } + } +}; +} // namespace kamayan::grid +#endif // GRID_REFINEMENT_OPERATIONS_HPP_ From 2c4316bcdd46967ecb05f7627e904b770b16f8c3 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Wed, 25 Mar 2026 21:26:54 -0700 Subject: [PATCH 14/55] tweaks --- src/grid/geometry.hpp | 52 +++++++++++++++++++++++++++++- src/grid/grid_types.hpp | 4 +++ src/grid/refinement_operations.hpp | 51 +++++++++++++++++------------ src/kamayan/fields.hpp | 3 +- src/physics/hydro/hydro.cpp | 11 ++++--- 5 files changed, 94 insertions(+), 27 deletions(-) diff --git a/src/grid/geometry.hpp b/src/grid/geometry.hpp index 10c3ea3..6823f31 100644 --- a/src/grid/geometry.hpp +++ b/src/grid/geometry.hpp @@ -33,6 +33,22 @@ struct Coordinates { const int block) : coords_(pack.GetCoordinates(block)) {} + template + KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const { + constexpr int dir = AxisToInt(ax); + if constexpr (ax == Axis::IAXIS && TopologicalOffsetI(el)) { + return coords_.X<1, el>(idx); + } else if constexpr (ax == Axis::JAXIS && TopologicalOffsetJ(el)) { + return coords_.X<2, el>(idx); + } else if constexpr (ax == Axis::JAXIS && TopologicalOffsetK(el)) { + return coords_.X<3, el>(idx); + } else { + return Xc(idx); + } + return 0; // This should never be reached, but w/o it some compilers generate + // warnings + } + // distance between cell centers template KOKKOS_FORCEINLINE_FUNCTION Real Dx() const { @@ -77,10 +93,20 @@ struct Coordinates { return coords_.Xf(idx); } - KOKKOS_FORCEINLINE_FUNCTION Real Xf(const Axis &ax, const int idx) const { + KOKKOS_FORCEINLINE_FUNCTION Real Xf(const Axis ax, const int idx) const { return AxisOverload([&, this]() { return Xf(idx); }, ax); } + // coordinate ax at center of element el + template + KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int idx) const { + if constexpr (static_cast(el) % 3 + 1 == AxisToInt(ax)) { + return Xf(idx); + } else { + return Xi(idx); + } + } + template KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int k, const int j, const int i) const { constexpr auto dir = AxisToInt(ax); @@ -125,6 +151,30 @@ struct Coordinates { return coords_.CellVolume(k, j, i); } + template + KOKKOS_FORCEINLINE_FUNCTION Real Volume(const int k, const int j, const int i) const { + using TE = TopologicalElement; + if constexpr (el == TE::CC) + return CellVolume(k, j, i); + else if constexpr (el == TE::F1) + return FaceArea(k, j, i); + else if constexpr (el == TE::F2) + return FaceArea(k, j, i); + else if constexpr (el == TE::F3) + return FaceArea(k, j, i); + else if constexpr (el == TE::E1) + return EdgeLength(k, j, i); + else if constexpr (el == TE::E2) + return EdgeLength(k, j, i); + else if constexpr (el == TE::E3) + return EdgeLength(k, j, i); + else if constexpr (el == TE::NN) + return 1.0; + PARTHENON_FAIL("ifyou reach this point, someone has added a new value to the the " + "TopologicalElement enum."); + return 0.0; + } + KOKKOS_FORCEINLINE_FUNCTION Real Volume(TopologicalElement el, const int k, const int j, const int i) const { using TE = TopologicalElement; diff --git a/src/grid/grid_types.hpp b/src/grid/grid_types.hpp index e2c0c1b..94555a9 100644 --- a/src/grid/grid_types.hpp +++ b/src/grid/grid_types.hpp @@ -26,6 +26,10 @@ constexpr Axis AxisFromTE(const TopologicalElement el) { return dir == 0 ? Axis::IAXIS : dir == 1 ? Axis::JAXIS : Axis::KAXIS; } +constexpr Axis AxisFromInt(const int dir) { + return dir == 1 ? Axis::IAXIS : dir == 2 ? Axis::JAXIS : Axis::KAXIS; +} + template concept EdgeElement = (edge >= TopologicalElement::E1 && edge <= TopologicalElement::E3); diff --git a/src/grid/refinement_operations.hpp b/src/grid/refinement_operations.hpp index c48a0d0..b24cb2d 100644 --- a/src/grid/refinement_operations.hpp +++ b/src/grid/refinement_operations.hpp @@ -10,24 +10,29 @@ #include "grid/grid_types.hpp" #include "interface/variable_state.hpp" namespace kamayan::grid { - +// Parthenon's implementation of non-cartesian coordinates assumes that +// it is determined at compile time for the entire library. In order +// to enable runtime geometry we re-implement the refinement operations upstream +// templated on the geometry. +// namespace util { -// compute distances from cell center to the nearest center in the + or - -// coordinate direction. Do so for both coarse and fine grids. -template + +// distances between element centroids, so that prolongation is always +// done in the volume coordinate +template KOKKOS_FORCEINLINE_FUNCTION void -GetGridSpacings(const parthenon::Coordinates_t &coords, - const parthenon::Coordinates_t &coarse_coords, +GetGridSpacings(const Coordinates &coords, const Coordinates &coarse_coords, const parthenon::IndexRange &cib, const parthenon::IndexRange &ib, int i, int fi, Real *dxm, Real *dxp, Real *dxfm, Real *dxfp) { // here "f" signifies the fine grid, not face locations. - const Real xm = coarse_coords.X(i - 1); - const Real xc = coarse_coords.X(i); - const Real xp = coarse_coords.X(i + 1); + constexpr auto ax = AxisFromInt(DIM); + const Real xm = coarse_coords.template X(i - 1); + const Real xc = coarse_coords.template X(i); + const Real xp = coarse_coords.template X(i + 1); *dxm = xc - xm; *dxp = xp - xc; - const Real fxm = coords.X(fi); - const Real fxp = coords.X(fi + 1); + const Real fxm = coords.template X(fi); + const Real fxp = coords.template X(fi + 1); *dxfm = xc - fxm; *dxfp = fxp - xc; } @@ -131,6 +136,7 @@ struct ProlongateSharedGeneral { const parthenon::Coordinates_t &pcoarse_coords, const parthenon::ParArrayND *pcoarse, const parthenon::ParArrayND *pfine) { + using util::GetGridSpacings; using util::GradMinMod; using TE = TopologicalElement; @@ -160,8 +166,8 @@ struct ProlongateSharedGeneral { [[maybe_unused]] Real gx1m = 0, gx1p = 0; if constexpr (INCLUDE_X1) { Real dx1m, dx1p; - GetGridSpacings<1, el>(pcoords, pcoarse_coords, cib, ib, i, fi, &dx1m, &dx1p, - &dx1fm, &dx1fp); + GetGridSpacings<1, el, geom>(coords, coarse_coords, cib, ib, i, fi, &dx1m, &dx1p, + &dx1fm, &dx1fp); Real gx1c = GradMinMod(fc, coarse(element_idx, l, m, n, k, j, i - 1), @@ -177,8 +183,8 @@ struct ProlongateSharedGeneral { [[maybe_unused]] Real gx2m = 0, gx2p = 0; if constexpr (INCLUDE_X2) { Real dx2m, dx2p; - GetGridSpacings<2, el>(pcoords, pcoarse_coords, cjb, jb, j, fj, &dx2m, &dx2p, - &dx2fm, &dx2fp); + GetGridSpacings<2, el, geom>(coords, coarse_coords, cjb, jb, j, fj, &dx2m, &dx2p, + &dx2fm, &dx2fp); Real gx2c = GradMinMod(fc, coarse(element_idx, l, m, n, k, j - 1, i), coarse(element_idx, l, m, n, k, j + 1, i), dx2m, dx2p, gx2m, gx2p); @@ -193,8 +199,8 @@ struct ProlongateSharedGeneral { [[maybe_unused]] Real gx3m = 0, gx3p = 0; if constexpr (INCLUDE_X3) { Real dx3m, dx3p; - GetGridSpacings<3, el>(pcoords, pcoarse_coords, ckb, kb, k, fk, &dx3m, &dx3p, - &dx3fm, &dx3fp); + GetGridSpacings<3, el, geom>(coords, coarse_coords, ckb, kb, k, fk, &dx3m, &dx3p, + &dx3fm, &dx3fp); Real gx3c = GradMinMod(fc, coarse(element_idx, l, m, n, k - 1, j, i), coarse(element_idx, l, m, n, k + 1, j, i), dx3m, dx3p, gx3m, gx3p); @@ -346,9 +352,9 @@ struct ProlongateInternalTothAndRoe { const int dir1 = element_idx + 1; const int dir2 = (element_idx + 1) % 3 + 1; const int dir3 = (element_idx + 2) % 3 + 1; - const auto dx2 = std::pow(coarse_coords.Dxc(dir1, k, j, i), 2); - const auto dy2 = std::pow(coarse_coords.Dxc(dir2, k, j, i), 2); - const auto dz2 = std::pow(coarse_coords.Dxc(dir3, k, j, i), 2); + const auto dx2 = std::pow(coarse_coords.Dx(Axis::IAXIS), 2); + const auto dy2 = std::pow(coarse_coords.Dx(Axis::JAXIS), 2); + const auto dz2 = std::pow(coarse_coords.Dx(Axis::KAXIS), 2); Vxyz *= 0.125 * dz2 / (dx2 + dz2); Wxyz *= 0.125 * dy2 / (dx2 + dy2); @@ -363,5 +369,10 @@ struct ProlongateInternalTothAndRoe { } } }; + +// TODO(acreyes): ProlongateInternalAverage +// In principle I guess we should weight the averaging based on the generalized volume of +// the elements + } // namespace kamayan::grid #endif // GRID_REFINEMENT_OPERATIONS_HPP_ diff --git a/src/kamayan/fields.hpp b/src/kamayan/fields.hpp index 42e77e2..5907356 100644 --- a/src/kamayan/fields.hpp +++ b/src/kamayan/fields.hpp @@ -99,7 +99,8 @@ void AddField(parthenon::StateDescriptor *pkg, std::vector m, std::vector shape = T::Shape()) { // can also add refinement ops here depending on the metadata if (T::variable_rank == VariableRank::vector) m.push_back(Metadata::Vector); - pkg->AddField(Metadata(m, shape)); + auto md = Metadata(m, shape); + pkg->AddField(md); } template diff --git a/src/physics/hydro/hydro.cpp b/src/physics/hydro/hydro.cpp index bebf2bc..da226be 100644 --- a/src/physics/hydro/hydro.cpp +++ b/src/physics/hydro/hydro.cpp @@ -11,6 +11,7 @@ #include "driver/kamayan_driver_types.hpp" #include "grid/geometry.hpp" #include "grid/grid.hpp" +#include "grid/refinement_operations.hpp" #include "grid/subpack.hpp" #include "hydro_types.hpp" #include "kamayan/config.hpp" @@ -82,9 +83,9 @@ void SetupParams(KamayanUnit *unit) { } struct InitializeHydro { - using options = OptTypeList; + using options = OptTypeList; using value = void; - template + template requires(NonTypeTemplateSpecialization) value dispatch(KamayanUnit *unit, Config *cfg) { // --8<-- [start:hydro_add_fields] @@ -97,9 +98,9 @@ struct InitializeHydro { if constexpr (hydro_vars::MHD == Mhd::ct) { auto m = Metadata(std::vector{ FACE_FLAGS(Metadata::Independent, Metadata::WithFluxes)}); - m.RegisterRefinementOps(); + m.RegisterRefinementOps, + grid::RestrictAverage, + grid::ProlongateInternalTothAndRoe>(); unit->AddField(m); } From 7603d7e6879370c40e275422a1c7090e73a756c1 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sat, 28 Mar 2026 18:57:02 -0700 Subject: [PATCH 15/55] refining vars for caching geometry --- src/dispatcher/options.hpp | 30 +++++++++-- src/grid/coordinates.hpp | 100 +++++++++++++++++++++++++++++++++++ src/grid/geometry.hpp | 60 ++++++++++++++++++--- src/grid/grid_update.hpp | 1 + src/grid/tests/test_grid.cpp | 9 ++-- src/grid/tests/test_grid.hpp | 3 +- src/kamayan/fields.hpp | 22 +++++--- src/kamayan/kamayan.cpp | 6 +++ src/kamayan/unit.hpp | 3 ++ src/physics/hydro/hydro.cpp | 16 ++++-- 10 files changed, 224 insertions(+), 26 deletions(-) create mode 100644 src/grid/coordinates.hpp diff --git a/src/dispatcher/options.hpp b/src/dispatcher/options.hpp index 68d2a32..8409f1f 100644 --- a/src/dispatcher/options.hpp +++ b/src/dispatcher/options.hpp @@ -17,19 +17,41 @@ namespace kamayan { -template -struct OptList {}; +// checks that a functor can be instantiated for EnumV +// used to describe functors that can be used with OptList::dispatch below +template +concept OptionFunctor = requires { + requires std::is_void_v().template operator()( + std::declval()...))>; +}; // used to enumerate the allowed values of a PolyOpt in a given // dispatch functor. By default this list is exactly what is written // unless cmake has been configured with -DOPT_enum_opt="enum_v,..." // in that case OptInf::isdef == true and the ParmList() // will just be those defined at configure time -template +template requires(PolyOpt) -struct OptList { +struct OptList { using type = enum_opt; static constexpr auto value = OptInfo::ParmList(); + + // used to dispatch to a runtime template specialization of a type's operator() + template + requires(OptionFunctor && ...) + static bool dispatch(const Functor &functor, const enum_opt parm, Args &&...args) { + bool found_parm = false; + (void)(( + [&]() { + if (parm == enum_vs) { + found_parm = true; + functor.template operator()(std::forward(args)...); + } + }(), + !found_parm) && + ...); + return found_parm; + } }; template diff --git a/src/grid/coordinates.hpp b/src/grid/coordinates.hpp new file mode 100644 index 0000000..b2db85c --- /dev/null +++ b/src/grid/coordinates.hpp @@ -0,0 +1,100 @@ +#ifndef GRID_COORDINATES_HPP_ +#define GRID_COORDINATES_HPP_ +#include +#include + +#include + +#include "grid/grid_types.hpp" +#include "utils/strings.hpp" +#include "utils/type_list.hpp" + +namespace kamayan::grid { + +// Idea for this is to have a variable that is stored per meshblock +// that is associated with a topological element, but might be only allocated for a +// particular kji dimension. For example radial face areas in cylindrical coordinates +// these face areas on a mesh block only depend on the radial index (i) and are +// located on the TE::F1 faces. In the case that there are no axes this is just a scalar +// +// using AreaF1 = CoordVar<"geom.area_f1"> +// +// The raw parthenon way of accessing this would be to pack it, then access it +// with the i-index +// pack(b, AreaF1())(i); +// +// It is not necessarily the case that an arbitrary coordinate system would +// index with i, it might even be the full (k,j,i), so we need some kind of helper +// to forward that +// +// coords.Area(pack, b, k, j, i); +// +// I think we define all the variables as CoordVars with just the var_name, +// we collect them into a TypeList, and then each coordinate system becomes +// responsible for exporting the shapes, so we can add as +// +// const auto shape = coords.shape(nx3, nx2, nx1); // should include ghosts +// pkg->AddField(Metadata({Metadata::None, Metadata::OneCopy}, shape)); +// +// this could also be wrapped in a type_for +template +struct CoordVar : parthenon::variable_names::base_t { + template + KOKKOS_INLINE_FUNCTION CoordVar(Args &&...args) + : parthenon::variable_names::base_t(std::forward(args)...) {} + + static std::string name() { return std::string(var_name.value); } + static constexpr auto label = var_name; +}; + +namespace impl { +template +constexpr auto AxisNamer(const strings::CompileTimeString base, Axis ax) { + if (ax == Axis::IAXIS) { + return strings::concat_cts(base, strings::make_cts("1")); + } else if (ax == Axis::JAXIS) { + return strings::concat_cts(base, strings::make_cts("2")); + } else if (ax == Axis::KAXIS) { + return strings::concat_cts(base, strings::make_cts("3")); + } + return strings::concat_cts(base, strings::make_cts("1")); +} +} // namespace impl + +template +using AxisCoord = CoordVar; + +template +using Dx = AxisCoord<"geom.Dx", ax>; + +template +using X = AxisCoord<"geom.X", ax>; + +template +using Xc = AxisCoord<"geom.Xc", ax>; + +template +using Xf = AxisCoord<"geom.Xf", ax>; + +template +using FaceArea = AxisCoord<"geom.FaceArea", ax>; + +template +using EdgeLength = AxisCoord<"geom.EdgeLength", ax>; + +using Volume = CoordVar<"geom.Volume">; + +template