Refactor of postgkyl

[Maxwell-Rosen]

Refactor of postgkyl

This pull request addresses pain points in postgkyl development. First, postgkyl currently has its own sympy-generated copy of DG kernels that Gkeyll already implements in C. Second, it is a CLI-first library and an awkward scripting interface. Third, equation-informed commands are mixed with basic operations, making structure unclear. This PR fixes all three at once: all DG algebra now runs through Gkeyll’s own compiled kernels via a shim that lives in the gkeyll tree (so drift is a compile error). A fluent Python API becomes the primary surface through GData, with the CLI as a thin shell over it. The new hierarchy separates core operations from equation-informed diagnostics, with the separation enforced by a unit test. Unit testing coverage rises from 33% to 99%. These improvements to the codebase will make postgkyl easier to use, simpler to maintain, and straightforward to extend.

Related GitHub issues

Closes #221. Rename write to save. Performed in this refactor.
Closes #218. [DR] Separating loading data from non-uniform mappings. A verb map is introduced so that loading simply loads data.
Closes #210. [DR] API for scripting interface. This PR implements the API, therefore closing the design review.
closes #198. [DR] Change command “pr” to “print”. All commands are given their full name. The Python side does not use abbreviations or shorthand (interp is now interpolate). The CLI always performs minimal string parsing so that pr resolves to print; unless prandtl is also a command, you have to specify pra/pri.
Closes #188. Plot saving also saves an unnamed .png file. This PR fixes this issue. I’ve fixed this many times in branches.
Closes #159. Plotting grids in 2D and 3D. This PR includes plotly and pyvista inside. This PR was approved by Ammar, but is still blocked by Mana.
Closes #141—hybrid basis specification in interp. We add support for the hybrid basis via the Gkeyll C library in this PR, closing this issue.
Closes #108. Allow specifying interpolation points in interp command. This is addressed using Gkeyll’s quadrature points and nodal evaluation of the interp command. This is a different statement that allows users to specify exactly the logical coordinates at which they want points evaluated.
Closes #52. Add support for hybrid basis used in Vlasov simulations. The Gkeyll C library supports this; therefore, closing this issue.
Closes #13. Weak operations in postgkeyll. This PR supports weak division, multiplication, and integration, with many more features likely to be developed.
Closes #7. Alternative interpolation points. This PR supports quadrature node interpolation and cell edges, resolving this issue.
Closes #6. Refactor the outdated DG part of the code. Oldest issue in the codebase. This issue is resolved by using the C layer of Gkeyll to perform DG operations, utilizing the maxima auto-generated code, and unifying the codebases.

Gkeyll backend support

A linker gpython.c is inside gkeyll/core, which exposes operations such as array_ops and integration to Python. Postgkeyll has a sparse git submodule for gkeyll/core, which builds a dependency-free version of postgkyl at pip install time. Gkeyll’s functionality is linked through an internal compiled file. This design circumvents the need for ctypes or cffi. This structure is advantageous because if the wrapper becomes out of sync, BOTH postgkyl and gkeyll will error, so gkeyll will fail to compile. These errors are loud and diagnosable, while cffi and ctypes can silently fail. This installation works out of the box with pip, without requiring users to install Gkeyll themselves first. For systems where Gkeyll cannot compile, a pure-Python reader is available. Similar to the Lua-C API rewrite, this hardens the compatibility layer and helps with maintainability.

Modal versus nodal operations

Gkeyll owns everything modal, as this functionality is precisely what the Gkeyll repository is built for. Gkeyll performs weak arithmetic, integration, interpolation, and DG evaluation. On the other hand, pointwise operations are exact, so NumPy is responsible for these operations. NumPy affords a rich ecosystem of commands, affording great flexibility. Only interpolation crosses the line between Gkeyll arrays and NumPy arrays. .gkyl files default to loading as modal, but sometimes we store nodal data. As a result, we must store whether an array is modal or nodal in the metadata of the files we save.

Unified script / CLI interface

The script mode and CLI interface now mirror each other. In script mode, one can type

import postgkyl as pg

pg.load('file').interpolate().select(z0=0.0).plot()

