-
Notifications
You must be signed in to change notification settings - Fork 26
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 uses a slightly-customized fork of Parthenon, hosted here. This is necessary mostly because KHARMA changes the default coordinate system used in compiling Parthenon, from the smaller, Newtonian/Cartesian UniformCartesian object to KHARMA's general-relativistic version, called GRCoordinates (described later on this page). Since the coordinate object must be included as a header in Parthenon code, there is not (or I do not know of) a good option for injecting user coordinate systems, short of patching the original code, which is most easily done by maintaining a fork. This has the pleasant side effect of letting me control the Parthenon update cadence for different KHARMA versions.
Parthenon is an open-source code developed primarily at Los Alamos National Laboratory. It has a broad community outside of KHARMA: notably you can find help in their matrix room
Parthenon is a framework, meaning that it controls the flow of the KHARMA executable program, calling into particular functions that KHARMA defines and registers with Parthenon. Most of these user-defined functions are part of classes called Drivers and Packages.
A Driver represents a fundamental capability of the code. KHARMA currently has one driver, HARMDriver, which implements the HARM algorithm (and all current supported variations on the scheme, e.g. electrons, magnetic field constraints, fluxes, floors, etc). It is likely this will be the code's only driver for the foreseeable future, as straightforward physics extensions do not merit writing a whole new driver. That is, algorithms incorporating e.g. M1 closure, semi-implicit time-stepping, etc. likely share too much code to warrant switching out this fundamental piece of the design.
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.
Parthenon (and by extension, KHARMA) is designed to accommodate many different physical processes, 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 retain any extra capabilities (e.g. Monte Carlo transport, implicit solvers, etc.) while paying little or no performance cost in simpler simulations without these features, and little maintenance cost when making changes to other packages.
This is an ongoing goal for KHARMA -- currently it is difficult to imagine the code functioning without the default "GRMHD" package enabled, but this is conceivable in the future. Note that packages will be referred to here by their capitalized string names, defined at the top of each package's Initialize function. Convention dictates each package's functions are defined in a namespace corresponding to this name (e.g. Electrons for package "Electrons"), and implementations are kept in a folder and in files with an indicative snake_case name (e.g. electrons/electrons.cpp). Supporting files for a package are kept in the same folder.
The main KHARMA package, "GRMHD", handles the density, internal energy, and velocity primitives and conserved variables as described in Gammie et al. (2003). It currently requires, but does not itself evolve, the magnetic field.
Currently only one working B field package is available for evolving the magnetic field in KHARMA: "B_FluxCT". It updates a cell-centered representation of the magnetic field, handling both the conversion of primitive to conserved values (i.e. flux <-> field) and updates to the cell-centered values in a way that preserves zero field divergence, using the Flux-CT method of Toth (2000).
A package "B_CT" is planned which will implement face-centered constrained transport with various upwinding schemes, once face-centered fields are better supported in Parthenon. This scheme will also include prolongation and restriction operations which preserve the field divergence, making possible simulations with both magnetic fields and adaptive or static mesh refinement schemes.
The package "B_CD" implements constraint damping as in Dedner et al (2002), but is useful only in non-relativistic contexts as it does not implement a covariant version of the scheme. Given KHARMA is a GRMHD code, that makes it fairly useless at present.
The Electrons package evolves a set of one or more entropy values corresponding to specifically the electrons in a fluid. This is fed by a source term which calculates total dissipation and assigns some portion of the value to electron heating. This is useful because accretion plasmas are likely strongly two-temperature, necessitating some way of determining the electron temperature or likely energy distribution.
The package uses the scheme of Ressler et al. (2015), tracking electrons heated according to any or all of the models from Kawazura (2019) and (2021), Werner (2016), Rowan (2017) and Sharma (2007).
Small features orthogonal to the rest of KHARMA's code have been split into a number of smaller packages. These are largely just anything which would clutter grmhd.{h,c}pp and fluxes.{h,c}pp.
- A "wind" source term with a slew of its own parameters
- Tools to calculate the 4-current, which is computed only for output as a part of dump files
- Tools to calculate other reductions used only for output, e.g. accretion rate, horizon B field flux, etc.
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".
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 (e.g. in AMR).
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). Note that zones are equally spaced in a logically Cartesian grid in "native" coordinates.
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 `Emplace
If fate smiles upon you, you will never need to touch the CoordinateEmbedding class except:
- to add new
BaseCoordsandTransformoptions, which can be added just like the existing alternatives to the functionEmplaceSystems - to add syntactic sugar functions described above and implemented in terms of existing functions The comments document how to make these sorts of additions in good detail, and writing the particular classes for new coordinate systems should be as trivial as transliterating the metric to C.
However, if for some reason you want to redesign the way that KHARMA allows coordinate system selection at runtime... don't. The std::variant system we use to fake virtual functions on the device side is delicate and can produce some awful template errors at compile time, and baffling segfaults after access at runtime. If the class needs to be substantially modified for some reason, the most compelling option may be to eliminate it entirely, and make the coordinate system a compile-time option.
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. Note that initialization is generally done with Kokkos loops, even though it is run only once. This is just what's proven convenient, frankly, as the bounds checking and flat loop structure make for readable code without manually shuffling data from CPU->GPU memory. If you would prefer more traditional initialization in designing your problem, consult iharm_restart.cpp, which demonstrates correct use of the functions GetHostMirrorAndCopy() and DeepCopy() in declaring host-accessible versions of Kokkos Views, and transferring the resulting data to the device correctly.
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.
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 ever again. Instead, you'll be dealing with Parthenon's tools and datatypes.
Remember I said we'd return to our example eventually? Here we are. Let's look at GRMHD::AddSource again.
TaskStatus GRMHD::AddSource(MeshData<Real> *md, MeshData<Real> *mdudt)
{
FLAG("Adding GRMHD source");
// Pointers
auto pmesh = md->GetMeshPointer();
auto pmb0 = md->GetBlockData(0)->GetBlockPointer();
// Options
const Real gam = pmb0->packages.Get("GRMHD")->Param<Real>("gamma");
// Pack variables
PackIndexMap prims_map, cons_map;
auto P = GRMHD::PackMHDPrims(md, prims_map);
auto dUdt = GRMHD::PackMHDCons(mdudt, cons_map);
const VarMap m_u(cons_map, true), m_p(prims_map, false);
// Get sizes
IndexDomain domain = IndexDomain::interior;
auto ib = md->GetBoundsI(domain);
auto jb = md->GetBoundsJ(domain);
auto kb = md->GetBoundsK(domain);
auto block = IndexRange{0, P.GetDim(5)-1};
const auto& G = dUdt.coords;
pmb0->par_for("grmhd_source", block.s, block.e, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
KOKKOS_LAMBDA_MESH_3D {
// Then calculate and add the GRMHD source term
FourVectors Dtmp;
GRMHD::calc_4vecs(G(b), P(b), m_p, k, j, i, Loci::center, Dtmp);
GRMHD::add_source(G(b), P(b), m_p, Dtmp, gam, k, j, i, dUdt(b), m_u);
}
);
FLAG("Added");
return TaskStatus::complete;
}There are many elements of this function you'll see only in Parthenon-based codes. Let's run through them in order.
- The function returns a
TaskStatusobject. This is a Parthenon enumerated value with possible valuesfail,complete,incomplete,iterate,skip. You will only commonly needcompleteandfail; the others are for running tasks concurrently, or repetition of function calls within a task list. - It is a member of the namespace
GRMHD, used to indicate it is part of the "GRMHD" package. - It takes two MeshData arguments, representing the fluid state (U) and the calculated full derivative dU/dt (that is, what you would multiply by dt and add to dU in a first-order scheme). I'll get more into what MeshData means below, but I will mention that "Real" is Parthenon's name for "double." It is designed to allow for calculations in single precision, but as there is no facility for mixed-precision calculations with e.g. fluid variables as single-precision and magnetic fields as doubles, it is useless to an MHD code.
- The
MeshandMeshBlockobjects 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. TheMeshBlockobject is also the host of thepar_forcall, so oneMeshBlockor another is pulled out in some form in nearly every host function. In this case, we just take the first block in the blocklist kept byMesh, since we'll be launching a kernel over the whole mesh. - Options in KHARMA are first parsed by Parthenon into a
ParameterInputobject, then selectively copied by packages into their ownStateDescriptorobjects during theirInitializefunctions. The specific heat ratio, for example, is stored as a part of the "GRMHD" package, and this line demonstrates accessing and assigning it.
- Generally, variables in KHARMA are added one by one, with separate names ("rho", "B") for flexibility. However, it's clearly more efficient to be able to iterate over a bunch of variables at once, for example when computing fluxes, exchanging boundaries, or in this case just adding a source term to a particular set of variables. Thus Parthenon allows packing variables
rho(k, j, i)of zones into a single arrayP(p, k, j, i)indexed by variablep, followed by zones. In fact, Parthenon has taken this one step further, allowing packs containing multiple blocks of data, indexedP(b, p, k, j, i)by block, then variable, then zones[^1]. - Each
VariablePackgenerates aVarMapindicating which integer index corresponds to which actual 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. Parthenon's type for representing this (PackIndexMap) is very heavy, so we generate a KHARMA type called aVarMap-- a struct with this information, described later in the documentation. - The Parthenon IndexRange object is just a struct of two integers, beginning and end. Note again that the end is inclusive, as described. The
MeshDataandMeshBlockDataobjects both carry boundary information. This info is consistent across allMeshDatabecause regardless of refinement level, each block is required to be the same shape. - A GRCoordinates object describing the location of its cells is carried with each MeshBlock. When packing blocks into a
MeshDataobject, coordinates are also packed into a list, available as shown.
[^1]: You may also see a block indexed separately, q = P(b); q(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.
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.