Parthenon

Note on this page: information on Parthenon can quickly fall out of date. For details & up-to-date info on Parthenon, always consult the Parthenon project documentation.

KHARMA and Parthenon

Source code

The current KHARMA release uses a slightly-customized fork of Parthenon, hosted here. The next major release will revert to mainline Parthenon, applying the few necessary patches at compile-time.

Support

Parthenon is an open-source code developed primarily at Los Alamos National Laboratory. It has a broad community outside of KHARMA: for example, you may want to join the development matrix room, and check the open pull requests and issues for a view of ongoing development.

KHARMA Structure as a Parthenon Code

Parthenon is a framework, meaning that it controls the flow of the KHARMA executable program, calling particular functions that KHARMA defines and registers with Parthenon. Most of these user-defined functions are part of C++ classes called Drivers and Packages.

Drivers

A Driver represents a fundamental capability of the code, and nearly all Parthenon codes define just one. KHARMA is not (any longer) an exception, defining just a single KHARMADriver (kharma_driver.cpp). What are colloquially still referred to as "drivers" in KHARMA are really different TaskLists defining what constitutes one step of a simulation (kharma_step.cpp, imex_step.cpp, simple_step.cpp). Very few applications or physics extensions will warrant writing a new TaskList, much less a new Driver -- each Package is given broad flexibility for adding variables, functions, and parameters as necessary.

However, one could imagine implementing an entirely separate Driver to e.g. compute an average over previous snapshot files, or create a ray-traced image -- anything that represents a totally separate mode of operation, with completely different ideas of what "initialization" or "execution" means.

The Driver is discussed in detail here.

Packages

Parthenon (and by extension, KHARMA) is designed to evolve many different physical processes at once, while retaining the ability to disable any subset without paying a performance or complexity penalty when they are not in use. This allows the code to include many extra capabilities (e.g. Monte Carlo transport, implicit solvers, etc.) in the form of packages, yet pay no performance or debugging cost when running simpler simulations with these packages disabled.

Note that packages will be referred to here by their capitalized string names, defined at the top of each package's Initialize function. Convention is to define each package's functions in a namespace corresponding to this name (e.g. Electrons for package "Electrons"), and to keep implementations in a folder with an indicative snake_case name, with a main file of the same name (e.g. electrons/electrons.cpp). Supporting files for a package are kept in the same folder.

Descriptions of each package can be found in Packages.

The KHARMA namespace

In addition to package-specific calls, Parthenon allows the user code to register some calls to happen during each step. KHARMA generally uses these calls to loop through calls of the same name in each package, as well as to update a few global variables, available as parameters of a package called "Globals".

Coordinates

The coordinates folder is not a package. Instead, it implements subclass of the class representing coordinates for Parthenon MeshBlock objects. An instance is created with each MeshBlock, and destroyed with a deleted MeshBlock.

The more important and more accessible piece of this folder is the GRCoordinate class in gr_coordinates.{h,c}pp. This class more or less takes the role of the Grid class in iharm3D and similar, with the member arrays of that class replaced by function calls in GRCoordinates.

The class includes a member called coords which is an instance of a CoordinateEmbedding. An embedding consists of a "base" coordinate system (e.g. Cartesian Minkowski, or Spherical Kerr-Schild coordinates) and a "transformation" to "native" coordinates (potentially nothing, but commonly exponentiating the radial coordinate or otherwise concentrating coordinates around areas of interest). The zones of a simulation (or of the base grid, in SMR) are evenly spaced in a logically Cartesian grid spanning the "native" coordinate domain.

The CoordinateEmbedding object handles metric calculations and transformations, and the GRCoordinates object caches useful values (metric, connection coefficients) in "native" coordinates, to be used throughout the code. Perhaps the most intense use of these functions is in the file prob/bondi.cpp, which also demonstrates that in the pursuit of generality they are somewhat unwieldy. More intense use cases might require amendments to particular coordinate systems, or writing syntactic sugar functions in GRcoordinates and CoordinateEmbedding.

One particularly important caveat for adding coordinate systems: there are quite a number of functions in CoordinateEmbedding which are special-cased: e.g. if (base_system == KerrSchild) {} else {}. If you are adding a new BaseCoords or Transform, you'll want to comb through the functions in CoordinateEmbedding and either add your new class to one of the existing implementations, or write your own suitable implementation of the special-cased function. The separation of responsibility between CoordinateEmbedding and the various BaseCoords and Transforms is a bit hazy, as a result of trying to keep as much as possible out of the latter and therefore off the GPU. The whole mess around coordinate systems and CoordinateEmbedding could use an overhaul, though.

