Pith. sign in

REVIEW 3 major objections 4 minor 20 references

Massively parallel numerical simulations with Julia

T0 review · 3 major / 4 minor · reviewed 2026-08-01 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read Julia, a high-level dynamic language, can run a production CFD solver on 61,440 CPU cores with close-to-ideal scaling.

desk verdict A useful and honest Julia scaling study that oversells its headline efficiency: the superlinear cache baseline inflates the 0.83 number, but the paper itself discloses the issue. read the letter →

arxiv 2607.25866 v1 pith:2Z6ZARBX submitted 2026-07-28 math.NA cs.DCcs.MScs.NAcs.PF

Simon Candelaresi , Benedict Geihe , Marco Artiano , Lars Christmann , Valentin Churavy , Andrés Rueda-Ramírez , Hendrik Ranocha , Gregor J. Gassner
show 1 more author
Michael Schlottke-Lakemper
This is my paper
classification math.NAcs.DCcs.MScs.NAcs.PF MSC 65Y0565M60
keywords Juliahigh-performancecomputingparallelscalabilityMPIdiscontinuousGalerkinTaylor-Greenvortexmagnetohydrodynamicssystemimage
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper aims to close a gap in the evidence for Julia as an HPC language: it asks whether Julia's just-in-time compilation and high-level syntax can survive traditional MPI-scale workloads. It answers with a scaling study of Trixi.jl, a Julia solver for conservation laws, on up to 61,440 CPU cores across several large clusters, compared against the established Fortran solver FLUXO. The paper shows that, given enough work per core, Trixi.jl retains a parallel efficiency of 0.83 at 61,440 ranks relative to a single node, and it identifies two practical obstacles—startup code loading and the global time-step reduction—plus a fix for the first. The result matters because it suggests scientist-developers can get developer productivity and Fortran-class scaling in one language.

What carries the argument

The load-bearing mechanism is Trixi.jl's MPI communication strategy: one non-blocking message exchange between neighbors per right-hand-side evaluation, overlapped with local computation, at the cost of recomputing some fluxes (whereas the Fortran code FLUXO uses a 'ping-pong' exchange that computes each flux once but needs two messages). The second piece of machinery is a custom system image built with Julia's package compiler, which precompiles all methods for a given simulation into a single large file; this removes both the startup compilation latency and the many-small-files I/O bottleneck that otherwise appears at thousands of MPI processes.

What would settle it

Re-run the Taylor-Green strong-scaling test with a cache-corrected baseline—e.g., using the weak-scaling efficiency curve or padding the per-rank working set so the single-node run is not cache-thrashing—and compare marginal efficiency per doubling of ranks. If the 61,440-core efficiency drops below about 0.5 against that baseline, the headline scaling result is an artifact of the cache-effect speedup rather than true HPC scaling.

Watch

Extended reading notes

Core claim

Julia's standard MPI library is sufficient for large parallel simulations. On the Taylor-Green vortex, Trixi.jl reaches 61,440 MPI ranks on 480 nodes with a parallel efficiency of 0.83 relative to a single-node run, and it outperforms FLUXO in raw throughput while showing slightly weaker scaling at extreme rank counts. A superlinear speedup, with efficiency peaking at 1.92 for 16,384 ranks, is attributed to cache effects; the marginal efficiency beyond that point is about 0.65 per doubling. The adaptive time-step calculation—a global MPI allreduce—becomes a bottleneck around 18,432 ranks, overtaking the right-hand-side computation. With a fixed time step, the code shows near-ideal scaling to

Load-bearing premise

The claim of 0.83 parallel efficiency at 61,440 cores rests on using the single-node run as the baseline, even though the reported superlinear speedup indicates the single-node baseline is partly a cache-thrashing artifact; a cache-corrected baseline would lower the headline efficiency.

Editorial extensions

