Keaton J. Burns

Polynomial spectral methods produce fast, accurate, and flexible solvers for a wide range of PDEs with one bounded dimension, where the incorporation of general boundary conditions is well understood. However, extending these methods to domains with multiple bounded dimensions is challenging because of difficulties imposing boundary conditions at shared edges and corners. Previous approaches have included various workarounds, such as the anisotropic inclusion of boundary data at shared edges or techniques restricted to specific boundary conditions. Here, we present a general system for imposing boundary conditions in tensor-product spectral methods for elliptic equations on hypercubes. This system accommodates general boundary conditions that commute on intersecting subsurfaces, including any combinations of Dirichlet, Neumann, and Robin conditions on all boundaries. Our approach builds on the generalized tau method, a framework that unifies many different types of spectral schemes by specifying the polynomial residuals that determine their exact algebraic solutions. This system encompasses traditional collocation, classical tau, and sparse Galerkin formulations. As an essential requirement, we add specific tau corrections to the boundary conditions, in addition to the bulk PDE, which produce a unique set of compatible boundary data at shared subsurfaces. This method allows us to systematically construct generic solvers that are fully equivariant to hyperoctahedral symmetries (rotations and reflections of the domain). We present the method explicitly for the Poisson equation in two and three dimensions and describe its extension to arbitrary elliptic equations in any dimension. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/kburns/corner_taus and in the supplementary materials ( corner_taus-main.zip [4.46KB]). [Formula: see text]

Astrophysical and geophysical fluids commonly generate organized magnetic fields, despite having enormous magnetic Reynolds numbers Rm and abundant small-scale turbulence. Flow-induced dynamo action produces these fields, with the “kinematic dynamo problem” devoted to determining the rate at which a flow exponentially amplifies weak magnetic fields. However, previous studies on high-Rm kinematic dynamos have generated flows via imposed volumetric forcing or oscillatory boundary conditions. In this article, we investigate a system with three important attributes: realistic flow conditions, fast dynamo action (operational for Rm → ∞ ), and a subharmonic spatiotemporal structure. We show that unsteady Taylor-vortex flow, a regime observed in laboratory experiments, gives rise to fast dynamos with timescales and length scales twice those of the flow at high Rm . By numerically integrating a Floquet system driven by periodic oscillations of Taylor vortices, we solve the kinematic dynamo problem up to Rm = 3.2 × 10 6 , calculating the dynamo's growth rate as a function of Rm and streamwise wave number. We find the onset of instability and compute finite-time Lyapunov exponents, which identify the regions of Lagrangian chaos required for fast dynamo action. To our knowledge, unsteady Taylor-vortex flow produces the most physically motivated fast dynamo to date.

We present an automated procedure for computing model gradients for partial differential equation (PDE) solvers built on sparse spectral methods, a broad class of numerical techniques widely used in the study of fluid dynamics, continuum mechanics, waves, and pattern formation across disciplines. Our approach applies reverse-mode automatic differentiation to symbolic graph representations of PDEs and efficiently constructs adjoint solvers that retain the speed and flexibility of modern Fourier and polynomial spectral methods. We demonstrate the advantages of this approach with a comprehensive implementation in the open-source Dedalus framework. This work uniquely provides a differentiable spectral solver that supports a broad class of equations, geometries, and boundary conditions, and runs efficiently in parallel. It enables users to compute gradients and perform PDE-based optimization for a wide range of time-dependent and nonlinear models with minimal additional code. We demonstrate this system’s capabilities using canonical problems from the literature, showing both strong performance and practical utility for a wide variety of optimization tasks. By integrating automatic adjoints into a flexible solver, our work enables researchers to perform sensitivity analyses in spectral simulations with ease and efficiency.

