Application-Defined Abstractions to Separate Numerical Intent from Computation Orchestration

Application-Defined Abstractions to Separate Numerical Intent from Computation Orchestration

arXiv:2611.11100v1 [cs.SE] 8 Nov 2026

ELENA ASCHMANN, Swiss Federal Institute of Technology, Switzerland ANSHU DUBEY, RIKEN Center for Computational Science, Japan MOHAMED WAHIB, RIKEN Center for Computational Science, Japan JENS DOMKE, RIKEN Center for Computational Science, Japan

Fig. 1. A methodology to achieve separation of the numerical intent of the solver from its surrounding code through application defined and controlled abstractions In High-performance Computing the holy grail has been to translate numerical algorithms expressed in mathematical notation into code that performs well across a range of platforms. A variety of portability frameworks have abstracted data layout, parallelization and orchestration from application’s view, however, the resulting format is still code that conforms to the idioms of the framework. Additionally, substantial refactoring is typically needed to achieve performance goals. Increasing sophistication of Large Language Models (LLMs) naturally raises the question whether they can be used to make application codes and/or solvers both portable and performant. The obstacle is the lack of correctness guarantee in LLM produced code. In this paper we pose the question differently – can we start with expressing numerical intent of the code into an immutable core arithmetic component isolated from the surrounding code such that the surrounding code can be safely generated by the LLMs? We answer this question by using Spark, a mature compressible hydrodynamics solver in Flash-X, a multiphysics software system. We use LLM assistance in devising the abstractions and embedding Authors’ Contact Information: Elena Aschmann, Swiss Federal Institute of Technology, Zurich, Switzerland, easchmann@ethz.ch; Anshu Dubey, RIKEN Center for Computational Science, Kobe, Japan, anshu.dubey@riken.jp; Mohamed Wahib, RIKEN Center for Computational Science, Tokyo, Japan, mohamed.attia@riken.jp; Jens Domke, RIKEN Center for Computational Science, Kobe, Japan, jens.domke@riken.jp. 2026. Manuscript submitted to ACM Manuscript submitted to ACM

1

2

Elena Aschmann, Anshu Dubey, Mohamed Wahib, and Jens Domke

them in the solver. Once properly verified, the resulting solver is lowered to concrete realization deterministically with tools also built with LLM assistance. The driver of the numerics and the code for interfacing with the application are generated by the LLMs. The resulting solver is exercised within Flash-X and with a separate environment we built for this experiment. Our results show that this is a viable approach for producing correct and performance portable code. CCS Concepts: • Software and its engineering → Software architectures; • Applied computing → Physics; Additional Key Words and Phrases: scientific software, application-defined abstractions, software architecture, data-layout abstraction, performance portability, AI-assisted refactoring, legacy code modernization

1

Introduction

For domain scientists and applied mathematicians working in high-performance computing (HPC) environments the holy grail has been to express their algorithm in purely mathematical notation that can be used to generate code that performs well across a range of platforms. The increasing sophistication of large language models (LLMs), and agents that drive them, have brought us close to achieving this goal. Figure 1 summarizes the premise and the process involved. Blue box shows how an algorithm might be presented in mathematical notation in a paper. The green box describes step-by-step arithmetic needed to implement the algorithm – a representation preferred by the domain scientists and computational mathematicians. This description is implementation independent. The yellow box shows how the representation in green box can be converted to a code with embedded abstractions that are solver or application specific. These abstractions, and the tools for lowering them to target language of the code are developed with LLM assistance. The pink box shows what corresponding traditional code might have looked like, highlighting portions that are core arithmetic versus surrounding code. The lavender box shows the core arithmetic that is translated using a verified tool and the variants of the surrounding code that are generated using LLMs. Our ultimate goal is to go from blue box to lavendar box, but for this work we begin from the green box. We start with a mature solver, extract its arithmetic into mathematical notation and use that as the basis for generating a “new” implementation. We take this approach because with a well tested solver we have a way to validate our hypothesis. We use Spark [6], a solver for compressible hydrodynamics, embedded in a large multiphysics software system, Flash-X [10] for the experiment. We first extract the description of the algorithm in a detailed mathematical notation similar to what is shown in the green box in figure 1. Then we convert this representation of the core arithmetic into a succinct and reusable form that is independent of the surrounding infrastructure necessary to facilitate the computation. Components of such infrastructure include data representation and movement, coordinating interactions among physics components, exposing parallelism, managing communication, and accommodating the constraints of different computing architectures. We refer to this infrastructure as surrounding code and distilled numerical component of the code as core arithmetic from here on for ease of explanation. With this exercise we wished to explore how far could we push code generation with LLMs to produce reliable code suitable for different platforms. Our finding is that it is possible to design, build and maintain application-specific abstractions that distance the core arithmetic from surrounding code up to the degree of intrusiveness acceptable to the developers. Such a separation has long been desirable for maintainability and portability, and several efforts over the years have achieved it to varying degrees [2, 3, 5, 8, 11, 14, 21, 22]. Fully realizing this separation has not been possible because application-specific abstractions were not viable as long as every line of code had to be written by humans. The breath and depth of expertise necessary to design, develop, and maintain abstractions is rare enough that groups engaged in these activities had to necessarily target a broad range of applications. The non-specificity resulting from this constraint Manuscript submitted to ACM

Application-Defined Abstractions to Separate Numerical Intent from Computation Orchestration

3

led naturally to libraries that are large and complex and require significant resources to maintain and keep up to date with evolution in technology. From the applications perspective also these solutions only provide moderate benefits. They have their own type system, data structures, and programming idioms which the application must adopt, and even then obtaining performance requires the application structure to be suitable for the target architecture. AI-assisted code development [12, 16, 17] could change this trade-off. However, even a fully validated specification cannot simply be handed over to an LLM to produce concrete code directly because of the fundamental limitations in the abilities of the models to guarantee correct outcome. The key, therefore, is to convert the mathematical notation to application specified abstractions and use LLMs to build tools that can lower the abstractions to a concrete realization in the target language. Once the abstractions and tools have been validated the translation is deterministic and by definition guaranteed to be correct up to the rigor of the verification applied to the translation. Because the abstraction boundary is controlled by the application rather than by an external framework, it can evolve with the application and can be introduced incrementally into an existing code base. In this paper, we develop and demonstrate this approach with Spark and use two different sets of abstractions to demonstrate the range and flexibility of the approach. Computational kernels are expressed using application-defined abstractions in both versions, though at different levels of intrusiveness into the original code. The contributions of this work are: • Defining abstraction boundary that guarantees reliable outcome with LLM assisted code generation. • A stable, infrastructure-agnostic core arithmetic component where explicit data references can be replaced by symbolic, parameterized references, creating a maintainable scientific artifact. • AI-assisted transformation of existing scientific software in a workflow where use of generative techniques is shifted from repeatedly generating scientifically sensitive code to performing constrained, verifiable source transformations. • Demonstration of correctness across independently developed infrastructures and performance on different architectures. 2

Premise