While in the command line, one types

pgkyl file interp sel --z0 0.0 pl

Having a unified syntax for both the script and the CLI simplifies working with this package and makes post-processing review easier.

Scripts can go even further, and here is the syntax that this PR aims to support and extend script mode.

import numpy as np
import postgkyl as pg

m0 = pg.load("sim_elc_M0_5.gkyl")
m1 = pg.load("sim_elc_M1_5.gkyl")
upar = m1 / m0              # weak division, inside Gkeyll's kernels
quad = upar.to_quad()       # values at quadrature points
root = np.sqrt(quad)        # NumPy ufuncs work on point-value data
modal = root.to_modal()     # exact projection back to DG coefficients
print(modal)                # prints like a NumPy array

To keep script and CLI in sync, both call upon the same API. The click interface is a thin layer on top of the API, making the repository more modular, while keeping all aspects in sync.

A necessary split

There is a fundamental issue with using this command structure. GData presently contains the grid and values, but introducing methods creates a paradox. The operations need to import GData, but GData has methods of operations, forming a cycle. To break this cycle, GData is separated into 2 classes, GData and GDataState. GData contains the methods, while GDataState contains the data itself. These objects live in different layers of the hierarchy.

Unit testing

postgkyl on main only has 33% coverage. This refactor was meticulous in its use of test-driven development. Now, the repo has 99% code coverage, including our most fundamental functionality. The layer hierarchy is enforced by the test_postgkyl.py test, which walks through every real import and checks it against an allowlist. The architecture diagram cannot silently rot because a test fails when a test fails.

The postgkyl hierarchy

Similar to Gkeyll, this PR implements this proposed hierarchy. Having layers in postgkyl keeps our code organized, modular, and readable. This layering mirrors the Gkeyll repository.

What are diagnostics

This PR separates commands into a set of operations/ and diagnostics/. Operations are fundamental low-layer manipulations of the data (e.g. interpolate, plot, select), without any knowledge about the data they are operating on. Diagnostics are equation-informed tools which help users understand the Gkeyll simulations (e.g. 10m, gyrokinetics, and vlasov). Sitting on top of the API, they can call the fundamental operations, just as a Python script would. This layer is highly flexible, enabling future users to implement and share their own post-processing tools without worrying about the lower levels of the codebase. While the CLI must group commands to preserve chaining, refactoring it into a thin wrapper affords greater modularity for organizing the underlying scripts.

File tree

Below is the file tree, in hierarchical order. This hierarchy is enforced by a unit test test_postgkyl.py.

