Skip to content

Grid Unit

Kamayan is built on parthenon's grid, which can support uniform, static and adaptively refined meshes as a runtime option. The domain is initialized as a uniform grid of nx1xnx2xnx3 cells. These can then be partitioned onto an initial grid of MeshBlocks, which can be further refined based on various conditions. See Parameters for more details on the runtime parameters that control this behavior.

Fields

Kamayan strongly prefers to reference field data stored on the MeshBlocks through the type based interface. Fields are associated with a struct, and instantiations of that type can be used to reference a desired element of that field if it is declared with some tensor shape. Kamayan exports all possible fields in kamayan/fields.hpp, but which fields have been registered to the simulation will depend on the configuration of the various units. Here fields are declared with the VARIABLE macro, which takes in the type name in all capital letters, which also serves as the string label associated with the field, and optionally a list of ints that will become the tensor shape of the field when registered with the kamayan::AddField[s] interfaces.

kamayan/fields.hpp:cons
// conserved variables
using DENS = VariableBase<"dens">;
// will register with shape=std::vector<int>{3}
using MOMENTUM = VariableBase<"momentum", VariableRank::vector, 3>;
using ENER = VariableBase<"ener">;
using MAG = VariableBase<"mag">;

Note

KamayanUnits will register fields through the Initialize callback interface that returns a StateDescriptor.

Packs

Field data can be accessed through the SparsePacks provided by parthenon that will pack variables by type into a single data structure that can be used to index into cell level data on MeshBlocks and MeshData partitions. Kamayan provides an interface to directly get these packs, either directly through listing explicitly the desired variables or by passing a TypeList of the variable structs you wish to access.

problems/isentropic_vortex.cpp:pack
physics/hydro/hydro_add_flux_tasks.cpp:pack
// explicitly specify the fields we want
  auto pack = grid::GetPack<DENS, VELOCITY, PRES, MAGC>(mb);
// packing with a decl'ed `TypeList` from the HydroTraits
    auto pack_recon = grid::GetPack(reconstruct_vars(), md);
    // include Xf for cylindrical geometry flux corrections
    using flux_vars = ConcatTypeLists_t<conserved_vars, grid::Xface>;
    auto pack_flux = grid::GetPack(flux_vars(), md, {PDOpt::WithFluxes});

Note

Additional pack options can be specified. In this example the PDOpt::WithFluxes option means that the associated fluxes if any for the packed variables will also be accessible with the pack_flux.flux interface.

Packs can then be indexed by block, variable and (k,j,i) index

problems/isentropic_vortex.cpp:index
        pack(0, DENS(), k, j, i) = state(DENS());
        pack(0, PRES(), k, j, i) = state(PRES());
        pack(0, VELOCITY(0), k, j, i) = state(VELOCITY(0));
        pack(0, VELOCITY(1), k, j, i) = state(VELOCITY(1));
        pack(0, VELOCITY(2), k, j, i) = state(VELOCITY(2));

Subpacks

Often we don't need full access to all block, and cell indices of the pack. Rather it is advantageous only worry about a relative offset along the dimensions of interest for a particular function. Some examples are the conversion of primitive to conservative variables, which only cares about indexing into different fields at the same cell, or the reconstruction of Riemann states, which need only to index into the same field at neighboring zones along a one dimensional stencil.

Kamayan provides the SubPack abstractions to wrap packs and Kokkos::Views to do just that.

physics/hydro/primconsflux.cpp:make-idx
          auto U = SubPack(pack, b, k, j, i);
          Prim2Cons<hydro_traits>(U, U);
physics/hydro/primconsflux.hpp:use-idx
  U(DENS()) = V(DENS());
  Real emag = 0.;
  Real ekin = 0.;
  for (int dir = 0; dir < 3; dir++) {
    U(MOMENTUM(dir)) = V(DENS()) * V(VELOCITY(dir));
    ekin += V(VELOCITY(dir)) * V(VELOCITY(dir));
    if constexpr (hydro_traits::MHD != Mhd::off) {
      U(MAGC(dir)) = V(MAGC(dir));
      emag += V(MAGC(dir)) * V(MAGC(dir));
    }
  }

physics/hydro/hydro_add_flux_tasks.cpp:make-stncl
                auto stencil = SubPack<Axis::IAXIS>(pack_recon, b, var, k, j, i);
                Reconstruct<reconstruction_traits>(stencil, vM(var, i), vP(var, i));