Large scientific simulation codes differ in important ways from conventional application software in that scientific knowledge is encoded in the implementations whose correctness has been established over many years of scientific use. This is the primary reason why LLM assisted development has not made as many inroads in scientific computing as it has elsewhere. However, the scientifically meaningful core arithmetic of the solvers on which the reliability of the obtained results hinges is typically a small fraction of the overall code. Therefore, if this core could be effectively isolated and made immutable then LLMs can be used for generating the surrounding code. Given the lack of guaranteed correctness of LLM generated code, any adoption of generative practices must be accompanied by rigorous verification. The challenge has been that the organization of data is usually embedded in numerical source code through explicit array subscripts and derived-type component selections. Consequently, a change in the data layout requires changes throughout the code, therefore, any LLM generated surrounding code would also entangle the arithmetic. The premise behind our experiment is that completely detaching the core arithmetic from the surrounding code is possible and the exercise of detachment will result in a precise elucidation of the scientific core of the code without being contaminated by details of the infrastructure that frames it. This can be done by defining a vocabulary of symbols that have specific meanings within the context of the application. Symbols can be as limited as data references only, or they can represent Manuscript submitted to ACM

4

Elena Aschmann, Anshu Dubey, Mohamed Wahib, and Jens Domke

arbitrarily sized code-snippets with parameterization. One then build tools to lower these symbols into the concrete source-level expression required by a particular environment. Once verified and validated, this representation becomes the enduring and stable scientific representation of the application. The resulting development model differs from both conventional abstraction libraries and unconstrained AI code generation. The scientifically significant arithmetic remains an explicit, maintained software artifact. The application itself determines the abstraction, AI assists in creating and manipulating this representation, while deterministic code generation lowers the abstractions into infrastructure-specific source code. This division provides the basis for the approach demonstrated in the remainder of this paper. 3

Instruments of Experimentation

We give a brief outline of the solver and the code used in this study to facilitate the description of the methodology. FlashX is a component-based software system for building multiphysics simulation applications where different permutations of code components generate specific instantiations capable of simulating the requisite physics. Components can have multiple alternative implementations to account for either different physical regimes or different fidelity requirements. It can be run with a uniformly discretized mesh (UG) or with adaptive mesh refinement (AMR) where higher resolution is used in the mesh only where necessary. Flash-X has a customized macro-expander [9] that became the starting point for the abstractions derived in this work. These custom macros can have multiple alternative definitions which allow them to unify code variants into one source. The expander has a built-in ability to arbitrate on which definition to use where to meet the requirements of various configurations. Spark solves the equations of ideal MHD using finite-volume discretization and explicit time integration, although in this work we disable the magnetic field component of the solver and use only the hydrodynamics. The solver begins by performing high-order spatial reconstruction on the cell-averaged primitive variables to estimate their values at the interface between cells. Running offset stencils for this reconstruction results in two approximations for the vector of primitives at each interface. These left and right states at each interface are used as input to a Riemann solver to compute the average flux across each interface between cells. Then the interface fluxes are used to compute conservative flux divergences for the update of the field variables. Spark uses second- or third-order strong stability preserving (SSP) Runge-Kutta (RK) time integration [13, 26]. This approach requires multiple sub-stages in which the vector of conserved quantities is updated requiring that the general procedure outlined above be repeated two or three times depending on whether second- or third-order time integration is used. The default for reconstruction is fifth-order weighted essentially non-oscillatory (WENO) [25] with the smoothness indicators of [4], and the HLLD approximate Riemann solver [20]. These are the calculations that form the core arithmetic of the solver that we extract and to which we apply the abstractions. We use three different applications in this study that are considered standard problems for exercising hydrodynamics solvers. They are, Sedov blast [23], Sod shock tube [28] and IsentropicVortex [24] shown in Figure 2. 1 . The Sedov blast problem initializes the domain with a pressure spike where the shock wave moves out radially. Sod produces all three types of nonlinear hydrodynamic waves (shock, contact, and rarefaction) associated with a shock, and isentropic vortex has a smooth profile with no shocks. All three applications use ideal gas equation of state. We used both commercially hosted and local institution hosted models for this work. The commercially hosted model were GPT-6.1-sol light with Codex as the agent and Opus-5.5 Medium. Institutionally hosted models included 1We reproduce the visualizations from Flash-X User’s Guide to give the reader an idea of what these problems are like, these are not produced from our

experiments Manuscript submitted to ACM

Application-Defined Abstractions to Separate Numerical Intent from Computation Orchestration

Sedov blast application is initialized with a pressure spike at the center of the domain. The shock travels out in a radial direction as shown above. This figure shows pressure after the solution has evolved for 0.05 time units.

Sod shock tube is initialized by establishing a planar interface at 45 degrees angle to the axes. The fluid is initially at rest on either side of the interface, with the density and pressure jumps. The shock, travels to the right while the rarefaction wave travels to the left.

5

Isentropic vortex problem has a smooth flow-field (there are no shocks or contact discontinuities) and contains no steep gradients. Initialization places the vortex at the center of the domain. This figure shows evolution after 100 time units.

Fig. 2. General visualizations of Sedov blast, Sod-shocktube and Isentropic-vortex problems after they have evolved for a while

DeepSeek-V4 and GLM-5.3. Our objective in this exercise is not to focus on any one model. We list all the models that we used for transparency and completeness. 4

Abstracting Data Layout

The most obvious abstraction, and the one embraced by various C++ based performance portability solutions such as Kokkos [11, 22, 30], is to eliminate explicit data references from the arithmetic. We use a macro-like mechanism explained in [9] for these simplest abstractions . Similar to macros, these symbols are strings with parameters that are replaced with their pre-defined expansions by the translator. For example, the code below can be translated into several different versions depending on how the symbols spc, spc_off, field, and field_off are defined. rho(@spc(i,j,k)@) = U(@field(DENS,i,j,k)@) flux(@spc(i,j,k)@) = & F(@field_off(DENS,i,+1,j,+0,k,+0)@) - & F(@field_off(DENS,i,-1,j,+0,k,+0)@) Table 1 lists symbol definitions and their corresponding expansions into Fortran code. The first three columns show the symbol, its arguments and how it is referenced in the code snippet. The last two columns show how changing the definition changes the code translation and what goes into the source code. The juxtaposed arguments are concatenated so spc_off(i,+1,j,+0,k,+0) yields i+1,j+0, k+0. These examples illustrate that the symbolic representation preserves the logical identity and relative offsets of data references, while the mapping from logical indices to their physical ordering is determined entirely by how the symbols are defined. The entire transformation can be verified by asking for a set of definitions that are compatible with the original data layout of the code. Lowered to the level of assembly, the original code and the translated code with embedded symbols match exactly. Formal methods can also be applied to verify the Manuscript submitted to ACM

6

Elena Aschmann, Anshu Dubey, Mohamed Wahib, and Jens Domke

semantic equivalence of any alternative set of definitions representing a different data layout [27], though we did not use them here because the variants are simple enough to produce bit-wise identical outputs. Table 1. Translation of application-defined symbols using two alternative definitions. For each symbol, the first row shows the original spatial index ordering and the second row shows the reversed ordering. Symbol

Arguments

Reference

spc

ib,jb,kb

@spc(i,j,k)@

spc_off

ib,io,jb,jo,kb,ko

@spc_off(i,+1,j,-1,k,+0)@