If this is right

  • Julia can be a drop-in language for large MPI simulations in scientific domains, at least for explicit high-order DG methods.
  • Problem size per rank controls scaling: with 8x more elements per rank, Trixi.jl scales to 61,440 ranks with about 83% efficiency; with small problems it flattens near 16k-18k ranks.
  • The global time-step reduction is a measurable, localizable bottleneck—at around 18,432 ranks it costs more than the RHS itself—so replacing that single allreduce is a concrete optimization target.
  • Startup code loading at scale can be neutralized with custom system images, removing a barrier that had made Julia's first-run behavior incompatible with short HPC jobs.
  • Trixi.jl's throughput advantage over FLUXO at moderate rank counts comes with a communication strategy that duplicates some flux computations; at very high rank counts this trade-off starts to hurt.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • If the cache-effect explanation is right, then the reported 0.83 efficiency at 61,440 cores is optimistic: against a baseline that removes cache-thrash gains on small rank counts, the true efficiency would be lower and the headline would change. A simple test would be to measure with per-rank working sets that stay within cache at the baseline.
  • The system-image fix points to a wider lesson: Julia's JIT startup cost is not an obstacle once deployment artifacts are precompiled per machine—this could translate into HPC-friendly container images or job templates.
  • The same scaling methodology could be applied to other Julia solvers with different communication patterns; a natural extension is to study whether implicit time integrators or AMR-heavy workloads show the same scaling limits.
  • Since the bottleneck is one global allreduce per step, replacing it with hierarchical or asynchronous reduction is an immediately testable change that should push the scaling plateau well beyond 18k ranks if communication is the issue.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 4 minor

Summary. The paper reports a strong-scaling study of the Julia-based DGSEM code Trixi.jl against the Fortran code FLUXO on three CPU clusters (JURECA, JUWELS Booster, MOGON NHR), using the 3D Taylor-Green vortex and an MHD Alfvén wave test. The central claim, stated in the abstract and conclusions, is that Trixi.jl scales well to 61,440 CPU cores on 480 nodes using standard MPI.jl, with a parallel efficiency of 0.83 relative to a single-node baseline. The paper also documents Julia startup and code-loading bottlenecks and shows that a custom system image alleviates them, and it investigates an MPI all-reduce bottleneck in adaptive time stepping.

Significance. If the scaling result is robust, the paper provides a valuable data point for the HPC community: a high-level, JIT-compiled language can be a viable alternative to Fortran for massively parallel CPU simulations. The comparison to an established Fortran solver on multiple clusters is a strength, and the paper is commendably explicit about startup costs, compilation, and the role of problem size. The authors' transparency about the superlinear speedup and its effect on efficiency is also a positive feature. However, the headline efficiency metric is sensitive to the choice of baseline, and the manuscript's own numbers indicate that the 'scales well up to 61,440 cores' formulation overstates the robustness of the scaling data.

major comments (3)
  1. [Results (Fig. 1) and Summary and conclusions] The central claim of 'good scaling up to 61,440 cores' rests on a parallel efficiency of 0.83 measured against a single-node baseline. In the same results paragraph, the authors note that the efficiency peaks at 1.92 at 16,384 ranks because of cache effects and that the marginal efficiency beyond that point is only about 0.65 per doubling. Using the 16,384-rank point as a cache-corrected baseline, the efficiency at 61,440 ranks is approximately 0.43. The abstract and conclusions do not carry this qualification. Please report a cache-corrected efficiency or revise the scaling claim to something like 'good scaling up to about 16,384 ranks, with continued but substantially degraded throughput beyond.'
  2. [Experimental setup (software versions)] The manuscript states that JURECA runs use Julia v1.8.3 with bounds checking disabled, while JUWELS and MOGON use Julia v1.10.6 with bounds checking enabled. Because the headline 61,440-core result comes from JURECA only, this inconsistency does not invalidate that specific result, but it is a confound when comparing scaling behavior across clusters and when generalizing the claim that 'Julia scales well.' Please state whether the Julia version or bounds-checking setting materially affects the JURECA results, and explicitly restrict cross-cluster generalizations if they are not controlled.
  3. [Results (all scaling figures)] Each simulation was run five times and only the minimum is reported. The paper gives no measure of run-to-run variability. Given the superlinear speedup and the observed MPI time-step bottleneck at high rank counts, the minimum can be sensitive to outliers such as OS jitter or MPI startup effects. At minimum, report min-max ranges or a standard deviation for the key data points, and indicate whether the observed trends are stable across the five runs.