physics/hydro/reconstruction.hpp:use-stncl
  const Real dvL = stencil(0) - stencil(-1);
  const Real dvR = stencil(1) - stencil(0);

Scratch Variables

Kamayan provides the ability to represent fields on the mesh that are temporary in nature using collections called ScratchVariableLists. These provide types that can be registered to a package that will alias common Metadata::derived fields across different regions where they are used. In this way temporary block storage is achieved without needing to allocate more memory than is minimally required to satisfy a single ScratchVariableList.

ScratchVarableLists are constructed from specializations of the ScratchVariable type that is used to describe the variable's TopologicalType and vector/tensor shape, and finally are used to index into the ScratchVariableList in order to get a type that can be used exactly as any other field through the type-based packing.

grid/grid_refinement.hpp:scratch
// first-order derivative at centers used in Loehner estimator
using FirstDer = RuntimeScratchVariable<"firstder", TopologicalType::Cell, 3>;

using RefinementScratch = RuntimeScratchVariableList<FirstDer>;
grid/grid.cpp:addscratch

In the above example a single scratch variable to represent the first-order derivatives in each 3D direction is registered as a cell centered variable as a 3 component vector. When it is registered it will define three fields, scratch_cell_n, that can be reused by other scratch variables, possibly of different shape. The alias FirstDer is pulled out of the ScratchVariableList can then be used with the pack overloads just like any other field.

grid/grid_refinement.cpp:FirstDer
        pack_der(b, 0, k, j, i) =
            0.5 * (pack(b, var, k, j, i + 1) - pack(b, var, k, j, i - 1)) /
            coords.Dxc<1>();

Note

If cmake is configured with -Dkamayan_DEBUG_SCRATCH then each scratch variable will be independently registered using the name string template parameter to the ScratchVariable. In the above example a shape {3} field scratch_firstder will be registered.

Coordinates

Kamayan builds on Parthenon's mesh and coordinate layout (currently parthenon::UniformCartesian), but wraps it behind a small geometry layer so that code can be written against a consistent coordinate/metric API while the geometry is selected at runtime.

At the moment the grid unit supports:

  • Geometry::cartesian: standard 1D/2D/3D Cartesian
  • Geometry::cylindrical: 2D axisymmetric r-z with an implicit azimuthal direction (the third direction uses dphi = 2*pi)

CoordinateSystem

The coordinate wrappers implement the CoordinateSystem concept (src/grid/geometry.hpp), which is the practical "contract" used throughout the code. It includes both compile-time and runtime axis dispatch:

  • Compile-time: template <Axis ax> Dx(), Xc(idx), Xf(idx), FaceArea(k,j,i), ...
  • Runtime: Dx(Axis), Xc(Axis, idx), Xf(Axis, idx), FaceArea(Axis,k,j,i), ...

The main coordinate wrapper is grid::Coordinates<geom> (src/grid/geometry.hpp). It is a thin layer over Parthenon's coordinates that provides geometry-specific fixes for metric factors:

  • X/Xi: coordinate at the Parthenon cell-center location
  • Xc: coordinate at the cell centroid (differs from Xi for cylindrical r)

  • Xf: coordinate at the face location

  • In Cartesian geometry, methods delegate to Parthenon.

  • In cylindrical geometry:
  • Dx<Axis::KAXIS>() returns 2*pi (implicit azimuthal direction)
  • Xc<Axis::IAXIS>(idx) computes the radial centroid of the cell (not just the midpoint); Xi<Axis::IAXIS>(idx) remains the "cell center" value from Parthenon
  • FaceArea, EdgeLength, and CellVolume use r-z metric factors (e.g. V = 0.5 * dphi * (r_{i+1/2}^2 - r_{i-1/2}^2) * dz)

CoordinatePacks

Parthenon's coordinate methods are largely inline calculations. For geometry-dependent metric factors (areas, volumes, centroids, etc.) it is often cheaper and simpler to precompute them once per MeshBlock and store them as normal fields.