field

var,ib,jb,kb

@field(DENS,i,j,k)@

field_off

var,ib,io,jb,jo,kb,ko

@field_off(DENS,i,-1,j,+0,k,+0)@

Definition ib,jb,kb kb,jb,ib ibio,jbjo,kbko kbko,jbjo,ibio var,@spc(ib,jb,kb)@ var,@spc(ib,jb,kb)@ var,@spc_off(ib,io,jb,jo,kb,ko)@ var,@spc_off(ib,io,jb,jo,kb,ko)@

Translation i,j,k k,j,i i+1,j-1,k+0 k+0,j-1,i+1 DENS,i,j,k DENS,k,j,i DENS,i-1,j+0,k+0 DENS,k+0,j+0,i-1

To apply this methodology to Flash-X we used the already existing macro expansion mechanism reported in [9]. We needed an alternative implementation of the mesh to test the mutability of the data layout handled in the solver. With AI assistance we generated an alternative implementation of Uniform Grid in the code which changes the order of indices from variable first to variable last (see section 6 for more details). We verified correctness by comparing the output against that produced by the unchanged solver running on the unchanged Uniform Grid. 5

Implementation from Mathematical Notation

In the next stage of exploration we pushed the abstractions deeper into the code. Here, we completely detach core arithmetic from the surrounding code. Figure 3 summarizes the steps in the approach. We start with a reference implementation of the solver in Flash-X and use an LLM to translate each mathematically relevant step of the algorithm into a corresponding mathematical nota-

Recover numerical intent

Legacy code → detailed mathematical algorithm Equations independent of data references

tion. Note that these are not the equations that one would

Equations -> symbolic numerical implementation Write tools for lowering to compilable code

encounter in an algorithm paper. They detail all the compu-

Instantiate application

tations that must take place in the algorithm down to the smallest detail in a step-by-step exposition of the arithmetic in the solver. We need this degree of detail to fully verify the new implementation against the reference implementation.

LLM Generated, human inspected

Define & lower abstractions

Generate surrounding code Integrate numerical kernel and generated infrastructure 03

Verify & validate

Compare against the original implementation Equivalence tests + scientific validation

Done by humans with LLM Assist

A direct implementation from this derivation can be verified against the reference implementation, and any mismatch in the derivations can be corrected in the specification. Starting from this description it becomes possible to iden-

Fig. 3. Steps involved in algorithmic decomposition for separation of concerns, including the roles apportioned to AI and human scientists.

tify reusable patterns in the solver and encode them into a condensed vocabulary of symbolic representations specific to the solver. One needs to create tools to lower the symbolized representation into a concrete form in the target language for compilation, and this is where the availability of highly capable LLMs changes the cost-benefit equation. As mentioned earlier, the existing portability solutions were constrained by the limitation of resources in the form of human developers to even entertain the idea of specialized symbolization most suitable for a single solver or application. They had to necessarily design their solutions to be Manuscript submitted to ACM

Application-Defined Abstractions to Separate Numerical Intent from Computation Orchestration

Invocation in code @hy_assembleFaceStates(work%mi nus,work%plus,block,i,j,k,e,ql,qr)@

Definitions hy_assembleFaceStates args=minus,plus,block,i,j,k,offset,ql,qr scope=block declarations= integer :: v

definition= do v=1,7 ql(v)=@cell_ref(plus,v,block,ioffset(1),j-offset(2),k -offset(3))@ qr(v)=@cell_ref(minus,v, block,i,j,k)@

end do

cell_ref args=data,var,block,ix,jy,kz definition=${data}(${var},${ix},${jy},${kz})

7

Expansion macro_scope_19: do associate(macro_scope_19_arg_minus => scr_uMinus, macro_scope_19_arg_plus => scr_uPlus, & macro_scope_19_arg_block => block, macro_scope_19_arg_i => i, macro_scope_19_arg_j => j, & macro_scope_19_arg_k => k, macro_scope_19_arg_offset => e, macro_scope_19_arg_ql => ql, & macro_scope_19_arg_qr => qr) do macro_scope_19_local_v=1,7 macro_scope_19_arg_ql(macro_scope_19_local_v)= & macro_scope_19_arg_plus(macro_scop_19_local_v, & macro_scope_19_arg_i-macro_scope_19_arg_offset(1),& macro_scope_19_arg_j- macro_scope_19_arg_offset(2), & macro_scope_19_arg_k-macro_scope_19_arg_offset(3)) macro_scope_19_arg_qr(macro_scope_19_local_v)=& macro_scope_19_arg_minus& (macro_scope_19_local_v,macro_scope_19_arg_i, & macro_scope_19_arg_j,macro_scope_19_arg_k)

end do end associate exit macro_scope_19

Fig. 4. An Example of a computational helper. The top left box shows how it is invoked in the code, the bottom left box shows the definitions of the macros used and the right box shows the expanded code.

general across a variety of solvers. The LLMs have largely eliminated this trade-off to the point where tools for lowering abstractions to a compilable language can be generated and verified even by developers with no compiler expertise. As with any concrete implementation of a numerical algorithm, there is no unique way to design the vocabulary and there is no binding to any specific language when equations are the starting point. We opted to use macros for symbolizing data references as with the first example described in Section 4, but we extended their role and also turned them into computational helpers. These computational helpers are templates that encode the algorithmic sequence of the solver with insertion markers for imports and declarations in the enclosing specification section. Helper definitions have a name, formal arguments, declarations when needed, and an executable body. Each of the helper also implies a scope that effectively becomes a unique namespace for the translator. Figure 4 shows an example of a helper that assembles face states. Selected color coded symbols in the figure show how the symbols are placed in the maintained abstracted code, how they are defined, and what do they look like when translated. Note that the translated symbols look very similar to how the linkers resolve symbols. The translator essentially performs a similar function, although at a much more limited scale. The expander hoists generated declarations to these markers, placing imports where appropriate in the generated code. Scratch storage is caller-owned and reusable across calls. Layout-specific allocation is performed by the caller through backend macros. Mesh services are not called directly: the whole-step (including all stages in a multi-stage method) driver receives a provider context and the procedures that are needed to perform these services. Manuscript submitted to ACM

8

Elena Aschmann, Anshu Dubey, Mohamed Wahib, and Jens Domke

6

Verification

We evaluate the abstracted code from two perspectives: verification – whether the arithmetic layer produces the same numerical results regardless of who owns the data and drives the time integration; and performance – whether the generated surrounding code can be tailored for a target platform to obtain good performance. 6.1

Verification of symbolic data-reference

