Packages

Introduction

A brief intro to packages in KHARMA is provided here; if you're interested in the details, we expound on KHARMA's packages here.

Packages are Parthenon's way of introducing modularity in the code in the way of physics and features. Each package has its own namespace within which we define its local variables, global constants (via the params object, an instance of parthenon::Params), its functions, and its tasks.

Unless marked private, packages can access variables declared in other packages. The ability to share variables eliminates the need to save multiple copies of the same variable. Probably the best example of this is the fluid adiabatic index gamma, a constant declared in the GRMHD package, which can be pulled like,

auto gam = pmb->packages.Get("GRMHD")->Param<Real>("gamma");

and is done so all across the code. Every package has an Initialize method declared in its namespace where its variables and constants are registered. The ProcessPackages function defined in the KHARMA namespace loads the required packages during startup. As explained in the section on the ImEx driver, the Implicit package is loaded by default if the driver is ImEx. Additionally, KHARMA initializes the GRMHD, Globals and floors package by default. The remainder of the packages are initialized only if specified during runtime. If you've written a new package, you must add its Initialize method to the packages object here.

In additional to global constants that are available through the params object, packages can also include variables that are defined over the entire mesh (also called fields). The variables contain two kinds of information - metadata and data. The metadata contains all the information about the variable, for eg. its name, size, datatype and any other additional info about the variable that facilitates accessing and manipulating the data. The metadata is specified through a vector of MetadataFlags. For example have a look at the following definition of a field in the Implicit package.

Metadata m_real = Metadata({Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::OneCopy});
pkg->AddField("solve_norm", m_real);

Here we have defined an object of Parthenon's Metadata class called m_real by called Metadata's constructor. We've provided a MetadataFlag vector that specifies four flags, Metadata::Real tells it that the field defined by this metadata will consist of Reals, Metadata::Cell clarifies that array of Reals will reside at the cell centers, Metadata::Derived explains that the data in this array can be computed on-the-fly from a more fundamental state of variables, and finally Metadata::OneCopy tells Parthenon that any field defined with this metadata will be shared between stages (recall that we have a multistage driver). The Parthenon docs have a page that lists all the metadata flags available. The second line in the code excerpt above defines a field solve_norm in the implicit package with m_real as its metadata.

We now go through the packages in KHARMA.

Globals

GRMHD

B field

Floors

Current

Electrons

Implicit

Overview of the solve

The implicit solver is essentially a multidimensional Newton-Raphson scheme. In its current form it finds the root for the residual of the system of balance laws, $\partial_t\boldsymbol{U}(\boldsymbol{P}) + \boldsymbol{\nabla}\boldsymbol{F}(\boldsymbol{P}) = \boldsymbol{S}(\boldsymbol{P})$ where $\boldsymbol{U}$ is the vector of conserved variables that are analytic functions of $\boldsymbol{P}$ - the vector of primitives, $\boldsymbol{F}$ are the corresponding fluxes and $\boldsymbol{S}$ are the source terms. To get a more detailed introduction to hyperbolic PDEs and balance laws, have a look at this page. The residual is then defined as $\boldsymbol{R}(\boldsymbol{P})\equiv\partial_t\boldsymbol{U}(\boldsymbol{P}) + \boldsymbol{\nabla}\boldsymbol{F}(\boldsymbol{P}) - \boldsymbol{S}(\boldsymbol{P})$. This is the cost/loss function that we'd like to minimize over the vector space of the implicitly evolved variables. This solver was necessitated by the introduction of extended MHD (EMHD) in KHARMA which contains implicit source terms.

As is the case for the multidimensional Newton-Raphson scheme, the solver computes a Jacobian $\boldsymbol{J}(\boldsymbol{P})$ and solves for the equation $\boldsymbol{J}(\boldsymbol{P})\cdot\boldsymbol{\delta P} = -\boldsymbol{R}(\boldsymbol{P})$ where $\boldsymbol{\delta P}$ is the update to the primitives in one iteration of the solve. The step length $\boldsymbol{\delta P}$ is the direction in which the update must take place; however, it is possible that when one is close to the solution, updating $\boldsymbol{P}$ with $\boldsymbol{\delta P}$ can overstep the solution. To avoid this, the solve includes a backtracking line search where we select a suitable step length $\lambda$ and update the primitives like $\boldsymbol{P} = \boldsymbol{P} + \lambda\boldsymbol{\delta P}$.

NOTE: This is a local solve ie., it minimizes the residual and computes the update to the primitives on a zone-by-zone basis.

The details of the scheme are present in section 3 of this paper.

Package variables and constants

The Implicit package defines a bunch of constants that are relevant to the implicit solve. These are,

• min_nonlinear_iter (default: 1) - The minimum number of iterations that must be performed by the solver.
• max_nonlinear_iter (default: 3) - The maximum number of iterations that must be performed by the solver.
• rootfind_tol (default: 1.e-12) - The tolerance that must be met during the solve to consider it a success.
• jacobian_delta (default: 4.e-8) - The small change considered for each primitive - $\Delta P$ - to compute the Jacobian.
• use_qr (default: true) - A boolean that dictates whether the QR decomposition should be considered while solving the system of linear equations. If false, the solver uses LU decomposition. While QR has certain stability benefits, we've found that it can lead to potentially inaccurate results when performing a >= 8x8 solve on GPUs.
• linesearch (default: true) - A boolean that decides if backtracking linesearch must be performed.

The following arguments are relevant only if linesearch has been enabled.

• max_linesearch_iter (default: 3) - The maximum number of iterations that must be performed by the linesearch.
• linesearch_eps (default: 1.e-4) -

EMHD