minor comments (4)
  1. [Adaptive Time Step Bottleneck] The sentence 'This shows the role that the problem size plays in the scalability. Keeping the problem size low increases the overhead...' appears in section 'Results' with a broken reference 'section .' and is repeated near the time-step analysis. The manuscript has several unresolved section cross-reference placeholders ('section .', 'As outlined in section') that must be fixed.
  2. [Figure 8 and text] The claim that 'at about the break-even point the scaling worsens' is based on visual coincidence. A quantitative statement (e.g., comparing the rate of increase of dt time to the rate of decrease of RHS time) would strengthen the causal argument.
  3. [Package loading / avoiding I/O bottleneck] The custom system image section is introduced as 'The code loading mechanism described in section poses challenges...' with another missing reference; also the text says 'The timing tests are preformed' (typo) and 'MPI paralellized' (typo) in the Introduction.
  4. [Figures 6 and 7] Figure 6's caption reads 'compressible MHD and Euler equations', but the text suggests the figure compares the Euler Taylor-Green problem with the MHD Alfvén wave for MOGON. Please clarify the exact configuration and ensure the caption matches the plotted data.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: scaling claim is an external benchmark against FLUXO; self-citations describe the software and do not carry the central result.

full rationale

The paper's central claim—that Trixi.jl scales to 61,440 cores on JURECA—is an empirical benchmark measured in this paper and compared against the external Fortran solver FLUXO. There is no fitted parameter that is later renamed as a prediction, and no quantity is defined in terms of the claim it is supposed to establish. Parallel efficiency is defined in the standard way as speedup relative to a single-node baseline, and the paper explicitly discloses the superlinear cache-effect artifact (efficiency 1.92 at 16,384 ranks and marginal efficiency ~0.65 per doubling beyond that), which weakens the headline interpretation but does not make the measurement circular. Self-citations to Trixi.jl, FLUXO, and related code-development papers are descriptive and are not the evidence for the scaling results; the scaling data are generated by the reported runs and compared to an independent code. The skeptical concern about the choice of baseline is a benchmarking-reasonableness issue, not a circularity issue: no reduction of a derived quantity to its own definition or to a self-citation chain is present. Therefore the paper receives a circularity score of 0.

Assumptions & free parameters 0 free parameters · 4 assumptions · 0 invented entities

The scaling claim rest on standard numerical-analysis assumptions (DGSEM stability, MPI correctness) and on methodological assumptions about benchmark comparability. No free parameters are fitted to data and no new entities are postulated.

assumptions (4)
  • domain assumption DGSEM discretization and low-storage Runge-Kutta time integration are representative of typical scientific HPC workloads.
    The paper generalizes from these two solvers and two test problems to a statement about Julia's HPC viability; the claim only holds for this class of workloads.
  • domain assumption The single-node baseline is a fair reference for parallel efficiency.
    Section Results, Fig. 1: the reported superlinear speedup (efficiency > 1) indicates the baseline is distorted by cache effects, undermining the 0.83 efficiency at 61,440 cores.
  • domain assumption Mixing Julia v1.8.3 with bounds checking disabled (JURECA) and v1.10.6 with bounds checking enabled (JUWELS/MOGON) does not affect the conclusions.
    Experimental setup states these different configurations, but no cross-check or sensitivity analysis is provided to show they are comparable.
  • domain assumption The minimum of five runs is a stable performance estimate.
    No variance or distribution is reported, so the break-even point in Fig. 8 and the efficiency numbers lack uncertainty information.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Massively parallel numerical simulations with Julia." pith.science (2026). https://pith.science/paper/2Z6ZARBX

@misc{pith2026260725866,
  author       = {Pith},
  title        = {Pith review of: Massively parallel numerical simulations with Julia},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/2Z6ZARBX}},
  note         = {Machine review of arXiv:2607.25866}
}
abstract

The Julia programming language aims to provide a modern approach to develop high-performance computing (HPC) applications. It tries to achieve this by combining a high-level, dynamic interface with just-in-time compilation to native machine code, thereby facilitating high developer productivity and native code performance at the same time. While this approach has already been shown to work well for serial applications, it is not clear if it readily translates to traditional, massively parallel HPC work loads. In this paper, we fill this gap by analyzing the parallel performance of the numerical computational fluid dynamics simulation code Trixi$.$jl, written in Julia, and compare it to the Fortran code FLUXO. We show some of the challenges of using Julia at scale and discuss possible solutions, specifically with respect to code loading and compilation at startup. Finally, we demonstrate the parallel scaling of our Julia code on up to 61440 CPU cores.

Figures

Figures reproduced from arXiv: 2607.25866 by the authors.