To verify the abstractions where only the explicit data-references are symbolized we use only the Sedov problem. Because the abstractions are shallow, verification with any one problem suffices. We use two alternative sets of definitions of macros that change the order of indices in the state data layout. The first set of definitions assumes the data layout to be 𝑈 (𝑓 𝑖𝑒𝑙𝑑_𝑣𝑎𝑟𝑖𝑎𝑏𝑙𝑒, 𝑖, 𝑗, 𝑘) where 𝑖, 𝑗, 𝑘 are space indices and 𝑓 𝑖𝑒𝑙𝑑_𝑣𝑎𝑟𝑖𝑎𝑏𝑙𝑒 represents physical variables (i.e. density, pressure etc). For the same set of indices the second set of definitions assumes the layout to be 𝑈 (𝑖, 𝑗, 𝑘, 𝑓 𝑖𝑒𝑙𝑑_𝑣𝑎𝑟𝑖𝑎𝑏𝑙𝑒). We refer to these definitions as variable-first and variable-last from here on. Note that Flash-X already has a built-in capability to change the order of indices during configuration because it supports two mesh packages, Paramesh [18], which assumes variable-first order, and AMReX [32], which assumes variable-last order. For the purpose of verifying the symbolic data referencing we disabled the native reordering capability in Spark and used the LLM to replace all explicit array references in the solver with appropriate symbolic representations. This set of verification runs used a Mac M5 laptop with Tahoe 26.5 operating system. The existing mechanism for

Version 1 Reference implementation Explicit array references Native Uniform Grid implementation Variable-first data layout

selection of alternative implementation in Flash-X was used to generate three versions of the problem shown in the figure 5. The first two versions have identical data layout. Their equivalence was verified by confirming that the assembly code

Version 2 Symbols with variable-first definitions Native Uniform Grid implementation Variable-first data layout

produced by the compiler was identical for all the relevant files. The equivalence of the second and the third versions was verified by confirming that the produced results were bit-wise Sedov shock application is initialized

Version 3 Symbols with variable-last definitions Reordered Uniform Grid implementation Variable-last data layout

with a pressure spike at the center of identical, and also through visualization of the results. the domain. The shock travels out in a

6.2

radial direction as shown above. This figure is shows pressure after the solution has evolved for 0.05 time Verification for Arithmetic Extraction units.

Verification by ensuring that compiler generated assembly is identical Verification through bitwise equivalence and visualization of results

This set of experiments probe the efficacy of separating the core arithmetic of the solver from the surrounding code. We

Fig. 5. Sedov blast problem run with three different configu-

refer to the version of the code with its core arithmetic ex- rations: the reference, abstracted with layout identical to the tracted as NewImpl for convenience. The questions we ask

reference, and abstracted with a different data layout

are: Q1 Does NewImpl reproduce the results of the reference implementation when it is inserted into Flash-X? Q2 Does NewImpl produce the same results when its data are owned and its time loop is driven by an independently developed infrastructure instead of Flash-X? Q3 How much code has to be written to attach an independent infrastructure, and does any of it touch the arithmetic layer? Q4 Does the arithmetic layer give the same results on a GPU? Manuscript submitted to ACM

Application-Defined Abstractions to Separate Numerical Intent from Computation Orchestration Config. 3 independent

Config. 4 bridge

Config. 5 GPU

Flash-X

standalone (no Flash-X code)

Flash-X

Flash-X

Spark

NewImpl

NewImpl

Spark

Config. 1 reference

Config. 2 native

host program

Flash-X

time-step driver

Spark Spark kernels

NewImpl kernels

NewImpl kernels

NewImpl kernels

NewImpl kernels (OpenMP target)

mapping

none (hard-coded)

wrappers + Flash-X layout

UMesh layout + grid provider

UMesh layout + grid provider

batching adapter + Paramesh layout

data owner

Flash-X

Flash-X

UMesh

UMesh (copied in/out)

Paramesh4 (packed in batches)

arithmetic

9

Q2: bitwise identical Q1: NewImpl vs. Spark on identical setups Q4: GPU vs. Spark on the same mesh and compiler Flash-X / Spark (original)

NewImpl (unchanged)

written to attach the mesh or GPU

independent mesh library

Fig. 6. Evaluated configurations. "Data owner" is the software that allocates and lays out the cell data; "time-step driver" is the code that sequences Runge–Kutta stages, guard-cell exchange, and flux synchronisation within a time step.

We compare five configurations described in Figure 6. They differ in who owns the cell data (and thus also the layout in memory), who drives the time step, and which hydrodynamics implementation performs the arithmetic. Here, the reference configuration, config 1, is the unmodified Spark solver running with Flash-X’s native UG. The second configuration, config 2, embeds the core arithmetic from NewImpl in Spark’s own driver through compatibility wrappers that we developed with LLM assistance. We run this configuration on Flash-X’s single UG block, like config 1. The third configuration, config 3, runs with core arithmetic of NewImpl with independent infrastructure. Here, we use a mesh developed as a C library with Fortran binding that divides a uniform Cartesian domain into equal blocks with halos, and whose memory layout, dimension, and halo width are compile-time macros. NewImpl’s whole-step driver advances the solution and UniformMesh supplies the data and the mesh services through a grid provider. The host program contains no Flash-X code. It reproduces Flash-X’s initial conditions, time-step control, and output formats from Flash-X’s sources, so that its results can be compared with config 1 and config 4 using Flash-X’s own tools. In the fourth configuration, config 4, we embed the whole of NewImpl and the independently developed mesh as a sub-unit of Spark in Flash-X. Here, Flash-X initializes the problem and controls the time step; in each step its cell data are copied into UniformMesh, advanced by NewImpl, and copied back. The configuration allows NewImpl and Spark to be verified with bit-wise identical results. The final configuration, config 5, uses OpenMP target offload to enable running the code on the GPU. In this configuration, NewImpl’s core arithmetic is embedded into the single stage face-flux calculation routine of Spark. Unlike all previous configurations, here the single stage driver for face-flux calculation needs to be modified to allow batching of blocks for offloading. This is because in AMR mode an individual block is too small for efficient use of the GPU. Several blocks must be computed in parallel on the GPU for it to be viable in production. To verify this version, the same generated source is also built for host execution (OpenMP on the CPU, "host kernels"). The mesh is a single refinement level of 4 × 4 blocks of 642 cells, so it has the same 2562 resolution as the other configurations, with the periodic boundaries and flux correction. We verify this configuration against the reference configuration using with the same compiler on the same mesh. For all of these verification runs NewImpl and Spark are taken from the same commit of Flash-X; NewImpl’s sources are not modified in any configuration, which each build verifies by recording checksums of all NewImpl inputs. All codes that have CPU only target are compiled with GNU Fortran 14.3.1 and Flash-X’s flags (-O2, 8-byte reals) and run Manuscript submitted to ACM

10

Elena Aschmann, Anshu Dubey, Mohamed Wahib, and Jens Domke

