ABSTRACT. Weather forecasting is critical for a range of human activities including transportation, agriculture, industry, as well as the safety of the general public. Over the last four years, machine learning models have shown that they have the potential to transform the complex weather prediction pipeline, but current approaches still rely on numerical weather prediction (NWP) systems, limiting forecast speed and accuracy. In this talk, some of the background on these developments will be given. A new generation of machine learning model will then be introduced which can replace the entire operational NWP pipeline. Aardvark Weather, an end-to-end data-driven weather prediction system, ingests raw observations and outputs global gridded forecasts and local station forecasts. Further, it can be optimised end-to-end to maximise performance over quantities of interest. It will be shown that the system outperforms an operational NWP baseline for multiple variables and lead times for gridded and station forecasts. Finally, the talk will end by discussing how these ideas might develop over the next few years, including their application to multiple parts of the Earth system on multiple time-scales, and their potential impact on climate modelling.

ABSTRACT. In-context operator networks (ICON) are a class of non-intrusive operator learning methods based on the novel architectures of foundation models. Trained on a diverse set of datasets of initial and boundary conditions paired with corresponding solutions to ordinary and partial differential equations (ODEs and PDEs), ICON learns to map example condition-solution pairs of a given differential equation to an approximation of its solution operator. In this talk, we present a probabilistic framework that reveals ICON as implicitly performing Bayesian inference, where it computes the mean of the posterior predictive distribution over solution operators conditioned on the provided context (example condition-solution pairs). By modeling the dependence of example pairs through random differential equations, we formalize how, given example condition-solution pairs, ICON is a point estimate of the posterior predictive distribution over operators. This probabilistic perspective provides a basis for extending ICON to generative settings, where one can sample from the posterior predictive distribution over solution operators. As a result, ICON is no longer limited to point prediction, as it can capture the underlying uncertainty in the solution operator. This enables principled uncertainty quantification in operator learning by using generative modeling to produce confidence intervals for predictive solutions.

ABSTRACT. In this talk, we present a new class of generative algorithms capable of efficiently learning arbitrary target distributions from possibly scarce, high-dimensional data and subsequently generating new samples. These particle-based generative algorithms are constructed as gradient flows of Lipschitz-regularised Kullback–Leibler or other f-divergences. In this framework, data from a source distribution can be stably transported as particles towards the vicinity of the target distribution. We demonstrate through numerical experiments that the proposed approach achieves high accuracy, even in small-sample regimes.

ABSTRACT. For multivariate normal distributions, a complete description of independence and conditional independence is encoded in the zeros of the covariance matrix and its inverse, the precision matrix. Thus, estimating independence reduces to matrix estimation, which is relatively straightforward to do from data. But when the distribution is non-Gaussian, the correspondence between independence and these matrix zeros no longer holds, and independence is, in general, much more difficult to estimate. In this talk, we'll give some recent results that bring together diagonal transport, analytic covariance formulas for non-Gaussian distributions, moment-generating functions, and a so-called generalized precision to yield estimates of independence for non-Gaussian distributions.

ABSTRACT. We study adversarial learning when the target distribution factorizes according to a known Bayesian network. Our focus is on interpolative divergences, such as $(f,\Gamma)$-divergences, which connect classical $f$-divergences and integral probability metrics through variational objectives with constrained discriminator classes. We establish a new \emph{infimal subadditivity} principle showing that, under suitable conditions, a graph-constrained global objective is bounded by an average of family-level objectives on the local neighborhoods of the graph, with equality in an additive regime. This provides a variational justification for replacing a graph-agnostic GAN with a monolithic discriminator by a graph-informed GAN with localized family-level discriminators. We also establish parallel results for integral probability metrics and proximal optimal transport divergences, and present numerical experiments illustrating improved stability and structural recovery in graph-informed training.

ABSTRACT. Autodiff is a core technique in Machine Learning. However the mathematics of autodiff is much less developed for non smooth activation functions such as ReLU and maxpool. It is a damaging for applications in SciML since classical Finite Element spaces start with piecewise affine basis functions which can be implemented with ReLU.

I will show recent progress in this direction. It allows the discussion of the minimal regularity of activation functions for DeepRitz and PINNs methods and explains a seemingly paradoxical numerical behavior. 2D numerical tests confirm the theory.

ABSTRACT. We present a preliminary study combining two approaches in the context of PDE solving: the classical finite element method (FEM) and more recent techniques based on neural networks. Indeed, in recent years, Physics-Informed Neural Networks (PINNs) have become particularly interesting for quickly solving such problems, especially in large dimensions. However, their lack of precision is a major drawback in this context, hence the interest in combining them with FEM, for which error estimators are already known. This combination will make it possible to correct and certify the prediction of neural networks in order to obtain a fast and accurate solution. The complete pipeline proposed here then consists of modifying the classical approximation spaces in FEM by taking the information of a prior, chosen here as the prediction of a PINNs. Current results show that pre-processing the problem using neural networks can achieve fixed error targets with coarser meshes than in standard finite element methods, thus saving time in computing the solution. Error estimates have been proven showing that enriched spaces outperform classical ones by a factor that depends only on the quality of the prior.

ABSTRACT. A significant limitation of applying neural networks to scientific computing is their reliance on large datasets. This challenge is particularly acute when data are expensive to obtain, for example when they come from experiments. A way to mitigate this is the use of structure-preserving methods that incorporate certain properties of the underlying physical system into the neural network. While the potential of structure-preserving properties to reduce data requirements has been acknowledged in the literature, a systematic study of diverse cases is currently lacking. In this presentation, we provide preliminary results of such a study, evaluating the performance regarding data efficiency of various structure-preserving architectures including Hamiltonian neural networks and volume-preserving neural networks.

ABSTRACT. Nonlinear shell models couple membrane, bending, and, for Reissner-Mindlin type theories, shear effects. Their discretization is challenging because standard finite element methods typically suffer from membrane and shear locking. At the same time, the bending energy involves curvature quantities that are only well-defined on differentiable surfaces. In this talk, we present a structure-preserving finite element framework for nonlinear shells that combines ideas from mixed finite elements and discrete differential geometry to obtain locking-free discretization methods.

For the bending energy, we introduce a generalized Weingarten tensor (shape operator) that yields a weak notion of curvature on discretized surfaces, including affine triangulations. The construction is based on a mixed formulation with the bending moment tensor as an additional unknown, thereby extending the Hellan-Herrmann-Johnson method from linear Kirchhoff-Love plates to nonlinear shells. This avoids the need for C1-conforming elements and provides a consistent discretization of shell bending on nonsmooth geometries.

To address membrane locking, we show that the membrane strains should be projected into the Regge finite element space, whose tangential-tangential continuity makes Regge elements the natural discrete setting for strain and metric quantities. For Reissner-Mindlin shells, shear locking is avoided by a hierarchical approach using tangential continuous Nedelec elements, analogous to the TDNNS method. Several examples implemented in the open-source finite element software NGSolve (www.ngsolve.org) demonstrate the excellent performance of the proposed shell elements.

ABSTRACT. We present a novel finite element discretisation of the Reissner-Mindlin plate problem based on a Hellinger-Reissner variational principle with square-integrable symmetric bending moments. Conforming Hu–Zhang elements are employed to approximate these bending moments, yielding highly accurate fields. To ensure optimal convergence across all variables in the Kirchhoff-Love limit, the formulation is extended with Raviart–Thomas type elements for the shear stresses. Existence and uniqueness are established in the continuous setting, with well-posedness inherited discretely via exact complexes. Implementation via polytopal templates allows an efficient construction of the finite element spaces, and in particular the Hu-Zhang element, on high-order mesh geometries.

Sky, A., Neunteufel, M., Hale, J. S., & Zilian, A. (2025). Formulae and transformations for simplicial tensorial finite elements via polytopal templates. Computers & Mathematics with Applications, 195, 322–348. https://doi.org/10.1016/j.camwa.2025.07.028

Sky, A., Neunteufel, M., Hale, J. S., & Zilian, A. (2023). A Reissner–Mindlin plate formulation using symmetric Hu-Zhang elements via polytopal transformations. Computer Methods in Applied Mechanics and Engineering, 416, 116291. https://doi.org/10.1016/j.cma.2023.116291

ABSTRACT. The elastic Cosserat micropolar continuum theory introduces independent rotational degrees of freedom~$\boldsymbol{\theta}$ alongside the displacement field $\boldsymbol{\varphi}$. These micro-rotations are described by a rotation tensor field $\bar{\boldsymbol{R}}(\boldsymbol{\theta})\in\mathrm{SO}(3)$. Their interaction with the macro-rotation field, described by the rotational part of the deformation gradient $\boldsymbol{F}$, is governed by the Cosserat couple modulus $\mu_{\mathrm{c}}$ via the Cosserat strain $\bar{\boldsymbol{R}}^{\mathsf{T}}\boldsymbol{F}-\boldsymbol{I}$. For large values of $\mu_{\mathrm{c}}$, the two rotations coincide.

In this work, we examine a structural compatibility issue that arises in finite element discretizations of finite-strain Cosserat micropolar continua for large $\mu_{\mathrm{c}}$. When using a formulation purely based on Lagrange elements for both the displacement and micro-rotation fields, the structural mismatch between the classical discretizations of $\boldsymbol{F}_h=\mathrm{D}\boldsymbol{\varphi}_h$ and $\bar{\boldsymbol{R}}_h$ can lead to locking phenomena.

We introduce a new method, dubbed \textbf{\textit{Geometric Structure-Preserving Interpolation ($\boldsymbol{\varGamma}$-SPIN)}} that relaxes the interaction between the rotation tensor and the polar part of the deformation gradient. This relaxation is achieved in two steps. First, the regularity of the Cosserat rotation tensor is reduced by interpolating it into the Nédélec space. Second, the resulting field is projected back onto the Lie-group of rotations. Together, these steps define a lower-regularity projection-based interpolation. We outline the formulation and present numerical benchmarks, establishing consistency and stability of the method and showing a significant reduction in locking across several test cases.

ABSTRACT. We present a consistent reduction of the relaxed micromorphic model to a two-dimensional planar formulation that preserves its ability to capture discontinuous dilatation fields. This reduction naturally reveals the need for $H^{\mathrm{dev}}(\mathrm{Curl})$-conforming finite element spaces that are intrinsic to the structure of the model.

To this end, we introduce two novel $H^{\mathrm{dev}}(\mathrm{Curl})$-conforming finite element constructions: one defined on regular triangulations and another given by a macro-element approach based on Clough–Tocher splits. We further derive compatible primal and mixed variational formulations of the planar relaxed micromorphic model.

The effectiveness of the proposed elements is demonstrated through two numerical examples, confirming their suitability for accurately resolving discontinuous dilatation fields and supporting the intrinsic nature of the discretization.