Problems

The prob folder is also not a package. It holds a main C++ file problem.cpp and a slew of particular initialization functions, defined in headers or source files as complexity warrants. Even though it's only run once, initialization is generally done on-device with Kokkos loops, because it's more convenient than filling and copying a host array. If you would prefer to initialize your problem on the host, consult iharm_restart.cpp, which demonstrates correct use of the functions GetHostMirrorAndCopy() and DeepCopy() to declare host-accessible versions of Kokkos Views, and to copy them to device.

The prob folder also contains a file post_initialize.cpp with code run after the whole mesh has been filled (it's run directly from main.cpp). This is used mostly to apply and normalize the magnetic field, but it's become home to

Parthenon Datatypes

Mostly, when writing new code for KHARMA, you won't be dealing with the overall structure above. Once you have stubs in place and compiling for some new feature, you likely won't need to consult the above again. Instead, you'll be dealing with Parthenon's utility functions and datatype classes.

Example function

Here is Flux::AddGeoSource from KHARMA 2024.8, defined in kharma/flux/flux.cpp, which you may recognize from the Kokkos example:

void Flux::AddGeoSource(MeshData<Real> *md, MeshData<Real> *mdudt, IndexDomain domain)
{
    // Pointers
    auto pmesh = md->GetMeshPointer();
    auto pmb0  = md->GetBlockData(0)->GetBlockPointer();
    auto pkgs = pmb0->packages;
    // Options
    const auto& pars = pkgs.Get("GRMHD")->AllParams();
    const Real gam   = pars.Get<Real>("gamma");

    // All connection coefficients are zero in Cartesian Minkowski space
    // TODO do we know this fully in init?
    if (pmb0->coords.coords.is_cart_minkowski()) return;

    // Pack variables
    PackIndexMap prims_map, cons_map;
    auto P    = md->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("Primitive")}, prims_map);
    auto dUdt = mdudt->PackVariables(std::vector<MetadataFlag>{Metadata::Conserved}, cons_map);
    const VarMap m_p(prims_map, false), m_u(cons_map, true);

    // EMHD params
    const EMHD::EMHD_parameters& emhd_params = EMHD::GetEMHDParameters(pmb0->packages);
    
    // Get sizes
    IndexRange3 bd = KDomain::GetRange(md, domain);
    auto block = IndexRange{0, P.GetDim(5)-1};

    pmb0->par_for("tmunu_source", block.s, block.e, bd.ks, bd.ke, bd.js, bd.je, bd.is, bd.ie,
        KOKKOS_LAMBDA (const int& b, const int &k, const int &j, const int &i) {
            const auto& G = dUdt.GetCoords(b);
            FourVectors D;
            GRMHD::calc_4vecs(G, P(b), m_p, k, j, i, Loci::center, D);
            // Call Flux::calc_tensor which will in turn call the right calc_tensor based on the number of primitives
            Real Tmu[GR_DIM]    = {0};
            Real new_du[GR_DIM] = {0};
            for (int mu = 0; mu < GR_DIM; ++mu) {
                Flux::calc_tensor(P(b), m_p, D, emhd_params, gam, k, j, i, mu, Tmu);
                for (int nu = 0; nu < GR_DIM; ++nu) {
                    // Contract mhd stress tensor with connection, and multiply by metric determinant
                    for (int lam = 0; lam < GR_DIM; ++lam) {
                        new_du[lam] += Tmu[nu] * G.gdet_conn(j, i, nu, lam, mu);
                    }
                }
            }

            dUdt(b, m_u.UU, k, j, i)           += new_du[0];
            VLOOP dUdt(b, m_u.U1 + v, k, j, i) += new_du[1 + v];
        }
    );
}

This function adds the geometric source term in the GRMHD equations, $\sqrt{-g} T_\lambda^\kappa \Gamma^\lambda_{\nu \kappa}$. It was chosen for being fairly standard, while demonstrating nearly everything you'll need for most new functions. There are many elements of this function you'll see only in Parthenon-based codes. Let's run through them in order.

Declaration