on a single core of an AMD EPYC 7302 node, one MPI rank. Config 5 and its Spark reference are compiled with NVHPC 26.3 (-O2 -r8, -mp=gpu -gpu=cc128) and run on an NVIDIA GeForce RTX5080 attached to an Intel Xeon Max 9468, with one MPI rank. Performance measurements on the GPU were conducted on a system with NVIDIA GB200 NVL4 with NVHPC 26.3 compiler suite. 6.2.1 Details of the runs. We use Sod shock tube and isentropic vortex problems for the CPU only runs. All variations of these problems run in two-dimensions with periodic boundaries on 256 × 256 cells. IsentropicVortex was setup on the physical domain sized at [0, 10] 2 , advected diagonally for one period (𝑡 = 10), where the exact solution is the translated initial condition. The Sod setup used a physical domain of size [0, 1] 2 with the interface normal to 𝑥, run to 𝑡 = 0.025. With periodic boundaries a second, mirror-image Riemann problem forms at 𝑥 = 0, so the exact solution is the superposition of two exact Riemann solutions until the wave systems meet, and it is symmetric about 𝑥 = 0.25. A second variant (Sod-hybrid) switches to the HLLE solver in shock-flagged cells. Two further variants switch flattening off in both solvers (Sod-noflat, Sod-hybrid-noflat), to identify differences at shocks. For the order of accuracy, the vortex is also run on 642 , 1282 , and 5122 cells with the remaining parameters left identical to the other runs. 1.84 × 10−5 / 4.0 × 10−4 Two vortex variants use a fixed Δ𝑡 = 5 × 10−3 after Flash-X’s start-up step (0.1 of the CFL limit). One runs a full period (vortex-fixed, 2001 steps), the other 100 steps (vortex-short). We run config 3 and config 4 with several decompositions of the domain into blocks of at least 8 × 8 cells, from strips (4 × 1, 1 × 16) to 32 × 32 blocks. Each problem uses two to six of them and each run is compared to the single block run of the same problem. 20 comparisons per configuration, 16 of them on the five main problems. The time per step is measured on vortex-short (100 steps) and Sod (30 steps), five times each on an otherwise idle node; we report the minimum. To determine which of the Flash-X behaviours reproduced by config 3 bit-wise identity depends on, configuration 3 is rerun with one of them replaced at a time and compared with configuration 4 on vortex-short, Sod, and the vortex: Flash-X’s equation of state by the grid provider’s simpler one, or without recomputing the internal energy; no start-up step; no limit on the growth of the time step; the last step shortened to the final time; and NewImpl’s own CFL formula instead of Spark’s. To exercise storage layouts NewImpl’s own example, a one-dimensional periodic Sod problem on its built-in grid provider, is built with each of NewImpl’s three storage mappings (variable-first, spatial-first, and per-cell records) and run with SSP-RK2 and SSP-RK3 on one and four blocks; the outputs are compared byte by byte. 6.2.2

Metrics. • bit-wise identity of checkpoints, compared with Flash-X’s sfocu tool at tolerance zero, and of the time series of integral quantities (which include the simulation time and thus the time-step sequence). • Difference between two configurations: the maximum and 𝐿1 norm of the cell-wise difference of all variables (from sfocu), and its ratio to the discretisation error. • Error against the exact solution: 𝐿1 and 𝐿∞ norms of the density error, with the exact solution evaluated at each checkpoint’s own time. • Mirror symmetry of the Sod solution: max𝑥 |𝑢 (𝑥) + 𝑢 (𝑥 ′ )|, where 𝑥 ′ is the mirror image of 𝑥 about 𝑥 = 0.25; it is zero for the exact solution. • Number of time steps and the time-step sequence. • Lines of code that attach the infrastructure, and lines of the arithmetic layer changed for it. • Conservation: the largest relative change of the total mass, momentum, and energy over a run.

Manuscript submitted to ACM

Application-Defined Abstractions to Separate Numerical Intent from Computation Orchestration

11

Table 2. NewImpl on UniformMesh in Flash-X (config 4) versus the original Spark solver (config 1); 256 × 256 cells, one rank. Differences are maxima over all cells and variables; errors are density errors against the exact solution (𝐿1 / 𝐿∞ ). vortex-short: 𝐿1 difference 2.5 × 10 −15 . Problem vortex-short vortex-fixed vortex (CFL) Sod Sod-hybrid

Steps

|NewImpl − Spark| max

NewImpl error

Spark error

100 2001 1877 30 30

1.9 × 10 −12

2.25 × 10 −6 / 8.9 × 10 −5

1.3 × 10 −8 2.5 × 10 −8 1.6 × 10 −13 1.6 × 10 −13

1.84 × 10 −5 / 4.0 × 10 −4 1.84 × 10 −5 / 4.0 × 10 −4 3.78 × 10 −3 / 8.47 × 10 −2 3.79 × 10 −3 / 8.54 × 10 −2

same same same same same

Table 3. Ablation of configuration 3: maximum difference to configuration 4 in the hydrodynamic variables when one reproduced Flash-X behaviour is replaced; – means bit-wise identical. Replaced behaviour

vortex-short

Sod

vortex

Equation of state of the grid provider No recomputation of the internal energy No start-up step No limit on time-step growth Last step shortened to final time NewImpl’s CFL formula

3.0 × 10 −12

1.2 × 10 −13

4.4 × 10 −13 3.3 × 10 −3 – – –

3.6 × 10 −14 – 1.1 × 10 −2 1.0 × 10 −1 –

2.7 × 10 −8 2.3 × 10 −8 7.4 × 10−4 1.6 × 10 −3 3.8 × 10 −3 2.8 × 10−8

• Time per step of the hydrodynamics: Flash-X’s Hydro timer for configurations 1 and 4 (for configuration 4 including the copies into and out of UniformMesh), and NewImpl’s whole-step driver for configuration 3; and the wall time of the whole run per step. 6.3

Verification Results

We begin by comparing the reference solver with core arithmetic of NewImpl (config 1 andconfig 4), summarized in Table 2. For smooth flow after 100 identical steps the two solvers differ by at most 1.9 × 10−12 in single cells, with an 𝐿1 difference of 2.5 × 10−15 , which can be considered as rounding errors. Over a full period the maximum difference grows to 1.3 × 10−8 with a fixed time step and 2.5 × 10−8 with the CFL-limited timestep. That is 0.003% and 0.006% of the maximum error against the exact solution. Even with the CFL-limited timestep, both solvers choose identical time-step series and have the same error against the exact solution. For smooth flow NewImpl reproduces Spark, including its time-step control, to within rounding effects (Q1). On the Sod problem, NewImpl reproduces Spark to rounding. The maximum difference is 1.6 × 10−13 (with and without the hybrid HLLE solver). Both solvers take the same 30 time steps, have the same error against the exact solution and preserve the mirror symmetry of the periodic problem to max |𝑢 (𝑥) + 𝑢 (𝑥 ′ )| = 10−13 . Shock detection, both Riemann solves, flattening, artificial viscosity and the update thus agree with Spark’s implementation. 6.3.1 Independent infrastructure (config 3 versus 4). Configuration 3 runs NewImpl on UniformMesh without any Flash-X code. On the same machine and compiler, its results are bit-wise identical to those of configuration 4 for all ten problems: every checkpoint (11 for each vortex run over a full period, compared with sfocu at tolerance zero), the time series of integral quantities including the simulation time, and therefore the entire time-step sequence. This answers Q2 – the result of the arithmetic does not depend on the host program. Mimicking Flash-X’s initialization, time-step control and output by an independent implementation produces bit-wise identical results. Manuscript submitted to ACM

12

Elena Aschmann, Anshu Dubey, Mohamed Wahib, and Jens Domke

Table 4. NewImpl on a GPU (configuration 5) versus the original Spark solver, both built with NVHPC 26.3 on the same Paramesh4 mesh (4 × 4 blocks of 642 cells, one level); one rank. Differences are maxima over all cells and variables of the final checkpoint; 𝐿1 is the mean cell-wise difference. ”Contraction off”: both codes built without fused multiply-add contraction. Last column: the GPU run compared with the same generated kernels executed on the host, both without contraction. Problem