ABSTRACT. High-order exponential integrators require computing linear combinations of exponential-like φ-functions of a large matrix acting on a vector. Krylov projection methods remain a highly efficient and general choice for evaluating these matrix-function-vector products, particularly when the matrix is too large to be explicitly stored or when obtaining spectral information is prohibitively expensive. However, the Krylov approximation relies on the Gram-Schmidt (GS) orthogonalization procedure to produce the orthonormal basis. In massively parallel environments, standard GS orthogonalization creates a significant bottleneck by requiring multiple global synchronizations for inner products and vector normalization. Reducing the number of global synchronizations is of paramount importance for algorithmic efficiency. In this talk, we improve the strong scaling properties and parallel efficiency of exponential integrators by addressing this linear algebra bottleneck via low-synchronization GS methods. The resulting algorithms maintain accuracy comparable to modified Gram-Schmidt but are better suited for distributed architectures, requiring only one global communication per orthogonalization step. We present numerical experiments to demonstrate that this reduction in global communication leads to better parallel scalability and a reduced time-to-solution for exponential integrators.

ABSTRACT. In this talk, exponential Runge-Kutta methods are constructed to solve semilinear parabolic problems with time-dependent and arbitrary state-dependent delay. The convergence properties and superconvergence results are analyzed. Numerical experiments are presented to illustrate our theoretical results.

ABSTRACT. In this talk, we present three new classes of exponential time integrators for stiff and highly oscillatory systems. For stiff evolution equations, including parabolic PDEs such as diffusion–reaction systems, we introduce two-derivative exponential integrators of Runge–Kutta and Rosenbrock types. For highly oscillatory problems, including hyperbolic PDEs such as wave equations, we introduce exponential Nyström integrators. Under reasonable regularity assumptions, we establish convergence results of up to fifth order for these methods in Banach and Hilbert spaces, with error bounds that remain independent of the stiffness or the high-frequency components of the problem. Numerical experiments illustrate the accuracy and efficiency of the proposed methods in comparison with existing schemes from the literature.

This is joint work with Hoang Nguyen and Huy Pham and was supported by NSF grants DMS-2309821 and DMS-2531805.

ABSTRACT. In this talk we will discuss several ways in which the ideas of exponential integration can be used to construct efficient schemes for stiff systems of differential equations. We will present a new framework to develop and to analyze new class of schemes we call stiffness resilient methods. Previously proposed exponential integrators are typically derived using either classical or stiff order conditions. These order conditions are complex and difficult to solve to construct high order schemes. Classically derived methods can also suffer from the order reduction phenomenon. The new $\varphi$-order conditions we propose allow greatly simplify construction of exponential methods with favorable properties. The structure of the error of these methods is designed to prevent order reduction for many important stiff problems. At the same time stiffness resilient schemes are easy to derive using our proposed approach. In addition, we will discuss applying exponential integration to problems in fields such as plasma physics and weather prediction. We will discuss how special exponential-type methods can be constructed to take advantage of the structure of the problem to further improve efficiency of the time integration.

ABSTRACT. Music audio signals are rich and complex, making the extraction of meaningful content (notes, instruments, structure, …) from complex audio mixtures a significant challenge in signal processing. The field of Music Information Retrieval (MIR) is dedicated to solving these tasks. However, the current prevalence of supervised deep learning in MIR has led to a heavy reliance on extensive, hard-to-obtain human annotations. Unsupervised or weakly-supervised methods are essential to bypass this limitation. In this context, low-rank factorization methods are well-positioned to excel: they provide powerful, elegant tools to blindly untangle complex audio mixtures and uncover interpretable musical information directly from the raw signal. This presentation explores the application of low-rank matrix and tensor factorization techniques (particularly Nonnegative Matrix Factorization (NMF) and Nonnegative Tucker Decomposition (NTD)) to audio signal analysis. The talk centers on two key contributions: the use of NTD for Music Structure Analysis (MSA) and the introduction of the nmf_audio_benchmark toolbox. First, we demonstrate the power of tensor-based factorization for music analysis. By structuring audio into time-frequency-bar tensors, low-rank NTD models naturally capture structural repetitions and isolate distinct musical motifs. We will highlight how these compressed, part-based representations can be leveraged for MSA, yielding an interpretable, blind pattern-extraction method. We also hypothesize that this method could be applied to musicological studies or to the development of controllable generative music models. Second, we introduce nmf_audio_benchmark, a comprehensive toolbox designed to bridge the gap between MIR and machine learning communities. This framework enables researchers to seamlessly evaluate their low-rank factorization models on music analysis tasks by integrating core factorization algorithms with standardized music datasets and domain-specific post-processing pipelines. Ultimately, this presentation demonstrates the strong synergy between complex audio signals and low-rank matrix and tensor methods. By providing both a concrete use-case in musical pattern extraction and a standardized benchmarking toolbox, this work aims to bridge the gap between MIR and tensor computation communities. We show that low-rank factorization remains a highly relevant, computationally efficient, and intrinsically interpretable approach for tackling the complex realities of modern audio analysis.

ABSTRACT. We introduce a new low-dimensional model of high-dimensional numerical simulation data based on low-rank tensor decompositions. Our new model aims to minimize differences between the model data and simulation data as well as functions of the model data and functions of the simulation data. This novel approach to dimensionality reduction of simulation data provides a means of directly incorporating quantities of interests and invariants associated with conservation principles associated with the simulation data into the low-dimensional model. Computational results of applying this approach to two standard low-rank tensor decompositions of data arising from simulation of combustion and plasma physics are presented, illustrating improved modeling of quantities of interestusing our approach versus the standard decompositions.