Figure 1
Figure 1. Degree of freedom updates per second for a flux differencing DG discretization of the 3D compressible Euler equations in FLUXO and Trixi.jl running on the JURECA cluster. Higher is better [PITH_FULL_IMAGE:figures/full_fig_p005_1.png] view at source ↗
Figure 4
Figure 4. Degree of freedom updates per second for a flux differencing DG discretization of the 3D compressible Euler equations in Trixi.jl running on the MOGON-NHR cluster. Higher is better. Next, Figures 4 and 5 show the results of the experiments on the MOGON NHR cluster. Here we run only the larger problem with 8× the number of degrees of freedom. Just as before, we get overall good scaling for up to 32 768 MPI ranks; bey… view at source ↗
Figure 5
Figure 5. Speedup for a flux differencing DG discretization of the 3D compressible Euler equations in Trixi.jl running on the MOGON-NHR cluster [PITH_FULL_IMAGE:figures/full_fig_p006_5.png] view at source ↗
Figures from the paper (2 more)
Figure 6
Figure 6. Figure 6: Degree of freedom updates per second for a flux differencing DG discretization of the 3D compressible MHD and Euler equations in Trixi.jl running on the MOGON-NHR cluster. Higher is better. Here we use the same resolution as for the larger Taylor￾Green vortex problem w…
Figure 9
Figure 9. Figure 9: Degree of freedom updates per second for a flux differencing DG discretization of the 3D compressible MHD equations in Trixi.jl running on the MOGON-NHR cluster for fixed time step size dt and different polynomial degrees p. Higher is better Fixed Time Step To test the…

Discussion (0). Sign in to comment.

Reference graph

Works this paper leans on

20 extracted references · 1 linked inside Pith

  1. [1]

    Beck A (2015) High order discontinuous galerkin methods for the simulation of multiscale problems

    2015

  2. [2]

    IEEE Transactions on Parallel and Distributed Systems doi:10.1109/TPDS.2018.2872064

    Besard T, Foket C and De Sutter B (2018) Effective extensible programming: Unleashing Julia on GPUs . IEEE Transactions on Parallel and Distributed Systems doi:10.1109/TPDS.2018.2872064

    arXiv 2018

  3. [3]

    SIAM Rev

    Bezanson J, Edelman A, Karpinski S and Shah VB (2017) Julia: A fresh approach to numerical computing . SIAM Rev. 59(1): 65--98. doi:10.1137/141000671

    doi:10.1137/141000671 2017

  4. [4]

    part i: Theory and numerical verification

    Bohm M, Winters AR, Gassner GJ, Derigs D, Hindenlang F and Saur J (2020) An entropy stable nodal discontinuous galerkin method for the resistive mhd equations. part i: Theory and numerical verification. J.\ Comput.\ Phys. 422: 108076

    2020

  5. [5]

    jl: Julia bindings for the message passing interface

    Byrne S, Wilcox LC and Churavy V (2021) Mpi. jl: Julia bindings for the message passing interface. In: Proceedings of the JuliaCon Conferences, volume 1. p. 68

    2021

  6. [6]

    SIAM J ournal on S cientific C omputing 33(3): 1103--1133

    C arsten B urstedde, L ucas C W ilcox and O mar G hattas (2011) p4est : S calable A lgorithms for P arallel A daptive M esh R efinement on F orests of O ctrees . SIAM J ournal on S cientific C omputing 33(3): 1103--1133. doi:10.1137/100791634

    doi:10.1137/100791634 2011

  7. [7]

    doi:10.48550/arXiv.2211.02740

    Churavy V, Godoy WF, Bauer C, Ranocha H, Schlottke-Lakemper M, R \"a ss L, Blaschke J, Giordano M, Schnetter E, Omlin S, Vetter JS and Edelman A (2022) Bridging HPC communities through the J ulia programming language. doi:10.48550/arXiv.2211.02740

    arXiv doi:10.48550/arxiv.2211.02740 2022

  8. [8]

    J.\ Comput.\ Phys

    Derigs D, Winters AR, Gassner GJ, Walch S and Bohm M (2018) Ideal glm-mhd: about the entropy consistent nine-wave magnetic field divergence diminishing ideal magnetohydrodynamics equations. J.\ Comput.\ Phys. 364: 420--467

    2018