The grid unit registers and fills these coordinate fields:

  • Field types live in src/grid/coordinates.hpp under kamayan::grid::coords::* and the full list is grid::CoordFields.
  • The fields are allocated with geometry-aware degenerate shapes via grid::CoordinateShape:
  • Cartesian: Dx*, FaceArea*, EdgeLength*, Volume are scalars; X*, Xc*, Xf* are 1D arrays per-axis.
  • Cylindrical: Dx* are scalars; r-dependent quantities (Volume, FaceArea*, EdgeLength*, and X*/Xc*/Xf* for r) are stored as 1D arrays in the radial direction.
  • grid::CalculateCoordinates (src/grid/coordinates.cpp) fills all CoordFields for each block (using grid::CoordinateIndexRanges to respect degenerate/face extents), and is called from the grid unit's InitMeshBlockData callback.

grid::CoordinatePack<geom, ...> (src/grid/coordinates.hpp) wraps a SparsePack holding these geom.* fields and exposes the same API as grid::Coordinates<geom> but indexed by (k,j,i). Internally it maps (k,j,i) onto the degenerate storage layout (scalar/1D), so call sites do not need to care about how each metric is stored.

Example usage (runtime geometry):

problems/isentropic_vortex.cpp:isen_coords_pack
        auto cp = grid::GenericCoordinatePack(geometry, cpack, 0);
        const Real x1 = cp.template Xc<Axis::IAXIS>(k, j, i);
        const Real x2 = cp.template Xc<Axis::JAXIS>(k, j, i);

Generic Coordinates (Packs)

When the geometry is only known at runtime, Kamayan provides variant-based wrappers:

  • grid::GenericCoordinate (src/grid/geometry.hpp) wraps grid::Coordinates<geom> in a variant and forwards calls via visit.
  • grid::GenericCoordinatePack (src/grid/coordinates.hpp) does the same for grid::CoordinatePack<geom, ...>.

In practice GenericCoordinatePack is constructed from a coordinate-field pack (e.g. grid::GetPack(grid::CoordFields(), ...)) and then used like a normal coordinate object.

These are convenient in code like problem generators, but they do add a small overhead compared to templating on Geometry. For performance-critical kernels prefer templating on geom and constructing the matching grid::Coordinates<geom>/grid::CoordinatePack<geom>.

Parameters

Paramter Type Default Allowed Description
<debug>
ndebug Integer 0 Number of debug fields to allocate at CC.
<geometry>
geometry String cartesian [ cartesian, cylindrical, ] grid geometry
<kamayan/refinement>
nref_vars Integer 1 Parameter determined at runtime for the number of registered refinement fields. Never any reason to be set.
<kamayan/refinement0>
derefine_tol Real 2.00000e-01 Error threshold for derefinement
field String dens Field to refine on.
filter Real 1.00000e-02 Noise filtering strength used in Loehner estimator.
max_level Integer 1 max refinement level for this field.
method String loehner [ loehner, derivative_order_1, derivative_order_2, ] Method to use for refinement
refine_tol Real 8.00000e-01 Error threshold for refinement
<parthenon/mesh>
ix1_bc String outflow [ periodic, outflow, reflect, user, axisymmetric, ] Inner boundary condition along x1.
ix2_bc String outflow [ periodic, outflow, reflect, user, ] Inner boundary condition along x2.
ix3_bc String outflow [ periodic, outflow, reflect, user, ] Inner boundary condition along x3.
nghost Integer 2 Number of ghost zones to use on each block.
numlevel Integer 1 Number of refinement levels.
nx1 Integer 32 Number of cells across the domain at level 0.
nx2 Integer 32 Number of cells across the domain at level 0. Set to 1 for 1D.
nx3 Integer 32 Number of cells across the domain at level 0. Set to 1 for 2D.
ox1_bc String outflow [ periodic, outflow, reflect, user, ] Outer boundary condition along x1.
ox2_bc String outflow [ periodic, outflow, reflect, user, ] Outer boundary condition along x2.
ox3_bc String outflow [ periodic, outflow, reflect, user, ] Outer boundary condition along x3.
refinement String adaptive [ adaptive, static, none, ] Mesh refinement strategy.
x1max Real 1.00000e+00 Maximum x1 value of domain.
x1min Real 0.00000e+00 Minimum x1 value of domain.
x2max Real 1.00000e+00 Maximum x2 value of domain.
x2min Real 0.00000e+00 Minimum x2 value of domain.
x3max Real 1.00000e+00 Maximum x3 value of domain.
x3min Real 0.00000e+00 Minimum x3 value of domain.
<parthenon/meshblock>
nx1 Integer 16 Size of meshblocks in x1.
nx2 Integer 16 Size of meshblocks in x2.
nx3 Integer 16 Size of meshblocks in x3.