ABSTRACT. The canonical polyadic decomposition (CPD) can be used as a low-rank representation of tensors appearing in various applications. However, obtaining a CPD can lead to a difficult optimization problem when, e.g., the rank is higher than some of the tensor dimensions, when the rank-1 terms are close to collinear, or when there are significant differences in magnitude of the terms. During the optimization, a degenerate subtensor can then be formed, i.e., an ill-conditioned tensor that tries to approximate a tensor of higher rank. This can lead to regions in the optimization landscape of very slow convergence, informally known as `swamps'. Efficient optimization frameworks are proposed to detect the degenerate terms and speed up the convergence of state-of-the-art methods for difficult-to-decompose tensors, while remaining as efficient on easier cases.

ABSTRACT. Tensor networks provide for machine learning flexible compressed representations that can address the curse of dimensionality and be used to discover patterns providing a versatile framework for both supervised and unsupervised learning. In the talk, it will first be highlighted how tensor networks can leverage multivariate polynomial regression for supervised learning. The talk will subsequently discuss how tensor networks can be used in unsupervised learning for the characterization of high-dimensional joint distributions useful for density estimation, outlier detection, and generative modeling. Finally, it will be demonstrated how uniqueness properties of certain tensor network representations can be used to discover patterns in functional neuroimaging datasets.

ABSTRACT. We study the Generalized Convolution Quadrature (gCQ) based on Runge--Kutta methods to approximate the solution of an important class of convolution equations. The gCQ generalizes Lubich's original Convolution Quadrature to variable steps. High order versions of the gCQ have been developed in the last decade in a rather general setting, which includes applications to hyperbolic problems, where the convolution kernels are typically non smooth. However, in the case of convolutions with smoother, sectorial kernels, recent developments show that the available stability and convergence theory was suboptimal, both in terms of convergence order and regularity requirements of the data. Optimal results for the gCQ of the first order are now available, for a special important class of sectorial problems. We address here the generalization of this theory to high order and prove the same order of convergence as for the original Runge--Kutta based CQ with fixed steps, under the same regularity hypotheses about the data, and for arbitrary time meshes. Moreover, for data with known singularities of algebraic type, we show how to choose optimally graded time meshes in order to achieve maximal order of convergence, overcoming the well-known order reduction of the original CQ in these situations. We also show that a fast and memory reduced implementation of the gCQ is possible for this class of problems and illustrate our theoretical results with several numerical experiments.

ABSTRACT. In this contribution, we are interested in the analysis and construction of low-cost, explicit, two-step peer methods [1] for the numerical solution of ordinary differential equations, where the second derivative of the solution is incorporated into the formula of the methods [2].

The main properties of this class of methods—such as consistency, convergence, and absolute stability—along with their order conditions, are analyzed. Several methods using only two evaluations of the first and second derivative per step are constructed; examples of such methods with three and four stages, up to order seven, are provided.

Finally, several numerical experiments demonstrate the efficiency of the constructed methods.

References:

[1] M. Calvo, J.I. Montijano, L. Rández, A. Saenz de la Torre Explicit two-step peer methods with reused stages Applied Numerical Mathematics 195, 2024, 75-88

[2] M. Sharifi, A. Abdi, M. Braś, G. Hojjati On implicit second derivative two-step peer methods with RK stability for ODEs Applied Numerical Mathematics 220, 2026, 329-345

ABSTRACT. In this contribution, we focus on the design and analysis of structure-preserving numerical methods for deterministic dynamical systems within the framework of Nonstandard Finite Difference methods (NSFD). Particular attention is devoted to the class of Patankar schemes and their extensions, including Modified Patankar methods, which have proven especially effective for production–destruction systems arising in applications. These schemes are characterized by their ability to enforce unconditional positivity and, in many cases, exact conservation, independently of the time step size.

We present recent advances in the development of Modified Patankar NSFD methods, with an emphasis on improving their consistency and accuracy while preserving their inherent structural properties. Numerical examples illustrate the behavior of the proposed approaches and highlight the good performance of the new methods, also through comparisons with well-known positivity-preserving schemes.

ABSTRACT. We introduce a numerical method for solving second-order stochastic differential equations describing a class of nonlinear oscillators with non-constant frequency, perturbed by white noise. The proposed scheme takes advantages of the Magnus approach to construct an integrator for this stochastic oscillator. Its convergence properties are rigorously analyzed and selected numerical experiments on relevant stochastic oscillators are carried out confirming the effectiveness of the proposed method.

ABSTRACT. Recent years have seen rapid progress in machine learning methods for partial differential equations, including physics-informed neural networks and neural operators, but challenges in reliability, accuracy, and out-of-distribution generalisation remain. The finite element method continues to serve as a robust workhorse of scientific computing, with accuracy governed by the representation capacity of the underlying discrete space and therefore in large part by the mesh. Adaptive meshing improves the resulting accuracy-cost trade-off by concentrating degrees of freedom where they are most effective. However, classical adaptive strategies incur substantial additional overhead through a posteriori error estimation, monitor-based mesh movement, adjoint-based optimisation, or the solution of auxiliary meshing PDEs such as Monge–Ampère formulations. Recent graph-based learning methods suggest an alternative paradigm in which the adaptation operator itself is learned within the numerical pipeline. G-Adaptivity uses graph neural networks to learn an $r$-adaptive mesh relocation policy by directly minimising finite element solution error, rather than approximating a classical meshing operator. HypeR Adaptivity extends this perspective to joint $hr$-adaptivity through a hypergraph-based multi-agent reinforcement learning formulation in which relocation and refinement are optimised simultaneously, overcoming the limitations of greedy heuristics. Both approaches deliver orders-of-magnitude speed-ups over classical adaptivity while improving the attainable accuracy at fixed mesh complexity.

ABSTRACT. Recently, there has been a lot of interest in machine learning based surrogate models for atmospheric fluid dynamics. Successful approaches such as GraphCast [Lam et al, Science, 382(6677), 2023] use an encoder-processor-decoder architecture where the processor describes the propagation of information on a Graph Neural Network (GNN) via message passing. Since the dynamics are approximated on a low-dimensional space, these models are much faster than traditional methods that solve the underlying PDE numerically on a fine mesh.

We explore a variation of this architecture where the processor is replaced by the solution of a time-dependent differential equation in a low-dimensional latent space. Compared to GNNs, this allows the application of standard techniques from numerical analysis, thus potentially improving stability and interpretability for example by exactly enforcing conservation laws.

To implement this, we use the recent Firedrake/PyTorch interface [Bouziani & Ham, 2023] to train and solve time-dependent PDE surrogate models on the sphere. The encoder in our model combines the interpolation of the initial condition to the latent space on a vertex-only-mesh with a learnable embedding; the decoder has a similar structure based on the adjoint of the interpolation. Our model is trained on numerical solutions of PDE which have been calculated a-priori using Firedrake. The aim of our work is to combine the reliability of Finite Elements with the efficiency of Neural Networks surrogates to produce a competitive model which has the potential to be applied to time-critical applications in weather forecasting.

ABSTRACT. Convection-diffusion equations arise in many physical applications, particularly in models with complex fluid flows. The numerical solution of discretized convection-diffusion equations is particularly difficult, even with the best-known preconditioners from the elliptic case, such as multigrid methods. In particular, as the problem becomes convection-dominated, standard geometric multigrid methods lose effectiveness, with poor relaxation leading to stagnating convergence, or even divergence with stationary iterations. We present a machine learning-enhanced multigrid method to find relaxation schemes that restore solver efficiency in these demanding regimes.

ABSTRACT. Deep learning offers a powerful approach to quantum many-body problems via neural network wavefunctions, but their optimization remains a key bottleneck, thereby limiting their potential. Existing techniques, such as natural gradient descent and its analog in quantum physics (stochastic reconfiguration), have slow convergence rates that critically depend on the Hamiltonian's spectral gap. This limits their applicability to open problems fraught with competing orders and nearly degenerate ground states, where predictive power hinges on high numerical accuracy, including frustrated magnets, strongly correlated electron materials, and multireference quantum chemistry. Here, we introduce Projected Inverse Iteration (PII), reframing the ground-state search as an eigenvalue problem. PII achieves rapid, gap-insensitive convergence while preserving the favorable polynomial computational scaling of stochastic reconfiguration. Demonstrated on challenging two-dimensional spin systems, including the highly frustrated $J_1$-$J_2$ model, PII outperforms standard natural gradient techniques. This provides a robust, scalable algorithmic strategy for discovering novel complex quantum states using deep learning, particularly in regimes where vanishing spectral gaps hinder optimization.

ABSTRACT. We consider numerical methods for the dynamical core in numerical weather prediction (NWP). Specifically, we are interested in schemes for the Euler equations in 3D at low Mach numbers with a gravity source term. As the trend in NWP is to increase the resolution, high order Discontinuous Galerkin (DG) methods are competitive, proving well suited for High Performance Computing due to their high arithmetic intensity. Due to stiffness, implicit discretizations are needed. The overall design goal is thus a suitable highly parallel stable and computational efficient DG implementation.

To deal with the gravitational source term, a well-balanced method is needed, suitable for large domains. Additionally, the scheme requires a numerical flux suitable for low Mach numbers, and some form of nonlinear stability such as entropy stability, suitable for large simulation times.

Environment and Climate Change Canada (ECCC) is currently redesigning its dynamical core along these lines within the WxFactory code. We present the implementation of a well-balanced DG method suitable for low Mach numbers. Specifically, we use a well-balanced Direct Flux Reconstruction implementation on GL nodes on quadrilateral meshes and a variant of the AUSM+up flux. For time integration, an exponential integrator is used. To achieve entropy stability, we incorporate an artificial diffusion term and a relaxation technique to stabilize the time integration. We provide numerical results showing the performance on atmospheric test cases in 2D and 3D.

ABSTRACT. FourCastNet 3 advances global weather modeling by implementing a scalable, geometric machine learning (ML) approach to probabilistic ensemble forecasting. The approach is designed to respect spherical geometry and to accurately model the spatially correlated probabilistic nature of the problem, resulting in stable spectra and realistic dynamics across multiple scales. FourCastNet 3 delivers forecasting accuracy that surpasses leading conventional ensemble models and rivals the best diffusion-based methods, while producing forecasts 8 to 60 times faster than these approaches. In contrast to other ML approaches, FourCastNet 3 demonstrates excellent probabilistic calibration and retains realistic spectra, even at extended lead times of up to 60 days. All of these advances are realized using a purely convolutional neural network architecture tailored for spherical geometry. Scalable and efficient large-scale training on 1024 GPUs and more is enabled by a novel training paradigm for combined model- and data-parallelism, inspired by domain decomposition methods in classical numerical models. Additionally, FourCastNet 3 enables rapid inference on a single GPU, producing a 60-day global forecast at 0.25°, 6-hourly resolution in under 4 minutes. Its computational efficiency, medium-range probabilistic skill, spectral fidelity, and rollout stability at subseasonal timescales make it a strong candidate for improving meteorological forecasting and early warning systems through large ensemble predictions.

ABSTRACT. Numerical algorithms designed for global weather prediction at horizontal resolutions of 10's to 100's of km are now being pushed towards km and sub-km resolution. Are they still fit for purpose in the very different flow regimes now being captured? I will discuss the performance of a semi-implicit semi-Lagrangian model in Implicit Large Eddy Simulations (ILES) of boundary layer flows. Away from the surface, the performance of the ILES approach is comparable to that of traditional LES. However, for shear driven flows, the ILES numerical diffusion does not adequately represent subgrid momentum fluxes near the surface, resulting in biases that are sensitive to numerical choices such as limiters. Sharp changes of gradient in the wind profile near the surface also result in significant momentum conservation errors. For convective flow driven by surface heating, on the other hand, ILES performs much better. New diagnostics are presented to diagnose energy dissipation by the numerics. The dissipation of kinetic energy has a realistic magnitude and is spatially correlated with the rate of strain. The model's transport scheme conserves entropy, neglecting sources due to mixing; this results in a spurious loss of energy by the model that is comparable in magnitude to the kinetic energy dissipation.

ABSTRACT. Currently, a new dynamical core for the established weather and climate forecast model ICON, based on the Discontinuous Galerkin (DG) method, is under development at the Deutscher Wetterdienst (DWD). The DG method combines conservation of the prognostic variable via the finite volume approach with higher order accuracy via the finite element approach. Additionally, it allows the use of explicit time integration schemes and is designed for massively-parallel computers due to its very compact discretization stencils.

One distinct aspect in the modeling of the atmosphere are large aspect ratios of the grid cells in particular near the lower boundary, which require implicit methods. Here, the horizontally explicit-vertically implicit (HEVI) treatment is used in combination with the collocation method (DG-SEM) and IMEX-Runge-Kutta schemes. Implications of the use of HEVI together with the boundary conditions are discussed. Two different versions of the Euler equations with regard to the thermodynamic variable are compared. They are written in covariant form to be applicable on arbitrary manifolds.

Since the combination of the Euler solver with several parameterizations is crucial, the combination with a one-equation TKE turbulence model is presented, in which the diffusion term is treated with HEVI, too.

ABSTRACT. Heating generated by high-intensity focused ultrasound waves is central to many emerging medical applications, including non-invasive cancer therapy and targeted drug delivery. In this talk, we present the analysis of a conforming finite element approximations of the underlying nonlinear models that describe ultrasound-heat interactions. These models are based on a coupling of a nonlinear Westervelt--Kuznetsov acoustic wave equation to the heat equation with a pressure-dependent source term. A particular challenging feature of the system is that the acoustic medium parameters may depend on the temperature. The core of our new arguments in the \emph{a priori} error analysis lies in devising energy estimates for the coupled semi-discrete system that can accommodate the nonlinearities present in the model. Theoretically obtained optimal convergence rates in the energy norm are confirmed by the numerical experiments.

ABSTRACT. Many mathematical physics problems can be formulated as well-conditioned integral equations. A key challenge for integral equation methods is the accurate discretization of singular and near-singular integral operators. While the literature offers myriad specialized quadrature methods, these techniques are often mathematically complex and difficult to implement. In this talk, we consider integral operators on a sphere or a deformed sphere. We propose a simple and efficient trapezoidal quadrature method on spherical meshes, which combines high-order singular correction derived from generalized Euler-Maclaurin formulas and fast trigonometric interpolation based on non-uniform FFT. The new method significantly simplifies the practical implementation of integral equation methods on deformed spheres for a variety of PDEs.

ABSTRACT. This work is focused on the propagation of electromagnetic waves in R^3 described by Maxwell's equation with large wave number and Silver-Muller radiation condition. The model problem is approximated by truncating the exact Dirichlet-to-Neumann (DtN) operator into a finite sum of vector spherical harmonics. We prove the well-posedness and wave-number-explicit H(curl)-stability of the solution to truncated problem by assuming that the truncation number N satisfies N≥akR for some a> 1, where k represents the wave number and R is the radius of the physical domain. Additionally, we demonstrate that the truncated solution is exponentially close, in terms of N, to the true scattering solution. Finally, we present the hp-finite element method (hp-FEM) for the truncated problem, along with its asymptotic error estimate. Some numerical experiments are provided to validate the theoretical findings.

ABSTRACT. In this talk we present an algorithm for solving nonlinear wave equations (NLW) based on an elliptic regularisation approach conjectured by Ennio De Giorgi. This approach guarantees existence and uniqueness to a solution of a regularised nonlinear wave equation which is elliptic and converges to a solution of NLW, even in cases where the NLW is not well posed. Using this framework we will propose a Discontinuous Galerkin symmetric interior penalty(SIP) method for the NLW which inherits desirable properties from the continuous problem with particular attention being paid to time discretisation, preconditioning and stability. Finally examples are illustrated with nonlinear reaction terms and degenerate diffusion.

ABSTRACT. In this talk, we consider a scattering problem in full space from an obstacle which is penetrable but governed by a nonlinear material law. Mathematically, this is modeled by coupling the free-space linear wave equation in the (unbounded) exterior domain with a semilinear wave equation inside of the (bounded) scatterer. The discretization is then carried out by coupling a finite element method with timestepping for the nonlinear part with a convolution quadrature based boundary integral discretization in the unbounded domain. We present a rigorous a-priori analysis of the fully discrete scheme and talk about the challenges inherent in this problem.

ABSTRACT. In this talk, we address the full discretization of Friedrichs’ systems with a two-field structure, such as Maxwell’s equations or the acoustic wave equation in div-grad form; cf. [W. Dörfler et al., Wave Phenomena: Mathematical Analysis and Numerical Approximation, Springer, Cham, 2023]. We focus on a discontinuous Galerkin space discretization applied to a locally refined mesh or a small region with high wave speed. This results in a stiff system of ordinary differential equations, where the stiffness is mainly caused by a small region of the spatial mesh. When using explicit time-integration schemes, the time stepsize is severely restricted by a few spatial elements, leading to a loss of efficiency. As a remedy, we propose and analyze a general leapfrog-based scheme which is motivated by [C. Carle and M. Hochbruck, SIAM J. Numer. Anal., 60 (2022), pp. 2897–2924]. The new, fully explicit, local time-integration method filters the stiff part of the system in such a way that its CFL condition is significantly weaker than that of the leapfrog scheme while its computational cost is only slightly larger. For this scheme, the filter function is a suitably scaled and shifted Chebyshev polynomial. While our main interest is in explicit local time-stepping schemes, the filter functions can be much more general, for instance, a certain rational function leads to the locally implicit method, proposed and analyzed in [M. Hochbruch and A. Sturm, SIAM J. Numer. Anal., 54 (2016), pp. 3167–3191]. Our analysis provides sufficient conditions on the filter function to ensure full order of convergence in space and second order in time for the whole class of local time-integration schemes.

ABSTRACT. This paper introduces a kernel-based framework for estimating the conditional mode in the context of right-censored survival data. By combining modal regression with inverse probability of censoring weighting within a Gaussian reproducing kernel Hilbert space, the proposed method demonstrates robustness against outliers and heavy-tailed noise. Theoretically, we establish fast learning rates that adapt to the intrinsic dimension of the covariate space, rather than the ambient dimension, resulting in enhanced generalization for data with low-dimensional structures. Numerical experiments validate that our approach outperforms existing robust alternatives and exhibits strong predictive performance on real-world datasets.

ABSTRACT. We propose sparse identification of port-Hamiltonian systems (SIPHy) enabling structure preserving symbolic regression from noisy trajectory observations. Port-Hamiltonian systems provides a general framework for describing dynamical systems in terms of energy exchange, dissipation and control. Here, we present an algorithm to jointly identify the Hamiltonian as well as the dissipation and input magnitudes. One of the main challenges of system identification for differential equations is to approximate derivatives of trajectory data that has been corrupted with noise or have missing points in time. To achieve this, we introduce Hamiltonian flow splines which assembles flows of piecewise polynomial Hamiltonians to produce a smooth, differentiable trajectory necessary for sparse regression. Combining flow splines with SIPHy provides a robust method for system identification that handles both high levels of noise and missing data for a range of dynamical systems with dissipation and control.

ABSTRACT. We study the localisation of gas leaks from line-of-sight concentration measurements using a convection–diffusion model with a measure-valued source term. By imposing sparsity-promoting regularisation on this Radon measure, we recover point sources—identifying both their locations and intensities—rather than diffuse approximations. We jointly estimate the underlying physical convection (wind) and diffusion parameters. We employ a grid-free optimisation approach for reconstructing the source measure. Numerical experiments demonstrate accurate localisation, highlighting the potential of the method for practical gas emission detection.

ABSTRACT. Hydrogen is a clean and high-energy fuel, yet its safe and efficient storage remains a key obstacle to widespread adoption. Metal–organic frameworks (MOFs), with their high surface area and tunable porosity, have emerged as promising candidates for solid-state hydrogen storage. In this work, we introduce a graph-based machine learning framework for predicting hydrogen uptake in MOFs by integrating spectral graph theory with data-driven modeling. Molecular structures are represented as weighted graphs from which we extract 20 graph-based descriptors─including Laplacian spectral features, degree statistics, and Zagreb indices─that capture both topological and geometric characteristics of the framework. These interpretable descriptors are used to train multiple regression models on a data set of 3300 MOFs from the Cambridge Structural Database. The XGBoost regressor achieved the highest performance in predicting hydrogen uptake, with a coefficient of determination (R2) of 0.737, RMSE of 0.850% wt, and MAE of 0.433% wt for gravimetric uptake (UG); and a coefficient of determination (R2) of 0.698, RMSE of 4.467 g H2/L, and MAE of 3.045 g H2/L for volumetric uptake (UV). Beyond accurate prediction, the framework enables inverse materials design by identifying graph-based motifs that contribute to improved storage capacity. This integration of chemical graph theory with machine learning provides a scalable, interpretable, and computationally efficient pathway for the discovery of next-generation MOFs tailored for hydrogen storage and other clean energy applications.

ABSTRACT. Accurate simulation of time-dependent PDEs typically requires fine discretizations, which may be computationally expensive. In practice, simplified or coarse solvers are often employed, leading to accumulated discretization errors. In this talk, we investigate a residual-based correction framework for numerical time-stepping solvers.Given coarse numerical solutions and high-fidelity reference solutions, we define the residual as their discrepancy. While residual correction has been explored in certain ODE settings, extending it to time-dependent PDE settings is challenging due to the high dimensionality of spatially discretized systems and the need for stable prediction beyond the calibration regime. In several problem settings, we observe that the residual snapshot matrix exhibits a rapidly decaying singular-value spectrum, suggesting an approximately low-rank spatial structure. Based on this observation, we construct a POD basis and represent the residual using reduced coefficients and discuss preliminary analysis of structural conditions under which such low-rank behavior may arise. We model the evolution of these reduced coefficients using a neural network trained on an offline calibration interval. The learned model is then rolled out forward in time to predict future residuals and to form an additive, non-intrusive correction of the baseline solver output. Numerical experiments for some time-dependent PDE benchmarks indicate improved accuracy over limited prediction horizons, while also illustrating challenges related to extrapolation and changes in solution behavior. Extensions to incremental updates of the residual basis and reduced-order dynamics in a streaming-data setting are also outlined.

ABSTRACT. We consider a 3D isentropic Euler system under the assumption of spherical symmetry and coupled with self gravity, using a modification of van der Waals real gas pressure law to incorporate realistic physical effects. This model provides a compact yet physically meaningful framework for describing stellar structure, gravitational collapse, outer stellar layers, and spherically symmetric accretion flows. It captures both hydrodynamic behaviour and realistic equation of state effects, making it suitable for applications in theoretical and computational astrophysics. The solution is presented in terms of classical wave patterns as well as through a numerical scheme. We illustrate the performance of the numerical method on several examples, highlighting the influence of adding a density dependent viscosity and showing how the system transitions from hyperbolic to parabolic behaviour when viscosity is introduced, all within a finite volume framework adapted to spherically symmetric coordinates.

ABSTRACT. Elongated filaments in slow viscous flows are difficult to resolve computationally. The separation of length scales enforces a high resolution in direct approaches, while the asymptotically the leading behaviour is that of an infinite cylinder which has no solution due to Stokes paradox. Classical slender-body theory approaches use matched asymptotics on an isolated fibre to overcome the paradox and reduce the problem to line integral. However, the addition of other bodies/boundaries significantly complicates the problem, typically preventing the classical slender-body theory approach and therefore again requiring full numerical resolution. Yet interacting fibres occur in many biological and industrial applications like the swarming of bacteria, organisation of cells and the manipulation of paper pulp. This presentation discusses how through careful manipulation and application of asymptotic approaches to the boundary integrals for Stokes flow, slender-body theory-like approaches for fibres by plane walls can be created.

ABSTRACT. The numerical analysis of closures and stabilizations for reduced-order models (ROMs) of convection-dominated buoyancy-driven flows is relatively scarce. In this work, we take a step in this direction and conduct numerical analysis of the evolve-filter-relax ROM (EFR-ROM), which uses spatial filtering to stabilize ROMs for convection-dominated flows. We establish stability and derive an a priori error bound for the EFR-ROM. Demonstrated in both 2D unsteady convection and 2D Rayleigh–Benard convection test problems, our numerical investigation shows that the theoretical convergence rates are recovered numerically, and that the EFR-ROM yields more accurate velocity and temperature solutions, as well as the kinetic energy and Nusselt number, compared to the standard Galerkin ROM (G-ROM).

ABSTRACT. The boundary condition of fluid-flow is one of the most important factors to determine its hydrodynamic behaviors. In microfluidic systems, boundary slip may have a significant effect on the performance of such system. As an effective technique to catch the boundary slip phenomenon, numerical method provides some theoretical guidance for related experimental research. In this talk, we analyze the lowest-order finite volume method for the Stokes equations with a nonlinear slip boundary condition of friction type, which is used to describe the flow in the blood vessel of arteriosclerosis, as well as the possible slip phenomena. Due to the subdifferentiability of such boundary condition, these models can be characterized by variational or hemivariational inequalities. We will design some stable and efficient finite volume schemes, and establish priori error analyses for such variational inequalities. Numerical tests are reported to verify the theoretical results.

ABSTRACT. A three-dimensional spherically symmetric model for the flow of a compressible real micropolar gas confined between two concentric spheres is studied. Micropolarity is employed to describe local micro-movements in fluids with pronounced particulate structure, such as biological fluids, polluted air, or combustion gases. In this framework, micropolar effects and real gas thermodynamics are unified within a single compressible fluid model.

The system couples the conservation of mass, linear momentum, angular momentum, and energy, and includes microrotation effects, couple stresses, and heat conduction. The pressure is described by a non-ideal, thermodynamically consistent equation of state depending on temperature and density, where the density dependence is nonlinear, allowing deviations from ideal gas behavior to be captured.

Working in a Lagrangian formulation and assuming strictly positive initial density and temperature, the global-in-time existence of spherically symmetric generalized solutions is established. The long-term behavior of solutions is analyzed, and convergence toward a stationary equilibrium state is proven. Furthermore, regularity is investigated, and it is shown that Hölder continuity of the initial data is preserved for all positive times.

ABSTRACT. This work presents communication-avoiding implementations of an asynchronous discontinuous Galerkin (ADG) method to solve high-dimensional (up to six dimensions) advection and Vlasov-Poisson equations in the matrix-free finite-element library hyper.deal. In such high-dimensional kinetic simulations, the exchange of five-dimensional surface data results in extremely large communication volumes, making inter-process data movement a dominant bottleneck, distinguishing them from typical low-dimensional, latency-bound problems. To address this challenge, we employ a communication-avoiding algorithm (CAA) based on the recently developed asynchronous discontinuous Galerkin (ADG) method in which communication and synchronization can be relaxed at a mathematical level. The method significantly reduces data movement among processing elements (processes) by skipping halo exchanges for a finite number of time steps and communicating only at fixed intervals. To retain accuracy despite delayed boundary data, we use asynchrony-tolerant (AT) fluxes that combine flux information from multiple past time levels at process-boundaries. We also present a new variant of CAA with improved stability properties. We first illustrate how standard synchronous DG implementations in hyper.deal become communication-bound at scale, with communication volume overwhelming computation as the number of processes increases. We then detail our CAA strategy and its integration into the matrix-free advection operator. Strong and weak scaling experiments for six-dimensional advection and Vlasov-Poisson simulations demonstrate that communication-avoiding implementations significantly reduce communication overhead and improve scalability compared to the synchronous baseline. These results highlight the effectiveness of the asynchronous DG method for extreme-scale high-dimensional simulations on modern distributed-memory systems.

ABSTRACT. We show that a penalty-free asymmetric Nitsche's method using Nédélec edge elements in problems with curl-curl operators for weak imposition of tangential Dirichlet boundary conditions is stable. The central result is an inf-sup stability estimate for the asymmetric bilinear form under an isolated patch condition on the tetrahedral mesh. Applications to a Maxwell-type problem and a magnetic advection-diffusion problem are discussed.

ABSTRACT. We investigate a coupled hyperbolic-parabolic system modeling thermoelastic diffusion (resp. thermo-poro\-elast\-icity) in plates, consisting of a fourth-order hyperbolic partial differential equation for plate deflection and two second-order parabolic partial differential equations for the first moments of temperature and chemical potential (resp. pore pressure). The unique solvability of the system is established via Galerkin approach, and the additional regularity of the solution is obtained under appropriately strengthened data. For numerical approximation, we employ the Newmark method for time discretization of the hyperbolic term and a continuous interior penalty scheme for the spatial discretization of displacement. For the parabolic equations that represent the first moments of temperature and chemical potential (resp. pore pressure), we use the Crank--Nicolson method for time discretization and conforming finite elements for spatial discretization. The convergence of the fully discrete scheme with quasi-optimal rates in space and time is established. The numerical experiments demonstrate the effectiveness of the 2D Kirchhoff--Love plate model in capturing thermoelastic diffusion and thermo-poroelastic behavior in specific materials. We illustrate that as plate thickness decreases, the two-dimensional simulations closely approximate the results of three-dimensional problem. Finally, the numerical experiments also validate the theoretical rates of convergence.

ABSTRACT. The nonlinear least-squares method addresses the approximation of smooth solutions to the Dirichlet Monge-Ampère problem by decoupling it into (i) a pointwise nonlinear minimization problem and (ii) a linear biharmonic variational problem. We derive a convergence framework for this iterative algorithm by reformulating the method as an alternating projections scheme in Sobolev spaces. Moreover, for the linear variational problem, we derive an equivalence to a biharmonic problem with Navier boundary conditions and solve it via mixed piecewise-linear finite elements. Reformulating this as a coupled second‐order system, we derive a priori and a posteriori P_1‐finite‐element error estimators and we design a robust adaptive‐mesh‐refinement strategy. Numerical tests confirm that the L^2 error scales as O(h^2), the H^1 error as O(h), and, using a post-processing strategy to approximate the Hessian, the H^2 error as O(h). Finally, we demonstrate the effectiveness of our a posteriori indicators in guiding mesh refinement. The talk is based on joint works with A. Caboussat, M. Picasso and M. Sorella.

ABSTRACT. In this report,I present a cell-centered Lagrangian discontinuous Galerkin method for solving resistive magneto-hydrodynamics. The equations are solved using an implicit-explicit method. The right-hand side of the equations is classifed as a hydrodynamic contribution, which is solved explicitly, and a magnetodynamic contribution which is solved implicitly.The conservative variables are discretized using Taylor basis functions within the reference element, and for magnetic feld discretization, we employ a locally divergence-free basis, and it is transformed to the reference element using the Piola transformation. Nodal velocities are determined through an approximate Riemann solver. Curvilinear mesh is achieved through basis function deformation in physical space. Numerical experiments demonstrate the accuracy and robustness of the method. (This work appears in Journal of Computational Physics. 551,15, 2026, 114698.)

ABSTRACT. We present a goal-oriented time-adaptive discretization method for linear index-1 port-Hamiltonian differential-algebraic equations (pH-DAEs). While structure-preserving integrators enforce energy balance exactly, they do not quantify discretisation error in other quantities of interest. We treat violations of the discrete energy balance as a computable goal functional and employ the Dual Weighted Residual (DWR) framework to derive a posteriori error representations.

The formulation uses a reduced variational system obtained by Schur-complement elimination of the algebraic variables. On this reduction, we define a differentiable goal functional from the squared energy-balance residual and formulate the associated adjoint problem, which exhibits impulsive jump conditions at discrete time nodes. Temporal discretisation is performed with a discontinuous Galerkin method of order zero (dG(0)), equivalent to the implicit Euler method.

Because the backward adjoint solve necessitates an additional solve backwards in time, we introduce a parallelisable Block-Jacobi sweep approximation and prove that the resulting iteration is a strict contraction in the E11​-norm under a coercivity assumption on the Schur complement. We verify that the contraction bound is sharp at the boundary of this assumption. Numerical experiments on a transmission-line pH-DAE benchmark with 199 differential variables confirm that the adaptive strategy achieves the expected O(N^{−3}) convergence rate for the quadratic QoI and reaches prescribed error tolerances using up to 91% fewer time intervals than uniform refinement, with bounded effectivity indices throughout the adaptive process.

ABSTRACT. Many physical laws, such as constitutive relations for flow in porous media, and the electromagnetic laws underlying Maxwell's equations, naturally arise in the framework of H(div), H^1 and L^2 function spaces. This talk presents a unified arbitrary-order Bernstein-Bézier finite element discretization of these spaces in 2D that can be applied across a wide range of such problems. The construction of the compatible bases is outlined, and key approximation and compatibility properties are highlighted. The implementation of these bases in code is described, emphasizing optimal complexity algorithms and their role as reusable building blocks for assembling mixed Poisson and Darcy flow formulations, Maxwell source problem, and Maxwell eigenvalue problem. Practical aspects and numerical examples with visualizations are presented to illustrate the flexibility, efficiency, and performance of the approach.

ABSTRACT. This abstract presents the work related to the paper [1] published in the \emph{Journal of Computational Physics} (Vol. 542, 2025). We propose an a leapfrog time stepping finite element scheme for simulating surface plasmon polaritons along complex graphene sheets.

An accurate simulation of SPPs presents several unique challenges, since SPPs often occur at complex interfaces between materials of different dielectric constants and appropriate boundary conditions at the graphene interfaces are crucial. To address this, our method develops a simplified graphene model and proposes a new finite element method by treating the graphene sheet as zero-thickness.

Numerical experiments are conducted on 2D benchmark problems. The new method is based on a reformulated system of governing equations for graphene with electric and magnetic fields as unknowns. Compared to our previous work [2], the new system does not explicitly contain the induced current. Hence, our new method is more efficient in memory and computational cost.

References: [1] Li, J., Neunteufel, M., Zhu, L., A novel finite element method for simulating surface plasmon polaritons along complex graphene sheets, Journal of Computational Physics, Vol. 542, 114372, 2025. [2] Li, J., Zhu, L., Arbogast, T., A new time-domain finite element method for simulating surface plasmon polaritons on graphene sheets., Computers & Mathematics with Applications, 142, 268-282(2023).

ABSTRACT. In this work, the geometric electromagnetic Particle-in-Cell (PIC) framework, GEMPICX, is extended to solve the quasineutral, fully kinetic Vlasov-Maxwell equations and also the quasineutral, hybrid drift-kinetic Vlasov-Maxwell model, with drift-kinetic electrons and fully kinetic ions. The discretization for both models is performed using structure-preserving finite differences on primal and dual grids. Unknowns are represented as point values, edge integrals, face integrals and volume integrals, on both grids. These mimetic discretizations employ operators that exactly imitate integrated versions of vector calculus identities on discrete spaces. A discrete action principle is derived for both models, taking into account the duality between the grids. The dynamical system for such models does not involve a temporal evolution term for the electric field. It can be obtained implicitly by solving a linear system of equations at each time-step, for both models. This also circumvents the need to obtain electric potentials. For the fully kinetic model, a discretized curl-curl equation is used to implicitly obtain the electric field at every time-step. A Lagrange multiplier is used to maintain the discretized divergence of the current density at machine zero. On the other hand, for the hybrid model, a similarly derived discretized curl-curl equation is used to implicitly obtain the electric field component that is parallel to the background magnetic field. The definition of the current in the drift kinetic model is used to obtain separate equations for the perpendicular component of the electric field. An explicit split time-stepping scheme is used for the fully kinetic model. For the hybrid model, a fully explicit scheme as well as two implicit-explicit (IMEX) schemes are tested. The two models are tested by verifying the various waves obtained from their dispersion relations.

ABSTRACT. Problems in plasma physics are rich in geometric structure. This is especially true when considering kinetic models (ones which describe the evolution of a probability density function for the charged particles, together with the electromagnetic fields, e.g. the Vlasov--Poisson or Vlasov--Maxwell equations) as they possess infinite-dimensional Hamiltonian / Poisson structure. Exploiting this geometric structure has led to the development of advanced numerical methods for solving (non-dissipative) kinetic problems very accurately. 

Issues arise, however, when one is interested in adding collisions to such models, as collisions are dissipative and do not fit into the Hamiltonian structure of the problem. This can be resolved through the use of the metriplectic formulation of Morrison (1986), which is a geometric way of describing Hamiltonian systems with dissipation. This formulation can also be utilised to construct accurate numerical discretisations when modelling such coupled systems, which is a much more recent topic of study. 

In this talk, I will introduce the relevant plasma physics models and their geometric structure. I will then present the metriplectic formulation, and detail how it can be used to describe Hamiltonian systems with dissipation. Finally, I will describe how the metriplectic structure of collision operators can be used to construct discretisations in a way that preserves the quantities of interest using relevant examples.

ABSTRACT. In magnetic confinement fusion, the plasma is held in place inside of a vacuum vessel, using a set of magnetic coils. In order to describe the magnetic field inside plasma, a steady-state solution of the ideal MHD equations, a so-called MHD equilibrium, must be computed numerically. For three-dimensional plasma configurations, finding an equilibrium solution typically requires an energy minimisation, starting from a given plasma boundary and some physical parameters, such as pressure and current profiles. However, the prescribed boundary is generally not consistent with the (external) field generated by the desired coils. To ensure the coil field and plasma boundary are self-consistent, the effect of the magnetic field generated by the plasma itself on the external coil field must be taken into account. This leads to the free-boundary problem, where the equilibrium code must be used in a inverse solve, and the plasma boundary is moved iteratively until the forces from the external and internal magnetic fields balance.

We discuss a method derived from the variation of the energy of the system. In this talk we will discuss the extension of the three-dimensional equilibrium solver GVEC (Galerkin Variational Equilibrium Code) from fixed to free-boundary. For the external field we use the Boundary Integral Equation Solver for Taylor states (BIEST). We also present some advantages of GVEC, such as its ability to use coordinate systems beyond those typically used in equilibrium codes, allowing new plasma configurations to be explored.

ABSTRACT. In this talk, we discuss the dynamics of self-similar solutions of the one-dimensional Schrödinger map and Landau–Lifshitz–Gilbert (LLG) equations from a numerical perspective. The LLG equation describes magnetization dynamics in ferromagnetic spin chains and interpolates between the Schrödinger map equation and the heat flow for harmonic maps via the Gilbert damping parameter $\alpha\in [0, 1]$.

In the first part, we revisit the self-similar solutions of the Schrödinger map, which form a one-parameter family developing corner singularities in finite time. We discuss the stereographic projection framework and boundary treatments that allow accurate approximation of the singularity formation on finite computational domains. In the second part, we present ongoing numerical work on the self-similar solutions of the LLG equation with $\alpha\in [0, 1]$, building on the existing analytical results. The introduction of damping breaks the time-reversibility of the Schrödinger map, leading to genuinely new phenomena that we explore numerically.

ABSTRACT. Hug is a recently proposed iterative mapping used to design efficient updates in Markov chain Monte Carlo (MCMC) methods. Hug generates proposals that remain very close to hypersurfaces (level sets) of constant probabilty density. We analyse a generalization of Hug from hypersurfaces to manifolds of arbitrary dimensions, not necessarily arising in a sampling context. The analysis is based on interpreting, in a nonstandard way, Hug as a consistent discretization of a system of differential equations with a rather complicated structure. The proof of convergence of this discretization includes a number of unusual features we explore fully. In particular there is a sumpraconvergence phenomenon, whereby second order convergence is achieved with first order consistency. We uncover and discuss an unexpected property of the solutions of the underlying dynamical system that manifest itself by the existence of Hug trajectories that fail to cover the manifold of interest. Joint work with Christophe Andrieu (Bristol).

ABSTRACT. We introduce a Nemytskii neural operator framework for nonlinear model reduction of parametrized steady-state partial differential equations. The method generalizes reduced basis approaches by replacing linear combinations of basis functions with a structured nonlinear mapping realized through a pointwise Nemytskii operator acting on fixed feature functions. Feature functions are learned offline via nonlinear dimension reduction from high-fidelity snapshots, and a hypernetwork maps model parameters to a lightweight reconstruction network, which is further refined online using physics-informed residual minimization. The Nemytskii structure preserves analytical regularity and enables efficient evaluation of spatial and parametric derivatives, leading to fast online adaptation. Numerical experiments demonstrate that the proposed method consistently outperforms linear model reduction techniques, particularly for complex solution manifolds.

ABSTRACT. Using an entropic formulation of the mean-field optimal control problem, we recast its convergence problem as the classical limit of mean-field Gibbs measures. This framework yields optimal convergence rates under path-dependent interactions and establishes sharp propagation of chaos estimates.

ABSTRACT. Stochastic unraveling schemes are powerful computational tools for simulating Lindblad equations, offering significant reductions in memory requirements. However, this advantage is accompanied by increased stochastic uncertainty, and the question of optimal unraveling remains open. Therefore, we investigate unraveling schemes driven by Brownian motion or Poisson processes and present a comprehensive parametric characterization of these approaches. For the case of a single Lindblad operator and one noise term, this parametric family provides a complete description for unraveling scheme with pathwise norm-preservation. We further analytically derive dynamically optimal quantum state diffusion (DO-QSD) and dynamically optimal quantum jump process (DO-QJP) that minimize the growth rate of the variance of an observable locally in time. Numerical results demonstrate that the proposed DO-QSD scheme may achieve substantial reductions in the variance of observables and the resulting simulation error.

ABSTRACT. This work presents a comprehensive analysis of coupled stochastic van der Pol oscillators, a paradigm for understanding synchronization, bifurcations, and chaos in nonlinear systems subject to random fluctuations. The system comprises two or more oscillators with nonlinear damping, linear diffusive coupling, and additive Gaussian white noise. We develop a unified framework that systematically connects global bifurcations, synchronization phenomena, and chaotic dynamics within a single coherent stochastic model. We explore the stochastic dynamics of coupled van der Pol oscillators by seamlessly blending theoretical principles with in-depth numerical simulations. This integrated approach forms a robust framework for analysis, with essential phenomena clearly depicted in the accompanying figures. We then extend this framework to a comprehensive investigation of large networks, focusing on their continuum limit, emergent pattern formation, the role of noise, and the onset of collective chaos.

ABSTRACT. In this talk, I will discuss the effect of coarse-graining (CG) on rare-event timescales in stochastic dynamics, especially in Langevin-type dynamics which are widely used as thermostats in molecular simulation. I will begin by presenting analytical results on how much rare event timescales change under coarse-graining within the framework of the Zwanzig projection formalism. Then, I will show how machine-learning methods to approximate the Koopman generator can be used to implicitly learn a predictive model for a coarse-grained dynamics. Finally, I will show how Koopman operator models can be practically used to assess and compare dynamical properties of CG models.

ABSTRACT. Phase-field models provide a fundamental continuum framework for describing phase separation and pattern formation in many physical and biological systems. Their predictive capability depends critically on constitutive quantities such as the bulk free-energy density and interfacial thickness parameter, which are often unknown and must be inferred from limited observations. In this work, we introduce an extended pseudo-spectral physics-informed neural network (ESPINN) framework for the inverse identification of phase-field models from transient snapshot data. The proposed method simultaneously reconstructs the bulk chemical potential and unknown gradient coefficients directly from dynamically evolving structures.

Numerical experiments show that ESPINN accurately recovers both the functional form of the free energy and the interfacial thickness parameter. Remarkably, substantial constitutive information can be extracted even from a single snapshot pair, while additional snapshots improve robustness and reduce variance across training runs. The framework remains stable in the presence of noise, with reconstruction accuracy improving as more observations are incorporated. These results highlight ESPINN as a data-efficient and physically consistent approach for learning constitutive structure in continuum models of phase separation.

ABSTRACT. Coarse-grained models have been a frequently used method for overcoming the computational bottleneck in traditional atomistic molecular dynamics simulations, at the cost of reduced transferability and accuracy. In recent times, however, developments in machine learning have enabled the construction of bottom-up coarse-grained models, with a hybrid physics-neural network architecture, that are transferable, data-efficient and quantitatively comparable to atomistic models.

In this talk I will give an overview of the state of the art in machine-learned coarse-graining (MLCG), showing that these models are not limited to globular proteins but can be extended to include metal ions, small molecules and even other biomolecules such as simple liquids and lipids. Furthermore, I will show that the neural network is not limited to the SchNet architecture but can be selected from a set of modern neural network potentials, each bringing some advantages and disadvantages.

I will conclude by discussing some of the future directions like the relationship between MLCG and generative models, and how to integrate the MLCG framework with experimental data.

ABSTRACT. Memory effects are a central feature of many-body open quantum dynamics and can be described through both continuous and discrete representations. In this talk, I will focus primarily on continuous memory kernels, which provide a natural operator-level description of non-Markovian dynamics, while using transfer-tensor representations as complementary viewpoints for comparison and interpretation.

I will discuss recent progress toward computing and characterizing continuous memory kernels from reduced dynamical information. The kernel construction combines reduced dynamical maps with a linked-diagram organization inspired by Inchworm formulations. This perspective provides access to memory objects beyond direct time propagation while retaining a close connection to the underlying reduced dynamics.

A main theme of the talk is the relation between continuous kernels and transfer-tensor representations of memory. In particular, I will discuss how transfer tensors relate to continuous memory kernels, and why this relation is subtler than a naive discretization picture may suggest. I will also comment on a continuous transfer-tensor viewpoint and its connection to the more familiar discrete formulation. The goal is to clarify the complementary roles of kernel-based and transfer-tensor descriptions of memory and to highlight possible directions for more effective numerical treatments of non-Markovian many-body open quantum systems.

ABSTRACT. Electronic excited states play a central role in determining the optical and reaction properties of atoms and molecules. Within variational approximate theories, state-specific methods estimate excited states as saddle points of associated energy landscapes defined on Riemannian manifolds. However, these methods often face convergence difficulties due to the intrinsic instability of saddle points, or explicitly involve the computationally expensive Riemannian Hessian. This talk introduces a constrained saddle dynamics on Riemannian manifolds that requires at most (estimates of) Riemannian Hessian-vector products. It couples a position governed by a reflected Riemannian gradient flow with a direction dynamics tracking the lowest invariant subspace of the Riemannian Hessian. By exploiting the underlying Riemannian geometry, we rigorously establish local theoretical results for the dynamics and the resulting discretized algorithm. The effectiveness of the algorithm is examined numerically on standard benchmark systems within restricted Hartree-Fock and complete active space self-consistent field theories.

ABSTRACT. Estimating thermal expectation values of quantum many-body systems is a central challenge in physics, chemistry, and materials science. Standard quantum Gibbs sampling protocols address this task by preparing the Gibbs state from scratch after every measurement, incurring a full mixing-time cost at each step. Recent advances in single-trajectory Gibbs sampling substantially reduce this overhead: once stationarity is reached, measurements can be collected along a single trajectory without re-thermalizing, provided the measurement channel preserves the Gibbs ensemble. However, explicit constructions of such non-destructive measurements have been limited primarily to observables that commute with the Hamiltonian. In this work, we fundamentally extend the single-trajectory framework to arbitrary, non-commuting observables. We provide two measurement constructions that extract measurement information without fully destroying the Gibbs state, thereby eliminating the need for full re-mixing between samples.

ABSTRACT. An efficient excited state method, named xCDFCI, in the configuration interaction framework, is proposed. xCDFCI extends the unconstrained nonconvex optimization problem in CDFCI to a multicolumn version, for low-lying excited states computation. The optimization problem is addressed via a tailored coordinate descent method. In each iteration, a determinant is selected based on an approximated gradient, and coefficients of all states associated with the selected determinant are updated. A deterministic compression is applied to limit memory usage. We test xCDFCI applied to H2O and N2 molecules under the cc-pVDZ basis set. For both systems, five low-lying excited states in the same symmetry sector are calculated together with the ground state. xCDFCI also produces accurate binding curves of carbon dimer in the cc-pVDZ basis with 10−2 mHa accuracy, where the ground state and four excited states in the same symmetry sector are benchmarked.

ABSTRACT. The numerical solution of kinetic equations is challenging due to the high dimensionality of the underlying phase space. In this work, we develop a dynamical low-rank method based on the projector-splitting integrator in tensor-train (TT) format. The key idea is to discretize the three-dimensional velocity variable using tensor trains while treating the spatial variable as a parameter, thereby exploiting the low-rank structure of the distribution function in velocity space. In contrast to the standard step-and-truncate approach, this method updates the tensor cores through a sweeping procedure, allowing the use of relatively small TT-ranks and leading to substantial reductions in memory usage and computational cost. We demonstrate the effectiveness of the proposed approach on several representative kinetic equations.

ABSTRACT. The radiative transfer equation (RTE) is a fundamental model for describing the transport, absorption, and scattering of radiation in various applications, ranging from medical imaging to nuclear physics. However, its dependence on both spatial and angular variables leads to high-dimensional problems that pose significant computational challenges. In this talk, I will present a low-rank tensor product framework for efficiently approximating stationary radiative transfer problems in plane-parallel geometry. The approach exploits the tensor product structure of the phase space to formulate the discrete RTE as a short sum of Kronecker products. This structure enables the use of a preconditioned and rank-controlled Richardson iteration in Hilbert spaces, allowing for rigorous control of both error and tensor rank. A key component of the method is a preconditioner constructed via exponential sum approximations, which is compatible with the low-rank tensor structure and transforms the problem into an equivalent formulation in a Euclidean metric. The resulting algorithm is able to identify quasi-optimal ranks during iteration automatically. We present numerical examples that indicate the computational performance.

ABSTRACT. Low-rank methods have emerged as a promising strategy for reducing the memory footprint and computational cost of discrete-ordinates discretizations of the radiative transfer equation (RTE). However, most existing rank-adaptive approaches rely on rank-proportional space augmentation, which can negate efficiency gains when the effective solution rank becomes moderately large. To overcome this limitation, we develop a rank-adaptive sweep-based source iteration with diffusion synthetic acceleration (SI–DSA) for the first-order steady-state RTE. The core of our method is a sweep-based inner-loop iterative low-rank solver that performs efficient rank adaptation by augmenting the basis space to a dimension slightly larger than the effective rank. In each inner iteration, the spatial basis is augmented with a small, rank-independent number of basis vectors without truncation, while a single truncation is performed only after the inner loop converges. Efficient rank adaptation is achieved through residual-based greedy angular subsampling strategy together with incremental updates of projection operators, enabling non-intrusive reuse of existing transport-sweep implementations. In the outer iteration, a DSA preconditioner is applied to accelerate convergence. Numerical experiments show that the proposed solver achieves accuracy and iteration counts comparable to those of full-rank SI–DSA while substantially reducing memory usage and runtime, even for challenging multiscale problems in which the effective rank reaches 30-45% of the full rank.

ABSTRACT. High-dimensional parametric partial differential equations (PDEs) often generate tensor-valued solution data whose storage and computational costs scale exponentially with dimension. Tensor-train (TT) decompositions provide efficient low-rank representations for such problems and have become increasingly important in large-scale scientific computing. Existing TT-based compression approaches, however, typically construct reduced bases from full tensor snapshots and subsequently compress them via TT-SVD, resulting in prohibitively expensive offline costs in high dimensions. Moreover, in many modern low-rank solvers, solution tensors are already represented in compressed formats such as TT, making reconstruction-based approaches both inefficient and unnecessary. In this work, we develop an incremental reduced basis construction algorithm that operates directly on sequential TT-formatted snapshots. Rather than reconstructing full tensors, the proposed method incrementally updates the reduced basis through core-wise projection, orthogonalization, and enrichment procedures acting entirely on TT cores. Incoming snapshots are projected onto the current TT subspace, while only complementary components are retained through local residual-based updates. As a result, the proposed approach enables reduced basis construction directly from low-rank tensor data while avoiding exponential scaling with respect to dimension. Theoretical analysis establishes approximation error bounds and demonstrates that the computational complexity scales linearly with the tensor dimension. Numerical experiments further demonstrate the accuracy, effectiveness, and computational efficiency of the proposed method.

ABSTRACT. Modern simulation-based design workflows often involve computationally intensive steps such as design optimization, uncertainty quantification or system-level coupled simulations. In order to reduce the computational burden of these steps for electric circuit design, we propose a machine learning approach that substitutes a computationally cheaper surrogate model for the circuit solver. The approach is based on modified nodal analysis (MNA), which is one of the most popular circuit formulations and forms the basis for many circuit solvers. A peculiarity of MNA is that it leads to systems of differential-algebraic equations (DAEs), i.e. systems that consist of both differential and purely algebraic equations. DAEs can in general be interpreted as constrained dynamical systems, such that their solutions do not evolve freely, but rather fulfill additional (hidden) constraints, which may or may not be directly visible from the algebraic equations. When it comes to MNA, its DAE structure is well-understood from a mathematical point of view and recently, we have leveraged this understanding to derive a topological decoupling of MNA. The decoupling separates the equations of MNA into their purely differential and algebraic parts, such that the differential part accounts for the dynamical evolution of the system, while the algebraic part explicitly describes all (hidden) constraints. Due to the topological nature of the decoupling, it can be implemented efficiently and therefore used as part of larger workflows. We leverage the decoupling to develop a surrogate model that provide physics-consistent, i.e. constraint fulfilling, predictions of circuit solutions, while simultaneously reducing the learning effort compared to conventional approaches.

ABSTRACT. Port-Hamiltonian (pH) systems provide a natural and unifying language for incorporating physical principles into machine learning models. Through energy-based representations, the port-Hamiltonian framework captures key properties such as interconnection, passivity, and stability, while also allowing the seamless integration of partial prior knowledge which is an essential feature when modeling complex systems with incomplete physical information. In this work, we build on our recent advances in learning nonlinear pH systems from input–state–output data and present a unified framework for identifying such systems using neural networks as structured function approximators. Incorporating physical structure into the hypothesis space improves data efficiency and yields models that are more accurate, physically consistent, and reliable for long-term prediction than purely data-driven approaches. In addition, we investigate how the performance of physics-informed system identification scales with available learning resources. Inspired by neural scaling laws observed in large-scale machine learning, we empirically study how identification accuracy depends on factors such as dataset size, model capacity, and computational budget. The resulting scaling laws provide practical guidance for forecasting achievable model performance, and designing data acquisition strategies.

ABSTRACT. In the numerical integration of differential equations, preserving the geometric structure and invariants of the original system is essential. Inspired by Galerkin variational integrators, we propose a new class of methods, nonlinear variational integrators—designed to enable large time-step integration using flexible, nonlinearly parameterized function representations. Unlike classical approaches that approximate generalized coordinates in the (space-)time domain using polynomials, our method employs more general nonlinear function spaces, such as those parameterized by neural networks or other expressive models. The degrees of freedom are determined by the corresponding discrete Euler–Lagrange (field) equations. Like other variational integrators, our approach ensures long-time stability, energy preservation, and a reduction in the number of degrees of freedom required over a given time interval.

ABSTRACT. The Met Office is developing a new weather and climate model, LFRic, to replace the current operational model, the Unified Model (UM) by the end of the decade. LFRic is designed to have better scalability and performance portability, whilst keeping the scientific skill the same as the UM [1]. The improved performance will enable kilometre-scale resolution, which is significant, as at this resolution deep convection is resolved and no longer requires a physical parameterisation scheme to represent its effects. The model is written in modern Fortran and uses a Domain Specific Language and Compiler, PSyclone, to generate the parallel code, using both MPI and OpenMP. The dynamical core, known as Gung Ho, contains most of the distributed memory communications and is the most time-consuming component.

We present performance scaling results for a purely dry dynamical core, as well as a performance analysis study across different supercomputers and compilers. We show the scaling of a global 1.5km resolution model up to 3,072 nodes (393,216 CPU cores). Further scaling comparisons and performance analyses of the model on different systems are presented, and the model is shown to have good scaling behaviour and overall computational performance. We present a detailed analysis of the communication and computational performance and find evidence that at very large scale, the scaling is eventually limited by synchronisation overheads in the solver.

[1] S.V. Adams, R.W. Ford, M. Hambley, J.M. Hobson, I. Kavcic, C.M. Maynard,

T. Melvin, E.H. Mueller, S. Mullerworth, A.R. Porter, M. Rezny, B.J. Shipway,

and R. Wong. LFRic: Meeting the challenges of scalability and performance porta-

bility in weather and climate models. J. of Parallel and Distributed Computing,

132:383–396, 2019.

ABSTRACT. Various formulations of the compressible Euler equations are currently being explored for the development of modern dynamical cores, differing in the choice of thermodynamic variables, such as potential temperature, entropy, or total energy, as well as in the form of the momentum equation, e.g. the conservative form or the vector-invariant form. In this talk, we present and compare novel structure-preserving methods for different formulations of the compressible Euler equations within the flux-differencing discontinuous Galerkin spectral element method (DGSEM) framework. For each formulation, we construct numerical fluxes ensuring conservation of key physical invariants, including entropy and total energy on general curvilinear meshes. Particular focus is given to the formulation with the potential temperature as the thermodynamic variable and the vector invariant form. We discuss the structural differences between the formulations, in particular the importance of a proper discretisation of non-conservative terms, which naturally leads to different non-conservative split forms, well-balanced schemes and increased robustness. Several atmospheric test cases are presented to assess the theoretical findings and to compare the different formulations.

ABSTRACT. We present a high-order entropy-stable discontinuous spectral-element method for the shallow water equations on the sphere. The method is formulated in a flux-differencing framework on tensor-product elements, using summation-by-parts operators and entropy-conserving fluxes to ensure discrete conservation and stability. A covariant formulation is employed to represent the spherical geometry exactly, avoiding the need to enforce metric identities. Numerical results demonstrate the method’s accuracy and robustness for standard benchmark problems.

ABSTRACT. In this work, we numerically study the bifurcation and stability properties at varying of the Reynolds number of a barotropic flow in a periodic channel. Although this model is highly idealized, it is able to provide relevant physical insights about large-scale phenomena, such as atmospheric blocking and sudden stratospheric warming. Reliable results are available for this model. However, they are based on a one-dimensional linear stability analysis. Here, we work in a two-dimensional framework and, after re-obtaining some findings through numerical simulations of the linearized equations, we perform numerical simulations of the fully non linear dynamic system for the purpose of investigating how the nonlinearities affect the spatio-temporal structure of the flow, especially in configurations where codimension-2 Hopf-Hopf bifurcations occur, i.e. when two pairs of linear eigenmodes become unstable simultaneously, marking transitions from steady to unsteady complex patterns.

ABSTRACT. In this study, we address a singularly perturbed parabolic turning point problem of convection-diffusion-reaction type involving a large space delay argument. When the convection term dominates over the diffusion term, the solution exhibits two strong boundary layers, while the presence of a delay argument in the reaction term produces a weak interior layer. An implicit finite-difference scheme based on a $\theta$-method is applied for time discretization on uniform mesh, and the streamline diffusion finite element method (SDFEM) is used for space discretization on various Shishkin-type layer-adapted meshes. The stability analysis of SDFEM is conducted by utilizing the discrete Green’s function estimates and a suitable stabilization parameter. A robust convergence of the proposed method has been proved in both spatial and temporal directions in the maximum norm. Some test problems are considered for numerical validation and to demonstrate the efficiency of the proposed numerical scheme. Additionally, a comparison of the convergence outcomes has been presented on various Shishkin-type meshes, and solution plots are provided to illustrate the impact of the delay term effectively.

ABSTRACT. This presentation addresses the numerical challenges inherent in solving the Fisher-KPP equation combining nonlocal diffusion with evolving free boundaries. The simultaneous presence of integral operators and a moving domain precludes standard parabolic solvers. We propose and rigorously compare two distinct numerical strategies: a Front-Tracking Runge-Kutta (FTRK) method and a Front-Fixing (FF) approach. The FTRK method operates on the physical domain, utilizing high-order time integration and dynamic grid adaptation to maintain spatial resolution during significant domain expansion. Conversely, the FF method employs a Landau-type transformation to map the moving boundary problem onto a fixed computational domain. This transformation converts the governing equation into a partial integro-differential equation (PIDE) with state-dependent advection terms. We validate these methods against theoretical results regarding the spreading-vanishing dichotomy and analyze asymptotic spreading rates, confirming the capability of both schemes to capture accelerated propagation driven by fat-tailed kernels. Finally, to address parametric uncertainty, we perform a variance-based global sensitivity analysis using Sobol’ indices. Results indicate that initial habitat size and boundary expansion rates are the dominant factors influencing domain evolution.

ABSTRACT. The aim of this study is to develop a numerical framework for simulating neural field dynamics on curved geometries, with a particular focus on multi-patch NURBS representations of complex surfaces such as the human cortex. We present a Computer-Aided Design (CAD) integrated approach, referred to as isogeometric collocation, for solving neural field integro-differential equations posed on multi-patch geometries. In this framework, each patch is described using NURBS and coupled through coefficient matching across patch interfaces to ensure global continuity.

The proposed methodology combines a precise geometric representation with efficient numerical techniques, including collocation-based spatial discretisation, Gaussian quadrature for integral evaluation, and semi-implicit schemes for time integration, alongside appropriate treatment of adaptation variables. In addition, efficient computation of geodesic distances on the surface is incorporated, enabling accurate modelling of spatial interactions between neural populations on curved geometries.

To demonstrate the effectiveness of the framework, we first implement and validate the method on a flat two-dimensional torus, before extending the approach to a curved three-dimensional torus geometry. This progression allows for systematic verification of the numerical scheme while highlighting its capability to handle increasingly complex curved surfaces. The aim is to provide a scalable foundation for simulating neural activity on physiologically realistic cortical surfaces derived from neuroimaging data.

ABSTRACT. The choice of the time step is important not only for stability, but also for efficiency and consistency with the equations under consideration. In the context of kinetic equations, these aspects are particularly critical. On the one hand, efficiency must be optimized because of the high dimensionality of phase space, which may range from 2 to 6 depending on the problem, leading to a very large computational cost. On the other hand, the simulation of models, especially in Plasma Physics, often involves very fast dynamics that require an appropriate time discretization in order to be accurately resolved. The aim of this talk is to present recent results on time-step selection strategies applied to a numerical scheme for the Vlasov–Poisson equation in the quasi-neutral limit. The scheme is based on a Hermite spectral decomposition in velocity and a finite-difference discretization in space, and was originally proposed in 2025 by Blaustein, Dimarco, Filbet, and Vignal. It relies on a fully implicit, L-stable DIRK time discretization, yielding an Asymptotic-Preserving method in the quasi-neutral limit, at the expense of a significant computational cost per time step. A key issue is that, although the L-stable scheme can capture the oscillating behavior in time of the Vlasov–Poisson system (particularly those of the electric field and current density) in the asymptotic regime, its dissipative nature may eventually damp these oscillations if the time step is not properly chosen. In this work, we instead adopt an IMEX approach, treating the nonlinear part implicitly and the linear part explicitly. This substantially improves the efficiency while maintaining stability, although the Asymptotic-Preserving property is lost. Combined with this semi-implicit treatment, we aim to exploit the accuracy control performed by the time-step selection technique and use it as main ingredient for recovering the fast and ample oscillations in time. This work was conducted in collaboration with prof. Francis Filbet (Institute of Mathematics of Toulouse)

ABSTRACT. An accurate prediction of the formation of composite polymer particles is vital for synthesis of high quality materials, but still not feasible due to its complexity. We present a set of Population Balance Equations (PBE) to model the kinetics of interest and attempt the prediction of particles synthesis. When experimentally grounded values of physical parameters are employed, a first difficulty is given by the great difference in magnitude of the involved variables. Thus, we propose a novel approach able to scale the PBE model to dimensionless variables, assuming a computationally tractable order of magnitude, even though experimental values of parameters are used. Then, with the aim of enhancing performance of the prospective PBE model, we investigate a possibility of reducing its complexity by neglecting the aggregation terms. In particular, we derive a quantitative criterion for locating regions of “slow” and “fast” aggregation and introduce a dimensionless PBE model of reduced complexity. When compared with the original PBE model, the resulting model demonstrates several orders of magnitude improved computational efficiency.

ABSTRACT. Modelling macroscopic thermal transport requires reliable upscaling from microscopic dynamics to continuum transport. This is particularly relevant for composites with macroscopic conductivities that vary widely due to variations in microscopic morphology, such as graphene-copper composites. Consequently, the description of heat transport in these composites is a multiscale modelling problem, ranging from microscopic models to continuum-scale partial differential equations.

Rigorous upscaling methods derive differential equations and transport coefficients describing macroscopic transport from microscopic dynamics, under explicit assumptions about scale separation, thermostats, and boundary conditions. While powerful, there exists a gap between model classes for which this analytical control is achievable and application-relevant models. We employ Nonequilibrium Molecular Dynamics (NEMD) as an alternative upscaling approach and apply it to a graphene-copper composite. NEMD provides a computational route to upscaling, regarding microscale models as n-body systems of ordinary differential equations with arbitrarily complex potentials. Under the prescribed macrostructure, the associated continuum coefficients of the partial differential equation can then be derived from computed nonequilibrium steady states, assuming local equilibrium to define ensemble-averaged observables. Features of the microscale dynamics compatible with the imposed Fourier closure are retained in the continuum model.

We use uncertainty quantification for NEMD to analyse sensitivity to the underlying high-dimensional parameter space describing model morphology and propagate these uncertainties to continuum-scale thermal transport coefficients, such as the effective conductivity and the thermal interface (Kapitza) resistance. This yields uncertainty estimates for MD-based upscaling, clarifies the limitations of the underlying thermodynamic assumptions, and provides a framework for comparing upscaling approaches.

ABSTRACT. We analyze a nonlinear beam network model for the elastic deformation of fiber-based materials such as paper. Individual fibers are modeled as Cosserat beams, allowing for large deformations and rotations in SO(3). The equilibrium configuration is obtained by minimizing the total energy under applied loads. We prove existence and regularity of a global minimizer. The continuous problem is discretized using a conforming finite-element approximation that respects the underlying geometric structure. We establish convergence of the method. Numerical results support the theoretical findings.

ABSTRACT. Quasicrystals are materials in which the atoms are ordered (in the sense of having a well-defined diffraction pattern) but not periodic, like in a Penrose tiling. The behaviour of electrons in a quasicrystal is governed by the Schrödinger equation with a quasiperiodic potential. As a model, we compute the spectrum of the one-dimensional Schrödinger operator −(d/dx)² + V(x) with a quasiperiodic potential like V(x) = cos(x) + cos(τx) where τ is an irrational number.

There are many methods for computing the spectrum of Schrödinger operators, but they assume that either the potential is either periodic or that the eigenstates are localized. Thus, one approach is to use a periodic approximant for the potential and then (for instance) expand in a Fourier series. We will compare this with another approach, which is to expand in a Fourier–Bohr series using a basis of quasiperiodic functions. Both methods show that the spectrum has an intricate Cantor-like structure as suggested by results from functional analysis. The relative performance of the two methods depends on which part of the spectrum is considered.

ABSTRACT. Accurate coarse-grained (CG) models depend on reliable estimates of the potential of mean force (PMF), but standard force matching approaches require expensive unbiased sampling from the Boltzmann distribution and often fail to accurately capture transition regions. We show that performing enhanced sampling directly in CG coordinate space, combined with unbiased recomputation of forces, speeds up exploration, improves convergence, and leads to better PMF estimation. Building on this, we introduce Coarse-Grained Boltzmann Generators, which combine flow-based generative models with importance sampling in CG space guided by a learned PMF. Our framework enables scalable and statistically exact sampling of complex molecular systems.

ABSTRACT. In many applications, systems of ordinary differential equations arise whose right-hand side can be split into N summands with respect to properties such as stiffness, nonlinearity, dimensionality, dynamical behavior, computational cost, or geometric structure. Splitting methods constitute a common choice for the numerical integration of such systems and are widely used in geometric numerical integration. For general N-split systems, state-of-the-art methods include the Lie—Trotter and the Strang splitting, as well as higher-order composition methods based on these schemes.

We introduce a hierarchical splitting approach that provides a design principle for splitting methods for N-split systems by iteratively applying splitting methods for two-split systems. The resulting class of hierarchical splitting methods contains the state-of-the-art approaches as special cases.

We further extend this approach with multiple time-stepping techniques, thereby turning it into a valuable framework at the intersection of geometric numerical integration and multirate integration.

Numerical results confirm the theoretical findings and demonstrate the computational efficiency of hierarchical splitting methods.

ABSTRACT. Port-Hamiltonian systems provide a powerful framework for the modeling of multiphysics systems with an inherent energy-based structure. Their geometric formulation naturally incorporates conservation laws, interconnection patterns, and dissipation mechanisms, making them an attractive class of systems for structure-preserving numerical simulation. Moreover, port-Hamiltonian systems offer a natural setting for coupling different subsystems, since the resulting interconnected system again inherits a port-Hamiltonian structure.

In this work, we study (higher-order) splitting methods for the time integration of coupled port-Hamiltonian systems. Splitting methods decompose the full system into simpler subproblems whose subflows can be computed more efficiently and often in a structure-preserving manner. Higher-order integrators are then obtained by suitably composing these subflows. A central difficulty, however, is that classical higher-order splitting schemes typically involve negative time steps. While such compositions are standard in the conservative setting, they may destroy the dissipative flow in the port-Hamiltonian context and thereby lead to instabilities.

To address this issue, we discuss splitting strategies tailored to coupled port-Hamiltonian systems and analyze how the underlying structural properties can be preserved at the discrete level. In particular, for decompositions into energy-conserving and dissipative parts, we consider commutator-based methods as an approach to achieving higher-order accuracy while preserving the structural properties. Based on the coupled port-Hamiltonian formulation, we compare different decompositions of the right-hand side with respect to computational efficiency, unconditional stability, and structure preservation.

ABSTRACT. In this talk, we present a numerical solver for Fokker–Planck equations based on linear combinations of Gaussian functions. This class of functions provides a natural representation for probability densities and therefore offers a convenient framework for approximating the evolution of solutions to the Fokker–Planck equation.

Our approach relies on an operator splitting strategy: the transport step is handled by an approximation of the action of the flow on Gaussians, while the diffusion step is treated analytically. The method adaptively adjusts the number of Gaussian components to during the computation to the given error tolerance.

We provide a convergence analysis of the proposed scheme and illustrate its performance through numerical experiments.

ABSTRACT. In this article we investigate the numerical solution of a scalar semilinear stochastic delay differential equation (SDDE) where the linear instantaneous feedback and nonlinear delayed feedback terms are perturbed by a pair of standard Brownian motions with correlation $\rho$. Such SDDEs may be naturally decomposed into two subsystems: a linear stochastic differential equation (SDE) without delay, and a nonlinear SDDE.

Splitting methods work by solving each subsystem separately and composing the results over a single step. Out main theoretical result provides a bound on the mean-square error of a particular strategy for doing this, known as Lie-Trotter splitting. This bound implies that the method is mean-square strongly convergent with order $1/2$ when $\rho=0$, so that the noises are uncorrelated, but assurances of convergence are lost when $\rho\neq 0$. Indeed we develop an upper bound on the global mean-square error with a term depends linearly on the magnitude of the correlation, and is independent of the stepsize.

While our theoretical error bound is an estimate from above, we conduct numerical experiments that confirm the order of mean-square strong convergence of Lie-Trotter splitting in the $\rho=0$ case, and demonstrate a rapid fall-off to effectively zero as $|\rho|$ increases. Similar numerical results are observed for an alternative commonly used strategy known as Strang splitting.

ABSTRACT. In this talk we present new space-time convergence results for numerical methods applied to differential Riccati equations posed on infinite-dimensional Hilbert spaces. The analysed scheme uses Lie-splitting in its temporal discretization, while the assumptions on the spatial discretizations allow for finite element methods. The scheme achieves first-order convergence in time and second-order convergence in space, up to logarithmic factors, under mild regularity conditions. Convergence results for the spatially semidiscretized problem date back to the 1990s, but as far as the authors are aware results extending to a fully discretized scheme have not previously been available in the given setting.

A typical application is the linear-quadratic regulator problem over a finite time interval, when the system dynamics are described by a parabolic partial differential equation. In this case the optimal control law may be expressed in terms of the solution to a differential Riccati equation, with parameters that depend on the system being controlled. Spatial discretization of this Riccati equation yields matrix-valued Riccati equations, which then also need to be discretized in time. The established results ensure that the computed approximative solutions will converge to the underlying operator-valued solution as both the time and space discretizations are refined.