Show all 20 references
  1. [9]

    In: 2022 IEEE International Conference on Cluster Computing (CLUSTER)

    Giordano M, Klöwer M and Churavy V (2022) Productivity meets performance: Julia on a64fx. In: 2022 IEEE International Conference on Cluster Computing (CLUSTER). pp. 549--555. doi:10.1109/CLUSTER51413.2022.00072

    2022

  2. [10]

    In: Proceedings of the SC '23 Workshops of The International Conference on High Performance Computing, Network, Storage, and Analysis, SC-W '23

    Godoy WF, Valero-Lara P, Anderson C, Lee KW, Gainaru A, Ferreira Da Silva R and Vetter JS (2023) Julia as a unifying end-to-end workflow language on the frontier exascale system. In: Proceedings of the SC '23 Workshops of The International Conference on High Performance Comput...

    2023

  3. [11]

    Hunold S and Steiner S (2020) Benchmarking julia’s communication performance: Is julia hpc ready or full hpc? In: 2020 IEEE/ACM Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems (PMBS). pp. 20--25. doi:10.1109/PMBS51919.2020.00008

    2020

  4. [12]

    Journal of large-scale research facilities 7(A183)

    J\" u lich Supercomputing Centre (2021) JUWELS Cluster and Booster: Exascale Pathfinder with Modular Supercomputing Architecture at Juelich Supercomputing Centre . Journal of large-scale research facilities 7(A183). doi:10.17815/jlsrf-7-183

    2021 doi

  5. [13]

    (2021) Flexi: A high order discontinuous galerkin framework for hyperbolic--parabolic conservation laws

    Krais N, Beck A, Bolemann T, Frank H, Flad D, Gassner G, Hindenlang F, Hoffmann M, Kuhn T, Sonntag M et al. (2021) Flexi: A high order discontinuous galerkin framework for hyperbolic--parabolic conservation laws. Comput.\ Math.\ Appl. 81: 186--219

    2021

  6. [14]

    In: Proceedings of the International Symposium on Code Generation and Optimization: Feedback-Directed and Runtime Optimization, CGO '04

    Lattner C and Adve V (2004) Llvm: A compilation framework for lifelong program analysis & transformation. In: Proceedings of the International Symposium on Code Generation and Optimization: Feedback-Directed and Runtime Optimization, CGO '04. USA: IEEE Computer Society, p. 75....

    2004

  7. [15]

    doi:10.1145/3625559

    Ranocha H, Schlottke-Lakemper M, Chan J, Rueda-Ram\' rez AM, Winters AR, Hindenlang F and Gassner GJ (2023) Efficient implementation of modern entropy stable and kinetic energy preserving discontinuous galerkin methods for conservation laws. doi:10.1145/3625559

    2023 doi

  8. [16]

    Proceedings of the JuliaCon Conferences 1(1): 77

    Ranocha H, Schlottke-Lakemper M, Winters AR, Faulhaber E, Chan J and Gassner G (2022) Adaptive numerical simulations with T rixi.jl: A case study of J ulia for scientific computing. Proceedings of the JuliaCon Conferences 1(1): 77. doi:10.21105/jcon.00077

    2022 doi

  9. [17]

    part ii: Subcell finite volume shock capturing

    Rueda-Ram \' rez AM, Hennemann S, Hindenlang FJ, Winters AR and Gassner GJ (2021) An entropy stable nodal discontinuous galerkin method for the resistive mhd equations. part ii: Subcell finite volume shock capturing. J.\ Comput.\ Phys. 444: 110580

    2021

  10. [18]

    https://github.com/trixi-framework/Trixi.jl

    Schlottke-Lakemper M, Gassner GJ, Ranocha H, Winters AR and Chan J (2021 a ) T rixi.jl: A daptive high-order numerical simulations of hyperbolic PDE s in J ulia. https://github.com/trixi-framework/Trixi.jl. doi:10.5281/zenodo.3996439

    2021 doi

  11. [19]

    J.\ Comput.\ Phys

    Schlottke-Lakemper M, Winters AR, Ranocha H and Gassner GJ (2021 b ) A purely hyperbolic discontinuous G alerkin approach for self-gravitating gas dynamics. J.\ Comput.\ Phys. 442: 110467. doi:10.1016/j.jcp.2021.110467

    2021

  12. [20]

    Strauss R, Bishnu S and Petersen M (2023) Julia for geophysical fluid dynamics: performance comparisons between cpu, gpu, and fortran-mpi doi:10.22541/essoar.167390522.28476012/v1

    2023

Pith tools

Reviewed August 1, 2026 · model on record in the stance chip above.