Steps

|GPU − Spark| max

𝐿1

contraction off GPU vs. host

vortex-short vortex-fixed vortex (CFL) Sod Sod-hybrid

100 2001 1877𝑎 30𝑎 30𝑎

5.2 × 10 −12 1.4 × 10 −7 6.9 × 10 −8 3.4 × 10 −13 2.3 × 10 −13

4.0 × 10 −15 2.3 × 10 −9 1.0 × 10 −9 7.1 × 10 −15 5.8 × 10 −15

3.8 × 10 −12 1.9 × 10 −7 8.2 × 10 −8 2.4 × 10 −13 9.3 × 10 −14

bit-wise bit-wise bit-wise bit-wise bit-wise

𝑎 Same number of steps as Spark; the CFL time steps differ only in the last digits (end times by at most 2.5 × 10 −12 for the vortex and 10 −17 for Sod).

6.3.2 Native environment (Config 2 vs config 1). Inside Spark’s own driver, NewImpl’s kernels reproduce Spark to rounding on every problem. After a full vortex period with a fixed time step the maximum difference is 1.7 × 10−12 , on the Sod problems 3.0 × 10−13 (2.4 × 10−13 with the hybrid solver), with identical step counts. The only larger difference, 3.0 × 10−8 for the CFL-limited vortex, arises because of rounding form the first steps on, so the time steps computed from them differ in the last bits and this error grows to the same order as for configuration 4 (2.5 × 10−8 ). With a fixed time step the difference stays at 1.7 × 10−12 The arithmetic layer thus behaves identically in its native environment and in the independently developed infrastructure. 6.3.3 Different architecture (Config 5 vs config 1). Table 4 compares the GPU configuration with Spark built with the same compiler on the same mesh. On every problem both take the same number of steps. With a fixed time step they take the identical time-step sequence. On Sod the maximum difference is 3.4 × 10−13 . Mirror symmetry is preserved to 3 × 10−13 , the same level as Spark itself. After 100 vortex steps the difference is 5.2×−12 . To separate the kernels from the hardware, the generated kernels were also run on the host. With contraction disabled, the GPU results are bit-wise identical to the host results on all five problems. The results are also independent of the batch size (1, 3 and 32). 6.4

Effort to attach the infrastructure

Table 5 lists the code written to attach UniformMesh to NewImpl (Q3). The arithmetic layer, 1102 lines of NewImpl templates, helper definitions, and the grid-provider contract, is used without any changes. Two files connect NewImpl to UniformMesh. The storage mapping (61 lines) defines the field type and the cell_ref and the face_ref references as a pointer into UniformMesh’s buffer, so no cell data is copied. Scratch fields remain ordinary arrays. The grid provider (194 lines) implements the four services of the contract: guard-cell exchange through UniformMesh, flux synchronisation at block borders, the global reductions (identity for one rank), and the equation of state. Everything else in Table 5 belongs to a host program, not to the attachment (the Flash-X sub-unit of configuration 4 and the standalone program of configuration 3, most of which reproduces Flash-X’s initialisation, time-step control, and output formats). Running NewImpl on a GPU required 732 lines of new code (Table 5), none of which touches the arithmetic layer. 7

Performance

We analyzed the performance of the code from two perspectives, that of overheads and portability. Table 6 lists the time per step. Expressing the arithmetic through application-defined abstractions has no measurable cost. Within Flash-X, NewImpl takes 0.127 s per step on the vortex where Spark takes 0.140 s. Configuration 3 takes 0.111 s. The difference to Manuscript submitted to ACM

Application-Defined Abstractions to Separate Numerical Intent from Computation Orchestration

13

Table 5. Lines of code (including comments) per component. Only the first two rows are needed to attach UniformMesh to NewImpl; only the rows marked 5 are needed to run NewImpl on a GPU. The GPU kernels themselves (841 lines) are generated from the unchanged NewImpl templates and are not maintained. Component

Lines

Used by

Storage mapping (layout definitions) Grid provider (contract implementation)

61 194

3, 4 3, 4

GPU kernel template and device layout definitions GPU kernel generator Batching adapter (packing, transfers, launches, flux publication)

140 200 392

5 5 5

Flash-X sub-unit and installer Standalone host program (incl. checkpoint writer) UniformMesh library (C, Fortran binding, layout macros)

588 1767 884

4 3 3, 4

NewImpl arithmetic layer, lines changed

0 of 1102

3, 4, 5

Table 6. Time per step in seconds (minimum of five runs), 256 × 256 cells, one core. “Hydro”: Flash-X’s Hydro timer configurations 1, 4) or NewImpl’s whole-step driver (configuration 3); “wall”: whole run per step.

Configuration

vortex-short Hydro wall

1: Spark 0.140 4: NewImpl in Flash-X, 1 block 0.127 4: NewImpl in Flash-X, 4 × 4 / 4 × 1 0.131 3: NewImpl standalone, 1 block 0.111 3: NewImpl standalone, 4 × 4 / 4 × 1 0.117

0.159 0.146 0.150 0.114 0.121

Sod Hydro

wall

0.140 0.124 0.136 0.109 0.112

0.198 0.183 0.195 0.113 0.117

configuration 4 is not in NewImpl: configuration 4’s timer additionally covers copying the solution between Flash-X’s array and UniformMesh in every step, which the bridge configuration needs because both own a copy of the data, and Flash-X’s general equation-of-state interface, which gathers the cell data into separate vectors; configuration 3 evaluates the same equation of state in place. Outside the hydrodynamics, Flash-X adds 0.02–0.06 s per step for time-step control and output, the standalone host almost nothing. These are single-core measurements on one machine with different timer scopes. On a GPU the surrounding code decides how much work each kernel launch gets. Figure 7 shows the time per step on Sedov for batch sizes 1 to 64. A batch costs 1.29 ms with 1 block and 1.48 ms with 64 (Figure 8), so each launch is dominated by launch overhead rather than computation. Grouping blocks therefore cuts the GPU time almost in proportion to the batch size. The time for the whole hydrodynamics step falls from 358 ms to 30ms (11.9×). At that point it is limited by host work that batching does not reduce. About 10 ms for guard-cell exchange and about 9 ms for packing, unpacking, the equation of state and fluxes. The arithmetic layer is the same for every batch size and gives bit-wise identical results. These measurements clearly show that the surrounding code can be adapted to an architecture’s execution model. 8

Related Work

This work lies at the boundary of performance-portability abstractions for HPC applications where code-generation approaches separate numerical intent from implementation, and AI-assisted transformation of existing software. Extant Manuscript submitted to ACM

14

Elena Aschmann, Anshu Dubey, Mohamed Wahib, and Jens Domke 400 GPU kernels and transfers host packing, unpacking, EOS, fluxes guard-cell fill reflux and other

358

Hydro time per step (ms)

350 300 250 194

200 150

113 100

69

50

45

35

30

16

32

64

0 1

2

4

8