src/postgkyl/
│
├─ __init__.py          facade · `import postgkyl as pg`              [SURFACE]
│
├─ cli/                 Thin Click shells: argv → gdata / diagnostics   [SURFACE]
│   ├─ app.py
│   └─ commands/        All the callable commands from the command line
│
├─ diagnostics/        equation-specific physics · one module        [COMPOSITION]
│                      per equation model
│
├─ gdata/                ★ THE FLUENT SURFACE  (sits ABOVE operations)  [FLUENT API]
│   ├─ gdata.py          class GData(GDataState) + .interpolate()/.plot()
│   ├─ gdatagroup.py     fluent GDataGroup: broadcasts verbs over its members
│   ├─ verbs.py          module-level fluent verbs with no single `self`
│   │                    (collect/evaluate/relchange/animate) — one-line
│   │                    delegations to `operations`, shared by GData and GDataGroup
│   └─ load.py           pg.load(...) → returns a GData
│
├─ operations/         one function per verb · the single seam        [VERBS]
│   ├─ interpolate.py    interpolate(d: GDataState) -> GDataState      (core verbs
│   ├─ select.py                                                       only —
│   ├─ plot.py                                                         equation-blind)
│   ├─ animate.py        terminal: sequence of datasets → render's animation engine
│   ├─ average.py        terminal-adjacent: weighted average over a dim subset,
│   │                    stays modal/gkyl-native (composes with further verbs)
│   ├─ eval_at_coord_proj.py  terminal-adjacent: eval at coords, project to the
│   │                    lower-dim basis for survivors, stays modal/gkyl-native
│   ├─ local_poly.py     modal coefficients → discontinuity-preserving plot mesh
│   └─ _materialize.py   shared modal → NumPy-shadow bridge used by plot/animate
│
├─ render/             matplotlib · plotly · pyvista → gdatastate/numerics [BACKEND]
│                      (below operations, which delegates plot() to it)
│
├─ gdatastate/         ★ THE CONTAINER  (state only, NO verbs)        [CONTAINER]
│   ├─ gdatastate.py          class GDataState: grid·values·ctx·_result·dunders
│   ├─ gdatastategroup.py     GDataStateGroup
│   └─ guards.py         shared field-domain guard (backend=="gkyl" -> raise);
│                        one home for the ".interpolate() first" check reused
│                        across operations/diagnostics instead of retyped per verb
│
├─ numerics/           pure NumPy math · 0 internal imports           [LEAF]
├─ dg/                 interpolation bridge + modal ops → gpython     [ENGINE]
├─ io/                 readers (C-native first) + writer → gpython    [ENGINE]
└─ gpython/            ★ THE FOREIGN FLOOR · compiled shim            [FLOOR]
    ├─ csrc/             _gpythonmodule.c — CPython extension over gkyl_gpython.h
    │                    (the shim itself lives in gkeyll/core/zero/)
    ├─ _gpython.so       built extension (scripts/build_gpython.sh)
    ├─ _lib.py           loads _gpython · GPYTHON_API_VERSION handshake
    ├─ array.py          GkylArray — capsule owner of a gkyl_array
    ├─ basis.py          basis cache + interpolation matrices via the shim
    ├─ rio.py            file loading via gkyl_array_rio
    └─ kernels.py        weak mul/div/inv · lincomb · reduce · integrate

A graphical user interface (GUI) can operate at the surface level, near the CLI, earning the abstraction of another interface.

Development notes

Version-breaking changes

Since this PR is an entire rewrite of the library, I anticipate it will have version-breaking changes. The old GData class has a fundamentally different structure. Please reach out if your old scripts no longer work and we can fix them. I am very open to feedback that improves the design of our codebase.

I paid special attention to porting the latest developments for gyrokinetics postprocessing, but this took several weeks, and the codebase is very active, so there is a non-negligible chance that I missed something. The best way to know is to get others' eyes on the code and use it.

Nodal arrays need a property value_forrm = "nodal" added to their metadata. Here is an example.

  struct gkyl_msgpack_map_elem meta[] = {
    { .key = "poly_order", .elem_type = GKYL_MP_UNSIGNED_INT, .uval = 0 },
    { .key = "basis_type", .elem_type = GKYL_MP_STRING, .cval = "serendipity" },
    { .key = "value_form", .elem_type = GKYL_MP_STRING, .cval = "nodal" },
  };

Agentic development

This refactor was carried out by Claude Code. I constructed a 15-staged refactoring plan with associated agents to orchestrate. I have attached the files for consideration. The agent was created in Fable 5, reviewed and then implemented with Sonnet 5. Correctness is enforced through a comprehensive layer of unit testing and human review.

Bug fixes

Along the way, the agents found several (5-10) genuine bugs in the codebase. They can be found in the agent's reports and the commit history. These bugs are fixed by the agents, and unit tests are added.

Known limitations

This refactor brings out some limitations of Gkeyll itself. Gkeyll does not have div/mul/integrate kernels for gkhybrid and hybrid basis functions. Integration for the GK distribution function is in gyrokinetic/ rather than core/. This is beyond the scope of this PR.

Modal integration is only supported in total. Per-axis integration was only supported in a nodal since in postgkyl, which is carried over. Gkeyll does not support per-axis integration in the kernels, so Gkeyll must be extended for this application.

Out of scope

In a perfect codebase, I'd like to have the website compiled directly from postgkyl, linked to tested examples and checked through CI. The codebase should also include linting and formatting, but I didn't implement them here. That's a separate discussion. I've included a few walkthrough scripts to guide new users through the codebase. We should prevent people from putting code in main that lacks unit tests. CI should check code coverage to keep our coverage high.