We consider the steady heat transfer between a collection of impermeable obstacles immersed in an incompressible two-dimensional (2-D) potential flow, when each obstacle has a prescribed boundary temperature distribution. Inside the fluid, the temperature satisfies a variable-coefficient elliptic partial differential equation (PDE), the solution of which usually requires expensive techniques. To solve this problem efficiently, we construct multiply connected conformal maps under which both the domain and governing equation are greatly simplified. In particular, each obstacle is mapped to a horizontal slit and the governing equation becomes a constant-coefficient elliptic PDE. We then develop a boundary integral approach in the mapped domain to solve for the temperature field when arbitrary Dirichlet temperature data are specified on the obstacles. The inverse conformal map is then used to compute the temperature field in the physical domain. We construct our multiply connected conformal maps by exploiting the flexible and highly accurate AAA-LS algorithm. In multiply connected domains and domains with non-constant boundary temperature data, we note similarities and key differences in the temperature fields and Nusselt number scalings as compared with the isothermal simply connected problem analysed by Choi et al. ( J. Fluid Mech. , vol. 536, 2005, pp. 155–184). In particular, we derive new asymptotic expressions for the Nusselt number in the case of arbitrary non-constant temperature data in singly connected domains at low Péclet number, and verify these scalings numerically. While our language focuses on the problem of conjugate heat transfer (the transfer of heat between objects in a flow), our methods and findings are equally applicable to the advection–diffusion of any passive scalar in a potential flow.

We describe the user interface, governing equations, and numerical methods underpinning the community ocean modeling software called "Oceananigans". Oceananigans development has been lead by the Climate Modeling Alliance to build a trainable climate model with quantifiable uncertainty. Oceananigans is written in the Julia programming language, which, like similar recent efforts based on modern programming languages, distinguishes it from usual software based on Fortran. Oceananigans can efficiently simulate all scales of ocean motion, ranging from millimeter-scale turbulence in a small box to planetary-scale ocean circulation. Oceananigans design combines (i) a basic structured finite volume algorithm (ii) optimized for high-resolution simulations on GPUs which is (iii) exposed behind a high-level, programmable user interface. This design negotiates a dual mandate for highest-possible performance (to support state-of-the-art applications) and enhanced accessibility (to facilitate adoption and development). The dual mandate aims ultimately to accelerate the progress of Earth system science. Achieving this aim, however, requires a substantial and sustained increase in the collective effort of Oceananigans development.

We discuss the use of systematic ‘a posteriori’ calibration in the development of complicated (but theory-based) parameterizations. With ‘a posteriori’ calibration, model error is assessed using the results of forward simulations, thereby incorporating numerical error, numerical stability, model-specific implementation details, and alleviating the need for explicit data for all parameterized model components. We show how calibration illuminates the parameterization development trade-off between reductions in model bias, producing better predictions, and increased parametric complexity, the latter which can decrease a model’s ability to extrapolate, increase both the data requirements and computational expense of the calibration. We illustrate the importance of a posteriori calibration by describing the iterative development of CATKE, a new parameterization we develop within CliMA for the fluxes associated with small- or "micro-scale" ocean turbulent mixing on scales between 1 and 100 meters. For calibration we use Ensemble Kalman Inversion to minimize the error between a set of large eddy simulations (="the truth") and predictions of the parameterization and this way find optimal values for the free parameters. Without systematic calibration we cannot make informed choices about parameterization development because we cannot distinguish between structural error and error due to non-optimal parameter values.

The magnetic dynamo cycle of the Sun features a distinct pattern: a propagating region of sunspot emergence appears around 30° latitude and vanishes near the equator every 11 years (ref. ¹). Moreover, longitudinal flows called torsional oscillations closely shadow sunspot migration, undoubtedly sharing a common cause². Contrary to theories suggesting deep origins of these phenomena, helioseismology pinpoints low-latitude torsional oscillations to the outer 5–10% of the Sun, the near-surface shear layer^3,4. Within this zone, inwardly increasing differential rotation coupled with a poloidal magnetic field strongly implicates the magneto-rotational instability^5,6, prominent in accretion-disk theory and observed in laboratory experiments⁷. Together, these two facts prompt the general question: whether the solar dynamo is possibly a near-surface instability. Here we report strong affirmative evidence in stark contrast to traditional models⁸ focusing on the deeper tachocline. Simple analytic estimates show that the near-surface magneto-rotational instability better explains the spatiotemporal scales of the torsional oscillations and inferred subsurface magnetic field amplitudes⁹. State-of-the-art numerical simulations corroborate these estimates and reproduce hemispherical magnetic current helicity laws¹⁰. The dynamo resulting from a well-understood near-surface phenomenon improves prospects for accurate predictions of full magnetic cycles and space weather, affecting the electromagnetic infrastructure of Earth.