• The function returns a TaskStatus object. This is a Parthenon enumerated value with possible values fail, complete, incomplete, iterate, skip. You will only commonly need complete and fail; the others are for running tasks concurrently, or repetition of function calls within a task list.
• It is a member of the namespace Flux, used to indicate it is part of the "Flux" package. Placing the function in "Flux" rather than "GRMHD" indicates that it is a theory-independent operation: GRHD, GRMHD, or Extended GRMHD simulations all use this same function.
• It takes two MeshData<Real> arguments: md represents the conserved fluid state, and mdudt the full derivative $dU/dt$ (that is, what you would multiply by $dt$ and add to $dU$ in a first-order scheme). Real is KHARMA's name for double when referring to normal variables on the grid (see conventions).
• Like many KHARMA-specific functions, it also takes an IndexDomain over which to operate. This particular function is always called with a domain of IndexDomain::interior, representing all physical zones. Some functions will need to be called over IndexDomain::entire in order to update the whole array including ghost zones, or e.g., IndexDomain::inner_x1 to update the ghost zones on a particular side of each block. So, it is best practice to write the function generally as we've done here.

Pointers & Options

• The Mesh and MeshBlock objects hold most important information about the block size, boundaries, coordinates & locations. Anything you would make global in a single-mesh code is in these objects. These are not to be confused with the MeshData and MeshBlockData objects, which hold all the variables themselves. Each MeshBlock has one MeshBlockData containing its variables. Multiple MeshBlockData objects are packed together to form a MeshData object (not every MeshBlockData though, just one or more). The Mesh object, by contrast, does contain every MeshBlock.
• Options in KHARMA are first parsed by Parthenon into a ParameterInput object, then selectively copied by packages into their own StateDescriptor objects during their Initialize functions. The adiabatic index, for example, is stored as a part of the "GRMHD" package, and this line demonstrates accessing and assigning it. Generally, packages should try not to access one another's parameters, but $\gamma$ is an exception.
• If the coordinate system is specifically Cartesian Minkowski, we skip this function, since all connection coefficients are zero.

VariablePacks

• In order to access the variables in a MeshBlockData or MeshData object, we generate a "pack." When packing variables from MeshBlockData objects, the resulting arrays are of the form P(p, k, j, i) indexed by variable p, followed by zone indices. Packs over MeshData objects are P(b, p, k, j, i) by block, then variable, then zones[^1]. Note that many device-side functions expect VariablePack<Real>, the former pack over one block -- for example, GRMHD::calc_4vecs or GRMHD::calc_tensor above. When operating over a MeshBlockPack where there's an additional index, we pass these functions the sub-pack for each block, e.g. P(b).
• Each VariablePack generates a VarMap indicating which integer index corresponds to which packed variable. In cases like this source term, when we need to add different terms to different variables, this map is important and must be accessible inside our device function later on (when applying the same operation to all variables, e.g. for boundary exchanges, it matters less which is which). Details on using the VarMap are on the Kokkos page.
• A GRCoordinates object is carried with each MeshBlock and packed into each VariablePack. When we're running a kernel over multiple meshblocks, as in this example, we can pull our particular block's coordinates out of the VariablePack, as we do in the first line of the device-side function in the example. When running over a single block, we can pull out the GRCoordinates object outside the kernel, since it won't change within the kernel.

IndexRanges

• The Parthenon IndexRange object is just a struct of two integers, beginning and end. Note again that the end is inclusive, as described. The MeshData and MeshBlockData objects both carry boundary information. This info is consistent across all MeshData because regardless of refinement level, each block is required to be the same shape.
• The IndexRange3 object is (currently) KHARMA-exclusive -- it's just a struct with six members defining three ranges. The functions in the KDomain namespace (see kharma/domain.hpp) provide an easy way to flexibly get different ranges. For example, the interior with a one-zone halo is GetRange(md, IndexDomain::interior, -1, 1), or the interior range of a face-centered field (one index greater in the face direction) is GetRange(md, IndexDomain::interior, F1), etc, etc. GetRange will never return an index outside the array, or a nonzero index in a trivial dimension.

---

[^1]: You may also see a block indexed separately, q = P(b); q(p, k, j, i), as in this example where P(b) is passed as an argument to a function expecting to index into a single MeshBlock of data.

Further info on Datatypes

For more complete and guaranteed up-to-date info on all of these datatypes, see Parthenon's documentation, notably the pages on the mesh types and ParArray.

For more info on Scratchpads (and the nested loop constructions generally), see the Kokkos documentation on hierarchical parallelism.