blocks per GPU batch (hy_gpuBatchSize) Fig. 7. Hydrodynamics time per step on the GPU for Sedov (256 blocks of 162 cells, 100 steps) versus the number of blocks per GPU batch, split into GPU kernels and transfers, host-side adapter work, guard-cell exchange, and the rest. (a)

(b) 56 £

ideal ( / batch size)

GPU kernels and transfers whole Hydro step

32 16

11:9 £ 8 4 2

1.8 1.6

GPU time per batch (ms)

speed-up over batch size 1

64

1.4

1.29

1.31

1

2

1.36

1.39

1.38

4

8

16

1.43

1.48

32

64

1.2 1.0 0.8 0.6 0.4 0.2

1

0.0 1

2

4

8

16

32

64

blocks per GPU batch

blocks per GPU batch

Fig. 8. (a) Speed-up over batch size 1 of the GPU kernels and transfers and of the whole hydrodynamics step. (b) GPU time per batch, almost independent of the number of blocks in it.

performance portability approaches have concerned themselves with providing abstractions for parallel execution, memory spaces, and data access while mapping these abstractions onto architecture-specific backends [29–31]. Directivebased approaches such as OpenMP [7] and OpenACC [15] target offloading allow architecture mappings to be expressed through annotations rather than by adopting a C++ abstraction library. The commonality among them is that the portability mechanism defines the abstractions available to the application. LLMs have introduced a new mechanism for automating software transformations that previously required substantial manual effort. However, the reliability of such transformations remains a central concern [1, 19]. This is particularly important for scientific software, where apparently small changes to arithmetic expressions, array indexing, or update ordering can silently alter numerical results. Manuscript submitted to ACM

Application-Defined Abstractions to Separate Numerical Intent from Computation Orchestration

15

Our study shows that a constrained use of LLMs with a robust verification framework can be used to shift the abstraction boundary to the application. The application defines a small symbolic vocabulary corresponding to its own existing computational and data-access constraints. Infrastructure-specific code generation then gives those symbols their concrete representation. Consequently, an application’s numerical implementation need not adopt the types, containers, or access syntax of a particular portability framework. It also differs from conventional source-to-source porting. Rather than transforming an application directly from one infrastructure-specific representation to another, there is an intermediate, application-controlled representation into the maintained source. Once this transformation has been performed and validated, subsequent infrastructure changes alter the lowering of the abstractions rather than the arithmetic that expresses the numerical method. 9

Conclusion

We have introduced a new way of thinking about abstractions where instead of conforming to the a general portability system’s requirements, the applications devise their own customized portability framework. Until the arrival of highly capable LLMs and agentic frameworks this wasn’t feasible. Each transformation had to be manually constructed and validated. The other extreme of using LLM to produce all the code is still not feasible because of reliability concerns. The central proposition of this work is not that scientific software should be regenerated whenever computing technology changes. Rather, the availability of reliable code-generation assistance allows us to reconsider where the stable boundary in scientific software should lie. By maintaining the scientifically meaningful arithmetic as a durable artifact and making its connection to infrastructure explicit through application-defined abstractions, it is possible to preserve validated scientific computation while allowing the surrounding software mechanisms to evolve with substantially less disruption. The GPU configuration tests this proposition directly. Moving the arithmetic to a GPU, and making it efficient there by batching blocks required only new surrounding code. The arithmetic layer stayed the same and reproduced the CPU results to rounding; with contraction disabled, the GPU and host executions of the generated kernels were bit-wise identical. We have exercised this proposition by introducing two levels of abstractions into an existing, scientifically validated implementation. One is a shallow abstraction that only affects array references leaving the remainder of the implementation unaltered, while the other is a fundamentally different way of implementing an algorithm from its detailed mathematical exposition. In the second approach the role of the reference implementation is limited. It is used to extract the detailed mathematical specification of the algorithm and as a baseline for verifying the generated code. Our ultimate objective is to use LLM assistance to go from the detailed mathematical description of the algorithm to generated executable code. Here we use extraction of mathematical notation from a mature solver because having a reference implementation allows us to validate the methodology. It is also worth noting that with this approach the degree of intrusiveness of the abstractions remains firmly in control of the application developer. References [1] 2026. Using LLMs to Enhance Code Quality: A Systematic Literature Review. Information and Software Technology 190 (2026), 107960. doi:10.1016/j. infsof.2025.107960 [2] Michael Bauer, Sean Treichler, Elliott Slaughter, and Alex Aiken. 2012. Legion: Expressing locality and independence with logical regions. In SC’12: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis. IEEE, 1–11. [3] Mauro Bianco, Lucas Benedicic, et al. 2020. GridTools. https://github.com/GridTools/gridtools [4] R. Borges, M. Carmona, B. Costa, and W. S. Don. 2008. An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws. J. Comput. Phys. 227 (March 2008), 3191–3211. Manuscript submitted to ACM

16

Elena Aschmann, Anshu Dubey, Mohamed Wahib, and Jens Domke

[5] Valentin Clement, Sylvaine Ferrachat, Oliver Fuhrer, Xavier Lapillonne, Carlos E Osuna, Robert Pincus, Jon Rood, and William Sawyer. 2018. The CLAW DSL: Abstractions for performance portable weather and climate models. In Proceedings of PASC. 1–10. [6] Sean M Couch, Jared Carlson, Michael Pajkos, Brian W O’Shea, Anshu Dubey, and Tom Klosterman. 2021. Towards performance portability in the Spark astrophysical magnetohydrodynamics solver in the Flash-X simulation framework. Parallel Comput. 108 (2021), 102830. [7] Leonardo Dagum and Ramesh Menon. 1998. OpenMP: An Industry-Standard API for Shared-Memory Programming. IEEE Comput. Sci. Eng. 5, 1 (jan 1998), 46–55. doi:10.1109/99.660313 [8] Zhe Dai, Liang Deng, YongGang Che, Ming Li, Jian Zhang, and Yueqing Wang. 2024. Evaluating Performance Portability of Five Shared-Memory Programming Models Using a High-Order Unstructured CFD Solver. J. Parallel and Distrib. Comput. 187 (2024), 104831. doi:10.1016/j.jpdc.2023.104831 [9] Anshu Dubey, Youngjun Lee, Tom Klosterman, and Emil Vatai. 2025. A tool and a methodology to use macros for abstracting variations in code for different computational demands. Future Generation Computer Systems 173 (2025), 106937. doi:10.1016/j.future.2023.07.014 [10] Anshu Dubey, Klaus Weide, Jared O’Neal, Akash Dhruv, Sean Couch, J. Austin Harris, Tom Klosterman, Rajeev Jain, Johann Rudi, Bronson Messer, Michael Pajkos, Jared Carlson, Ran Chu, Mohamed Wahib, Saurabh Chawdhary, Paul M. Ricker, Dongwook Lee, Katie Antypas, Katherine M. Riley, Christopher Daley, Murali Ganapathy, Francis X. Timmes, Dean M. Townsley, Marcos Vanella, John Bachan, Paul M. Rich, Shravan Kumar, Eirik Endeve, W. Raphael Hix, Anthony Mezzacappa, and Thomas Papatheodore. 2022. Flash-X: A Multiphysics Simulation Software Instrument. SoftwareX 19 (2022), 101168. doi:10.1016/j.softx.2022.101168 [11] H. Carter Edwards, Christian R. Trott, and Daniel Sunderland. 2014. Kokkos: Enabling Manycore Performance Portability through Polymorphic Memory Access Patterns. J. Parallel and Distrib. Comput. 74, 12 (2014), 3202–3216. doi:10.1016/j.jpdc.2014.07.003 [12] Angela Fan, Beliz Gokkaya, Mark Harman, Mitya Lyubarskiy, Shubho Sengupta, Shin Yoo, and Jie M. Zhang. 2023. Large Language Models for Software Engineering: Survey and Open Problems. In 2023 IEEE/ACM International Conference on Software Engineering: Future of Software Engineering (ICSE-FoSE). 31–53. [13] Sigal Gottlieb, David I. Ketcheson, and Chi-Wang Shu. 2009. High Order Strong Stability Preserving Time Discretizations. Journal of Scientific Computing 38, 3 (Mar 2009), 251–289. [14] Tobias Gysi, Carlos Osuna, Oliver Fuhrer, Mauro Bianco, and Thomas C. Schulthess. 2015. STELLA: A domain-specific tool for structured grid methods in weather and climate models. In SC ’15: Proceedings of the international conference for high performance computing, networking, storage and analysis. 1–12. doi:10.1145/2807591.2807627 [15] J. A. Herdman, W. P. Gaudin, O. Perks, D. A. Beckingsale, A. C. Mallinson, and S. A. Jarvis. 2014. Achieving Portability and Performance through OpenACC. In 2014 First Workshop on Accelerator Programming using Directives. 19–26. doi:10.1109/WACCPD.2014.10 [16] Xinyi Hou, Yanjie Zhao, Yue Liu, Zhou Yang, et al. 2024. Large Language Models for Software Engineering: A Systematic Literature Review. ACM Transactions on Software Engineering and Methodology 33, 8, Article 220 (2024), 79 pages. doi:10.1145/3695988 [17] Juyong Jiang, Fan Wang, Jiasi Shen, Sungju Kim, and Sunghun Kim. 2026. A Survey on Large Language Models for Code Generation. ACM Transactions on Software Engineering and Methodology 35, 2 (2026), 1–72. doi:10.1145/3747588 [18] Peter MacNeice, Kevin M Olson, Clark Mobarry, Rosalinda De Fainchtein, and Charles Packer. 2000. PARAMESH: A parallel adaptive mesh refinement community toolkit. Computer physics communications 126, 3 (2000), 330–354. [19] Sofia Martinez, Luo Xu, Mariam Elnaggar, and Eman Abdullah Alomar. 2025. Software Refactoring Research with Large Language Models: A Systematic Literature Review. Journal of Systems and Software (2025). doi:10.1016/j.jss.2025.112762 [20] T. Miyoshi and K. Kusano. 2005. A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics. J. Comput. Phys. 208 (Sept. 2005), 315–344. [21] Jonathan Ragan-Kelley, Connelly Barnes, Andrew Adams, Sylvain Paris, Frédo Durand, and Saman Amarasinghe. 2013. Halide: a language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines. ACM Sigplan Notices 48, 6 (2013), 519–530. [22] Xuanzhengbo Ren, Tetsuya Hoshino, Daichi Mukunoki, and Takahiro Katagiri. 2026. Porting NICAM Microphysics to Kokkos: A Case Study in Performance Portability for Atmospheric Models. Journal of Information Processing 34 (2026), 503–513. doi:10.2197/ipsjjip.34.503 [23] Leonid Ivanovich Sedov. 1993. Similarity and Dimensional Methods in Mechanics (10 ed.). CRC Press. 496 pages. doi:10.1201/9780203739730 [24] Chi-Wang Shu. 1997. Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws. Technical Report ICASE Report 97-65. Institute for Computer Applications in Science and Engineering, NASA Langley Research Center, Hampton, Virginia. NASA/CR-97-206253. [25] C.-W. Shu. 2009. High Order Weighted Essentially Nonoscillatory Schemes for Convection Dominated Problems. SIAM Rev. 51 (Jan. 2009), 82–126. [26] C.-W. Shu and S. Osher. 1988. Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes. J. Comput. Phys. 77 (Aug. 1988), 439–471. [27] Stephen F Siegel, Manchun Zheng, Ziqing Luo, Timothy K Zirkel, Andre V Marianiello, John G Edenhofner, Matthew B Dwyer, and Michael S Rogers. 2015. CIVL: the concurrency intermediate verification language. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. Association for Computing Machinery, New York, NY, USA, 1–12. doi:10.1145/2807591.2807635 [28] Gary A. Sod. 1978. A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws. J. Comput. Phys. 27, 1 (1978), 1–31. doi:10.1016/0021-9991(78)90023-2 [29] Jan Stephan, Sergei Bastrakov, Andrea Di Pilato, Sebastian Ehrig, Benjamin M. Gruber, Jan Vyskočil, René Widera, and Michael Bussmann. 2026. Performance Portability with Alpaka. Journal of Physics: Conference Series 3206, 1 (2026), 012054. doi:10.1088/1742-6596/3206/1/012054 Manuscript submitted to ACM

Application-Defined Abstractions to Separate Numerical Intent from Computation Orchestration

17

[30] Christian R. Trott, Luc Berger-Vergiat, David Poliakoff, Sivasankaran Rajamanickam, Damien Lebrun-Grandie, Jonathan Madsen, Nader Al Awar, Milos Gligoric, Galen Shipman, and Geoff Womeldorff. 2021. The Kokkos Ecosystem: Comprehensive Performance Portability for High Performance Computing. Computing in Science & Engineering 23, 5 (2021), 10–18. doi:10.1109/MCSE.2021.3098509 [31] Erik Zenker, Benjamin Worpitz, René Widera, Axel Huebl, Guido Juckeland, Andreas Knüpfer, Wolfgang E. Nagel, and Michael Bussmann. 2016. Alpaka – An Abstraction Library for Parallel Kernel Acceleration. In 2016 IEEE International Parallel and Distributed Processing Symposium Workshops. [32] Weiqun Zhang, Ann Almgren, Vince Beckner, John Bell, Johannes Blaschke, Cy Chan, Marcus Day, Brian Friesen, Kevin Gott, Daniel Graves, Max Katz, Andrew Myers, Tan Nguyen, Andrew Nonaka, Michele Rosso, Samuel Williams, and Michael Zingale. 2019. AMReX: a framework for block-structured adaptive mesh refinement. Journal of Open Source Software 4, 37 (May 2019), 1370. doi:10.21105/joss.01370

Manuscript submitted to ACM

Related documents

• HEL: A Heterogeneous Extensible Layer for GPU-Resident Runtime Orchestration#1389776
• HEL: A Heterogeneous Extensible Layer for GPU-Resident Runtime Orchestration#1381210
• Deterministic vs. LLM-Controlled Orchestration for COBOL-to-Python Modernization#175332
• Migrating a Legacy Scientific Application with an AI Coding Agent: NICAM from Fortran to Kokkos#1534494
• Migrating a Legacy Scientific Application with an AI Coding Agent: NICAM from Fortran to Kokkos#1514890