Society for Industrial and Applied Mathematics: SIAM Journal on Scientific Computing: Table of Contents
Table of Contents for SIAM Journal on Scientific Computing. List of articles from both the latest and ahead of print issues.
https://epubs.siam.org/loi/sjoce3?af=R
Society for Industrial and Applied Mathematics: SIAM Journal on Scientific Computing: Table of Contents
Society for Industrial and Applied Mathematics
enUS
SIAM Journal on Scientific Computing
SIAM Journal on Scientific Computing
https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg
https://epubs.siam.org/loi/sjoce3?af=R

Inverse Source Problem of the Biharmonic Equation from Multifrequency Phaseless Data
https://epubs.siam.org/doi/abs/10.1137/24M162889X?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page A2799A2818, October 2024. <br/> Abstract. This work deals with an inverse source problem for the biharmonic wave equation. A twostage numerical method is proposed to identify the unknown source from the multifrequency phaseless data. In the first stage, we introduce some artificial auxiliary point sources to the inverse source system and establish a phase retrieval formula. Theoretically, we point out that the phase can be uniquely determined and estimate the stability of this phase retrieval approach. Once the phase information is retrieved, the Fourier method is adopted to reconstruct the source function from the phased multifrequency data. The proposed method is easy to implement and there is no forward solver involved in the reconstruction. Numerical experiments are conducted to verify the performance of the proposed method.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page A2799A2818, October 2024. <br/> Abstract. This work deals with an inverse source problem for the biharmonic wave equation. A twostage numerical method is proposed to identify the unknown source from the multifrequency phaseless data. In the first stage, we introduce some artificial auxiliary point sources to the inverse source system and establish a phase retrieval formula. Theoretically, we point out that the phase can be uniquely determined and estimate the stability of this phase retrieval approach. Once the phase information is retrieved, the Fourier method is adopted to reconstruct the source function from the phased multifrequency data. The proposed method is easy to implement and there is no forward solver involved in the reconstruction. Numerical experiments are conducted to verify the performance of the proposed method. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Inverse Source Problem of the Biharmonic Equation from Multifrequency Phaseless Data
10.1137/24M162889X
SIAM Journal on Scientific Computing
20240903T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Yan Chang
Yukun Guo
Yue Zhao
Inverse Source Problem of the Biharmonic Equation from Multifrequency Phaseless Data
46
5
A2799
A2818
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/24M162889X
https://epubs.siam.org/doi/abs/10.1137/24M162889X?af=R
© 2024 Society for Industrial and Applied Mathematics

Generalizing Lloyd’s Algorithm for Graph Clustering
https://epubs.siam.org/doi/abs/10.1137/23M1556800?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page A2819A2847, October 2024. <br/> Abstract. Clustering is a commonplace problem in many areas of data science, with applications in biology and bioinformatics, understanding chemical structure, image segmentation, building recommender systems, and many more fields. While there are many different clustering variants (based on given distance or graph structure, probability distributions, or data density), we consider here the problem of clustering nodes in a graph, motivated by the problem of aggregating discrete degrees of freedom in multigrid and domain decomposition methods for solving sparse linear systems. Specifically, we consider the challenge of forming balanced clusters in the graph of a sparse matrix for use in algebraic multigrid, although the algorithm has general applicability. Based on an extension of the Bellman–Ford algorithm, we generalize Lloyd’s algorithm for partitioning subsets of [math] to balance the number of nodes in each cluster; this is accompanied by a rebalancing algorithm that reduces the overall energy in the system. The algorithm provides control over the number of clusters and leads to “well centered” partitions of the graph. Theoretical results are provided to establish linear complexity and numerical results in the context of algebraic multigrid highlight the benefits of improved clustering. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/lukeolson/paperlloyddata and in the supplementary materials (paperlloyddata.zip [88.1MB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page A2819A2847, October 2024. <br/> Abstract. Clustering is a commonplace problem in many areas of data science, with applications in biology and bioinformatics, understanding chemical structure, image segmentation, building recommender systems, and many more fields. While there are many different clustering variants (based on given distance or graph structure, probability distributions, or data density), we consider here the problem of clustering nodes in a graph, motivated by the problem of aggregating discrete degrees of freedom in multigrid and domain decomposition methods for solving sparse linear systems. Specifically, we consider the challenge of forming balanced clusters in the graph of a sparse matrix for use in algebraic multigrid, although the algorithm has general applicability. Based on an extension of the Bellman–Ford algorithm, we generalize Lloyd’s algorithm for partitioning subsets of [math] to balance the number of nodes in each cluster; this is accompanied by a rebalancing algorithm that reduces the overall energy in the system. The algorithm provides control over the number of clusters and leads to “well centered” partitions of the graph. Theoretical results are provided to establish linear complexity and numerical results in the context of algebraic multigrid highlight the benefits of improved clustering. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/lukeolson/paperlloyddata and in the supplementary materials (paperlloyddata.zip [88.1MB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Generalizing Lloyd’s Algorithm for Graph Clustering
10.1137/23M1556800
SIAM Journal on Scientific Computing
20240903T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Tareq Zaman
Nicolas Nytko
Ali Taghibakhshi
Scott MacLachlan
Luke Olson
Matthew West
Generalizing Lloyd’s Algorithm for Graph Clustering
46
5
A2819
A2847
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1556800
https://epubs.siam.org/doi/abs/10.1137/23M1556800?af=R
© 2024 Society for Industrial and Applied Mathematics

A New Provably Stable Weighted State Redistribution Algorithm
https://epubs.siam.org/doi/abs/10.1137/23M1597484?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page A2848A2873, October 2024. <br/> Abstract. We propose a practical finite volume method on cut cells using state redistribution. Our algorithm is provably monotone, total variation diminishing, and GKS (Gustafsson, Kreiss, Sundström) stable in many situations, and shuts off continuously as the cut cell size approaches a target value. Our analysis reveals why original state redistribution works so well: it results in a monotone scheme for most configurations, though at times subject to a slightly smaller CFL condition. Our analysis also explains why a premerging step is beneficial. We show computational experiments in two and three dimensions.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page A2848A2873, October 2024. <br/> Abstract. We propose a practical finite volume method on cut cells using state redistribution. Our algorithm is provably monotone, total variation diminishing, and GKS (Gustafsson, Kreiss, Sundström) stable in many situations, and shuts off continuously as the cut cell size approaches a target value. Our analysis reveals why original state redistribution works so well: it results in a monotone scheme for most configurations, though at times subject to a slightly smaller CFL condition. Our analysis also explains why a premerging step is beneficial. We show computational experiments in two and three dimensions. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A New Provably Stable Weighted State Redistribution Algorithm
10.1137/23M1597484
SIAM Journal on Scientific Computing
20240904T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Marsha Berger
Andrew Giuliani
A New Provably Stable Weighted State Redistribution Algorithm
46
5
A2848
A2873
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1597484
https://epubs.siam.org/doi/abs/10.1137/23M1597484?af=R
© 2024 Society for Industrial and Applied Mathematics

Stable Backward Differentiation Formula Time Discretization of BGNBased Parametric Finite Element Methods for Geometric Flows
https://epubs.siam.org/doi/abs/10.1137/23M1625597?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page A2874A2898, October 2024. <br/> Abstract. We propose a novel class of temporal highorder parametric finite element methods for solving a wide range of geometric flows of curves and surfaces. By incorporating the backward differentiation formula (BDF) for time discretization into the BGN formulation, originally proposed by Barrett, Garcke, and Nürnberg (J. Comput. Phys., 222 (2007), pp. 441–467), we successfully develop highorder BGN/BDF[math] schemes. The proposed BGN/BDF[math] schemes not only retain almost all the advantages of the classical firstorder BGN scheme such as computational efficiency and good mesh quality, but also exhibit the desired [math]thorder temporal accuracy in terms of shape metrics, ranging from secondorder to fourthorder accuracy. Furthermore, we validate the performance of our proposed BGN/BDF[math] schemes through extensive numerical examples, demonstrating their highorder temporal accuracy for various types of geometric flows while maintaining good mesh quality throughout the evolution.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page A2874A2898, October 2024. <br/> Abstract. We propose a novel class of temporal highorder parametric finite element methods for solving a wide range of geometric flows of curves and surfaces. By incorporating the backward differentiation formula (BDF) for time discretization into the BGN formulation, originally proposed by Barrett, Garcke, and Nürnberg (J. Comput. Phys., 222 (2007), pp. 441–467), we successfully develop highorder BGN/BDF[math] schemes. The proposed BGN/BDF[math] schemes not only retain almost all the advantages of the classical firstorder BGN scheme such as computational efficiency and good mesh quality, but also exhibit the desired [math]thorder temporal accuracy in terms of shape metrics, ranging from secondorder to fourthorder accuracy. Furthermore, we validate the performance of our proposed BGN/BDF[math] schemes through extensive numerical examples, demonstrating their highorder temporal accuracy for various types of geometric flows while maintaining good mesh quality throughout the evolution. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Stable Backward Differentiation Formula Time Discretization of BGNBased Parametric Finite Element Methods for Geometric Flows
10.1137/23M1625597
SIAM Journal on Scientific Computing
20240906T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Wei Jiang
Chunmei Su
Ganghui Zhang
Stable Backward Differentiation Formula Time Discretization of BGNBased Parametric Finite Element Methods for Geometric Flows
46
5
A2874
A2898
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1625597
https://epubs.siam.org/doi/abs/10.1137/23M1625597?af=R
© 2024 Society for Industrial and Applied Mathematics

BoundPreserving Framework for CentralUpwind Schemes for General Hyperbolic Conservation Laws
https://epubs.siam.org/doi/abs/10.1137/23M1628024?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page A2899A2924, October 2024. <br/> Abstract. Centralupwind (CU) schemes are Riemannproblemsolverfree finitevolume methods widely applied to a variety of hyperbolic systems of PDEs. Exact solutions of these systems typically satisfy certain bounds, and it is highly desirable and even crucial for the numerical schemes to preserve these bounds. In this paper, we develop and analyze boundpreserving (BP) CU schemes for general hyperbolic systems of conservation laws. Unlike many other Godunovtype methods, CU schemes cannot, in general, be recast as convex combinations of firstorder BP schemes. Consequently, standard BP analysis techniques are invalidated. We address these challenges by establishing a novel framework for analyzing the BP property of CU schemes. To this end, we discover that the CU schemes can be decomposed as a convex combination of several intermediate solution states. Thanks to this key finding, the goal of designing BPCU schemes is simplified to the enforcement of four more accessible BP conditions, each of which can be achieved with the help of a minor modification of the CU schemes. We employ the proposed approach to construct provably BPCU schemes for the Euler equations of gas dynamics. The robustness and effectiveness of the BPCU schemes are validated by several demanding numerical examples, including highspeed jet problems, flow past a forwardfacing step, and a shock diffraction problem.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page A2899A2924, October 2024. <br/> Abstract. Centralupwind (CU) schemes are Riemannproblemsolverfree finitevolume methods widely applied to a variety of hyperbolic systems of PDEs. Exact solutions of these systems typically satisfy certain bounds, and it is highly desirable and even crucial for the numerical schemes to preserve these bounds. In this paper, we develop and analyze boundpreserving (BP) CU schemes for general hyperbolic systems of conservation laws. Unlike many other Godunovtype methods, CU schemes cannot, in general, be recast as convex combinations of firstorder BP schemes. Consequently, standard BP analysis techniques are invalidated. We address these challenges by establishing a novel framework for analyzing the BP property of CU schemes. To this end, we discover that the CU schemes can be decomposed as a convex combination of several intermediate solution states. Thanks to this key finding, the goal of designing BPCU schemes is simplified to the enforcement of four more accessible BP conditions, each of which can be achieved with the help of a minor modification of the CU schemes. We employ the proposed approach to construct provably BPCU schemes for the Euler equations of gas dynamics. The robustness and effectiveness of the BPCU schemes are validated by several demanding numerical examples, including highspeed jet problems, flow past a forwardfacing step, and a shock diffraction problem. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
BoundPreserving Framework for CentralUpwind Schemes for General Hyperbolic Conservation Laws
10.1137/23M1628024
SIAM Journal on Scientific Computing
20240910T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Shumo Cui
Alexander Kurganov
Kailiang Wu
BoundPreserving Framework for CentralUpwind Schemes for General Hyperbolic Conservation Laws
46
5
A2899
A2924
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1628024
https://epubs.siam.org/doi/abs/10.1137/23M1628024?af=R
© 2024 Society for Industrial and Applied Mathematics

The SparseGridBased Adaptive Spectral Koopman Method
https://epubs.siam.org/doi/abs/10.1137/23M1578292?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page A2925A2950, October 2024. <br/> Abstract. The adaptive spectral Koopman (ASK) method was introduced to numerically solve autonomous dynamical systems that laid the foundation for numerous applications across different fields in science and engineering. Although ASK achieves high accuracy, it is computationally more expensive for multidimensional systems compared with conventional time integration schemes like Runge–Kutta. In this work, we combine the sparse grid and ASK to accelerate the computation for multidimensional systems. This sparsegridbased ASK (SASK) method uses the Smolyak structure to construct multidimensional collocation points as well as associated polynomials that are used to approximate eigenfunctions of the Koopman operator of the system. In this way, the number of collocation points is reduced compared with using the tensor product rule. We demonstrate that SASK can be used to solve ordinary differential equations (ODEs) and partial differential equations (PDEs) based on their semidiscrete forms. Numerical experiments are illustrated to compare the performance of SASK and stateoftheart ODE solvers.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page A2925A2950, October 2024. <br/> Abstract. The adaptive spectral Koopman (ASK) method was introduced to numerically solve autonomous dynamical systems that laid the foundation for numerous applications across different fields in science and engineering. Although ASK achieves high accuracy, it is computationally more expensive for multidimensional systems compared with conventional time integration schemes like Runge–Kutta. In this work, we combine the sparse grid and ASK to accelerate the computation for multidimensional systems. This sparsegridbased ASK (SASK) method uses the Smolyak structure to construct multidimensional collocation points as well as associated polynomials that are used to approximate eigenfunctions of the Koopman operator of the system. In this way, the number of collocation points is reduced compared with using the tensor product rule. We demonstrate that SASK can be used to solve ordinary differential equations (ODEs) and partial differential equations (PDEs) based on their semidiscrete forms. Numerical experiments are illustrated to compare the performance of SASK and stateoftheart ODE solvers. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
The SparseGridBased Adaptive Spectral Koopman Method
10.1137/23M1578292
SIAM Journal on Scientific Computing
20240911T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Bian Li
Yue Yu
Xiu Yang
The SparseGridBased Adaptive Spectral Koopman Method
46
5
A2925
A2950
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1578292
https://epubs.siam.org/doi/abs/10.1137/23M1578292?af=R
© 2024 Society for Industrial and Applied Mathematics

Robust Iterative Method for Symmetric Quantum Signal Processing in All Parameter Regimes
https://epubs.siam.org/doi/abs/10.1137/23M1598192?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page A2951A2971, October 2024. <br/> Abstract. This paper addresses the problem of solving nonlinear systems in the context of symmetric quantum signal processing (QSP), a powerful technique for implementing matrix functions on quantum computers. Symmetric QSP focuses on representing target polynomials as products of matrices in SU(2) that possess symmetry properties. We present a novel Newton’s method tailored for efficiently solving the nonlinear system involved in determining the phase factors within the symmetric QSP framework. Our method demonstrates rapid and robust convergence in all parameter regimes, including the challenging scenario with illconditioned Jacobian matrices, using standard double precision arithmetic operations. For instance, solving symmetric QSP for a highly oscillatory target function [math] (polynomial degree [math]) takes 6 iterations to converge to machine precision when [math], and the number of iterations only increases to 18 iterations when [math] with a highly illconditioned Jacobian matrix. Leveraging the matrix product state structure of symmetric QSP, the computation of the Jacobian matrix incurs a computational cost comparable to a single function evaluation. Moreover, we introduce a reformulation of symmetric QSP using realnumber arithmetics, further enhancing the method’s efficiency. Extensive numerical tests validate the effectiveness and robustness of our approach, which has been implemented in the QSPPACK software package.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page A2951A2971, October 2024. <br/> Abstract. This paper addresses the problem of solving nonlinear systems in the context of symmetric quantum signal processing (QSP), a powerful technique for implementing matrix functions on quantum computers. Symmetric QSP focuses on representing target polynomials as products of matrices in SU(2) that possess symmetry properties. We present a novel Newton’s method tailored for efficiently solving the nonlinear system involved in determining the phase factors within the symmetric QSP framework. Our method demonstrates rapid and robust convergence in all parameter regimes, including the challenging scenario with illconditioned Jacobian matrices, using standard double precision arithmetic operations. For instance, solving symmetric QSP for a highly oscillatory target function [math] (polynomial degree [math]) takes 6 iterations to converge to machine precision when [math], and the number of iterations only increases to 18 iterations when [math] with a highly illconditioned Jacobian matrix. Leveraging the matrix product state structure of symmetric QSP, the computation of the Jacobian matrix incurs a computational cost comparable to a single function evaluation. Moreover, we introduce a reformulation of symmetric QSP using realnumber arithmetics, further enhancing the method’s efficiency. Extensive numerical tests validate the effectiveness and robustness of our approach, which has been implemented in the QSPPACK software package. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Robust Iterative Method for Symmetric Quantum Signal Processing in All Parameter Regimes
10.1137/23M1598192
SIAM Journal on Scientific Computing
20240913T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Yulong Dong
Lin Lin
Hongkang Ni
Jiasu Wang
Robust Iterative Method for Symmetric Quantum Signal Processing in All Parameter Regimes
46
5
A2951
A2971
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1598192
https://epubs.siam.org/doi/abs/10.1137/23M1598192?af=R
© 2024 Society for Industrial and Applied Mathematics

Finite Element Approximation for the Delayed Generalized Burgers–Huxley Equation with Weakly Singular Kernel: Part II Nonconforming and DG Approximation
https://epubs.siam.org/doi/abs/10.1137/23M1612196?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page A2972A2998, October 2024. <br/> Abstract. In this paper, the numerical approximation of the generalized Burgers–Huxley equation (GBHE) with weakly singular kernels using nonconforming methods will be presented. Specifically, we discuss two new formulations. The first formulation is based on the nonconforming finite element method. The other formulation is based on discontinuous Galerkin finite element methods. The wellposedness results for both formulations are proved. Then, a priori error estimates for both the semidiscrete and fully discrete schemes are derived. Specific numerical examples, including some applications for the GBHE with a weakly singular model, are discussed to validate the theoretical results.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page A2972A2998, October 2024. <br/> Abstract. In this paper, the numerical approximation of the generalized Burgers–Huxley equation (GBHE) with weakly singular kernels using nonconforming methods will be presented. Specifically, we discuss two new formulations. The first formulation is based on the nonconforming finite element method. The other formulation is based on discontinuous Galerkin finite element methods. The wellposedness results for both formulations are proved. Then, a priori error estimates for both the semidiscrete and fully discrete schemes are derived. Specific numerical examples, including some applications for the GBHE with a weakly singular model, are discussed to validate the theoretical results. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Finite Element Approximation for the Delayed Generalized Burgers–Huxley Equation with Weakly Singular Kernel: Part II Nonconforming and DG Approximation
10.1137/23M1612196
SIAM Journal on Scientific Computing
20240919T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Sumit Mahahjan
Arbaz Khan
Finite Element Approximation for the Delayed Generalized Burgers–Huxley Equation with Weakly Singular Kernel: Part II Nonconforming and DG Approximation
46
5
A2972
A2998
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1612196
https://epubs.siam.org/doi/abs/10.1137/23M1612196?af=R
© 2024 Society for Industrial and Applied Mathematics

Conditional QuasiMonte Carlo with Constrained Active Subspaces
https://epubs.siam.org/doi/abs/10.1137/23M1548918?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page A2999A3021, October 2024. <br/> Abstract. Conditional Monte Carlo or preintegration is a powerful tool for reducing variance and improving the regularity of integrands when using Monte Carlo and quasiMonte Carlo (QMC) methods. To select the variable to preintegrate, one must consider both the variable’s importance and the tractability of the conditional expectation. For integrals over a Gaussian distribution, any linear combination of variables can potentially be preintegrated. Liu and Owen [SIAM J. Numer. Anal., 61 (2023), pp. 495–514] propose to select the linear combination based on an active subspace decomposition of the integrand. However, preintegrating the selected direction might be intractable. In this work, we address this issue by finding the active subspace subject to constraints such that preintegration can be easily carried out. The proposed algorithm also provides a computationally efficient alternative to dimension reduction for preintegrated functions. The method is applied to examples from computational finance, density estimation, and computational chemistry and is shown to achieve smaller errors than previous methods.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page A2999A3021, October 2024. <br/> Abstract. Conditional Monte Carlo or preintegration is a powerful tool for reducing variance and improving the regularity of integrands when using Monte Carlo and quasiMonte Carlo (QMC) methods. To select the variable to preintegrate, one must consider both the variable’s importance and the tractability of the conditional expectation. For integrals over a Gaussian distribution, any linear combination of variables can potentially be preintegrated. Liu and Owen [SIAM J. Numer. Anal., 61 (2023), pp. 495–514] propose to select the linear combination based on an active subspace decomposition of the integrand. However, preintegrating the selected direction might be intractable. In this work, we address this issue by finding the active subspace subject to constraints such that preintegration can be easily carried out. The proposed algorithm also provides a computationally efficient alternative to dimension reduction for preintegrated functions. The method is applied to examples from computational finance, density estimation, and computational chemistry and is shown to achieve smaller errors than previous methods. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Conditional QuasiMonte Carlo with Constrained Active Subspaces
10.1137/23M1548918
SIAM Journal on Scientific Computing
20240923T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Sifan Liu
Conditional QuasiMonte Carlo with Constrained Active Subspaces
46
5
A2999
A3021
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1548918
https://epubs.siam.org/doi/abs/10.1137/23M1548918?af=R
© 2024 Society for Industrial and Applied Mathematics

Randomized Sketching of Nonlinear Eigenvalue Problems
https://epubs.siam.org/doi/abs/10.1137/22M153656X?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page A3022A3043, October 2024. <br/> Abstract. Rational approximation is a powerful tool to obtain accurate surrogates for nonlinear functions that are easy to evaluate and linearize. The interpolatory adaptive Antoulas–Anderson (AAA) method is one approach to construct such approximants numerically. For largescale vector and matrixvalued functions, however, the direct application of the setvalued variant of AAA becomes inefficient. We propose and analyze a new sketching approach for such functions called sketchAAA that, with high probability, leads to much better approximants than previously suggested approaches while retaining efficiency. The sketching approach works in a blackbox fashion where only evaluations of the nonlinear function at sampling points are needed. Numerical tests with nonlinear eigenvalue problems illustrate the efficacy of our approach, with speedups over 200 for sampling largescale blackbox functions without sacrificing accuracy. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://zenodo.org/doi/10.5281/zenodo.11634528.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page A3022A3043, October 2024. <br/> Abstract. Rational approximation is a powerful tool to obtain accurate surrogates for nonlinear functions that are easy to evaluate and linearize. The interpolatory adaptive Antoulas–Anderson (AAA) method is one approach to construct such approximants numerically. For largescale vector and matrixvalued functions, however, the direct application of the setvalued variant of AAA becomes inefficient. We propose and analyze a new sketching approach for such functions called sketchAAA that, with high probability, leads to much better approximants than previously suggested approaches while retaining efficiency. The sketching approach works in a blackbox fashion where only evaluations of the nonlinear function at sampling points are needed. Numerical tests with nonlinear eigenvalue problems illustrate the efficacy of our approach, with speedups over 200 for sampling largescale blackbox functions without sacrificing accuracy. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://zenodo.org/doi/10.5281/zenodo.11634528. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Randomized Sketching of Nonlinear Eigenvalue Problems
10.1137/22M153656X
SIAM Journal on Scientific Computing
20240924T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Stefan Güttel
Daniel Kressner
Bart Vandereycken
Randomized Sketching of Nonlinear Eigenvalue Problems
46
5
A3022
A3043
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/22M153656X
https://epubs.siam.org/doi/abs/10.1137/22M153656X?af=R
© 2024 Society for Industrial and Applied Mathematics

Shape Optimization by Constrained FirstOrder System Least Mean Approximation
https://epubs.siam.org/doi/abs/10.1137/23M1605570?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page A3044A3066, October 2024. <br/> Abstract. In this work, the problem of shape optimization, subject to PDE constraints, is reformulated as an [math] best approximation problem under divergence constraints to the shape tensor introduced by Laurain and Sturm [ESAIM Math. Model. Numer. Anal., 50 (2016), pp. 1241–1267]. More precisely, the main result of this paper states that the [math] distance of the above approximation problem is equal to the dual norm of the shape derivative considered as a functional on [math] (where [math]). This implies that for any given shape, one can evaluate its distance from being stationary with respect to the shape derivative by simply solving the associated [math]type least mean approximation problem. Moreover, the Lagrange multiplier for the divergence constraint turns out to be the shape deformation of steepest descent. This provides, as an alternative to the approach by Deckelnick, Herbert, and Hinze, a way to compute shape gradients in [math] for [math]. The discretization of the least mean approximation problem is done with (lowestorder) matrixvalued Raviart–Thomas finite element spaces leading to piecewise constant approximations of the shape deformation acting as a Lagrange multiplier. Admissible deformations in [math] to be used in a shape gradient iteration are reconstructed locally. Our computational results confirm that the [math] distance of the best approximation does indeed measure the distance of the considered shape to optimality. Also confirmed by our computational tests is the observation of Deckelnick, Herbert, and Hinze [ESAIM Control Optim. Calc. Var., 28 (2022), 2] that choosing [math] (much) larger than 2 (which means that [math] must be close to 1 in our best approximation problem) decreases the chance of encountering mesh degeneracy during the shape gradient iteration. Reproducibility of computational results.This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available in the supplemental material (2DshapeclmSISC160557R2.zip [34.6KB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page A3044A3066, October 2024. <br/> Abstract. In this work, the problem of shape optimization, subject to PDE constraints, is reformulated as an [math] best approximation problem under divergence constraints to the shape tensor introduced by Laurain and Sturm [ESAIM Math. Model. Numer. Anal., 50 (2016), pp. 1241–1267]. More precisely, the main result of this paper states that the [math] distance of the above approximation problem is equal to the dual norm of the shape derivative considered as a functional on [math] (where [math]). This implies that for any given shape, one can evaluate its distance from being stationary with respect to the shape derivative by simply solving the associated [math]type least mean approximation problem. Moreover, the Lagrange multiplier for the divergence constraint turns out to be the shape deformation of steepest descent. This provides, as an alternative to the approach by Deckelnick, Herbert, and Hinze, a way to compute shape gradients in [math] for [math]. The discretization of the least mean approximation problem is done with (lowestorder) matrixvalued Raviart–Thomas finite element spaces leading to piecewise constant approximations of the shape deformation acting as a Lagrange multiplier. Admissible deformations in [math] to be used in a shape gradient iteration are reconstructed locally. Our computational results confirm that the [math] distance of the best approximation does indeed measure the distance of the considered shape to optimality. Also confirmed by our computational tests is the observation of Deckelnick, Herbert, and Hinze [ESAIM Control Optim. Calc. Var., 28 (2022), 2] that choosing [math] (much) larger than 2 (which means that [math] must be close to 1 in our best approximation problem) decreases the chance of encountering mesh degeneracy during the shape gradient iteration. Reproducibility of computational results.This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available in the supplemental material (2DshapeclmSISC160557R2.zip [34.6KB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Shape Optimization by Constrained FirstOrder System Least Mean Approximation
10.1137/23M1605570
SIAM Journal on Scientific Computing
20241001T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Gerhard Starke
Shape Optimization by Constrained FirstOrder System Least Mean Approximation
46
5
A3044
A3066
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1605570
https://epubs.siam.org/doi/abs/10.1137/23M1605570?af=R
© 2024 Society for Industrial and Applied Mathematics

A Scalable TwoLevel Domain Decomposition Eigensolver for Periodic Schrödinger Eigenstates in Anisotropically Expanding Domains
https://epubs.siam.org/doi/abs/10.1137/23M161848X?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page A3067A3093, October 2024. <br/> Abstract. Accelerating iterative eigenvalue algorithms is often achieved by employing a spectral shifting strategy. Unfortunately, improved shifting typically leads to a smaller eigenvalue for the resulting shifted operator, which in turn results in a high condition number of the underlying solution matrix, posing a major challenge for iterative linear solvers. This paper introduces a twolevel domain decomposition preconditioner that addresses this issue for the linear Schrödinger eigenvalue problem, even in the presence of a vanishing eigenvalue gap in nonuniform, expanding domains. Since the quasioptimal shift, which is already available as the solution to a spectral cell problem, is required for the eigenvalue solver, it is logical to also use its associated eigenfunction as a generator to construct a coarse space. We analyze the resulting twolevel additive Schwarz preconditioner and obtain a condition number bound that is independent of the domain’s anisotropy, despite the need for only one basis function per subdomain for the coarse solver. Several numerical examples are presented to illustrate its flexibility and efficiency. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://doi.org/10.5281/zenodo.11072863.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page A3067A3093, October 2024. <br/> Abstract. Accelerating iterative eigenvalue algorithms is often achieved by employing a spectral shifting strategy. Unfortunately, improved shifting typically leads to a smaller eigenvalue for the resulting shifted operator, which in turn results in a high condition number of the underlying solution matrix, posing a major challenge for iterative linear solvers. This paper introduces a twolevel domain decomposition preconditioner that addresses this issue for the linear Schrödinger eigenvalue problem, even in the presence of a vanishing eigenvalue gap in nonuniform, expanding domains. Since the quasioptimal shift, which is already available as the solution to a spectral cell problem, is required for the eigenvalue solver, it is logical to also use its associated eigenfunction as a generator to construct a coarse space. We analyze the resulting twolevel additive Schwarz preconditioner and obtain a condition number bound that is independent of the domain’s anisotropy, despite the need for only one basis function per subdomain for the coarse solver. Several numerical examples are presented to illustrate its flexibility and efficiency. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://doi.org/10.5281/zenodo.11072863. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A Scalable TwoLevel Domain Decomposition Eigensolver for Periodic Schrödinger Eigenstates in Anisotropically Expanding Domains
10.1137/23M161848X
SIAM Journal on Scientific Computing
20241001T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Lambert Theisen
Benjamin Stamm
A Scalable TwoLevel Domain Decomposition Eigensolver for Periodic Schrödinger Eigenstates in Anisotropically Expanding Domains
46
5
A3067
A3093
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M161848X
https://epubs.siam.org/doi/abs/10.1137/23M161848X?af=R
© 2024 Society for Industrial and Applied Mathematics

Bounds on Nonlinear Errors for Variance Computation with Stochastic Rounding
https://epubs.siam.org/doi/abs/10.1137/23M1563001?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page B579B599, October 2024. <br/> Abstract. The main objective of this work is to investigate nonlinear errors and pairwise summation using stochastic rounding (SR) in variance computation algorithms. We estimate the forward error of computations under SR through two methods: the first is based on a bound of the variance and the Bienaymé–Chebyshev inequality, while the second is based on martingales and the Azuma–Hoeffding inequality. The study shows that for pairwise summation, using SR results in a probabilistic bound of the forward error proportional to [math] rather than the deterministic bound in [math] when using the default rounding mode. We examine two algorithms that compute the variance, one called “textbook” and the other “twopass,” which both exhibit nonlinear errors. Using the two methods mentioned above, we show that the forward errors of these algorithms have probabilistic bounds under SR in [math] instead of [math] for the deterministic bounds. We show that this advantage holds using pairwise summation for both textbook and twopass, with probabilistic bounds of the forward error proportional to [math]. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow the reader to reproduce the results in this paper are available at https://github.com/verificarlo/srnonlinearbounds and in the supplementary material (srnonlinearboundsmain.zip [8.62KB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page B579B599, October 2024. <br/> Abstract. The main objective of this work is to investigate nonlinear errors and pairwise summation using stochastic rounding (SR) in variance computation algorithms. We estimate the forward error of computations under SR through two methods: the first is based on a bound of the variance and the Bienaymé–Chebyshev inequality, while the second is based on martingales and the Azuma–Hoeffding inequality. The study shows that for pairwise summation, using SR results in a probabilistic bound of the forward error proportional to [math] rather than the deterministic bound in [math] when using the default rounding mode. We examine two algorithms that compute the variance, one called “textbook” and the other “twopass,” which both exhibit nonlinear errors. Using the two methods mentioned above, we show that the forward errors of these algorithms have probabilistic bounds under SR in [math] instead of [math] for the deterministic bounds. We show that this advantage holds using pairwise summation for both textbook and twopass, with probabilistic bounds of the forward error proportional to [math]. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow the reader to reproduce the results in this paper are available at https://github.com/verificarlo/srnonlinearbounds and in the supplementary material (srnonlinearboundsmain.zip [8.62KB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Bounds on Nonlinear Errors for Variance Computation with Stochastic Rounding
10.1137/23M1563001
SIAM Journal on Scientific Computing
20240903T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
EM. El Arar
D. Sohier
P. de Oliveira Castro
E. Petit
Bounds on Nonlinear Errors for Variance Computation with Stochastic Rounding
46
5
B579
B599
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1563001
https://epubs.siam.org/doi/abs/10.1137/23M1563001?af=R
© 2024 Society for Industrial and Applied Mathematics

Interpolating Parametrized Quantum Circuits Using Blackbox Queries
https://epubs.siam.org/doi/abs/10.1137/23M1609543?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page B600B620, October 2024. <br/> Abstract. This article focuses on developing classical surrogates for parametrized quantum circuits using interpolation via (trigonometric) polynomials. We develop two algorithms for the construction of such surrogates and prove performance guarantees. The constructions are based on circuit evaluations which are blackbox in the sense that no structural specifics of the circuits are exploited. While acknowledging the limitations of the blackbox approach compared to whitebox evaluations, which exploit specific circuit properties, we demonstrate scenarios in which the blackbox approach might prove beneficial. Sample applications include but are not restricted to the approximation of variational quantum eigensolvers and the alleviaton of the barren plateau problem.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page B600B620, October 2024. <br/> Abstract. This article focuses on developing classical surrogates for parametrized quantum circuits using interpolation via (trigonometric) polynomials. We develop two algorithms for the construction of such surrogates and prove performance guarantees. The constructions are based on circuit evaluations which are blackbox in the sense that no structural specifics of the circuits are exploited. While acknowledging the limitations of the blackbox approach compared to whitebox evaluations, which exploit specific circuit properties, we demonstrate scenarios in which the blackbox approach might prove beneficial. Sample applications include but are not restricted to the approximation of variational quantum eigensolvers and the alleviaton of the barren plateau problem. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Interpolating Parametrized Quantum Circuits Using Blackbox Queries
10.1137/23M1609543
SIAM Journal on Scientific Computing
20240913T07:00:00Z
© 2024 Bundesdruckerei GmbH
Lars Simon
Holger Eble
HagenHenrik Kowalski
Manuel Radons
Interpolating Parametrized Quantum Circuits Using Blackbox Queries
46
5
B600
B620
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1609543
https://epubs.siam.org/doi/abs/10.1137/23M1609543?af=R
© 2024 Bundesdruckerei GmbH

A Neural Network Approach for Stochastic Optimal Control
https://epubs.siam.org/doi/abs/10.1137/23M155832X?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page C535C556, October 2024. <br/> Abstract. We present a neural network approach for approximating the value function of highdimensional stochastic control problems. Our training process simultaneously updates our value function estimate and identifies the part of the state space likely to be visited by optimal trajectories. Our approach leverages insights from optimal control theory and the fundamental relation between semilinear parabolic partial differential equations and forwardbackward stochastic differential equations. To focus the sampling on relevant states during neural network training, we use the stochastic Pontryagin maximum principle (PMP) to obtain the optimal controls for the current value function estimate. By design, our approach coincides with the method of characteristics for the nonviscous Hamilton–Jacobi–Bellman equation arising in deterministic control problems. Our training loss consists of a weighted sum of the objective functional of the control problem and penalty terms that enforce the HJB equations along the sampled trajectories. Importantly, training is unsupervised in that it does not require solutions of the control problem. Our numerical experiments highlight our scheme’s ability to identify the relevant parts of the state space and produce meaningful value estimates. Using a twodimensional model problem, we demonstrate the importance of the stochastic PMP to inform the sampling and compare it to a finite element approach. With a nonlinear control affine quadcopter example, we illustrate that our approach can handle complicated dynamics. For a 100dimensional benchmark problem, we demonstrate that our approach improves accuracy and timetosolution, and, via a modification, we show the wider applicability of our scheme. Reproducibility of computational results.This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/EmoryMLIP/NeuralSOC and in the supplementary material (NeuralSOCmain.zip [ 29.9MB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page C535C556, October 2024. <br/> Abstract. We present a neural network approach for approximating the value function of highdimensional stochastic control problems. Our training process simultaneously updates our value function estimate and identifies the part of the state space likely to be visited by optimal trajectories. Our approach leverages insights from optimal control theory and the fundamental relation between semilinear parabolic partial differential equations and forwardbackward stochastic differential equations. To focus the sampling on relevant states during neural network training, we use the stochastic Pontryagin maximum principle (PMP) to obtain the optimal controls for the current value function estimate. By design, our approach coincides with the method of characteristics for the nonviscous Hamilton–Jacobi–Bellman equation arising in deterministic control problems. Our training loss consists of a weighted sum of the objective functional of the control problem and penalty terms that enforce the HJB equations along the sampled trajectories. Importantly, training is unsupervised in that it does not require solutions of the control problem. Our numerical experiments highlight our scheme’s ability to identify the relevant parts of the state space and produce meaningful value estimates. Using a twodimensional model problem, we demonstrate the importance of the stochastic PMP to inform the sampling and compare it to a finite element approach. With a nonlinear control affine quadcopter example, we illustrate that our approach can handle complicated dynamics. For a 100dimensional benchmark problem, we demonstrate that our approach improves accuracy and timetosolution, and, via a modification, we show the wider applicability of our scheme. Reproducibility of computational results.This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/EmoryMLIP/NeuralSOC and in the supplementary material (NeuralSOCmain.zip [ 29.9MB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A Neural Network Approach for Stochastic Optimal Control
10.1137/23M155832X
SIAM Journal on Scientific Computing
20240903T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Xingjian Li
Deepanshu Verma
Lars Ruthotto
A Neural Network Approach for Stochastic Optimal Control
46
5
C535
C556
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M155832X
https://epubs.siam.org/doi/abs/10.1137/23M155832X?af=R
© 2024 Society for Industrial and Applied Mathematics

A Domain Decomposition–Based CNNDNN Architecture for Model Parallel Training Applied to Image Recognition Problems
https://epubs.siam.org/doi/abs/10.1137/23M1562202?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page C557C582, October 2024. <br/> Abstract. Deep neural networks (DNNs) and, in particular, convolutional neural networks (CNNs) have brought significant advances in a wide range of modern computer application problems. However, the increasing availability of large numbers of datasets and the increasing available computational power of modern computers have led to steady growth in the complexity and size of DNN and CNN models, respectively, and thus, to longer training times. Hence, various methods and attempts have been developed to accelerate and parallelize the training of complex network architectures. In this work, a novel CNNDNN architecture is proposed that naturally supports a model parallel training strategy and that is loosely inspired by twolevel domain decomposition methods (DDMs). First, local CNN models, that is, subnetworks, are defined that operate on overlapping or nonoverlapping parts of the input data, for example, subimages. The subnetworks can be trained completely in parallel and independently of each other. Each subnetwork then outputs a local decision for the given machine learning problem which is exclusively based on the respective local input data. Subsequently, in a second step, an additional DNN model is trained which evaluates the local decisions of the local subnetworks and generates a final, global decision. With respect to the analogy to DDMs, the DNN models can be loosely interpreted as a coarse problem and hence, the new approach can be interpreted as a twolevel domain decomposition. In this paper, we apply the proposed architecture to image classification problems using CNNs. Experimental results for different twodimensional image classification problems are provided, as well as a face recognition problem and a classification problem for threedimensional computed tomography (CT) scans. Therefore, classical Residual Network (ResNet) and VGG architectures are considered. More modern architectures, such as, e.g., MobileNet2, are left for future work. The results show that the proposed approach can significantly accelerate the required training time compared to the global model and, additionally, can also help to improve the accuracy of the underlying classification problem.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page C557C582, October 2024. <br/> Abstract. Deep neural networks (DNNs) and, in particular, convolutional neural networks (CNNs) have brought significant advances in a wide range of modern computer application problems. However, the increasing availability of large numbers of datasets and the increasing available computational power of modern computers have led to steady growth in the complexity and size of DNN and CNN models, respectively, and thus, to longer training times. Hence, various methods and attempts have been developed to accelerate and parallelize the training of complex network architectures. In this work, a novel CNNDNN architecture is proposed that naturally supports a model parallel training strategy and that is loosely inspired by twolevel domain decomposition methods (DDMs). First, local CNN models, that is, subnetworks, are defined that operate on overlapping or nonoverlapping parts of the input data, for example, subimages. The subnetworks can be trained completely in parallel and independently of each other. Each subnetwork then outputs a local decision for the given machine learning problem which is exclusively based on the respective local input data. Subsequently, in a second step, an additional DNN model is trained which evaluates the local decisions of the local subnetworks and generates a final, global decision. With respect to the analogy to DDMs, the DNN models can be loosely interpreted as a coarse problem and hence, the new approach can be interpreted as a twolevel domain decomposition. In this paper, we apply the proposed architecture to image classification problems using CNNs. Experimental results for different twodimensional image classification problems are provided, as well as a face recognition problem and a classification problem for threedimensional computed tomography (CT) scans. Therefore, classical Residual Network (ResNet) and VGG architectures are considered. More modern architectures, such as, e.g., MobileNet2, are left for future work. The results show that the proposed approach can significantly accelerate the required training time compared to the global model and, additionally, can also help to improve the accuracy of the underlying classification problem. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A Domain Decomposition–Based CNNDNN Architecture for Model Parallel Training Applied to Image Recognition Problems
10.1137/23M1562202
SIAM Journal on Scientific Computing
20240904T07:00:00Z
© 2024 Axel Klawonn, Martin Lanser, and Janine Weber
Axel Klawonn
Martin Lanser
Janine Weber
A Domain Decomposition–Based CNNDNN Architecture for Model Parallel Training Applied to Image Recognition Problems
46
5
C557
C582
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1562202
https://epubs.siam.org/doi/abs/10.1137/23M1562202?af=R
© 2024 Axel Klawonn, Martin Lanser, and Janine Weber

Nonlinear Embeddings for Conserving Hamiltonians and Other Quantities with Neural Galerkin Schemes
https://epubs.siam.org/doi/abs/10.1137/23M1607799?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page C583C607, October 2024. <br/> Abstract. This work focuses on the conservation of quantities such as Hamiltonians, mass, and momentum when solution fields of partial differential equations are approximated with nonlinear parametrizations such as deep networks. The proposed approach builds on Neural Galerkin schemes that are based on the Dirac–Frenkel variational principle to train nonlinear parametrizations sequentially in time. We first show that only adding constraints that aim to conserve quantities in continuous time can be insufficient because the nonlinear dependence on the parameters implies that even quantities that are linear in the solution fields become nonlinear in the parameters and thus are challenging to discretize in time. Instead, we propose Neural Galerkin schemes that compute at each time step an explicit embedding onto the manifold of nonlinearly parametrized solution fields to guarantee conservation of quantities. The embeddings can be combined with standard explicit and implicit time integration schemes. Numerical experiments demonstrate that the proposed approach conserves quantities up to machine precision.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page C583C607, October 2024. <br/> Abstract. This work focuses on the conservation of quantities such as Hamiltonians, mass, and momentum when solution fields of partial differential equations are approximated with nonlinear parametrizations such as deep networks. The proposed approach builds on Neural Galerkin schemes that are based on the Dirac–Frenkel variational principle to train nonlinear parametrizations sequentially in time. We first show that only adding constraints that aim to conserve quantities in continuous time can be insufficient because the nonlinear dependence on the parameters implies that even quantities that are linear in the solution fields become nonlinear in the parameters and thus are challenging to discretize in time. Instead, we propose Neural Galerkin schemes that compute at each time step an explicit embedding onto the manifold of nonlinearly parametrized solution fields to guarantee conservation of quantities. The embeddings can be combined with standard explicit and implicit time integration schemes. Numerical experiments demonstrate that the proposed approach conserves quantities up to machine precision. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Nonlinear Embeddings for Conserving Hamiltonians and Other Quantities with Neural Galerkin Schemes
10.1137/23M1607799
SIAM Journal on Scientific Computing
20240925T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Paul Schwerdtner
Philipp Schulze
Jules Berman
Benjamin Peherstorfer
Nonlinear Embeddings for Conserving Hamiltonians and Other Quantities with Neural Galerkin Schemes
46
5
C583
C607
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1607799
https://epubs.siam.org/doi/abs/10.1137/23M1607799?af=R
© 2024 Society for Industrial and Applied Mathematics

Learning InBetween Imagery Dynamics via Physical Latent Spaces
https://epubs.siam.org/doi/abs/10.1137/23M1609440?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/5">Volume 46, Issue 5</a>, Page C608C632, October 2024. <br/> Abstract. We present a framework designed to learn the underlying dynamics between two images observed at consecutive time steps. The complex nature of image data and the lack of temporal information pose significant challenges in capturing the unique evolving patterns. Our proposed method focuses on estimating the intermediary stages of image evolution, allowing for interpretability through latent dynamics while preserving spatial correlations with the image. By incorporating a latent variable that follows a physical model expressed in partial differential equations, our approach ensures the interpretability of the learned model and provides insight into corresponding image dynamics. We demonstrate the robustness and effectiveness of our learning framework through a series of numerical tests using geoscientific imagery data.
SIAM Journal on Scientific Computing, Volume 46, Issue 5, Page C608C632, October 2024. <br/> Abstract. We present a framework designed to learn the underlying dynamics between two images observed at consecutive time steps. The complex nature of image data and the lack of temporal information pose significant challenges in capturing the unique evolving patterns. Our proposed method focuses on estimating the intermediary stages of image evolution, allowing for interpretability through latent dynamics while preserving spatial correlations with the image. By incorporating a latent variable that follows a physical model expressed in partial differential equations, our approach ensures the interpretability of the learned model and provides insight into corresponding image dynamics. We demonstrate the robustness and effectiveness of our learning framework through a series of numerical tests using geoscientific imagery data. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Learning InBetween Imagery Dynamics via Physical Latent Spaces
10.1137/23M1609440
SIAM Journal on Scientific Computing
20240927T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Jihun Han
Yoonsang Lee
Anne Gelb
Learning InBetween Imagery Dynamics via Physical Latent Spaces
46
5
C608
C632
20241031T07:00:00Z
20241031T07:00:00Z
10.1137/23M1609440
https://epubs.siam.org/doi/abs/10.1137/23M1609440?af=R
© 2024 Society for Industrial and Applied Mathematics

Global Minimization of Polynomial Integral Functionals
https://epubs.siam.org/doi/abs/10.1137/23M1592584?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2123A2149, August 2024. <br/> Abstract. We describe a “discretizethenrelax” strategy to globally minimize integral functionals over functions [math] in a Sobolev space subject to Dirichlet boundary conditions. The strategy applies whenever the integral functional depends polynomially on [math] and its derivatives, even if it is nonconvex. The “discretize” step uses a bounded finite element scheme to approximate the integral minimization problem with a convergent hierarchy of polynomial optimization problems over a compact feasible set, indexed by the decreasing size [math] of the finite element mesh. The “relax” step employs sparse momentsumofsquares relaxations to approximate each polynomial optimization problem with a hierarchy of convex semidefinite programs, indexed by an increasing relaxation order [math]. We prove that, as [math] and [math], solutions of such semidefinite programs provide approximate minimizers that converge in a suitable sense (including in certain [math] norms) to the global minimizer of the original integral functional if it is unique. We also report computational experiments showing that our numerical strategy works well even when technical conditions required by our theoretical analysis are not satisfied.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2123A2149, August 2024. <br/> Abstract. We describe a “discretizethenrelax” strategy to globally minimize integral functionals over functions [math] in a Sobolev space subject to Dirichlet boundary conditions. The strategy applies whenever the integral functional depends polynomially on [math] and its derivatives, even if it is nonconvex. The “discretize” step uses a bounded finite element scheme to approximate the integral minimization problem with a convergent hierarchy of polynomial optimization problems over a compact feasible set, indexed by the decreasing size [math] of the finite element mesh. The “relax” step employs sparse momentsumofsquares relaxations to approximate each polynomial optimization problem with a hierarchy of convex semidefinite programs, indexed by an increasing relaxation order [math]. We prove that, as [math] and [math], solutions of such semidefinite programs provide approximate minimizers that converge in a suitable sense (including in certain [math] norms) to the global minimizer of the original integral functional if it is unique. We also report computational experiments showing that our numerical strategy works well even when technical conditions required by our theoretical analysis are not satisfied. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Global Minimization of Polynomial Integral Functionals
10.1137/23M1592584
SIAM Journal on Scientific Computing
20240701T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Giovanni Fantuzzi
Federico Fuentes
Global Minimization of Polynomial Integral Functionals
46
4
A2123
A2149
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1592584
https://epubs.siam.org/doi/abs/10.1137/23M1592584?af=R
© 2024 Society for Industrial and Applied Mathematics

A Highly Efficient and Accurate DivergenceFree Spectral Method for the CurlCurl Equation in Two and Three Dimensions
https://epubs.siam.org/doi/abs/10.1137/23M1587038?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2150A2177, August 2024. <br/> Abstract. In this paper, we present a fast divergencefree spectral algorithm for the curlcurl problem. Divergencefree bases in two and three dimensions are constructed by using the generalized Jacobi polynomials. An accurate spectral method with exact preservation of the divergencefree constraint pointwisely is then proposed, and its corresponding error estimate is established. We then present a highly efficient solution algorithm based on a combination of the matrixfree preconditioned Krylov subspace iterative method and a fully diagonalizable auxiliary problem, which is derived from the spectral discretizations of generalized eigenvalue problems of Laplace and biharmonic operators. We rigorously prove that the dimensions of the invariant subspace of the preconditioned linear system resulting from the divergencefree spectral method with respect to the dominant eigenvalue 1 are [math] and [math] for two and threedimensional problems with [math] and [math] unknowns, respectively. Thus, the proposed method usually takes only several iterations to converge, and, astonishingly, as the problem size (polynomial order) increases, the number of iterations will decrease, even for highly indefinite system and oscillatory solutions. As a result, the computational cost of the solution algorithm is only a small multiple of [math] and [math] floating number operations for two and threedimensional problems, respectively. Plenty of numerical examples for solving the curlcurl problem with both constant and variable coefficients in two and three dimensions are presented to demonstrate the accuracy and efficiency of the proposed method.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2150A2177, August 2024. <br/> Abstract. In this paper, we present a fast divergencefree spectral algorithm for the curlcurl problem. Divergencefree bases in two and three dimensions are constructed by using the generalized Jacobi polynomials. An accurate spectral method with exact preservation of the divergencefree constraint pointwisely is then proposed, and its corresponding error estimate is established. We then present a highly efficient solution algorithm based on a combination of the matrixfree preconditioned Krylov subspace iterative method and a fully diagonalizable auxiliary problem, which is derived from the spectral discretizations of generalized eigenvalue problems of Laplace and biharmonic operators. We rigorously prove that the dimensions of the invariant subspace of the preconditioned linear system resulting from the divergencefree spectral method with respect to the dominant eigenvalue 1 are [math] and [math] for two and threedimensional problems with [math] and [math] unknowns, respectively. Thus, the proposed method usually takes only several iterations to converge, and, astonishingly, as the problem size (polynomial order) increases, the number of iterations will decrease, even for highly indefinite system and oscillatory solutions. As a result, the computational cost of the solution algorithm is only a small multiple of [math] and [math] floating number operations for two and threedimensional problems, respectively. Plenty of numerical examples for solving the curlcurl problem with both constant and variable coefficients in two and three dimensions are presented to demonstrate the accuracy and efficiency of the proposed method. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A Highly Efficient and Accurate DivergenceFree Spectral Method for the CurlCurl Equation in Two and Three Dimensions
10.1137/23M1587038
SIAM Journal on Scientific Computing
20240701T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Lechang Qin
Changtao Sheng
Zhiguo Yang
A Highly Efficient and Accurate DivergenceFree Spectral Method for the CurlCurl Equation in Two and Three Dimensions
46
4
A2150
A2177
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1587038
https://epubs.siam.org/doi/abs/10.1137/23M1587038?af=R
© 2024 Society for Industrial and Applied Mathematics

Transformed Model Reduction for Partial Differential Equations with Sharp Inner Layers
https://epubs.siam.org/doi/abs/10.1137/23M1589980?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2178A2201, August 2024. <br/> Abstract. Small parameters in partial differential equations can give rise to solutions with sharp inner layers that evolve over time. However, the standard model reduction method becomes inefficient when applied to these problems due to the slow decaying Kolmogorov [math]width of the solution manifold. To address this issue, a natural approach is to transform the equation in such a way that the transformed solution manifold exhibits a fast decaying Kolmogorov [math]width. In this paper, we focus on the Allen–Cahn equation as a model problem. We employ asymptotic analysis to identify slow variables and perform a transformation of the partial differential equations accordingly. Subsequently, we apply the proper orthogonal decomposition method and a QR discrete empirical interpolation method (qDEIM) technique to the transformed equation with the slow variables. Numerical experiments demonstrate that the new model reduction method yields significantly improved results compared to direct model reduction applied to the original equation. Furthermore, this approach can be extended to other equations, such as the convection equation and the Burgers equation. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/toreanony/TransformedModelReduction and in the supplementary materials (TransformedModelReductionmaster.zip [19.1KB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2178A2201, August 2024. <br/> Abstract. Small parameters in partial differential equations can give rise to solutions with sharp inner layers that evolve over time. However, the standard model reduction method becomes inefficient when applied to these problems due to the slow decaying Kolmogorov [math]width of the solution manifold. To address this issue, a natural approach is to transform the equation in such a way that the transformed solution manifold exhibits a fast decaying Kolmogorov [math]width. In this paper, we focus on the Allen–Cahn equation as a model problem. We employ asymptotic analysis to identify slow variables and perform a transformation of the partial differential equations accordingly. Subsequently, we apply the proper orthogonal decomposition method and a QR discrete empirical interpolation method (qDEIM) technique to the transformed equation with the slow variables. Numerical experiments demonstrate that the new model reduction method yields significantly improved results compared to direct model reduction applied to the original equation. Furthermore, this approach can be extended to other equations, such as the convection equation and the Burgers equation. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/toreanony/TransformedModelReduction and in the supplementary materials (TransformedModelReductionmaster.zip [19.1KB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Transformed Model Reduction for Partial Differential Equations with Sharp Inner Layers
10.1137/23M1589980
SIAM Journal on Scientific Computing
20240704T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Tianyou Tang
Xianmin Xu
Transformed Model Reduction for Partial Differential Equations with Sharp Inner Layers
46
4
A2178
A2201
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1589980
https://epubs.siam.org/doi/abs/10.1137/23M1589980?af=R
© 2024 Society for Industrial and Applied Mathematics

Modeling Diffusion in One Dimensional Discontinuous Media under Generalized Permeable Interface Conditions: Theory and Algorithms
https://epubs.siam.org/doi/abs/10.1137/23M1590846?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2202A2223, August 2024. <br/> Abstract. Diffusive transport in media with discontinuous properties is a challenging problem that arises in many applications. This paper focuses on one dimensional discontinuous media with generalized permeable boundary conditions at the discontinuity interface. It presents novel analytical expressions from the method of images to simulate diffusive processes, such as mass or thermal transport. The analytical expressions are used to formulate a generalization of the existing Skew Brownian Motion, HYMLA, and Uffink’s method, here named as GSBM, GHYMLA, and GUM, respectively, to handle generic interface conditions. The algorithms rely upon the random walk method and are tested by simulating transport in a bimaterial and in a multilayered medium with piecewise constant properties. The results indicate that the GUM algorithm provides the best performance in terms of accuracy and computational cost. The methods proposed can be applied for simulation of a wide range of differential problems.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2202A2223, August 2024. <br/> Abstract. Diffusive transport in media with discontinuous properties is a challenging problem that arises in many applications. This paper focuses on one dimensional discontinuous media with generalized permeable boundary conditions at the discontinuity interface. It presents novel analytical expressions from the method of images to simulate diffusive processes, such as mass or thermal transport. The analytical expressions are used to formulate a generalization of the existing Skew Brownian Motion, HYMLA, and Uffink’s method, here named as GSBM, GHYMLA, and GUM, respectively, to handle generic interface conditions. The algorithms rely upon the random walk method and are tested by simulating transport in a bimaterial and in a multilayered medium with piecewise constant properties. The results indicate that the GUM algorithm provides the best performance in terms of accuracy and computational cost. The methods proposed can be applied for simulation of a wide range of differential problems. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Modeling Diffusion in One Dimensional Discontinuous Media under Generalized Permeable Interface Conditions: Theory and Algorithms
10.1137/23M1590846
SIAM Journal on Scientific Computing
20240704T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Elisa Baioni
Antoine Lejay
Géraldine Pichot
Giovanni Michele Porta
Modeling Diffusion in One Dimensional Discontinuous Media under Generalized Permeable Interface Conditions: Theory and Algorithms
46
4
A2202
A2223
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1590846
https://epubs.siam.org/doi/abs/10.1137/23M1590846?af=R
© 2024 Society for Industrial and Applied Mathematics

A New Thermodynamically Compatible Finite Volume Scheme for Lagrangian Gas Dynamics
https://epubs.siam.org/doi/abs/10.1137/23M1580863?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2224A2247, August 2024. <br/> Abstract. The equations of Lagrangian gas dynamics fall into the larger class of overdetermined hyperbolic and thermodynamically compatible (HTC) systems of partial differential equations. They satisfy an entropy inequality (second principle of thermodynamics) and conserve total energy (first principle of thermodynamics). The aim of this work is to construct a novel thermodynamically compatible cellcentered Lagrangian finite volume scheme on unstructured meshes. Unlike in existing schemes, we choose to directly discretize the entropy inequality, hence obtaining total energy conservation as a consequence of the new thermodynamically compatible discretization of the other equations. First, the governing equations are written in fluctuation form. Next, the noncompatible centered numerical fluxes are corrected according to the approach recently introduced by Abgrall et al. using a scalar correction factor that is defined at the nodes of the grid. This perfectly fits into the formalism of nodal solvers which is typically adopted in cellcentered Lagrangian finite volume methods. Semidiscrete entropy conservative and entropy stable Lagrangian schemes are devised, and they are adequately blended together via a convex combination based on either a priori or a posteriori detectors of discontinuous solutions. The nonlinear stability in the energy norm is rigorously demonstrated, and the new schemes are provably positivity preserving for density and pressure. Furthermore, they exhibit zero numerical diffusion for isentropic flows while still being nonlinearly stable. The new schemes are tested against classical benchmarks for Lagrangian hydrodynamics, assessing their convergence and robustness and comparing their numerical dissipation with classical Lagrangian finite volume methods.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2224A2247, August 2024. <br/> Abstract. The equations of Lagrangian gas dynamics fall into the larger class of overdetermined hyperbolic and thermodynamically compatible (HTC) systems of partial differential equations. They satisfy an entropy inequality (second principle of thermodynamics) and conserve total energy (first principle of thermodynamics). The aim of this work is to construct a novel thermodynamically compatible cellcentered Lagrangian finite volume scheme on unstructured meshes. Unlike in existing schemes, we choose to directly discretize the entropy inequality, hence obtaining total energy conservation as a consequence of the new thermodynamically compatible discretization of the other equations. First, the governing equations are written in fluctuation form. Next, the noncompatible centered numerical fluxes are corrected according to the approach recently introduced by Abgrall et al. using a scalar correction factor that is defined at the nodes of the grid. This perfectly fits into the formalism of nodal solvers which is typically adopted in cellcentered Lagrangian finite volume methods. Semidiscrete entropy conservative and entropy stable Lagrangian schemes are devised, and they are adequately blended together via a convex combination based on either a priori or a posteriori detectors of discontinuous solutions. The nonlinear stability in the energy norm is rigorously demonstrated, and the new schemes are provably positivity preserving for density and pressure. Furthermore, they exhibit zero numerical diffusion for isentropic flows while still being nonlinearly stable. The new schemes are tested against classical benchmarks for Lagrangian hydrodynamics, assessing their convergence and robustness and comparing their numerical dissipation with classical Lagrangian finite volume methods. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A New Thermodynamically Compatible Finite Volume Scheme for Lagrangian Gas Dynamics
10.1137/23M1580863
SIAM Journal on Scientific Computing
20240710T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Walter Boscheri
Michael Dumbser
PierreHenri Maire
A New Thermodynamically Compatible Finite Volume Scheme for Lagrangian Gas Dynamics
46
4
A2224
A2247
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1580863
https://epubs.siam.org/doi/abs/10.1137/23M1580863?af=R
© 2024 Society for Industrial and Applied Mathematics

Local StructurePreserving Relaxation Method for Equilibrium of Charged Systems on Unstructured Meshes
https://epubs.siam.org/doi/abs/10.1137/23M1607234?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2248A2269, August 2024. <br/> Abstract. This work considers charged systems described by the modified Poisson–Nernst–Planck (PNP) equations, which incorporate ionic steric effects and the Born solvation energy for dielectric inhomogeneity. Solving the equilibrium of modified PNP equations poses numerical challenges due to the emergence of sharp boundary layers caused by small Debye lengths, particularly when local ionic concentrations reach saturation. To address this, we first reformulate the problem as a constraint optimization, where the ionic concentrations on unstructured Delaunay nodes are treated as fractional particles moving along edges between nodes. The electric fields are then updated to minimize the objective free energy while satisfying the discrete Gauss law. We develop a local relaxation method on unstructured meshes that inherently respects the discrete Gauss law, ensuring curlfree electric fields. Numerical analysis demonstrates that the optimal mass of the moving fractional particles guarantees the positivity of both ionic and solvent concentrations. Additionally, the free energy of the charged system consistently decreases during successive updates of ionic concentrations and electric fields. We conduct numerical tests to validate the expected numerical accuracy, positivity, freeenergy dissipation, and robustness of our method in simulating charged systems with sharp boundary layers.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2248A2269, August 2024. <br/> Abstract. This work considers charged systems described by the modified Poisson–Nernst–Planck (PNP) equations, which incorporate ionic steric effects and the Born solvation energy for dielectric inhomogeneity. Solving the equilibrium of modified PNP equations poses numerical challenges due to the emergence of sharp boundary layers caused by small Debye lengths, particularly when local ionic concentrations reach saturation. To address this, we first reformulate the problem as a constraint optimization, where the ionic concentrations on unstructured Delaunay nodes are treated as fractional particles moving along edges between nodes. The electric fields are then updated to minimize the objective free energy while satisfying the discrete Gauss law. We develop a local relaxation method on unstructured meshes that inherently respects the discrete Gauss law, ensuring curlfree electric fields. Numerical analysis demonstrates that the optimal mass of the moving fractional particles guarantees the positivity of both ionic and solvent concentrations. Additionally, the free energy of the charged system consistently decreases during successive updates of ionic concentrations and electric fields. We conduct numerical tests to validate the expected numerical accuracy, positivity, freeenergy dissipation, and robustness of our method in simulating charged systems with sharp boundary layers. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Local StructurePreserving Relaxation Method for Equilibrium of Charged Systems on Unstructured Meshes
10.1137/23M1607234
SIAM Journal on Scientific Computing
20240711T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Zhonghua Qiao
Zhenli Xu
Qian Yin
Shenggao Zhou
Local StructurePreserving Relaxation Method for Equilibrium of Charged Systems on Unstructured Meshes
46
4
A2248
A2269
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1607234
https://epubs.siam.org/doi/abs/10.1137/23M1607234?af=R
© 2024 Society for Industrial and Applied Mathematics

Efficient TensorProduct SpectralElement Operators with the SummationbyParts Property on Curved Triangles and Tetrahedra
https://epubs.siam.org/doi/abs/10.1137/23M1573963?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2270A2297, August 2024. <br/> Abstract. We present an extension of the summationbyparts (SBP) framework to tensorproduct spectralelement operators in collapsed coordinates. The proposed approach enables the construction of provably stable discretizations of arbitrary order which combine the geometric flexibility of unstructured triangular and tetrahedral meshes with the efficiency of sumfactorization algorithms. Specifically, a methodology is developed for constructing triangular and tetrahedral spectralelement operators of any order which possess the SBP property (i.e., satisfying a discrete analogue of integration by parts) as well as a tensorproduct decomposition. Such operators are then employed within the context of discontinuous spectralelement methods based on nodal expansions collocated at the tensorproduct quadrature nodes as well as modal expansions employing Proriol–Koornwinder–Dubiner polynomials, the latter approach resolving the time step limitation associated with the singularity of the collapsed coordinate transformation. Energystable formulations for curvilinear meshes are obtained using a skewsymmetric splitting of the metric terms, and a weightadjusted approximation is used to efficiently invert the curvilinear modal mass matrix. The proposed schemes are compared to those using nontensorial multidimensional SBP operators and are found to offer comparable accuracy to such schemes in the context of smooth linear advection problems on curved meshes, but at a reduced computational cost for higher polynomial degrees. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and Data Available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/tristanmontoya/ReproduceSBPSimplex and in the supplementary materials (reproducibility_repository.zip [35.7MB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2270A2297, August 2024. <br/> Abstract. We present an extension of the summationbyparts (SBP) framework to tensorproduct spectralelement operators in collapsed coordinates. The proposed approach enables the construction of provably stable discretizations of arbitrary order which combine the geometric flexibility of unstructured triangular and tetrahedral meshes with the efficiency of sumfactorization algorithms. Specifically, a methodology is developed for constructing triangular and tetrahedral spectralelement operators of any order which possess the SBP property (i.e., satisfying a discrete analogue of integration by parts) as well as a tensorproduct decomposition. Such operators are then employed within the context of discontinuous spectralelement methods based on nodal expansions collocated at the tensorproduct quadrature nodes as well as modal expansions employing Proriol–Koornwinder–Dubiner polynomials, the latter approach resolving the time step limitation associated with the singularity of the collapsed coordinate transformation. Energystable formulations for curvilinear meshes are obtained using a skewsymmetric splitting of the metric terms, and a weightadjusted approximation is used to efficiently invert the curvilinear modal mass matrix. The proposed schemes are compared to those using nontensorial multidimensional SBP operators and are found to offer comparable accuracy to such schemes in the context of smooth linear advection problems on curved meshes, but at a reduced computational cost for higher polynomial degrees. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and Data Available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/tristanmontoya/ReproduceSBPSimplex and in the supplementary materials (reproducibility_repository.zip [35.7MB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Efficient TensorProduct SpectralElement Operators with the SummationbyParts Property on Curved Triangles and Tetrahedra
10.1137/23M1573963
SIAM Journal on Scientific Computing
20240712T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Tristan Montoya
David W. Zingg
Efficient TensorProduct SpectralElement Operators with the SummationbyParts Property on Curved Triangles and Tetrahedra
46
4
A2270
A2297
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1573963
https://epubs.siam.org/doi/abs/10.1137/23M1573963?af=R
© 2024 Society for Industrial and Applied Mathematics

On Generalized Preconditioners for TimeParallel Parabolic Optimal Control
https://epubs.siam.org/doi/abs/10.1137/23M1553194?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2298A2323, August 2024. <br/> Abstract. The ParaDiag family of algorithms solves differential equations by using preconditioners that can be inverted in parallel through diagonalization. In the context of optimal control of linear parabolic PDEs, the stateoftheart ParaDiag method is limited to solving selfadjoint problems with a tracking objective. We propose three improvements to the ParaDiag method: the use of alphacirculant matrices to construct an alternative preconditioner, a generalization of the algorithm for solving nonselfadjoint equations, and the formulation of an algorithm for terminalcost objectives. We present novel analytic results about the eigenvalues of the preconditioned systems for all discussed ParaDiag algorithms in the case of selfadjoint equations, which proves the favorable properties of the alphacirculant preconditioner. We use these results to perform a theoretical parallelscaling analysis of ParaDiag for selfadjoint problems. Numerical tests confirm our findings and suggest that the selfadjoint behavior, which is backed by theory, generalizes to the nonselfadjoint case. We provide a sequential, opensource reference solver in MATLAB for all discussed algorithms. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available,” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://gitlab.kuleuven.be/numa/public/pintopt or in the supplementary materials (reprogeneralizedparadiag.zip [86.4KB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2298A2323, August 2024. <br/> Abstract. The ParaDiag family of algorithms solves differential equations by using preconditioners that can be inverted in parallel through diagonalization. In the context of optimal control of linear parabolic PDEs, the stateoftheart ParaDiag method is limited to solving selfadjoint problems with a tracking objective. We propose three improvements to the ParaDiag method: the use of alphacirculant matrices to construct an alternative preconditioner, a generalization of the algorithm for solving nonselfadjoint equations, and the formulation of an algorithm for terminalcost objectives. We present novel analytic results about the eigenvalues of the preconditioned systems for all discussed ParaDiag algorithms in the case of selfadjoint equations, which proves the favorable properties of the alphacirculant preconditioner. We use these results to perform a theoretical parallelscaling analysis of ParaDiag for selfadjoint problems. Numerical tests confirm our findings and suggest that the selfadjoint behavior, which is backed by theory, generalizes to the nonselfadjoint case. We provide a sequential, opensource reference solver in MATLAB for all discussed algorithms. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available,” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://gitlab.kuleuven.be/numa/public/pintopt or in the supplementary materials (reprogeneralizedparadiag.zip [86.4KB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
On Generalized Preconditioners for TimeParallel Parabolic Optimal Control
10.1137/23M1553194
SIAM Journal on Scientific Computing
20240712T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Arne Bouillon
Giovanni Samaey
Karl Meerbergen
On Generalized Preconditioners for TimeParallel Parabolic Optimal Control
46
4
A2298
A2323
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1553194
https://epubs.siam.org/doi/abs/10.1137/23M1553194?af=R
© 2024 Society for Industrial and Applied Mathematics

An Efficient FrequencyIndependent Numerical Method for Computing the FarField Pattern Induced by Polygonal Obstacles
https://epubs.siam.org/doi/abs/10.1137/23M1612160?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2324A2350, August 2024. <br/> Abstract. For problems of timeharmonic scattering by rational polygonal obstacles, embedding formulae express the farfield pattern induced by any incident plane wave in terms of the farfield patterns for a relatively small (frequencyindependent) set of canonical incident angles. Although these remarkable formulae are exact in theory, here we demonstrate that (i) they are highly sensitive to numerical errors in practice, and (ii) direct calculation of the coefficients in these formulae may be impossible for particular sets of canonical incident angles, even in exact arithmetic. Only by overcoming these practical issues can embedding formulae provide a highly efficient approach to computing the farfield pattern induced by a large number of incident angles. Here we address challenges (i) and (ii), supporting our theory with numerical experiments. Challenge (i) is solved using techniques from computational complex analysis: we reformulate the embedding formula as a complex contour integral and prove that this is much less sensitive to numerical errors. In practice, this contour integral can be efficiently evaluated by residue calculus. Challenge (ii) is addressed using techniques from numerical linear algebra: we oversample, considering more canonical incident angles than are necessary, thus expanding the set of valid coefficient vectors. The coefficient vector can then be selected using either a least squares approach or column subset selection.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2324A2350, August 2024. <br/> Abstract. For problems of timeharmonic scattering by rational polygonal obstacles, embedding formulae express the farfield pattern induced by any incident plane wave in terms of the farfield patterns for a relatively small (frequencyindependent) set of canonical incident angles. Although these remarkable formulae are exact in theory, here we demonstrate that (i) they are highly sensitive to numerical errors in practice, and (ii) direct calculation of the coefficients in these formulae may be impossible for particular sets of canonical incident angles, even in exact arithmetic. Only by overcoming these practical issues can embedding formulae provide a highly efficient approach to computing the farfield pattern induced by a large number of incident angles. Here we address challenges (i) and (ii), supporting our theory with numerical experiments. Challenge (i) is solved using techniques from computational complex analysis: we reformulate the embedding formula as a complex contour integral and prove that this is much less sensitive to numerical errors. In practice, this contour integral can be efficiently evaluated by residue calculus. Challenge (ii) is addressed using techniques from numerical linear algebra: we oversample, considering more canonical incident angles than are necessary, thus expanding the set of valid coefficient vectors. The coefficient vector can then be selected using either a least squares approach or column subset selection. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
An Efficient FrequencyIndependent Numerical Method for Computing the FarField Pattern Induced by Polygonal Obstacles
10.1137/23M1612160
SIAM Journal on Scientific Computing
20240717T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Andrew Gibbs
Stephen Langdon
An Efficient FrequencyIndependent Numerical Method for Computing the FarField Pattern Induced by Polygonal Obstacles
46
4
A2324
A2350
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1612160
https://epubs.siam.org/doi/abs/10.1137/23M1612160?af=R
© 2024 Society for Industrial and Applied Mathematics

An Adaptive Factorized Nyström Preconditioner for Regularized Kernel Matrices
https://epubs.siam.org/doi/abs/10.1137/23M1565139?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2351A2376, August 2024. <br/> Abstract. The spectrum of a kernel matrix significantly depends on the parameter values of the kernel function used to define the kernel matrix. This makes it challenging to design a preconditioner for a regularized kernel matrix that is robust across different parameter values. This paper proposes the adaptive factorized Nyström (AFN) preconditioner. The preconditioner is designed for the case where the rank [math] of the Nyström approximation is large, i.e., for kernel function parameters that lead to kernel matrices with eigenvalues that decay slowly. AFN deliberately chooses a wellconditioned submatrix to solve with and corrects a Nyström approximation with a factorized sparse approximate matrix inverse. This makes AFN efficient for kernel matrices with large numerical ranks. AFN also adaptively chooses the size of this submatrix to balance accuracy and cost. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/scalablematrix/H2Pack/tree/AFN_precond and in the supplementary materials (H2Pack.zip [3.99MB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2351A2376, August 2024. <br/> Abstract. The spectrum of a kernel matrix significantly depends on the parameter values of the kernel function used to define the kernel matrix. This makes it challenging to design a preconditioner for a regularized kernel matrix that is robust across different parameter values. This paper proposes the adaptive factorized Nyström (AFN) preconditioner. The preconditioner is designed for the case where the rank [math] of the Nyström approximation is large, i.e., for kernel function parameters that lead to kernel matrices with eigenvalues that decay slowly. AFN deliberately chooses a wellconditioned submatrix to solve with and corrects a Nyström approximation with a factorized sparse approximate matrix inverse. This makes AFN efficient for kernel matrices with large numerical ranks. AFN also adaptively chooses the size of this submatrix to balance accuracy and cost. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/scalablematrix/H2Pack/tree/AFN_precond and in the supplementary materials (H2Pack.zip [3.99MB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
An Adaptive Factorized Nyström Preconditioner for Regularized Kernel Matrices
10.1137/23M1565139
SIAM Journal on Scientific Computing
20240717T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Shifan Zhao
Tianshi Xu
Hua Huang
Edmond Chow
Yuanzhe Xi
An Adaptive Factorized Nyström Preconditioner for Regularized Kernel Matrices
46
4
A2351
A2376
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1565139
https://epubs.siam.org/doi/abs/10.1137/23M1565139?af=R
© 2024 Society for Industrial and Applied Mathematics

SuperLocalized Orthogonal Decomposition for HighFrequency Helmholtz Problems
https://epubs.siam.org/doi/abs/10.1137/21M1465950?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2377A2397, August 2024. <br/> Abstract. We propose a novel variant of the Localized Orthogonal Decomposition (LOD) method for timeharmonic scattering problems of Helmholtz type with high wavenumber [math]. On a coarse mesh of width [math], the proposed method identifies local finite element source terms that yield rapidly decaying responses under the solution operator. They can be constructed to high accuracy from independent local snapshot solutions on patches of width [math] and are used as problemadapted basis functions in the method. In contrast to the classical LOD and other stateoftheart multiscale methods, two and threedimensional numerical computations show that the localization error decays superexponentially as the oversampling parameter [math] is increased. This suggests that optimal convergence is observed under the substantially relaxed oversampling condition [math] with [math] denoting the spatial dimension. Numerical experiments demonstrate the significantly improved offline and online performance of the method also in the case of heterogeneous media and perfectly matched layers.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2377A2397, August 2024. <br/> Abstract. We propose a novel variant of the Localized Orthogonal Decomposition (LOD) method for timeharmonic scattering problems of Helmholtz type with high wavenumber [math]. On a coarse mesh of width [math], the proposed method identifies local finite element source terms that yield rapidly decaying responses under the solution operator. They can be constructed to high accuracy from independent local snapshot solutions on patches of width [math] and are used as problemadapted basis functions in the method. In contrast to the classical LOD and other stateoftheart multiscale methods, two and threedimensional numerical computations show that the localization error decays superexponentially as the oversampling parameter [math] is increased. This suggests that optimal convergence is observed under the substantially relaxed oversampling condition [math] with [math] denoting the spatial dimension. Numerical experiments demonstrate the significantly improved offline and online performance of the method also in the case of heterogeneous media and perfectly matched layers. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
SuperLocalized Orthogonal Decomposition for HighFrequency Helmholtz Problems
10.1137/21M1465950
SIAM Journal on Scientific Computing
20240718T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Philip Freese
Moritz Hauck
Daniel Peterseim
SuperLocalized Orthogonal Decomposition for HighFrequency Helmholtz Problems
46
4
A2377
A2397
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/21M1465950
https://epubs.siam.org/doi/abs/10.1137/21M1465950?af=R
© 2024 Society for Industrial and Applied Mathematics

Computing [math]Conforming Finite Element Approximations Without Having to Implement [math]Elements
https://epubs.siam.org/doi/abs/10.1137/23M1615486?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2398A2420, August 2024. <br/> Abstract. We develop a method to compute the [math]conforming finite element approximation to planar fourth order elliptic problems without having to implement [math] elements. The algorithm consists of replacing the original [math]conforming scheme with preprocessing and postprocessing steps that require only an [math]conforming Poisson type solve and an inner Stokeslike problem that again only requires at most [math]conformity. We then demonstrate the method applied to the Morgan–Scott elements with three numerical examples. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://doi.org/10.5281/zenodo.10070565.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2398A2420, August 2024. <br/> Abstract. We develop a method to compute the [math]conforming finite element approximation to planar fourth order elliptic problems without having to implement [math] elements. The algorithm consists of replacing the original [math]conforming scheme with preprocessing and postprocessing steps that require only an [math]conforming Poisson type solve and an inner Stokeslike problem that again only requires at most [math]conformity. We then demonstrate the method applied to the Morgan–Scott elements with three numerical examples. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://doi.org/10.5281/zenodo.10070565. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Computing [math]Conforming Finite Element Approximations Without Having to Implement [math]Elements
10.1137/23M1615486
SIAM Journal on Scientific Computing
20240722T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Mark Ainsworth
Charles Parker
Computing [math]Conforming Finite Element Approximations Without Having to Implement [math]Elements
46
4
A2398
A2420
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1615486
https://epubs.siam.org/doi/abs/10.1137/23M1615486?af=R
© 2024 Society for Industrial and Applied Mathematics

A Full Approximation Scheme Multilevel Method for Nonlinear Variational Inequalities
https://epubs.siam.org/doi/abs/10.1137/23M1594200?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2421A2444, August 2024. <br/> Abstract. We present the full approximation scheme constraint decomposition (FASCD) multilevel method for solving variational inequalities (VIs). FASCD is a joint extension of both the full approximation scheme multigrid technique for nonlinear partial differential equations, due to A. Brandt, and the constraint decomposition (CD) method introduced by X.C. Tai for VIs arising in optimization. We extend the CD idea by exploiting the telescoping nature of certain subset decompositions arising from multilevel mesh hierarchies. When a reducedspace (active set) Newton method is applied as a smoother, with work proportional to the number of unknowns on a given mesh level, FASCD Vcycles exhibit nearly meshindependent convergence rates. The full multigrid cycle version is an optimal solver. The example problems include differential operators which are symmetric linear, nonsymmetric linear, and nonlinear, in unilateral and bilateral VI problems. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://bitbucket.org/pefarrell/fascd/, where the software used to produce the results in section 8 is archived at tag v1.0, and at https://doi.org/10.5281/zenodo.10476845 or in the supplementary materials (pefarrellfascd6407e9f547d6.zip [225KB]). The authors used Firedrake master revision c5e939dde.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2421A2444, August 2024. <br/> Abstract. We present the full approximation scheme constraint decomposition (FASCD) multilevel method for solving variational inequalities (VIs). FASCD is a joint extension of both the full approximation scheme multigrid technique for nonlinear partial differential equations, due to A. Brandt, and the constraint decomposition (CD) method introduced by X.C. Tai for VIs arising in optimization. We extend the CD idea by exploiting the telescoping nature of certain subset decompositions arising from multilevel mesh hierarchies. When a reducedspace (active set) Newton method is applied as a smoother, with work proportional to the number of unknowns on a given mesh level, FASCD Vcycles exhibit nearly meshindependent convergence rates. The full multigrid cycle version is an optimal solver. The example problems include differential operators which are symmetric linear, nonsymmetric linear, and nonlinear, in unilateral and bilateral VI problems. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://bitbucket.org/pefarrell/fascd/, where the software used to produce the results in section 8 is archived at tag v1.0, and at https://doi.org/10.5281/zenodo.10476845 or in the supplementary materials (pefarrellfascd6407e9f547d6.zip [225KB]). The authors used Firedrake master revision c5e939dde. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A Full Approximation Scheme Multilevel Method for Nonlinear Variational Inequalities
10.1137/23M1594200
SIAM Journal on Scientific Computing
20240723T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Ed Bueler
Patrick E. Farrell
A Full Approximation Scheme Multilevel Method for Nonlinear Variational Inequalities
46
4
A2421
A2444
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1594200
https://epubs.siam.org/doi/abs/10.1137/23M1594200?af=R
© 2024 Society for Industrial and Applied Mathematics

Domain Decomposition Learning Methods for Solving Elliptic Problems
https://epubs.siam.org/doi/abs/10.1137/22M1515392?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2445A2474, August 2024. <br/> Abstract. With the aid of hardware and software developments, there has been a surge of interest in solving PDEs by deep learning techniques, and the integration with domain decomposition strategies has recently attracted considerable attention due to its enhanced representation and parallelization capacity of the network solution. While there are already several works that substitute the numerical solver of overlapping Schwarz methods with the deep learning approach, the nonoverlapping counterpart has not been thoroughly studied yet because of the inevitable interface overfitting problem that would propagate the errors to neighboring subdomains and eventually hamper the convergence of outer iteration. In this work, a novel learning approach, i.e., the compensated deep Ritz method using neural network extension operators, is proposed to enable the flux transmission across subregion interfaces with guaranteed accuracy, thereby allowing us to construct effective learning algorithms for realizing the more general nonoverlapping domain decomposition methods in the presence of overfitted interface conditions. Numerical experiments on a series of elliptic boundary value problems, including the regular and irregular interfaces, low and high dimensions, and smooth and highcontrast coefficients on multidomains, are carried out to validate the effectiveness of our proposed domain decomposition learning algorithms. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available" as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available in https://github.com/AI4SCTJU or in the supplementary materials.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2445A2474, August 2024. <br/> Abstract. With the aid of hardware and software developments, there has been a surge of interest in solving PDEs by deep learning techniques, and the integration with domain decomposition strategies has recently attracted considerable attention due to its enhanced representation and parallelization capacity of the network solution. While there are already several works that substitute the numerical solver of overlapping Schwarz methods with the deep learning approach, the nonoverlapping counterpart has not been thoroughly studied yet because of the inevitable interface overfitting problem that would propagate the errors to neighboring subdomains and eventually hamper the convergence of outer iteration. In this work, a novel learning approach, i.e., the compensated deep Ritz method using neural network extension operators, is proposed to enable the flux transmission across subregion interfaces with guaranteed accuracy, thereby allowing us to construct effective learning algorithms for realizing the more general nonoverlapping domain decomposition methods in the presence of overfitted interface conditions. Numerical experiments on a series of elliptic boundary value problems, including the regular and irregular interfaces, low and high dimensions, and smooth and highcontrast coefficients on multidomains, are carried out to validate the effectiveness of our proposed domain decomposition learning algorithms. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available" as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available in https://github.com/AI4SCTJU or in the supplementary materials. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Domain Decomposition Learning Methods for Solving Elliptic Problems
10.1137/22M1515392
SIAM Journal on Scientific Computing
20240723T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Qi Sun
Xuejun Xu
Haotian Yi
Domain Decomposition Learning Methods for Solving Elliptic Problems
46
4
A2445
A2474
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/22M1515392
https://epubs.siam.org/doi/abs/10.1137/22M1515392?af=R
© 2024 Society for Industrial and Applied Mathematics

Multilevel Particle Filters for a Class of Partially Observed Piecewise Deterministic Markov Processes
https://epubs.siam.org/doi/abs/10.1137/23M1600505?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2475A2502, August 2024. <br/> Abstract. In this paper we consider the filtering of a class of partially observed piecewise deterministic Markov processes. In particular, we assume that an ordinary differential equation (ODE) drives the deterministic element and can only be solved numerically via a time discretization. We develop, based upon the approach in Lemaire, Thieullen, and Thomas [Adv. Appl. Probab., 52 (2020), pp. 138–172], a new particle and multilevel particle filter (MLPF) in order to approximate the filter associated to the discretized ODE. We provide a bound on the mean square error associated to the MLPF which provides guidance on setting the simulation parameters of the algorithm and implies that significant computational gains can be obtained versus using a particle filter. Our theoretical claims are confirmed in several numerical examples.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2475A2502, August 2024. <br/> Abstract. In this paper we consider the filtering of a class of partially observed piecewise deterministic Markov processes. In particular, we assume that an ordinary differential equation (ODE) drives the deterministic element and can only be solved numerically via a time discretization. We develop, based upon the approach in Lemaire, Thieullen, and Thomas [Adv. Appl. Probab., 52 (2020), pp. 138–172], a new particle and multilevel particle filter (MLPF) in order to approximate the filter associated to the discretized ODE. We provide a bound on the mean square error associated to the MLPF which provides guidance on setting the simulation parameters of the algorithm and implies that significant computational gains can be obtained versus using a particle filter. Our theoretical claims are confirmed in several numerical examples. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Multilevel Particle Filters for a Class of Partially Observed Piecewise Deterministic Markov Processes
10.1137/23M1600505
SIAM Journal on Scientific Computing
20240726T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Ajay Jasra
Kengo Kamatani
Mohamed Maama
Multilevel Particle Filters for a Class of Partially Observed Piecewise Deterministic Markov Processes
46
4
A2475
A2502
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1600505
https://epubs.siam.org/doi/abs/10.1137/23M1600505?af=R
© 2024 Society for Industrial and Applied Mathematics

A SemiImplicit Fully Exactly WellBalanced Relaxation Scheme for the Shallow Water System
https://epubs.siam.org/doi/abs/10.1137/23M1621289?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2503A2527, August 2024. <br/> Abstract. This article focuses on the design of semiimplicit schemes that are fully wellbalanced for the onedimensional shallow water equations, that is, schemes that preserve all smooth steady states of the system and not just wateratrest equilibria. The proposed methods outperform standard explicit schemes in the lowFroude regime, where the celerity is much larger than the fluid velocity, eliminating the need for a large number of iterations on large time intervals. In this work, splitting and relaxation techniques are combined in order to obtain fully wellbalanced semiimplicit first and second order schemes. In contrast to recent Lagrangianbased approaches, this one allows the preservation of all the steady states while avoiding the complexities associated with Lagrangian formalism.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2503A2527, August 2024. <br/> Abstract. This article focuses on the design of semiimplicit schemes that are fully wellbalanced for the onedimensional shallow water equations, that is, schemes that preserve all smooth steady states of the system and not just wateratrest equilibria. The proposed methods outperform standard explicit schemes in the lowFroude regime, where the celerity is much larger than the fluid velocity, eliminating the need for a large number of iterations on large time intervals. In this work, splitting and relaxation techniques are combined in order to obtain fully wellbalanced semiimplicit first and second order schemes. In contrast to recent Lagrangianbased approaches, this one allows the preservation of all the steady states while avoiding the complexities associated with Lagrangian formalism. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A SemiImplicit Fully Exactly WellBalanced Relaxation Scheme for the Shallow Water System
10.1137/23M1621289
SIAM Journal on Scientific Computing
20240731T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Celia CaballeroCárdenas
Manuel Jesús Castro
Christophe Chalons
Tomás Morales de Luna
María Luz MuñozRuiz
A SemiImplicit Fully Exactly WellBalanced Relaxation Scheme for the Shallow Water System
46
4
A2503
A2527
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1621289
https://epubs.siam.org/doi/abs/10.1137/23M1621289?af=R
© 2024 Society for Industrial and Applied Mathematics

Energetic Variational Neural Network Discretizations of Gradient Flows
https://epubs.siam.org/doi/abs/10.1137/22M1529427?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2528A2556, August 2024. <br/> Abstract. We present a structurepreserving Eulerian algorithm for solving [math]gradient flows and a structurepreserving Lagrangian algorithm for solving generalized diffusions. Both algorithms employ neural networks as tools for spatial discretization. Unlike most existing methods that construct numerical discretizations based on the strong or weak form of the underlying PDE, the proposed schemes are constructed based on the energydissipation law directly. This guarantees the monotonic decay of the system’s free energy, which avoids unphysical states of solutions and is crucial for the longterm stability of numerical computations. To address challenges arising from nonlinear neural network discretization, we perform temporal discretizations on these variational systems before spatial discretizations. This approach is computationally memoryefficient when implementing neural networkbased algorithms. The proposed neural networkbased schemes are meshfree, allowing us to solve gradient flows in high dimensions. Various numerical experiments are presented to demonstrate the accuracy and energy stability of the proposed numerical schemes.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2528A2556, August 2024. <br/> Abstract. We present a structurepreserving Eulerian algorithm for solving [math]gradient flows and a structurepreserving Lagrangian algorithm for solving generalized diffusions. Both algorithms employ neural networks as tools for spatial discretization. Unlike most existing methods that construct numerical discretizations based on the strong or weak form of the underlying PDE, the proposed schemes are constructed based on the energydissipation law directly. This guarantees the monotonic decay of the system’s free energy, which avoids unphysical states of solutions and is crucial for the longterm stability of numerical computations. To address challenges arising from nonlinear neural network discretization, we perform temporal discretizations on these variational systems before spatial discretizations. This approach is computationally memoryefficient when implementing neural networkbased algorithms. The proposed neural networkbased schemes are meshfree, allowing us to solve gradient flows in high dimensions. Various numerical experiments are presented to demonstrate the accuracy and energy stability of the proposed numerical schemes. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Energetic Variational Neural Network Discretizations of Gradient Flows
10.1137/22M1529427
SIAM Journal on Scientific Computing
20240805T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Ziqing Hu
Chun Liu
Yiwei Wang
Zhiliang Xu
Energetic Variational Neural Network Discretizations of Gradient Flows
46
4
A2528
A2556
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/22M1529427
https://epubs.siam.org/doi/abs/10.1137/22M1529427?af=R
© 2024 Society for Industrial and Applied Mathematics

A Levenberg–Marquardt Method for Nonsmooth Regularized Least Squares
https://epubs.siam.org/doi/abs/10.1137/22M1538971?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2557A2581, August 2024. <br/> Abstract. We develop a Levenberg–Marquardt method for minimizing the sum of a smooth nonlinear leastsquares term [math] and a nonsmooth term [math]. Both [math] and [math] may be nonconvex. Steps are computed by minimizing the sum of a regularized linear leastsquares model and a model of [math] using a firstorder method such as the proximal gradient method. We establish global convergence to a firstorder stationary point under the assumptions that [math] and its Jacobian are Lipschitz continuous and [math] is proper and lower semicontinuous. In the worst case, our method performs [math] iterations to bring a measure of stationarity below [math]. We also derive a trustregion variant that enjoys similar asymptotic worstcase iteration complexity as a special case of the trustregion algorithm of Aravkin, Baraldi, and Orban [SIAM J. Optim., 32 (2022), pp. 900–929]. We report numerical results on three examples: a grouplasso basispursuit denoise example, a nonlinear support vector machine, and parameter estimation in a neuroscience application. To implement those examples, we describe in detail how to evaluate proximal operators for separable [math] and for the group lasso with trustregion constraint. In all cases, the Levenberg–Marquardt methods perform fewer outer iterations than either a proximal gradient method with adaptive step length or a quasiNewton trustregion method, neither of which exploit the leastsquares structure of the problem. Our results also highlight the need for more sophisticated subproblem solvers than simple firstorder methods.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2557A2581, August 2024. <br/> Abstract. We develop a Levenberg–Marquardt method for minimizing the sum of a smooth nonlinear leastsquares term [math] and a nonsmooth term [math]. Both [math] and [math] may be nonconvex. Steps are computed by minimizing the sum of a regularized linear leastsquares model and a model of [math] using a firstorder method such as the proximal gradient method. We establish global convergence to a firstorder stationary point under the assumptions that [math] and its Jacobian are Lipschitz continuous and [math] is proper and lower semicontinuous. In the worst case, our method performs [math] iterations to bring a measure of stationarity below [math]. We also derive a trustregion variant that enjoys similar asymptotic worstcase iteration complexity as a special case of the trustregion algorithm of Aravkin, Baraldi, and Orban [SIAM J. Optim., 32 (2022), pp. 900–929]. We report numerical results on three examples: a grouplasso basispursuit denoise example, a nonlinear support vector machine, and parameter estimation in a neuroscience application. To implement those examples, we describe in detail how to evaluate proximal operators for separable [math] and for the group lasso with trustregion constraint. In all cases, the Levenberg–Marquardt methods perform fewer outer iterations than either a proximal gradient method with adaptive step length or a quasiNewton trustregion method, neither of which exploit the leastsquares structure of the problem. Our results also highlight the need for more sophisticated subproblem solvers than simple firstorder methods. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A Levenberg–Marquardt Method for Nonsmooth Regularized Least Squares
10.1137/22M1538971
SIAM Journal on Scientific Computing
20240812T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Aleksandr Y. Aravkin
Robert Baraldi
Dominique Orban
A Levenberg–Marquardt Method for Nonsmooth Regularized Least Squares
46
4
A2557
A2581
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/22M1538971
https://epubs.siam.org/doi/abs/10.1137/22M1538971?af=R
© 2024 Society for Industrial and Applied Mathematics

A HighOrder Fast Direct Solver for Surface PDEs
https://epubs.siam.org/doi/abs/10.1137/22M1525259?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2582A2606, August 2024. <br/> Abstract. We introduce a fast direct solver for variablecoefficient elliptic PDEs on surfaces based on the hierarchical Poincaré–Steklov method. The method takes as input an unstructured, highorder quadrilateral mesh of a surface and discretizes surface differential operators on each element using a highorder spectral collocation scheme. Elemental solution operators and DirichlettoNeumann maps tangent to the surface are precomputed and merged in a pairwise fashion to yield a hierarchy of solution operators that may be applied in [math] operations for a mesh with [math] degrees of freedom. The resulting fast direct solver may be used to accelerate highorder implicit timestepping schemes, as the precomputed operators can be reused for fast elliptic solves on surfaces. On a standard laptop, precomputation for a 12thorder surface mesh with over 1 million degrees of freedom takes 10 seconds, while subsequent solves take only 0.25 seconds. We apply the method to a range of problems on both smooth surfaces and surfaces with sharp corners and edges, including the static Laplace–Beltrami problem, the Hodge decomposition of a tangential vector field, and some timedependent nonlinear reactiondiffusion systems. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: code and data available”, as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/danfortunato/surfacehpssisc.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2582A2606, August 2024. <br/> Abstract. We introduce a fast direct solver for variablecoefficient elliptic PDEs on surfaces based on the hierarchical Poincaré–Steklov method. The method takes as input an unstructured, highorder quadrilateral mesh of a surface and discretizes surface differential operators on each element using a highorder spectral collocation scheme. Elemental solution operators and DirichlettoNeumann maps tangent to the surface are precomputed and merged in a pairwise fashion to yield a hierarchy of solution operators that may be applied in [math] operations for a mesh with [math] degrees of freedom. The resulting fast direct solver may be used to accelerate highorder implicit timestepping schemes, as the precomputed operators can be reused for fast elliptic solves on surfaces. On a standard laptop, precomputation for a 12thorder surface mesh with over 1 million degrees of freedom takes 10 seconds, while subsequent solves take only 0.25 seconds. We apply the method to a range of problems on both smooth surfaces and surfaces with sharp corners and edges, including the static Laplace–Beltrami problem, the Hodge decomposition of a tangential vector field, and some timedependent nonlinear reactiondiffusion systems. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: code and data available”, as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/danfortunato/surfacehpssisc. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A HighOrder Fast Direct Solver for Surface PDEs
10.1137/22M1525259
SIAM Journal on Scientific Computing
20240813T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Daniel Fortunato
A HighOrder Fast Direct Solver for Surface PDEs
46
4
A2582
A2606
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/22M1525259
https://epubs.siam.org/doi/abs/10.1137/22M1525259?af=R
© 2024 Society for Industrial and Applied Mathematics

A Preconditioned Krylov Subspace Method for Linear Inverse Problems with GeneralForm Tikhonov Regularization
https://epubs.siam.org/doi/abs/10.1137/23M1593802?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2607A2633, August 2024. <br/> Abstract. Tikhonov regularization is a widely used technique in solving inverse problems that can enforce prior properties on the desired solution. In this paper, we propose a Krylov subspace based iterative method for solving linear inverse problems with generalform Tikhonov regularization term [math], where [math] is a positive semidefinite matrix. An iterative process called the preconditioned Golub–Kahan bidiagonalization (pGKB) is designed, which implicitly utilizes a proper preconditioner to generate a series of solution subspaces with desirable properties encoded by the regularizer [math]. Based on the pGKB process, we propose an iterative regularization algorithm via projecting the original problem onto small dimensional solution subspaces. We analyze the regularization properties of this algorithm, including the incorporation of prior properties of the desired solution into the solution subspace and the semiconvergence behavior of the regularized solution. To overcome instabilities caused by semiconvergence, we further propose two pGKB based hybrid regularization algorithms. All the proposed algorithms are tested on both smallscale and largescale linear inverse problems. Numerical results demonstrate that these iterative algorithms exhibit excellent performance, outperforming other stateoftheart algorithms in some cases.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2607A2633, August 2024. <br/> Abstract. Tikhonov regularization is a widely used technique in solving inverse problems that can enforce prior properties on the desired solution. In this paper, we propose a Krylov subspace based iterative method for solving linear inverse problems with generalform Tikhonov regularization term [math], where [math] is a positive semidefinite matrix. An iterative process called the preconditioned Golub–Kahan bidiagonalization (pGKB) is designed, which implicitly utilizes a proper preconditioner to generate a series of solution subspaces with desirable properties encoded by the regularizer [math]. Based on the pGKB process, we propose an iterative regularization algorithm via projecting the original problem onto small dimensional solution subspaces. We analyze the regularization properties of this algorithm, including the incorporation of prior properties of the desired solution into the solution subspace and the semiconvergence behavior of the regularized solution. To overcome instabilities caused by semiconvergence, we further propose two pGKB based hybrid regularization algorithms. All the proposed algorithms are tested on both smallscale and largescale linear inverse problems. Numerical results demonstrate that these iterative algorithms exhibit excellent performance, outperforming other stateoftheart algorithms in some cases. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A Preconditioned Krylov Subspace Method for Linear Inverse Problems with GeneralForm Tikhonov Regularization
10.1137/23M1593802
SIAM Journal on Scientific Computing
20240814T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Haibo Li
A Preconditioned Krylov Subspace Method for Linear Inverse Problems with GeneralForm Tikhonov Regularization
46
4
A2607
A2633
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1593802
https://epubs.siam.org/doi/abs/10.1137/23M1593802?af=R
© 2024 Society for Industrial and Applied Mathematics

The Effect of Approximate CoarsestLevel Solves on the Convergence of Multigrid VCycle Methods
https://epubs.siam.org/doi/abs/10.1137/23M1578255?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2634A2659, August 2024. <br/> Abstract. The multigrid Vcycle method is a popular method for solving systems of linear equations. It computes an approximate solution by using smoothing on fine levels and solving a system of linear equations on the coarsest level. Solving on the coarsest level depends on the size and difficulty of the problem. If the size permits, it is typical to use a direct method based on LU or Cholesky decomposition. In settings with large coarsestlevel problems, approximate solvers such as iterative Krylov subspace methods, or direct methods based on lowrank approximation, are often used. The accuracy of the coarsestlevel solver is typically determined based on the experience of the users with the concrete problems and methods. In this paper, we present an approach to analyzing the effects of approximate coarsestlevel solves on the convergence of the Vcycle method for symmetric positive definite problems. Using these results, we derive coarsestlevel stopping criterion through which we may control the difference between the approximation computed by a Vcycle method with approximate coarsestlevel solver and the approximation which would be computed if the coarsestlevel problems were solved exactly. The coarsestlevel stopping criterion may thus be set up such that the Vcycle method converges to a chosen finestlevel accuracy in (nearly) the same number of Vcycle iterations as the Vcycle method with exact coarsestlevel solver. We also utilize the theoretical results to discuss how the convergence of the Vcycle method may be affected by the choice of a tolerance in a coarsestlevel stopping criterion based on the relative residual norm. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://doi.org/10.5281/zenodo.11178544.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2634A2659, August 2024. <br/> Abstract. The multigrid Vcycle method is a popular method for solving systems of linear equations. It computes an approximate solution by using smoothing on fine levels and solving a system of linear equations on the coarsest level. Solving on the coarsest level depends on the size and difficulty of the problem. If the size permits, it is typical to use a direct method based on LU or Cholesky decomposition. In settings with large coarsestlevel problems, approximate solvers such as iterative Krylov subspace methods, or direct methods based on lowrank approximation, are often used. The accuracy of the coarsestlevel solver is typically determined based on the experience of the users with the concrete problems and methods. In this paper, we present an approach to analyzing the effects of approximate coarsestlevel solves on the convergence of the Vcycle method for symmetric positive definite problems. Using these results, we derive coarsestlevel stopping criterion through which we may control the difference between the approximation computed by a Vcycle method with approximate coarsestlevel solver and the approximation which would be computed if the coarsestlevel problems were solved exactly. The coarsestlevel stopping criterion may thus be set up such that the Vcycle method converges to a chosen finestlevel accuracy in (nearly) the same number of Vcycle iterations as the Vcycle method with exact coarsestlevel solver. We also utilize the theoretical results to discuss how the convergence of the Vcycle method may be affected by the choice of a tolerance in a coarsestlevel stopping criterion based on the relative residual norm. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://doi.org/10.5281/zenodo.11178544. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
The Effect of Approximate CoarsestLevel Solves on the Convergence of Multigrid VCycle Methods
10.1137/23M1578255
SIAM Journal on Scientific Computing
20240814T07:00:00Z
© 2024 SIAM. Published by SIAM under the terms of the Creative Commons 4.0 license
Petr Vacek
Erin Carson
Kirk M. Soodhalter
The Effect of Approximate CoarsestLevel Solves on the Convergence of Multigrid VCycle Methods
46
4
A2634
A2659
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1578255
https://epubs.siam.org/doi/abs/10.1137/23M1578255?af=R
© 2024 SIAM. Published by SIAM under the terms of the Creative Commons 4.0 license

Implicitexplicit Schemes for Incompressible Flow Problems with Variable Viscosity
https://epubs.siam.org/doi/abs/10.1137/23M1606526?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2660A2682, August 2024. <br/> Abstract. This article investigates different implicitexplicit (IMEX) methods for incompressible flows with variable viscosity. The viscosity field may depend on space and time alone or, for example, on velocity gradients. Unlike most previous works on IMEX schemes, which focus on the convective term, we propose also treating parts of the diffusive term explicitly, which can reduce the coupling between the velocity components. We present different IMEX alternatives for the variableviscosity Navier–Stokes system, analyzing their theoretical and algorithmic properties. Temporal stability is proven for all the methods presented, including monolithic and fractionalstep variants. These results are unconditional except for one of the fractionalstep discretizations, whose stability is shown for timestep sizes under an upper bound that depends solely on the problem data. The key finding of this work is a class of IMEX schemes whose steps decouple the velocity components and are fully linearized (even if the viscosity depends nonlinearly on the velocity) without requiring any CFL condition for stability. Moreover, in the presence of Neumann boundaries, some of our formulations lead naturally to conditions involving normal pseudotractions. This generalizes to the variableviscosity case what happens for the standard Laplacian form with constant viscosity. Our analysis is supported by a series of numerical experiments.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2660A2682, August 2024. <br/> Abstract. This article investigates different implicitexplicit (IMEX) methods for incompressible flows with variable viscosity. The viscosity field may depend on space and time alone or, for example, on velocity gradients. Unlike most previous works on IMEX schemes, which focus on the convective term, we propose also treating parts of the diffusive term explicitly, which can reduce the coupling between the velocity components. We present different IMEX alternatives for the variableviscosity Navier–Stokes system, analyzing their theoretical and algorithmic properties. Temporal stability is proven for all the methods presented, including monolithic and fractionalstep variants. These results are unconditional except for one of the fractionalstep discretizations, whose stability is shown for timestep sizes under an upper bound that depends solely on the problem data. The key finding of this work is a class of IMEX schemes whose steps decouple the velocity components and are fully linearized (even if the viscosity depends nonlinearly on the velocity) without requiring any CFL condition for stability. Moreover, in the presence of Neumann boundaries, some of our formulations lead naturally to conditions involving normal pseudotractions. This generalizes to the variableviscosity case what happens for the standard Laplacian form with constant viscosity. Our analysis is supported by a series of numerical experiments. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Implicitexplicit Schemes for Incompressible Flow Problems with Variable Viscosity
10.1137/23M1606526
SIAM Journal on Scientific Computing
20240819T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Gabriel Barrenechea
Ernesto Castillo
Douglas Pacheco
Implicitexplicit Schemes for Incompressible Flow Problems with Variable Viscosity
46
4
A2660
A2682
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1606526
https://epubs.siam.org/doi/abs/10.1137/23M1606526?af=R
© 2024 Society for Industrial and Applied Mathematics

Adaptivity in Local Kernel Based Methods for Approximating the Action of Linear Operators
https://epubs.siam.org/doi/abs/10.1137/23M1598052?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2683A2708, August 2024. <br/> Abstract. Building on the successes of local kernel methods for approximating the solutions to partial differential equations (PDEs) and the evaluation of definite integrals (quadrature/cubature), a local estimate of the error in such approximations is developed. This estimate is useful for determining locations in the solution domain where increased node density (equivalently, reduction in the spacing between nodes) can decrease the error in the solution. An adaptive procedure for adding nodes to the domain for both the approximation of derivatives and the approximate evaluation of definite integrals is described. This method efficiently computes the error estimate at a set of prescribed points and adds new nodes for approximation where the error is too large. Computational experiments demonstrate close agreement between the error estimate and actual absolute error in the approximation. Such methods are necessary or desirable when approximating solutions to PDEs (or in the case of quadrature/cubature), where the initial data and subsequent solution (or integrand) exhibit localized features that require significant refinement to resolve and where uniform increases in the density of nodes across the entire computational domain is not possible or too burdensome.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2683A2708, August 2024. <br/> Abstract. Building on the successes of local kernel methods for approximating the solutions to partial differential equations (PDEs) and the evaluation of definite integrals (quadrature/cubature), a local estimate of the error in such approximations is developed. This estimate is useful for determining locations in the solution domain where increased node density (equivalently, reduction in the spacing between nodes) can decrease the error in the solution. An adaptive procedure for adding nodes to the domain for both the approximation of derivatives and the approximate evaluation of definite integrals is described. This method efficiently computes the error estimate at a set of prescribed points and adds new nodes for approximation where the error is too large. Computational experiments demonstrate close agreement between the error estimate and actual absolute error in the approximation. Such methods are necessary or desirable when approximating solutions to PDEs (or in the case of quadrature/cubature), where the initial data and subsequent solution (or integrand) exhibit localized features that require significant refinement to resolve and where uniform increases in the density of nodes across the entire computational domain is not possible or too burdensome. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Adaptivity in Local Kernel Based Methods for Approximating the Action of Linear Operators
10.1137/23M1598052
SIAM Journal on Scientific Computing
20240819T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Jonah A. Reeger
Adaptivity in Local Kernel Based Methods for Approximating the Action of Linear Operators
46
4
A2683
A2708
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1598052
https://epubs.siam.org/doi/abs/10.1137/23M1598052?af=R
© 2024 Society for Industrial and Applied Mathematics

Multilevel Parareal Algorithm with Averaging for Oscillatory Problems
https://epubs.siam.org/doi/abs/10.1137/23M1547123?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2709A2736, August 2024. <br/> Abstract. The present study is an extension of the work done by Peddle, Haut, and Wingate [SIAM J. Sci. Comput., 41 (2019), pp. A3476–A3497] and Haut and Wingate [SIAM J. Sci. Comput., 36 (2014), pp. A693–A713], where a twolevel Parareal method with mapping and averaging is examined. The method proposed in this paper is a multilevel Parareal method with arbitrarily many levels, which is not restricted to the twolevel case. We give an asymptotic error estimate which reduces to the twolevel estimate for the case when only two levels are considered. Introducing more than two levels has important consequences for the averaging procedure, as we choose separate averaging windows for each of the different levels, which is an additional new feature of the present study. The different averaging windows make the proposed method especially appropriate for nonlinear multiscale problems, because we can introduce a level for each intrinsic scale of the problem and adapt the averaging procedure such that we reproduce the behavior of the model on the particular scale resolved by the level. The method is applied to nonlinear differential equations. The nonlinearities can generate a range of frequencies in the problem. The computational cost of the new method is investigated and studied on several examples.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2709A2736, August 2024. <br/> Abstract. The present study is an extension of the work done by Peddle, Haut, and Wingate [SIAM J. Sci. Comput., 41 (2019), pp. A3476–A3497] and Haut and Wingate [SIAM J. Sci. Comput., 36 (2014), pp. A693–A713], where a twolevel Parareal method with mapping and averaging is examined. The method proposed in this paper is a multilevel Parareal method with arbitrarily many levels, which is not restricted to the twolevel case. We give an asymptotic error estimate which reduces to the twolevel estimate for the case when only two levels are considered. Introducing more than two levels has important consequences for the averaging procedure, as we choose separate averaging windows for each of the different levels, which is an additional new feature of the present study. The different averaging windows make the proposed method especially appropriate for nonlinear multiscale problems, because we can introduce a level for each intrinsic scale of the problem and adapt the averaging procedure such that we reproduce the behavior of the model on the particular scale resolved by the level. The method is applied to nonlinear differential equations. The nonlinearities can generate a range of frequencies in the problem. The computational cost of the new method is investigated and studied on several examples. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Multilevel Parareal Algorithm with Averaging for Oscillatory Problems
10.1137/23M1547123
SIAM Journal on Scientific Computing
20240820T07:00:00Z
© 2024 SIAM. Published by SIAM under the terms of the Creative Commons 4.0 license
Juliane Rosemeier
Terry Haut
Beth Wingate
Multilevel Parareal Algorithm with Averaging for Oscillatory Problems
46
4
A2709
A2736
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1547123
https://epubs.siam.org/doi/abs/10.1137/23M1547123?af=R
© 2024 SIAM. Published by SIAM under the terms of the Creative Commons 4.0 license

A Fast Iterative PDEBased Algorithm for Feedback Controls of Nonsmooth MeanField Control Problems
https://epubs.siam.org/doi/abs/10.1137/21M1441158?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2737A2773, August 2024. <br/> Abstract. We propose a PDEbased accelerated gradient algorithm for optimal feedback controls of McKean–Vlasov dynamics that involve meanfield interactions both in the state and action. The method exploits a forwardbackward splitting approach and iteratively refines the approximate controls based on the gradients of smooth costs, the proximal maps of nonsmooth costs, and dynamically updated momentum parameters. At each step, the state dynamics is approximated via a particle system, and the required gradient is evaluated through a coupled system of nonlocal linear PDEs. The latter is solved by finite difference approximation or neural networkbased residual approximation, depending on the state dimension. We present exhaustive numerical experiments for low and highdimensional meanfield control problems, including sparse stabilization of stochastic Cucker–Smale models, which reveal that our algorithm captures important structures of the optimal feedback control and achieves a robust performance with respect to parameter perturbation.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2737A2773, August 2024. <br/> Abstract. We propose a PDEbased accelerated gradient algorithm for optimal feedback controls of McKean–Vlasov dynamics that involve meanfield interactions both in the state and action. The method exploits a forwardbackward splitting approach and iteratively refines the approximate controls based on the gradients of smooth costs, the proximal maps of nonsmooth costs, and dynamically updated momentum parameters. At each step, the state dynamics is approximated via a particle system, and the required gradient is evaluated through a coupled system of nonlocal linear PDEs. The latter is solved by finite difference approximation or neural networkbased residual approximation, depending on the state dimension. We present exhaustive numerical experiments for low and highdimensional meanfield control problems, including sparse stabilization of stochastic Cucker–Smale models, which reveal that our algorithm captures important structures of the optimal feedback control and achieves a robust performance with respect to parameter perturbation. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A Fast Iterative PDEBased Algorithm for Feedback Controls of Nonsmooth MeanField Control Problems
10.1137/21M1441158
SIAM Journal on Scientific Computing
20240820T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Christoph Reisinger
Wolfgang Stockinger
Yufei Zhang
A Fast Iterative PDEBased Algorithm for Feedback Controls of Nonsmooth MeanField Control Problems
46
4
A2737
A2773
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/21M1441158
https://epubs.siam.org/doi/abs/10.1137/21M1441158?af=R
© 2024 Society for Industrial and Applied Mathematics

A SketchandSelect Arnoldi Process
https://epubs.siam.org/doi/abs/10.1137/23M1588007?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page A2774A2797, August 2024. <br/> Abstract. A sketchandselect Arnoldi process to generate a wellconditioned basis of a Krylov space at low cost is proposed. At each iteration the procedure utilizes randomized sketching to select a limited number of previously computed basis vectors to project out of the current basis vector. The computational cost grows linearly with the dimension of the Krylov space. The subset selection problem for the projection step is approximately solved with a number of heuristic algorithms and greedy methods used in statistical learning and compressive sensing. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/simunec/sketchselectarnoldi and in the supplementary materials (sketchselectarnoldimain.zip [2.21MB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2774A2797, August 2024. <br/> Abstract. A sketchandselect Arnoldi process to generate a wellconditioned basis of a Krylov space at low cost is proposed. At each iteration the procedure utilizes randomized sketching to select a limited number of previously computed basis vectors to project out of the current basis vector. The computational cost grows linearly with the dimension of the Krylov space. The subset selection problem for the projection step is approximately solved with a number of heuristic algorithms and greedy methods used in statistical learning and compressive sensing. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/simunec/sketchselectarnoldi and in the supplementary materials (sketchselectarnoldimain.zip [2.21MB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A SketchandSelect Arnoldi Process
10.1137/23M1588007
SIAM Journal on Scientific Computing
20240821T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Stefan Güttel
Igor Simunec
A SketchandSelect Arnoldi Process
46
4
A2774
A2797
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1588007
https://epubs.siam.org/doi/abs/10.1137/23M1588007?af=R
© 2024 Society for Industrial and Applied Mathematics

Shape Optimization of Optical Microscale Inclusions
https://epubs.siam.org/doi/abs/10.1137/23M158262X?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page B377B402, August 2024. <br/> Abstract. This paper describes a class of shape optimization problems for optical metamaterials comprised of periodic microscale inclusions composed of a dielectric, lowdimensional material suspended in a nonmagnetic bulk dielectric. The shape optimization approach is based on a homogenization theory for timeharmonic Maxwell’s equations that describes effective material parameters for the propagation of electromagnetic waves through the metamaterial. The control parameter of the optimization is a deformation field representing the deviation of the microscale geometry from a reference configuration of the cell problem. This allows for describing the homogenized effective permittivity tensor as a function of the deformation field. We show that the underlying deformed cell problem is wellposed and regular. This, in turn, proves that the shape optimization problem is wellposed. In addition, a numerical scheme is formulated that utilizes an adjoint formulation with either gradient descent or BFGS as optimization algorithms. The developed algorithm is tested numerically on a number of prototypical shape optimization problems with a prescribed effective permittivity tensor as the target. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://zenodo.org/records/10459309.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page B377B402, August 2024. <br/> Abstract. This paper describes a class of shape optimization problems for optical metamaterials comprised of periodic microscale inclusions composed of a dielectric, lowdimensional material suspended in a nonmagnetic bulk dielectric. The shape optimization approach is based on a homogenization theory for timeharmonic Maxwell’s equations that describes effective material parameters for the propagation of electromagnetic waves through the metamaterial. The control parameter of the optimization is a deformation field representing the deviation of the microscale geometry from a reference configuration of the cell problem. This allows for describing the homogenized effective permittivity tensor as a function of the deformation field. We show that the underlying deformed cell problem is wellposed and regular. This, in turn, proves that the shape optimization problem is wellposed. In addition, a numerical scheme is formulated that utilizes an adjoint formulation with either gradient descent or BFGS as optimization algorithms. The developed algorithm is tested numerically on a number of prototypical shape optimization problems with a prescribed effective permittivity tensor as the target. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://zenodo.org/records/10459309. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Shape Optimization of Optical Microscale Inclusions
10.1137/23M158262X
SIAM Journal on Scientific Computing
20240703T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Manaswinee Bezbaruah
Matthias Maier
Winnifried Wollner
Shape Optimization of Optical Microscale Inclusions
46
4
B377
B402
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M158262X
https://epubs.siam.org/doi/abs/10.1137/23M158262X?af=R
© 2024 Society for Industrial and Applied Mathematics

Quantum Simulation for Quantum Dynamics with Artificial Boundary Conditions
https://epubs.siam.org/doi/abs/10.1137/23M1563451?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page B403B421, August 2024. <br/> Abstract. Quantum dynamics, typically expressed in the form of a timedependent Schrödinger equation with a Hermitian Hamiltonian, is a natural application for quantum computing. However, when simulating quantum dynamics that involves the emission of electrons, it is necessary to use artificial boundary conditions (ABCs) to confine the computation within a fixed domain. The introduction of ABCs alters the Hamiltonian structure of the dynamics, and existing quantum algorithms cannot be directly applied since the evolution is no longer unitary. The current paper utilizes a recently introduced Schrödingerization method that converts nonHermitian dynamics into a Schrödinger form for the artificial boundary problems [S. Jin, N. Liu, and Y. Yu, Quantum Simulation of Partial Differential Equations via Schrödingerisation, preprint, arXiv:2212.13969, 2022], [S. Jin, N. Liu, and Y. Yu, Phys. Rev. A, 108 (2023), 032603]. We implement this method for three types of ABCs, including the complex absorbing potential technique, perfectly matched layer methods, and DirichlettoNeumann approach. We analyze the query complexity of these algorithms and perform numerical experiments to demonstrate the validity of this approach. This helps to bridge the gap between available quantum algorithms and computational models for quantum dynamics in unbounded domains.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page B403B421, August 2024. <br/> Abstract. Quantum dynamics, typically expressed in the form of a timedependent Schrödinger equation with a Hermitian Hamiltonian, is a natural application for quantum computing. However, when simulating quantum dynamics that involves the emission of electrons, it is necessary to use artificial boundary conditions (ABCs) to confine the computation within a fixed domain. The introduction of ABCs alters the Hamiltonian structure of the dynamics, and existing quantum algorithms cannot be directly applied since the evolution is no longer unitary. The current paper utilizes a recently introduced Schrödingerization method that converts nonHermitian dynamics into a Schrödinger form for the artificial boundary problems [S. Jin, N. Liu, and Y. Yu, Quantum Simulation of Partial Differential Equations via Schrödingerisation, preprint, arXiv:2212.13969, 2022], [S. Jin, N. Liu, and Y. Yu, Phys. Rev. A, 108 (2023), 032603]. We implement this method for three types of ABCs, including the complex absorbing potential technique, perfectly matched layer methods, and DirichlettoNeumann approach. We analyze the query complexity of these algorithms and perform numerical experiments to demonstrate the validity of this approach. This helps to bridge the gap between available quantum algorithms and computational models for quantum dynamics in unbounded domains. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Quantum Simulation for Quantum Dynamics with Artificial Boundary Conditions
10.1137/23M1563451
SIAM Journal on Scientific Computing
20240710T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Shi Jin
Xiantao Li
Nana Liu
Yue Yu
Quantum Simulation for Quantum Dynamics with Artificial Boundary Conditions
46
4
B403
B421
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1563451
https://epubs.siam.org/doi/abs/10.1137/23M1563451?af=R
© 2024 Society for Industrial and Applied Mathematics

Bayesian Parameter Identification in Impedance Boundary Conditions for Helmholtz Problems
https://epubs.siam.org/doi/abs/10.1137/23M1591517?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page B422B447, August 2024. <br/> Abstract. We consider the problem of identifying the acoustic impedance of a wall surface from noisy pressure measurements in a closed room using a Bayesian approach. The room acoustics are modeled by the interior Helmholtz equation with impedance boundary conditions. The aim is to compute moments of the acoustic impedance to estimate a suitable density function of the impedance coefficient. For the computation of moments we use ratio estimators and Monte Carlo sampling. We consider two different experimental scenarios. In the first scenario, the noisy measurements correspond to a wall modeled by impedance boundary conditions. In this case, the Bayesian algorithm uses a model that is (up to the noise) consistent with the measurements and our algorithm is able to identify acoustic impedance with high accuracy. In the second scenario, the noisy measurements come from a coupled acousticstructural problem, modeling a wall made of glass, whereas the Bayesian algorithm still uses a model with impedance boundary conditions. In this case, the parameter identification model is inconsistent with the measurements and therefore is not capable to represent them well. Nonetheless, for particular frequency bands the Bayesian algorithm identifies estimates with high likelihood. Outside these frequency bands the algorithm fails. We discuss the results of both examples and possible reasons for the failure of the latter case for particular frequency values.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page B422B447, August 2024. <br/> Abstract. We consider the problem of identifying the acoustic impedance of a wall surface from noisy pressure measurements in a closed room using a Bayesian approach. The room acoustics are modeled by the interior Helmholtz equation with impedance boundary conditions. The aim is to compute moments of the acoustic impedance to estimate a suitable density function of the impedance coefficient. For the computation of moments we use ratio estimators and Monte Carlo sampling. We consider two different experimental scenarios. In the first scenario, the noisy measurements correspond to a wall modeled by impedance boundary conditions. In this case, the Bayesian algorithm uses a model that is (up to the noise) consistent with the measurements and our algorithm is able to identify acoustic impedance with high accuracy. In the second scenario, the noisy measurements come from a coupled acousticstructural problem, modeling a wall made of glass, whereas the Bayesian algorithm still uses a model with impedance boundary conditions. In this case, the parameter identification model is inconsistent with the measurements and therefore is not capable to represent them well. Nonetheless, for particular frequency bands the Bayesian algorithm identifies estimates with high likelihood. Outside these frequency bands the algorithm fails. We discuss the results of both examples and possible reasons for the failure of the latter case for particular frequency values. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Bayesian Parameter Identification in Impedance Boundary Conditions for Helmholtz Problems
10.1137/23M1591517
SIAM Journal on Scientific Computing
20240710T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Nick Wulbusch
Reinhild Roden
Matthias Blau
Alexey Chernov
Bayesian Parameter Identification in Impedance Boundary Conditions for Helmholtz Problems
46
4
B422
B447
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1591517
https://epubs.siam.org/doi/abs/10.1137/23M1591517?af=R
© 2024 Society for Industrial and Applied Mathematics

Integrated Nested Laplace Approximations for LargeScale Spatiotemporal Bayesian Modeling
https://epubs.siam.org/doi/abs/10.1137/23M1561531?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page B448B473, August 2024. <br/> Abstract. Bayesian inference tasks continue to pose a computational challenge. This especially holds for spatiotemporal modeling, where highdimensional latent parameter spaces are ubiquitous. The methodology of integrated nested Laplace approximations (INLA) provides a framework for performing Bayesian inference applicable to a large subclass of additive Bayesian hierarchical models. In combination with the stochastic partial differential equation (SPDE) approach, it gives rise to an efficient method for spatiotemporal modeling. In this work, we build on the INLASPDE approach by putting forward a performant distributed memory variant, INLADIST, for largescale applications. To perform the arising computational kernel operations, consisting of Cholesky factorizations, solving linear systems, and selected matrix inversions, we present two numerical solver options: a sparse CPUbased library and a novel blocked GPUaccelerated approach which we propose. We leverage the recurring nonzero block structure in the arising precision (inverse covariance) matrices, which allows us to employ dense subroutines within a sparse setting. Both versions of INLADIST are highly scalable, capable of performing inference on models with millions of latent parameters. We demonstrate their accuracy and performance on synthetic as well as realworld climate dataset applications.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page B448B473, August 2024. <br/> Abstract. Bayesian inference tasks continue to pose a computational challenge. This especially holds for spatiotemporal modeling, where highdimensional latent parameter spaces are ubiquitous. The methodology of integrated nested Laplace approximations (INLA) provides a framework for performing Bayesian inference applicable to a large subclass of additive Bayesian hierarchical models. In combination with the stochastic partial differential equation (SPDE) approach, it gives rise to an efficient method for spatiotemporal modeling. In this work, we build on the INLASPDE approach by putting forward a performant distributed memory variant, INLADIST, for largescale applications. To perform the arising computational kernel operations, consisting of Cholesky factorizations, solving linear systems, and selected matrix inversions, we present two numerical solver options: a sparse CPUbased library and a novel blocked GPUaccelerated approach which we propose. We leverage the recurring nonzero block structure in the arising precision (inverse covariance) matrices, which allows us to employ dense subroutines within a sparse setting. Both versions of INLADIST are highly scalable, capable of performing inference on models with millions of latent parameters. We demonstrate their accuracy and performance on synthetic as well as realworld climate dataset applications. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Integrated Nested Laplace Approximations for LargeScale Spatiotemporal Bayesian Modeling
10.1137/23M1561531
SIAM Journal on Scientific Computing
20240712T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Lisa GaedkeMerzhäuser
Elias Krainski
Radim Janalik
Håvard Rue
Olaf Schenk
Integrated Nested Laplace Approximations for LargeScale Spatiotemporal Bayesian Modeling
46
4
B448
B473
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1561531
https://epubs.siam.org/doi/abs/10.1137/23M1561531?af=R
© 2024 Society for Industrial and Applied Mathematics

SOPT: A Points Selection Algorithm for HyperReduction in Reduced Order Models
https://epubs.siam.org/doi/abs/10.1137/22M1484018?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page B474B501, August 2024. <br/> Abstract. While projectionbased reduced order models can reduce the dimension of full order solutions, the resulting reduced models may still contain terms that scale with the full order dimension. Hyperreduction techniques are samplingbased methods that further reduce this computational complexity by approximating such terms with a much smaller dimension. The goal of this work is to introduce the points selection algorithm developed by Shin and Xiu [SIAM J. Sci. Comput., 38 (2016), pp. A385–A411] as a hyperreduction method. The selection algorithm was originally proposed as a stochastic collocation method for uncertainty quantification. Since the algorithm aims at maximizing a quantity [math] that measures both the column orthogonality and the determinant, we refer to the algorithm as SOPT. Numerical examples are provided to demonstrate the performance of SOPT and to compare its performance with a gappy proper orthogonal decomposition (POD) algorithm. We found that using the SOPT algorithm is shown to predict the full order solutions with higher accuracy than gappy POD especially when the number of sampling points is small, although we note that SOPT shows slow asymptotic convergence with respect to the number of samples for some applications, e.g., Lagrangian hydrodynamics.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page B474B501, August 2024. <br/> Abstract. While projectionbased reduced order models can reduce the dimension of full order solutions, the resulting reduced models may still contain terms that scale with the full order dimension. Hyperreduction techniques are samplingbased methods that further reduce this computational complexity by approximating such terms with a much smaller dimension. The goal of this work is to introduce the points selection algorithm developed by Shin and Xiu [SIAM J. Sci. Comput., 38 (2016), pp. A385–A411] as a hyperreduction method. The selection algorithm was originally proposed as a stochastic collocation method for uncertainty quantification. Since the algorithm aims at maximizing a quantity [math] that measures both the column orthogonality and the determinant, we refer to the algorithm as SOPT. Numerical examples are provided to demonstrate the performance of SOPT and to compare its performance with a gappy proper orthogonal decomposition (POD) algorithm. We found that using the SOPT algorithm is shown to predict the full order solutions with higher accuracy than gappy POD especially when the number of sampling points is small, although we note that SOPT shows slow asymptotic convergence with respect to the number of samples for some applications, e.g., Lagrangian hydrodynamics. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
SOPT: A Points Selection Algorithm for HyperReduction in Reduced Order Models
10.1137/22M1484018
SIAM Journal on Scientific Computing
20240716T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Jessica T. Lauzon
Siu Wun Cheung
Yeonjong Shin
Youngsoo Choi
Dylan M. Copeland
Kevin Huynh
SOPT: A Points Selection Algorithm for HyperReduction in Reduced Order Models
46
4
B474
B501
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/22M1484018
https://epubs.siam.org/doi/abs/10.1137/22M1484018?af=R
© 2024 Society for Industrial and Applied Mathematics

Rounding Error Using Low Precision Approximate Random Variables
https://epubs.siam.org/doi/abs/10.1137/23M1552814?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page B502B526, August 2024. <br/> Abstract. For numerical approximations to stochastic differential equations using the Euler–Maruyama scheme, we propose incorporating approximate random variables computed using low precisions, such as single and half precision. We propose and justify a model for the rounding error incurred, and produce an average case error bound for two and four way differences, appropriate for regular and nested multilevel Monte Carlo estimations. Our rounding error model recovers and extends the statistical model by Arciniega and Allen [Stoch. Anal. Appl., 21 (2003), pp. 281–300], while bounding the size that systematic and biased rounding errors are permitted to be. By considering the variance structure of multilevel Monte Carlo correction terms in various precisions with and without a Kahan compensated summation, we compute the potential speed ups offered from the various precisions. We find single precision offers the potential for approximate speed improvements by a factor of 7 across a wide span of discretization levels. Half precision offers comparable improvements for several levels of coarse simulations, and even offers improvements by a factor of 10–12 for the very coarsest few levels, which are likely to dominate higher order methods such as the Milstein scheme.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page B502B526, August 2024. <br/> Abstract. For numerical approximations to stochastic differential equations using the Euler–Maruyama scheme, we propose incorporating approximate random variables computed using low precisions, such as single and half precision. We propose and justify a model for the rounding error incurred, and produce an average case error bound for two and four way differences, appropriate for regular and nested multilevel Monte Carlo estimations. Our rounding error model recovers and extends the statistical model by Arciniega and Allen [Stoch. Anal. Appl., 21 (2003), pp. 281–300], while bounding the size that systematic and biased rounding errors are permitted to be. By considering the variance structure of multilevel Monte Carlo correction terms in various precisions with and without a Kahan compensated summation, we compute the potential speed ups offered from the various precisions. We find single precision offers the potential for approximate speed improvements by a factor of 7 across a wide span of discretization levels. Half precision offers comparable improvements for several levels of coarse simulations, and even offers improvements by a factor of 10–12 for the very coarsest few levels, which are likely to dominate higher order methods such as the Milstein scheme. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Rounding Error Using Low Precision Approximate Random Variables
10.1137/23M1552814
SIAM Journal on Scientific Computing
20240716T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Michael B. Giles
Oliver SheridanMethven
Rounding Error Using Low Precision Approximate Random Variables
46
4
B502
B526
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1552814
https://epubs.siam.org/doi/abs/10.1137/23M1552814?af=R
© 2024 Society for Industrial and Applied Mathematics

Cut Finite Element Discretizations of CellbyCell EMI Electrophysiology Models
https://epubs.siam.org/doi/abs/10.1137/23M1580632?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page B527B553, August 2024. <br/> Abstract. The EMI (extracellularmembraneintracellular) model describes electrical activity in excitable tissue, where the extracellular and intracellular spaces and cellular membrane are explicitly represented. The model couples a system of partial differential equations (PDEs) in the intracellular and extracellular spaces with a system of ordinary differential equations (ODEs) on the membrane. A key challenge for the EMI model is the generation of highquality meshes conforming to the complex geometries of brain cells. To overcome this challenge, we propose a novel cut finite element method (CutFEM) where the membrane geometry can be represented independently of a structured and easytogenerate background mesh for the remaining computational domain. Starting from a Godunov splitting scheme, the EMI model is split into separate PDE and ODE parts. The resulting PDE part is a nonstandard elliptic interface problem, for which we devise two different CutFEM formulations: one singledimensional formulation with the intra/extracellular electrical potentials as unknowns, and a multidimensional formulation that also introduces the electrical current over the membrane as an additional unknown leading to a penalized saddle point problem. Both formulations are augmented by suitably designed ghost penalties to ensure stability and convergence properties that are insensitive to how the membrane surface mesh cuts the background mesh. For the ODE part, we introduce a new unfitted discretization to solve the membrane bound ODEs on a membrane interface that is not aligned with the background mesh. Finally, we perform extensive numerical experiments to demonstrate that CutFEM is a promising approach to efficiently simulate electrical activity in geometrically resolved brain cells. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://zenodo.org/record/8068506.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page B527B553, August 2024. <br/> Abstract. The EMI (extracellularmembraneintracellular) model describes electrical activity in excitable tissue, where the extracellular and intracellular spaces and cellular membrane are explicitly represented. The model couples a system of partial differential equations (PDEs) in the intracellular and extracellular spaces with a system of ordinary differential equations (ODEs) on the membrane. A key challenge for the EMI model is the generation of highquality meshes conforming to the complex geometries of brain cells. To overcome this challenge, we propose a novel cut finite element method (CutFEM) where the membrane geometry can be represented independently of a structured and easytogenerate background mesh for the remaining computational domain. Starting from a Godunov splitting scheme, the EMI model is split into separate PDE and ODE parts. The resulting PDE part is a nonstandard elliptic interface problem, for which we devise two different CutFEM formulations: one singledimensional formulation with the intra/extracellular electrical potentials as unknowns, and a multidimensional formulation that also introduces the electrical current over the membrane as an additional unknown leading to a penalized saddle point problem. Both formulations are augmented by suitably designed ghost penalties to ensure stability and convergence properties that are insensitive to how the membrane surface mesh cuts the background mesh. For the ODE part, we introduce a new unfitted discretization to solve the membrane bound ODEs on a membrane interface that is not aligned with the background mesh. Finally, we perform extensive numerical experiments to demonstrate that CutFEM is a promising approach to efficiently simulate electrical activity in geometrically resolved brain cells. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://zenodo.org/record/8068506. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Cut Finite Element Discretizations of CellbyCell EMI Electrophysiology Models
10.1137/23M1580632
SIAM Journal on Scientific Computing
20240806T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Nanna Berre
Marie E. Rognes
André Massing
Cut Finite Element Discretizations of CellbyCell EMI Electrophysiology Models
46
4
B527
B553
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1580632
https://epubs.siam.org/doi/abs/10.1137/23M1580632?af=R
© 2024 Society for Industrial and Applied Mathematics

Implicit Adaptive Mesh Refinement for Dispersive Tsunami Propagation
https://epubs.siam.org/doi/abs/10.1137/23M1585210?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page B554B578, August 2024. <br/> Abstract. We present an algorithm to solve the dispersive depthaveraged Serre–Green–Naghdi equations using patchbased adaptive mesh refinement. These equations require adding additional higher derivative terms to the nonlinear shallow water equations. This has been implemented as a new component of the open source GeoClaw software that is widely used for modeling tsunamis, storm surge, and related hazards, improving its accuracy on shorter wavelength phenomena. We use a formulation that requires solving an elliptic system of equations at each time step, making the method implicit. The adaptive algorithm allows different time steps on different refinement levels and solves the implicit equations level by level. Computational examples are presented to illustrate the stability and accuracy on a radially symmetric test case and two realistic tsunami modeling problems, including a hypothetical asteroid impact creating a short wavelength tsunami for which dispersive terms are necessary. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/rjleveque/ImplicitAMRpaper and in the supplementary materials (ImplicitAMRpaper.zip [174KB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page B554B578, August 2024. <br/> Abstract. We present an algorithm to solve the dispersive depthaveraged Serre–Green–Naghdi equations using patchbased adaptive mesh refinement. These equations require adding additional higher derivative terms to the nonlinear shallow water equations. This has been implemented as a new component of the open source GeoClaw software that is widely used for modeling tsunamis, storm surge, and related hazards, improving its accuracy on shorter wavelength phenomena. We use a formulation that requires solving an elliptic system of equations at each time step, making the method implicit. The adaptive algorithm allows different time steps on different refinement levels and solves the implicit equations level by level. Computational examples are presented to illustrate the stability and accuracy on a radially symmetric test case and two realistic tsunami modeling problems, including a hypothetical asteroid impact creating a short wavelength tsunami for which dispersive terms are necessary. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/rjleveque/ImplicitAMRpaper and in the supplementary materials (ImplicitAMRpaper.zip [174KB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Implicit Adaptive Mesh Refinement for Dispersive Tsunami Propagation
10.1137/23M1585210
SIAM Journal on Scientific Computing
20240814T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Marsha J. Berger
Randall J. LeVeque
Implicit Adaptive Mesh Refinement for Dispersive Tsunami Propagation
46
4
B554
B578
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1585210
https://epubs.siam.org/doi/abs/10.1137/23M1585210?af=R
© 2024 Society for Industrial and Applied Mathematics

On the Training and Generalization of Deep Operator Networks
https://epubs.siam.org/doi/abs/10.1137/23M1598751?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page C273C296, August 2024. <br/> Abstract. We present a novel training method for deep operator networks (DeepONets), one of the most popular neural network models for operators. DeepONets are constructed by two subnetworks, namely the branch and trunk networks. Typically, the two subnetworks are trained simultaneously, which amounts to solving a complex optimization problem in a high dimensional space. In addition, the nonconvex and nonlinear nature makes training very challenging. To tackle such a challenge, we propose a twostep training method that trains the trunk network first and then sequentially trains the branch network. The core mechanism is motivated by the divideandconquer paradigm and is the decomposition of the entire complex training task into two subtasks with reduced complexity. Therein the Gram–Schmidt orthonormalization process is introduced which significantly improves stability and generalization ability. On the theoretical side, we establish a generalization error estimate in terms of the number of training data, the width of DeepONets, and the number of input and output sensors. Numerical examples are presented to demonstrate the effectiveness of the twostep training method, including Darcy flow in heterogeneous porous media.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C273C296, August 2024. <br/> Abstract. We present a novel training method for deep operator networks (DeepONets), one of the most popular neural network models for operators. DeepONets are constructed by two subnetworks, namely the branch and trunk networks. Typically, the two subnetworks are trained simultaneously, which amounts to solving a complex optimization problem in a high dimensional space. In addition, the nonconvex and nonlinear nature makes training very challenging. To tackle such a challenge, we propose a twostep training method that trains the trunk network first and then sequentially trains the branch network. The core mechanism is motivated by the divideandconquer paradigm and is the decomposition of the entire complex training task into two subtasks with reduced complexity. Therein the Gram–Schmidt orthonormalization process is introduced which significantly improves stability and generalization ability. On the theoretical side, we establish a generalization error estimate in terms of the number of training data, the width of DeepONets, and the number of input and output sensors. Numerical examples are presented to demonstrate the effectiveness of the twostep training method, including Darcy flow in heterogeneous porous media. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
On the Training and Generalization of Deep Operator Networks
10.1137/23M1598751
SIAM Journal on Scientific Computing
20240701T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Sanghyun Lee
Yeonjong Shin
On the Training and Generalization of Deep Operator Networks
46
4
C273
C296
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1598751
https://epubs.siam.org/doi/abs/10.1137/23M1598751?af=R
© 2024 Society for Industrial and Applied Mathematics

Slow Invariant Manifolds of Singularly Perturbed Systems via PhysicsInformed Machine Learning
https://epubs.siam.org/doi/abs/10.1137/23M1602991?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page C297C322, August 2024. <br/> Abstract. We present a physicsinformed machine learning (PIML) approach for the approximation of slow invariant manifolds of singularly perturbed systems, providing functionals in an explicit form that facilitates the construction and numerical integration of reducedorder models (ROMs). The proposed scheme solves the partial differential equation corresponding to the invariance equation (IE) within the geometric singular perturbation theory (GSPT) framework. For the solution of the IE, we used two neural network structures, namely, feedforward neural networks and random projection neural networks, with symbolic differentiation for the computation of the gradients required for the learning process. The efficiency of our PIML method is assessed via three benchmark problems, namely, the Michaelis–Menten, the targetmediated drug disposition reaction mechanism, and the 3D Sel’kov model. We show that the proposed PIML scheme provides approximations of equivalent or even higher accuracy than those provided by other traditional GSPTbased methods, and importantly, for any practical purposes, it is not affected by the magnitude of the perturbation parameter. This is of particular importance because there are many systems for which the gap between the fast and slow timescales is not that big, but still, ROMs can be constructed. A comparison of the computational costs between symbolic, automatic, and numerical approximation of the required derivatives in the learning process is also provided. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://epubs.siam.org/doi/suppl/10.1137/23M1602991/suppl_file/131735_1_supp_551502_s5k7wy_sc.pdf and https://epubs.siam.org/doi/suppl/10.1137/23M1602991/suppl_file/SISC_PIML_SIMs_SPmain.zip.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C297C322, August 2024. <br/> Abstract. We present a physicsinformed machine learning (PIML) approach for the approximation of slow invariant manifolds of singularly perturbed systems, providing functionals in an explicit form that facilitates the construction and numerical integration of reducedorder models (ROMs). The proposed scheme solves the partial differential equation corresponding to the invariance equation (IE) within the geometric singular perturbation theory (GSPT) framework. For the solution of the IE, we used two neural network structures, namely, feedforward neural networks and random projection neural networks, with symbolic differentiation for the computation of the gradients required for the learning process. The efficiency of our PIML method is assessed via three benchmark problems, namely, the Michaelis–Menten, the targetmediated drug disposition reaction mechanism, and the 3D Sel’kov model. We show that the proposed PIML scheme provides approximations of equivalent or even higher accuracy than those provided by other traditional GSPTbased methods, and importantly, for any practical purposes, it is not affected by the magnitude of the perturbation parameter. This is of particular importance because there are many systems for which the gap between the fast and slow timescales is not that big, but still, ROMs can be constructed. A comparison of the computational costs between symbolic, automatic, and numerical approximation of the required derivatives in the learning process is also provided. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://epubs.siam.org/doi/suppl/10.1137/23M1602991/suppl_file/131735_1_supp_551502_s5k7wy_sc.pdf and https://epubs.siam.org/doi/suppl/10.1137/23M1602991/suppl_file/SISC_PIML_SIMs_SPmain.zip. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Slow Invariant Manifolds of Singularly Perturbed Systems via PhysicsInformed Machine Learning
10.1137/23M1602991
SIAM Journal on Scientific Computing
20240705T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Dimitrios Patsatzis
Gianluca Fabiani
Lucia Russo
Constantinos Siettos
Slow Invariant Manifolds of Singularly Perturbed Systems via PhysicsInformed Machine Learning
46
4
C297
C322
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1602991
https://epubs.siam.org/doi/abs/10.1137/23M1602991?af=R
© 2024 Society for Industrial and Applied Mathematics

Render unto Numerics: Orthogonal Polynomial Neural Operator for PDEs with Nonperiodic Boundary Conditions
https://epubs.siam.org/doi/abs/10.1137/23M1556320?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page C323C348, August 2024. <br/> Abstract. By learning the mappings between infinite function spaces using carefully designed neural networks, the operator learning methodology has exhibited significantly more efficiency than traditional methods in solving differential equations, but faces concerns about their accuracy and reliability. To overcome these limitations through robustly enforcing boundary conditions (BCs), a general neural architecture named spectral operator learning is introduced by combining with the structures of the spectral numerical method. One variant called the orthogonal polynomial neural operator (OPNO) is proposed later, aiming at PDEs with Dirichlet, Neumann, and Robin BCs. The strict BC satisfaction properties and the universal approximation capacity of the OPNO are theoretically proven. A variety of numerical experiments with physical backgrounds demonstrate that the OPNO outperforms other existing deep learning methodologies, showcasing potential of comparable accuracy with the traditional secondorder finite difference method that employs a considerably fine mesh (with relative errors on the order of [math]), and is up to almost 5 magnitudes faster over the traditional method. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/liuziyuanmath/spectral_operator_learning/tree/main/OPNO/Reproduce and in the supplementary materials (spectral_operator_learningmain.zip [669KB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C323C348, August 2024. <br/> Abstract. By learning the mappings between infinite function spaces using carefully designed neural networks, the operator learning methodology has exhibited significantly more efficiency than traditional methods in solving differential equations, but faces concerns about their accuracy and reliability. To overcome these limitations through robustly enforcing boundary conditions (BCs), a general neural architecture named spectral operator learning is introduced by combining with the structures of the spectral numerical method. One variant called the orthogonal polynomial neural operator (OPNO) is proposed later, aiming at PDEs with Dirichlet, Neumann, and Robin BCs. The strict BC satisfaction properties and the universal approximation capacity of the OPNO are theoretically proven. A variety of numerical experiments with physical backgrounds demonstrate that the OPNO outperforms other existing deep learning methodologies, showcasing potential of comparable accuracy with the traditional secondorder finite difference method that employs a considerably fine mesh (with relative errors on the order of [math]), and is up to almost 5 magnitudes faster over the traditional method. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/liuziyuanmath/spectral_operator_learning/tree/main/OPNO/Reproduce and in the supplementary materials (spectral_operator_learningmain.zip [669KB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Render unto Numerics: Orthogonal Polynomial Neural Operator for PDEs with Nonperiodic Boundary Conditions
10.1137/23M1556320
SIAM Journal on Scientific Computing
20240712T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Ziyuan Liu
Haifeng Wang
Hong Zhang
Kaijun Bao
Xu Qian
Songhe Song
Render unto Numerics: Orthogonal Polynomial Neural Operator for PDEs with Nonperiodic Boundary Conditions
46
4
C323
C348
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1556320
https://epubs.siam.org/doi/abs/10.1137/23M1556320?af=R
© 2024 Society for Industrial and Applied Mathematics

Temporal Difference Learning for HighDimensional PIDEs with Jumps
https://epubs.siam.org/doi/abs/10.1137/23M1584538?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page C349C368, August 2024. <br/> Abstract. In this paper, we propose a deep learning framework for solving highdimensional partial integrodifferential equations (PIDEs) based on the temporal difference learning. We introduce a set of Lévy processes and construct a corresponding reinforcement learning model. To simulate the entire process, we use deep neural networks to represent the solutions and nonlocal terms of the equations. Subsequently, we train the networks using the temporal difference error, the termination condition, and properties of the nonlocal terms as the loss function. The relative error of the method reaches [math] in 100dimensional experiments and [math] in onedimensional pure jump problems. Additionally, our method demonstrates the advantages of low computational cost and robustness, making it wellsuited for addressing problems with different forms and intensities of jumps.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C349C368, August 2024. <br/> Abstract. In this paper, we propose a deep learning framework for solving highdimensional partial integrodifferential equations (PIDEs) based on the temporal difference learning. We introduce a set of Lévy processes and construct a corresponding reinforcement learning model. To simulate the entire process, we use deep neural networks to represent the solutions and nonlocal terms of the equations. Subsequently, we train the networks using the temporal difference error, the termination condition, and properties of the nonlocal terms as the loss function. The relative error of the method reaches [math] in 100dimensional experiments and [math] in onedimensional pure jump problems. Additionally, our method demonstrates the advantages of low computational cost and robustness, making it wellsuited for addressing problems with different forms and intensities of jumps. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Temporal Difference Learning for HighDimensional PIDEs with Jumps
10.1137/23M1584538
SIAM Journal on Scientific Computing
20240718T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Liwei Lu
Hailong Guo
Xu Yang
Yi Zhu
Temporal Difference Learning for HighDimensional PIDEs with Jumps
46
4
C349
C368
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1584538
https://epubs.siam.org/doi/abs/10.1137/23M1584538?af=R
© 2024 Society for Industrial and Applied Mathematics

Solving Poisson Problems in Polygonal Domains with Singularity Enriched Physics Informed Neural Networks
https://epubs.siam.org/doi/abs/10.1137/23M1601195?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page C369C398, August 2024. <br/> Abstract. Physicsinformed neural networks (PINNs) are a powerful class of numerical solvers for partial differential equations, employing deep neural networks with successful applications across a diverse set of problems. However, their effectiveness is somewhat diminished when addressing issues involving singularities, such as point sources or geometric irregularities, where the approximations they provide often suffer from reduced accuracy due to the limited regularity of the exact solution. In this work, we investigate PINNs for solving Poisson equations in polygonal domains with geometric singularities and mixed boundary conditions. We propose a novel singularity enriched PINN, by explicitly incorporating the singularity behavior of the analytic solution, e.g., corner singularity, mixed boundary condition, and edge singularities, into the ansatz space, and present a convergence analysis of the scheme. We present extensive numerical simulations in two and three dimensions to illustrate the efficiency of the method, and also a comparative study with several existing neural network based approaches. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/hhjcweb/SEPINN.git and in the supplementary materials (M160119_SuppMat.pdf [399KB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C369C398, August 2024. <br/> Abstract. Physicsinformed neural networks (PINNs) are a powerful class of numerical solvers for partial differential equations, employing deep neural networks with successful applications across a diverse set of problems. However, their effectiveness is somewhat diminished when addressing issues involving singularities, such as point sources or geometric irregularities, where the approximations they provide often suffer from reduced accuracy due to the limited regularity of the exact solution. In this work, we investigate PINNs for solving Poisson equations in polygonal domains with geometric singularities and mixed boundary conditions. We propose a novel singularity enriched PINN, by explicitly incorporating the singularity behavior of the analytic solution, e.g., corner singularity, mixed boundary condition, and edge singularities, into the ansatz space, and present a convergence analysis of the scheme. We present extensive numerical simulations in two and three dimensions to illustrate the efficiency of the method, and also a comparative study with several existing neural network based approaches. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/hhjcweb/SEPINN.git and in the supplementary materials (M160119_SuppMat.pdf [399KB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Solving Poisson Problems in Polygonal Domains with Singularity Enriched Physics Informed Neural Networks
10.1137/23M1601195
SIAM Journal on Scientific Computing
20240730T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Tianhao Hu
Bangti Jin
Zhi Zhou
Solving Poisson Problems in Polygonal Domains with Singularity Enriched Physics Informed Neural Networks
46
4
C369
C398
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1601195
https://epubs.siam.org/doi/abs/10.1137/23M1601195?af=R
© 2024 Society for Industrial and Applied Mathematics

Graph Neural Reaction Diffusion Models
https://epubs.siam.org/doi/abs/10.1137/23M1576700?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page C399C420, August 2024. <br/> Abstract. The integration of graph neural networks (GNNs) and neural ordinary and partial differential equations has been extensively studied in recent years. GNN architectures powered by neural differential equations allow us to reason about their behavior, and develop GNNs with desired properties such as controlled smoothing or energy conservation. In this paper we take inspiration from Turing instabilities in a reaction diffusion (RD) system of partial differential equations, and propose a novel family of GNNs based on neural RD systems, called RDGNN. We show that our RDGNN is powerful for the modeling of various data types, from homophilic, to heterophilic, and spatiotemporal datasets. We discuss the theoretical properties of our RDGNN, its implementation, and show that it improves or offers competitive performance to stateoftheart methods.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C399C420, August 2024. <br/> Abstract. The integration of graph neural networks (GNNs) and neural ordinary and partial differential equations has been extensively studied in recent years. GNN architectures powered by neural differential equations allow us to reason about their behavior, and develop GNNs with desired properties such as controlled smoothing or energy conservation. In this paper we take inspiration from Turing instabilities in a reaction diffusion (RD) system of partial differential equations, and propose a novel family of GNNs based on neural RD systems, called RDGNN. We show that our RDGNN is powerful for the modeling of various data types, from homophilic, to heterophilic, and spatiotemporal datasets. We discuss the theoretical properties of our RDGNN, its implementation, and show that it improves or offers competitive performance to stateoftheart methods. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Graph Neural Reaction Diffusion Models
10.1137/23M1576700
SIAM Journal on Scientific Computing
20240801T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Moshe Eliasof
Eldad Haber
Eran Treister
Graph Neural Reaction Diffusion Models
46
4
C399
C420
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1576700
https://epubs.siam.org/doi/abs/10.1137/23M1576700?af=R
© 2024 Society for Industrial and Applied Mathematics

An Alternating Flux Learning Method for Multidimensional Nonlinear Conservation Laws
https://epubs.siam.org/doi/abs/10.1137/23M1556605?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page C421C447, August 2024. <br/> Abstract. In a recent work [Q. Li and S. Evje, Netw. Heterog. Media, 18 (2023), pp. 48–79], it was explored how to identify the unknown flux function in a onedimensional scalar conservation law. Key ingredients are symbolic neural networks to represent the candidate flux functions, entropysatisfying numerical schemes, and a proper combination of initial data. The purpose of this work is to extend this methodology to a twodimensional scalar conservation law ([math]) [math]. Straightforward extension of the method from the 1D to the 2D problem results in poor identification of the unknown [math] and [math]. Relying on ideas from joint and alternating equations training, a learning strategy is designed that enables accurate identification of the flux functions, even when 2D observations are sparse. It involves an alternating flux training approach where a first set of candidate flux functions obtained from joint training is improved through an alternating directiondependent training strategy. Numerical investigations demonstrate that the method can effectively identify the true underlying flux functions [math] and [math] in the general case when they are nonconvex and unequal.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C421C447, August 2024. <br/> Abstract. In a recent work [Q. Li and S. Evje, Netw. Heterog. Media, 18 (2023), pp. 48–79], it was explored how to identify the unknown flux function in a onedimensional scalar conservation law. Key ingredients are symbolic neural networks to represent the candidate flux functions, entropysatisfying numerical schemes, and a proper combination of initial data. The purpose of this work is to extend this methodology to a twodimensional scalar conservation law ([math]) [math]. Straightforward extension of the method from the 1D to the 2D problem results in poor identification of the unknown [math] and [math]. Relying on ideas from joint and alternating equations training, a learning strategy is designed that enables accurate identification of the flux functions, even when 2D observations are sparse. It involves an alternating flux training approach where a first set of candidate flux functions obtained from joint training is improved through an alternating directiondependent training strategy. Numerical investigations demonstrate that the method can effectively identify the true underlying flux functions [math] and [math] in the general case when they are nonconvex and unequal. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
An Alternating Flux Learning Method for Multidimensional Nonlinear Conservation Laws
10.1137/23M1556605
SIAM Journal on Scientific Computing
20240801T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Qing Li
Steinar Evje
An Alternating Flux Learning Method for Multidimensional Nonlinear Conservation Laws
46
4
C421
C447
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1556605
https://epubs.siam.org/doi/abs/10.1137/23M1556605?af=R
© 2024 Society for Industrial and Applied Mathematics

LeastSquares Neural Network (LSNN) Method for Linear AdvectionReaction Equation: Discontinuity Interface
https://epubs.siam.org/doi/abs/10.1137/23M1568107?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page C448C478, August 2024. <br/> Abstract. We studied the leastsquares ReLU neural network (LSNN) method for solving a linear advectionreaction equation with discontinuous solution in [Z. Cai et al., J. Comput. Phys., 443 (2021), 110514]. The method is based on a leastsquares formulation and uses a new class of approximating functions: ReLU neural network (NN) functions. A critical and additional component of the LSNN method, differing from other NNbased methods, is the introduction of a properly designed and physics preserved discrete differential operator. In this paper, we study the LSNN method for problems with discontinuity interfaces. First, we show that ReLU NN functions with depth [math] can approximate any [math]dimensional step function on a discontinuity interface generated by a vector field as streamlines with any prescribed accuracy. By decomposing the solution into continuous and discontinuous parts, we prove theoretically that the discretization error of the LSNN method using ReLU NN functions with depth [math] is mainly determined by the continuous part of the solution provided that the solution jump is constant. Numerical results for both two and threedimensional test problems with various discontinuity interfaces show that the LSNN method with enough layers is accurate and does not exhibit the common Gibbs phenomena along discontinuity interfaces.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C448C478, August 2024. <br/> Abstract. We studied the leastsquares ReLU neural network (LSNN) method for solving a linear advectionreaction equation with discontinuous solution in [Z. Cai et al., J. Comput. Phys., 443 (2021), 110514]. The method is based on a leastsquares formulation and uses a new class of approximating functions: ReLU neural network (NN) functions. A critical and additional component of the LSNN method, differing from other NNbased methods, is the introduction of a properly designed and physics preserved discrete differential operator. In this paper, we study the LSNN method for problems with discontinuity interfaces. First, we show that ReLU NN functions with depth [math] can approximate any [math]dimensional step function on a discontinuity interface generated by a vector field as streamlines with any prescribed accuracy. By decomposing the solution into continuous and discontinuous parts, we prove theoretically that the discretization error of the LSNN method using ReLU NN functions with depth [math] is mainly determined by the continuous part of the solution provided that the solution jump is constant. Numerical results for both two and threedimensional test problems with various discontinuity interfaces show that the LSNN method with enough layers is accurate and does not exhibit the common Gibbs phenomena along discontinuity interfaces. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
LeastSquares Neural Network (LSNN) Method for Linear AdvectionReaction Equation: Discontinuity Interface
10.1137/23M1568107
SIAM Journal on Scientific Computing
20240812T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Zhiqiang Cai
Junpyo Choi
Min Liu
LeastSquares Neural Network (LSNN) Method for Linear AdvectionReaction Equation: Discontinuity Interface
46
4
C448
C478
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1568107
https://epubs.siam.org/doi/abs/10.1137/23M1568107?af=R
© 2024 Society for Industrial and Applied Mathematics

A Meshless Solver for Blood Flow Simulations in Elastic Vessels Using a PhysicsInformed Neural Network
https://epubs.siam.org/doi/abs/10.1137/23M1622696?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page C479C507, August 2024. <br/> Abstract. Investigating blood flow in the cardiovascular system is crucial for assessing cardiovascular health. Computational approaches offer some noninvasive alternatives to measure blood flow dynamics. Numerical simulations based on traditional methods such as finiteelement and other numerical discretizations have been extensively studied and have yielded excellent results. However, adapting these methods to reallife simulations remains a complex task. In this paper, we propose a method that offers flexibility and can efficiently handle reallife simulations. We suggest utilizing the physicsinformed neural network to solve the Navier–Stokes equation in a deformable domain, specifically addressing the simulation of blood flow in elastic vessels. Our approach models blood flow using an incompressible, viscous Navier–Stokes equation in an arbitrary Lagrangian–Eulerian form. The mechanical model for the vessel wall structure is formulated by an equation of Newton’s second law of momentum and linear elasticity to the force exerted by the fluid flow. Our method is a meshfree approach that eliminates the need for discretization and meshing of the computational domain. This makes it highly efficient in solving simulations involving complex geometries. Additionally, with the availability of welldeveloped opensource machine learning framework packages and parallel modules, our method can easily be accelerated through GPU computing and parallel computing. To evaluate our approach, we conducted experiments on regular cylinder vessels as well as vessels with plaque on their walls. We compared our results to a solution calculated by finite element methods using a dense grid and small time steps, which we considered as the ground truth solution. We report the relative error and the time consumed to solve the problem, highlighting the advantages of our method.
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C479C507, August 2024. <br/> Abstract. Investigating blood flow in the cardiovascular system is crucial for assessing cardiovascular health. Computational approaches offer some noninvasive alternatives to measure blood flow dynamics. Numerical simulations based on traditional methods such as finiteelement and other numerical discretizations have been extensively studied and have yielded excellent results. However, adapting these methods to reallife simulations remains a complex task. In this paper, we propose a method that offers flexibility and can efficiently handle reallife simulations. We suggest utilizing the physicsinformed neural network to solve the Navier–Stokes equation in a deformable domain, specifically addressing the simulation of blood flow in elastic vessels. Our approach models blood flow using an incompressible, viscous Navier–Stokes equation in an arbitrary Lagrangian–Eulerian form. The mechanical model for the vessel wall structure is formulated by an equation of Newton’s second law of momentum and linear elasticity to the force exerted by the fluid flow. Our method is a meshfree approach that eliminates the need for discretization and meshing of the computational domain. This makes it highly efficient in solving simulations involving complex geometries. Additionally, with the availability of welldeveloped opensource machine learning framework packages and parallel modules, our method can easily be accelerated through GPU computing and parallel computing. To evaluate our approach, we conducted experiments on regular cylinder vessels as well as vessels with plaque on their walls. We compared our results to a solution calculated by finite element methods using a dense grid and small time steps, which we considered as the ground truth solution. We report the relative error and the time consumed to solve the problem, highlighting the advantages of our method. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A Meshless Solver for Blood Flow Simulations in Elastic Vessels Using a PhysicsInformed Neural Network
10.1137/23M1622696
SIAM Journal on Scientific Computing
20240819T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Han Zhang
Raymond H. Chan
XueCheng Tai
A Meshless Solver for Blood Flow Simulations in Elastic Vessels Using a PhysicsInformed Neural Network
46
4
C479
C507
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1622696
https://epubs.siam.org/doi/abs/10.1137/23M1622696?af=R
© 2024 Society for Industrial and Applied Mathematics

A Pseudoreversible Normalizing Flow for Stochastic Dynamical Systems with Various Initial Distributions
https://epubs.siam.org/doi/abs/10.1137/23M1585635?af=R
SIAM Journal on Scientific Computing, <a href="https://epubs.siam.org/toc/sjoce3/46/4">Volume 46, Issue 4</a>, Page C508C533, August 2024. <br/> Abstract. We present a pseudoreversible normalizing flow method for efficiently generating samples of the state of a stochastic differential equation (SDE) with various initial distributions. The primary objective is to construct an accurate and efficient sampler that can be used as a surrogate model for computationally expensive numerical integration of SDEs, such as those employed in particle simulation. After training, the normalizing flow model can directly generate samples of the SDE’s final state without simulating trajectories. The existing normalizing flow model for SDEs depends on the initial distribution, meaning the model needs to be retrained when the initial distribution changes. The main novelty of our normalizing flow model is that it can learn the conditional distribution of the state, i.e., the distribution of the final state conditional on any initial state, such that the model only needs to be trained once and the trained model can be used to handle various initial distributions. This feature can provide a significant computational saving in studies of how the final state varies with the initial distribution. Additionally, we propose to use a pseudoreversible network architecture to define the normalizing flow model, which has sufficient expressive power and training efficiency for a variety of SDEs in science and engineering, e.g., in particle physics. We provide a rigorous convergence analysis of the pseudoreversible normalizing flow model to the target probability density function in the Kullback–Leibler divergence metric. Numerical experiments are provided to demonstrate the effectiveness of the proposed normalizing flow model. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/mlmathphy/PRNF and in the supplementary materials (PRNFmain.zip [27.4MB]).
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C508C533, August 2024. <br/> Abstract. We present a pseudoreversible normalizing flow method for efficiently generating samples of the state of a stochastic differential equation (SDE) with various initial distributions. The primary objective is to construct an accurate and efficient sampler that can be used as a surrogate model for computationally expensive numerical integration of SDEs, such as those employed in particle simulation. After training, the normalizing flow model can directly generate samples of the SDE’s final state without simulating trajectories. The existing normalizing flow model for SDEs depends on the initial distribution, meaning the model needs to be retrained when the initial distribution changes. The main novelty of our normalizing flow model is that it can learn the conditional distribution of the state, i.e., the distribution of the final state conditional on any initial state, such that the model only needs to be trained once and the trained model can be used to handle various initial distributions. This feature can provide a significant computational saving in studies of how the final state varies with the initial distribution. Additionally, we propose to use a pseudoreversible network architecture to define the normalizing flow model, which has sufficient expressive power and training efficiency for a variety of SDEs in science and engineering, e.g., in particle physics. We provide a rigorous convergence analysis of the pseudoreversible normalizing flow model to the target probability density function in the Kullback–Leibler divergence metric. Numerical experiments are provided to demonstrate the effectiveness of the proposed normalizing flow model. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/mlmathphy/PRNF and in the supplementary materials (PRNFmain.zip [27.4MB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A Pseudoreversible Normalizing Flow for Stochastic Dynamical Systems with Various Initial Distributions
10.1137/23M1585635
SIAM Journal on Scientific Computing
20240821T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Minglei Yang
Pengjun Wang
Diego delCastilloNegrete
Yanzhao Cao
Guannan Zhang
A Pseudoreversible Normalizing Flow for Stochastic Dynamical Systems with Various Initial Distributions
46
4
C508
C533
20240831T07:00:00Z
20240831T07:00:00Z
10.1137/23M1585635
https://epubs.siam.org/doi/abs/10.1137/23M1585635?af=R
© 2024 Society for Industrial and Applied Mathematics

Reducing Operator Complexity of Galerkin Coarsegrid Operators with Machine Learning
https://epubs.siam.org/doi/abs/10.1137/23M1583533?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. We propose a datadriven and machinelearningbased approach to compute nonGalerkin coarsegrid operators in multigrid (MG) methods, addressing the wellknown issue of increasing operator complexity. Guided by the MG theory on spectrally equivalent coarsegrid operators, we have developed novel machine learning algorithms that utilize neural networks combined with smooth test vectors from multigrid eigenvalue problems. The proposed method demonstrates promise in reducing the complexity of coarsegrid operators while maintaining overall MG convergence for solving parametric partial differential equation problems. Numerical experiments on anisotropic rotated Laplacian and linear elasticity problems are provided to showcase the performance and comparison with existing methods for computing nonGalerkin coarsegrid operators. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/liruipeng/SparseCoarseOperator.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. We propose a datadriven and machinelearningbased approach to compute nonGalerkin coarsegrid operators in multigrid (MG) methods, addressing the wellknown issue of increasing operator complexity. Guided by the MG theory on spectrally equivalent coarsegrid operators, we have developed novel machine learning algorithms that utilize neural networks combined with smooth test vectors from multigrid eigenvalue problems. The proposed method demonstrates promise in reducing the complexity of coarsegrid operators while maintaining overall MG convergence for solving parametric partial differential equation problems. Numerical experiments on anisotropic rotated Laplacian and linear elasticity problems are provided to showcase the performance and comparison with existing methods for computing nonGalerkin coarsegrid operators. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/liruipeng/SparseCoarseOperator. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Reducing Operator Complexity of Galerkin Coarsegrid Operators with Machine Learning
10.1137/23M1583533
SIAM Journal on Scientific Computing
20240906T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Ru Huang
Kai Chang
Huan He
Ruipeng Li
Yuanzhe Xi
Reducing Operator Complexity of Galerkin Coarsegrid Operators with Machine Learning
S296
S316
10.1137/23M1583533
https://epubs.siam.org/doi/abs/10.1137/23M1583533?af=R
© 2024 Society for Industrial and Applied Mathematics

Multigrid Preconditioning for Regularized LeastSquares Problems
https://epubs.siam.org/doi/abs/10.1137/23M1583417?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. In this paper, we are concerned with efficiently solving the sequences of regularized linear leastsquares problems associated with employing Tikhonovtype regularization with regularization operators designed to enforce edge recovery. An optimal regularization parameter, which balances the fidelity to the data with the edgeenforcing constraint term, is typically not known a priori. This adds to the total number of regularized linear leastsquares problems that must be solved before the final image can be recovered. Therefore, in this paper, we determine effective multigrid preconditioners for these sequences of systems. We focus our approach on the sequences that arise as a result of the edgepreserving method introduced in [S. Gazzola et al., Inverse Problems, 36 (2020), 124004], where we can exploit an interpretation of the regularization term as a diffusion operator; however, our methods are also applicable in other edgepreserving settings, such as iteratively reweighted leastsquares problems. Particular attention is paid to the selection of components of the multigrid preconditioner in order to achieve robustness for different ranges of the regularization parameter value. In addition, we present a parameter trimming approach that, when used with the Lcurve heuristic, reduces the total number of solves required. We demonstrate our preconditioning and parameter trimming routines on examples in computed tomography and image deblurring.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. In this paper, we are concerned with efficiently solving the sequences of regularized linear leastsquares problems associated with employing Tikhonovtype regularization with regularization operators designed to enforce edge recovery. An optimal regularization parameter, which balances the fidelity to the data with the edgeenforcing constraint term, is typically not known a priori. This adds to the total number of regularized linear leastsquares problems that must be solved before the final image can be recovered. Therefore, in this paper, we determine effective multigrid preconditioners for these sequences of systems. We focus our approach on the sequences that arise as a result of the edgepreserving method introduced in [S. Gazzola et al., Inverse Problems, 36 (2020), 124004], where we can exploit an interpretation of the regularization term as a diffusion operator; however, our methods are also applicable in other edgepreserving settings, such as iteratively reweighted leastsquares problems. Particular attention is paid to the selection of components of the multigrid preconditioner in order to achieve robustness for different ranges of the regularization parameter value. In addition, we present a parameter trimming approach that, when used with the Lcurve heuristic, reduces the total number of solves required. We demonstrate our preconditioning and parameter trimming routines on examples in computed tomography and image deblurring. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Multigrid Preconditioning for Regularized LeastSquares Problems
10.1137/23M1583417
SIAM Journal on Scientific Computing
20240905T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Matthias Bolten
Misha E. Kilmer
Scott MacLachlan
Multigrid Preconditioning for Regularized LeastSquares Problems
S271
S295
10.1137/23M1583417
https://epubs.siam.org/doi/abs/10.1137/23M1583417?af=R
© 2024 Society for Industrial and Applied Mathematics

Efficient GMRES+AMG on GPUs: Composite Smoothers And Mixed [math]Cycles
https://epubs.siam.org/doi/abs/10.1137/23M1578632?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. In this study, we introduce algorithms optimized for GPU architectures, aimed at efficiently solving large sparse linear systems, a central challenge in Navier–Stokes pressure projection problems. Our approach includes an adaptation of the GMRES algorithm, drawing inspiration from the merged vector operations first proposed by Bielich et al. [Parallel Comput., 112 (2022), 102940]. This adaptation increases computational intensity on GPU platforms through optimized vector update strategies. The algorithm incorporates modified and classical Gram–Schmidt methods with an algebraic multigrid (AMG) preconditioner, each tailored for GPU performance. A key innovation in our work is the development of a Gram–Schmidt projector [math] employing a rank1 perturbation of the identity matrix. Designed to maximize the high memory bandwidth utilization of the AMD MI250X GPU, this approach includes a strategy for treating the unit diagonal that minimizes memory reads, leading to a 25% increase in computational efficiency. The application of perturbation theory further ensures that orthogonality loss is limited to [math], where [math] is the number of iterations. Additionally, we introduce a mixed AMG [math]cycle strategy combining ILU(0) and [math]Jacobi smoothers, which achieves a 30–50% reduction in GPU compute times compared to conventional methods, while maintaining low backward error. This strategy, alongside our novel treatment of the diagonal in triangular matrices, marks a substantial increase in AMG efficicency for GPU systems. We believe that these contributions represent a significant advance in optimizing GMRES+AMG algorithms for GPU computations. The empirical results demonstrate notable speed increments and maintain rigorous backward error bounds, underscoring the potential of our methods to substantially increase computational efficiency in largescale scientific applications.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. In this study, we introduce algorithms optimized for GPU architectures, aimed at efficiently solving large sparse linear systems, a central challenge in Navier–Stokes pressure projection problems. Our approach includes an adaptation of the GMRES algorithm, drawing inspiration from the merged vector operations first proposed by Bielich et al. [Parallel Comput., 112 (2022), 102940]. This adaptation increases computational intensity on GPU platforms through optimized vector update strategies. The algorithm incorporates modified and classical Gram–Schmidt methods with an algebraic multigrid (AMG) preconditioner, each tailored for GPU performance. A key innovation in our work is the development of a Gram–Schmidt projector [math] employing a rank1 perturbation of the identity matrix. Designed to maximize the high memory bandwidth utilization of the AMD MI250X GPU, this approach includes a strategy for treating the unit diagonal that minimizes memory reads, leading to a 25% increase in computational efficiency. The application of perturbation theory further ensures that orthogonality loss is limited to [math], where [math] is the number of iterations. Additionally, we introduce a mixed AMG [math]cycle strategy combining ILU(0) and [math]Jacobi smoothers, which achieves a 30–50% reduction in GPU compute times compared to conventional methods, while maintaining low backward error. This strategy, alongside our novel treatment of the diagonal in triangular matrices, marks a substantial increase in AMG efficicency for GPU systems. We believe that these contributions represent a significant advance in optimizing GMRES+AMG algorithms for GPU computations. The empirical results demonstrate notable speed increments and maintain rigorous backward error bounds, underscoring the potential of our methods to substantially increase computational efficiency in largescale scientific applications. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Efficient GMRES+AMG on GPUs: Composite Smoothers And Mixed [math]Cycles
10.1137/23M1578632
SIAM Journal on Scientific Computing
20240903T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Stephen Thomas
Allison H. Baker
Efficient GMRES+AMG on GPUs: Composite Smoothers And Mixed [math]Cycles
S246
S270
10.1137/23M1578632
https://epubs.siam.org/doi/abs/10.1137/23M1578632?af=R
© 2024 Society for Industrial and Applied Mathematics

LFATuned MatrixFree Multigrid Method for the Elastic Helmholtz Equation
https://epubs.siam.org/doi/abs/10.1137/23M1583466?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. We present an efficient matrixfree geometric multigrid method for the elastic Helmholtz equation, and a suitable discretization. Many discretization methods had been considered in the literature for the Helmholtz equations, as well as many solvers and preconditioners, some of which are adapted for the elastic version of the equation. However, there is very little work considering the reciprocity of discretization and a solver. In this work, we aim to bridge this gap. By choosing an appropriate stencil for rediscretization of the equation on the coarse grid, we develop a multigrid method that can be easily implemented as matrixfree, relying on stencils rather than sparse matrices. This is crucial for efficient implementation on modern hardware. Using twogrid local Fourier analysis, we validate the compatibility of our discretization with our solver, and tune a choice of weights for the stencil for which the convergence rate of the multigrid cycle is optimal. It results in a scalable multigrid preconditioner that can tackle large realworld threedimensional scenarios.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. We present an efficient matrixfree geometric multigrid method for the elastic Helmholtz equation, and a suitable discretization. Many discretization methods had been considered in the literature for the Helmholtz equations, as well as many solvers and preconditioners, some of which are adapted for the elastic version of the equation. However, there is very little work considering the reciprocity of discretization and a solver. In this work, we aim to bridge this gap. By choosing an appropriate stencil for rediscretization of the equation on the coarse grid, we develop a multigrid method that can be easily implemented as matrixfree, relying on stencils rather than sparse matrices. This is crucial for efficient implementation on modern hardware. Using twogrid local Fourier analysis, we validate the compatibility of our discretization with our solver, and tune a choice of weights for the stencil for which the convergence rate of the multigrid cycle is optimal. It results in a scalable multigrid preconditioner that can tackle large realworld threedimensional scenarios. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
LFATuned MatrixFree Multigrid Method for the Elastic Helmholtz Equation
10.1137/23M1583466
SIAM Journal on Scientific Computing
20240305T08:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Rachel Yovel
Eran Treister
LFATuned MatrixFree Multigrid Method for the Elastic Helmholtz Equation
S1
S21
10.1137/23M1583466
https://epubs.siam.org/doi/abs/10.1137/23M1583466?af=R
© 2024 Society for Industrial and Applied Mathematics

On the Convergence of Monolithic Multigrid for Implicit Runge–Kutta Time Stepping of Finite Element Problems
https://epubs.siam.org/doi/abs/10.1137/23M1569344?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. Finite element discretizations of timedependent problems also require effective timestepping schemes. While implicit Runge–Kutta methods provide favorable accuracy and stability properties, they give rise to large and complicated systems of equations to solve for each time step. These algebraic systems couple all Runge–Kutta stages together, giving a much larger system than for singlestage methods. We consider an approach to these systems based on monolithic smoothing. If stagecoupled smoothers possess a certain kind of structure, then the question of convergence of a twogrid or multigrid iteration reduces to convergence of a related strategy for a singlestage system with a complexvalued time step. In addition to providing a general theoretical approach to the convergence of monolithic multigrid methods, several numerical examples are given to illustrate the theory and show how higherorder Runge–Kutta methods can be made effective in practice. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and Data Available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/rckirby/CodeForMMGPaper as well as in the supplemental material.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. Finite element discretizations of timedependent problems also require effective timestepping schemes. While implicit Runge–Kutta methods provide favorable accuracy and stability properties, they give rise to large and complicated systems of equations to solve for each time step. These algebraic systems couple all Runge–Kutta stages together, giving a much larger system than for singlestage methods. We consider an approach to these systems based on monolithic smoothing. If stagecoupled smoothers possess a certain kind of structure, then the question of convergence of a twogrid or multigrid iteration reduces to convergence of a related strategy for a singlestage system with a complexvalued time step. In addition to providing a general theoretical approach to the convergence of monolithic multigrid methods, several numerical examples are given to illustrate the theory and show how higherorder Runge–Kutta methods can be made effective in practice. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and Data Available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/rckirby/CodeForMMGPaper as well as in the supplemental material. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
On the Convergence of Monolithic Multigrid for Implicit Runge–Kutta Time Stepping of Finite Element Problems
10.1137/23M1569344
SIAM Journal on Scientific Computing
20240313T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Robert C. Kirby
On the Convergence of Monolithic Multigrid for Implicit Runge–Kutta Time Stepping of Finite Element Problems
S22
S45
10.1137/23M1569344
https://epubs.siam.org/doi/abs/10.1137/23M1569344?af=R
© 2024 Society for Industrial and Applied Mathematics

Enhancing Training of PhysicsInformed Neural Networks Using Domain Decomposition–Based Preconditioning Strategies
https://epubs.siam.org/doi/abs/10.1137/23M1583375?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. We propose to enhance the training of physicsinformed neural networks. To this aim, we introduce nonlinear additive and multiplicative preconditioning strategies for the widely used LBFGS optimizer. The nonlinear preconditioners are constructed by utilizing the Schwarz domain decomposition framework, where the parameters of the network are decomposed in a layerwise manner. Through a series of numerical experiments, we demonstrate that both additive and multiplicative preconditioners significantly improve the convergence of the standard LBFGS optimizer while providing more accurate solutions of the underlying PDEs. Moreover, the additive preconditioner is inherently parallel, thus giving rise to a novel approach to model parallelism.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. We propose to enhance the training of physicsinformed neural networks. To this aim, we introduce nonlinear additive and multiplicative preconditioning strategies for the widely used LBFGS optimizer. The nonlinear preconditioners are constructed by utilizing the Schwarz domain decomposition framework, where the parameters of the network are decomposed in a layerwise manner. Through a series of numerical experiments, we demonstrate that both additive and multiplicative preconditioners significantly improve the convergence of the standard LBFGS optimizer while providing more accurate solutions of the underlying PDEs. Moreover, the additive preconditioner is inherently parallel, thus giving rise to a novel approach to model parallelism. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Enhancing Training of PhysicsInformed Neural Networks Using Domain Decomposition–Based Preconditioning Strategies
10.1137/23M1583375
SIAM Journal on Scientific Computing
20240411T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Alena Kopaničáková
Hardik Kothari
George E. Karniadakis
Rolf Krause
Enhancing Training of PhysicsInformed Neural Networks Using Domain Decomposition–Based Preconditioning Strategies
S46
S67
10.1137/23M1583375
https://epubs.siam.org/doi/abs/10.1137/23M1583375?af=R
© 2024 Society for Industrial and Applied Mathematics

A Reduced Conjugate Gradient Basis Method for Fractional Diffusion
https://epubs.siam.org/doi/abs/10.1137/23M1575913?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. This work is on a fast and accurate reduced basis method for solving discretized fractional elliptic partial differential equations (PDEs) of the form [math] by rational approximation. A direct computation of the action of such an approximation would require solving multiple (20[math]30) largescale sparse linear systems. Our method constructs the reduced basis using the first few directions obtained from the preconditioned conjugate gradient method applied to one of the linear systems. As shown in the theory and experiments, only a small number of directions (5[math]10) are needed to approximately solve all largescale systems on the reduced basis subspace. This reduces the computational cost dramatically because: (1) We only use one of the largescale problems to construct the basis; and (2) all largescale problems restricted to the subspace have much smaller sizes. We test our algorithms for fractional PDEs on a 3d Euclidean domain, a 2d surface, and random combinatorial graphs. We also use a novel approach to construct the rational approximation for the fractional power function by the orthogonal greedy algorithm (OGA). Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/yuwenli925/RCGBM and in the supplementary materials (RCGBMmain.zip [19.8KB]).
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. This work is on a fast and accurate reduced basis method for solving discretized fractional elliptic partial differential equations (PDEs) of the form [math] by rational approximation. A direct computation of the action of such an approximation would require solving multiple (20[math]30) largescale sparse linear systems. Our method constructs the reduced basis using the first few directions obtained from the preconditioned conjugate gradient method applied to one of the linear systems. As shown in the theory and experiments, only a small number of directions (5[math]10) are needed to approximately solve all largescale systems on the reduced basis subspace. This reduces the computational cost dramatically because: (1) We only use one of the largescale problems to construct the basis; and (2) all largescale problems restricted to the subspace have much smaller sizes. We test our algorithms for fractional PDEs on a 3d Euclidean domain, a 2d surface, and random combinatorial graphs. We also use a novel approach to construct the rational approximation for the fractional power function by the orthogonal greedy algorithm (OGA). Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/yuwenli925/RCGBM and in the supplementary materials (RCGBMmain.zip [19.8KB]). <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A Reduced Conjugate Gradient Basis Method for Fractional Diffusion
10.1137/23M1575913
SIAM Journal on Scientific Computing
20240418T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Yuwen Li
Ludmil Zikatanov
Cheng Zuo
A Reduced Conjugate Gradient Basis Method for Fractional Diffusion
S68
S87
10.1137/23M1575913
https://epubs.siam.org/doi/abs/10.1137/23M1575913?af=R
© 2024 Society for Industrial and Applied Mathematics

RoundingError Analysis of Multigrid [math]Cycles
https://epubs.siam.org/doi/abs/10.1137/23M1582898?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. Earlier work on roundingerror analysis of multigrid was restricted to cycles that used one relaxation step before coarsening and none afterwards. The present paper extends this analysis to twogrid methods that use one relaxation step both before and after coarsening. The analysis is based on floating point arithmetic and focuses on a twogrid scheme that is perturbed on the coarse grid to allow for an approximate coarsegrid solve. Leveraging previously published results, this twogrid theory can then be extended to general [math]cycles, as well as full multigrid. It can also be extended to mixedprecision iterative refinement based on these cycles. An added benefit of the theory here over previous work is that it is obtained in a more organized, transparent, and simpler way.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. Earlier work on roundingerror analysis of multigrid was restricted to cycles that used one relaxation step before coarsening and none afterwards. The present paper extends this analysis to twogrid methods that use one relaxation step both before and after coarsening. The analysis is based on floating point arithmetic and focuses on a twogrid scheme that is perturbed on the coarse grid to allow for an approximate coarsegrid solve. Leveraging previously published results, this twogrid theory can then be extended to general [math]cycles, as well as full multigrid. It can also be extended to mixedprecision iterative refinement based on these cycles. An added benefit of the theory here over previous work is that it is obtained in a more organized, transparent, and simpler way. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
RoundingError Analysis of Multigrid [math]Cycles
10.1137/23M1582898
SIAM Journal on Scientific Computing
20240419T07:00:00Z
© 2024 Disney
Stephen F. McCormick
Rasmus Tamstorf
RoundingError Analysis of Multigrid [math]Cycles
S88
S95
10.1137/23M1582898
https://epubs.siam.org/doi/abs/10.1137/23M1582898?af=R
© 2024 Disney

Constrained Local Approximate Ideal Restriction for AdvectionDiffusion Problems
https://epubs.siam.org/doi/abs/10.1137/23M1583442?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. This paper focuses on developing a reductionbased algebraic multigrid (AMG) method that is suitable for solving general (non)symmetric linear systems and is naturally robust from pure advection to pure diffusion. Initial motivation comes from a new reductionbased AMG approach, [math] (local approximate ideal restriction), that was developed for solving advectiondominated problems. Though this new solver is very effective in the advectiondominated regime, its performance degrades in cases where diffusion becomes dominant. This is consistent with the fact that in general, reductionbased AMG methods tend to suffer from growth in complexity and/or convergence rates as the problem size is increased, especially for diffusiondominated problems in two or three dimensions. Motivated by the success of [math] in the advective regime, our aim in this paper is to generalize the AIR framework with the goal of improving the performance of the solver in diffusiondominated regimes. To do so, we propose a novel way to combine mode constraints as used commonly in energyminimization AMG methods with the local approximation of ideal operators used in [math]. The resulting constrained [math] algorithm is able to achieve fast scalable convergence on advective and diffusive problems. In addition, it is able to achieve standard low complexity hierarchies in the diffusive regime through aggressive coarsening, something that was previously difficult for reductionbased methods.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. This paper focuses on developing a reductionbased algebraic multigrid (AMG) method that is suitable for solving general (non)symmetric linear systems and is naturally robust from pure advection to pure diffusion. Initial motivation comes from a new reductionbased AMG approach, [math] (local approximate ideal restriction), that was developed for solving advectiondominated problems. Though this new solver is very effective in the advectiondominated regime, its performance degrades in cases where diffusion becomes dominant. This is consistent with the fact that in general, reductionbased AMG methods tend to suffer from growth in complexity and/or convergence rates as the problem size is increased, especially for diffusiondominated problems in two or three dimensions. Motivated by the success of [math] in the advective regime, our aim in this paper is to generalize the AIR framework with the goal of improving the performance of the solver in diffusiondominated regimes. To do so, we propose a novel way to combine mode constraints as used commonly in energyminimization AMG methods with the local approximation of ideal operators used in [math]. The resulting constrained [math] algorithm is able to achieve fast scalable convergence on advective and diffusive problems. In addition, it is able to achieve standard low complexity hierarchies in the diffusive regime through aggressive coarsening, something that was previously difficult for reductionbased methods. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Constrained Local Approximate Ideal Restriction for AdvectionDiffusion Problems
10.1137/23M1583442
SIAM Journal on Scientific Computing
20240502T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Ahsan Ali
James J. Brannick
Karsten Kahl
Oliver A. Krzysik
Jacob B. Schroder
Ben S. Southworth
Constrained Local Approximate Ideal Restriction for AdvectionDiffusion Problems
S96
S122
10.1137/23M1583442
https://epubs.siam.org/doi/abs/10.1137/23M1583442?af=R
© 2024 Society for Industrial and Applied Mathematics

MultigridAugmented Deep Learning Preconditioners for the Helmholtz Equation Using Compact Implicit Layers
https://epubs.siam.org/doi/abs/10.1137/23M1583302?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. We present a deep learning–based iterative approach to solve the discrete heterogeneous Helmholtz equation for high wavenumbers. Combining classical iterative multigrid solvers and convolutional neural networks (CNNs) via preconditioning, we obtain a faster, learned neural solver that scales better than a standard multigrid solver. Our approach offers three main contributions over previous neural methods of this kind. First, we construct a multilevel UNetlike encodersolver CNN with an implicit layer on the coarsest grid of the UNet, where convolution kernels are inverted. This alleviates the field of view problem in CNNs and allows better scalability. Second, we improve upon the previous CNN preconditioner in terms of the number of parameters, computation time, and convergence rates. Third, we propose a multiscale training approach that enables the network to scale to problems of previously unseen dimensions while still maintaining a reasonable training procedure. Our encodersolver architecture can be used to generalize over different slowness models of various difficulties and is efficient at solving for many righthand sides per slowness model. We demonstrate the benefits of our novel architecture with numerical experiments on various heterogeneous twodimensional problems at high wavenumbers.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. We present a deep learning–based iterative approach to solve the discrete heterogeneous Helmholtz equation for high wavenumbers. Combining classical iterative multigrid solvers and convolutional neural networks (CNNs) via preconditioning, we obtain a faster, learned neural solver that scales better than a standard multigrid solver. Our approach offers three main contributions over previous neural methods of this kind. First, we construct a multilevel UNetlike encodersolver CNN with an implicit layer on the coarsest grid of the UNet, where convolution kernels are inverted. This alleviates the field of view problem in CNNs and allows better scalability. Second, we improve upon the previous CNN preconditioner in terms of the number of parameters, computation time, and convergence rates. Third, we propose a multiscale training approach that enables the network to scale to problems of previously unseen dimensions while still maintaining a reasonable training procedure. Our encodersolver architecture can be used to generalize over different slowness models of various difficulties and is efficient at solving for many righthand sides per slowness model. We demonstrate the benefits of our novel architecture with numerical experiments on various heterogeneous twodimensional problems at high wavenumbers. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
MultigridAugmented Deep Learning Preconditioners for the Helmholtz Equation Using Compact Implicit Layers
10.1137/23M1583302
SIAM Journal on Scientific Computing
20240603T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Bar Lerer
Ido BenYair
Eran Treister
MultigridAugmented Deep Learning Preconditioners for the Helmholtz Equation Using Compact Implicit Layers
S123
S144
10.1137/23M1583302
https://epubs.siam.org/doi/abs/10.1137/23M1583302?af=R
© 2024 Society for Industrial and Applied Mathematics

FaultTolerant Parallel Multigrid Method on Unstructured Adaptive Mesh
https://epubs.siam.org/doi/abs/10.1137/23M1582904?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. As the generation of exascale highperformance clusters begins, it has become evident that numerical algorithms will greatly benefit from builtin resilience features that can handle system faults. Prior studies of faulttolerant multigrid methods have focused on structured grids. In this work, however, we study the resilience of multigrid solvers on unstructured grids with adaptive refinement. The challenge lies in the fact that unstructured grids distributed across multiple processors may manifest as local hierarchical grids with unaligned boundaries. Our numerical experiments highlight that this disparity can result in divergence when employing standard local multigrid for fault recovery. We analyze this phenomenon by using an energy control condition. To tackle the divergence issue, we propose a simple variation of the multigrid Vcycle that scales the coarse problem. We present a convergence proof for the new algorithm. By implementing this new method for local recovery, our numerical experiments confirm that convergence can be recovered on unstructured grids while the algorithm agrees with the standard multigrid Vcycle on grids with aligned boundaries. More importantly, the impact of a fault can be mitigated and delays in the global multigrid iterations can be reduced. Finally, we investigate how local regions within the adaptive mesh, associated with different faulty processors, affect the effectiveness of fault recovery.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. As the generation of exascale highperformance clusters begins, it has become evident that numerical algorithms will greatly benefit from builtin resilience features that can handle system faults. Prior studies of faulttolerant multigrid methods have focused on structured grids. In this work, however, we study the resilience of multigrid solvers on unstructured grids with adaptive refinement. The challenge lies in the fact that unstructured grids distributed across multiple processors may manifest as local hierarchical grids with unaligned boundaries. Our numerical experiments highlight that this disparity can result in divergence when employing standard local multigrid for fault recovery. We analyze this phenomenon by using an energy control condition. To tackle the divergence issue, we propose a simple variation of the multigrid Vcycle that scales the coarse problem. We present a convergence proof for the new algorithm. By implementing this new method for local recovery, our numerical experiments confirm that convergence can be recovered on unstructured grids while the algorithm agrees with the standard multigrid Vcycle on grids with aligned boundaries. More importantly, the impact of a fault can be mitigated and delays in the global multigrid iterations can be reduced. Finally, we investigate how local regions within the adaptive mesh, associated with different faulty processors, affect the effectiveness of fault recovery. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
FaultTolerant Parallel Multigrid Method on Unstructured Adaptive Mesh
10.1137/23M1582904
SIAM Journal on Scientific Computing
20240606T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Frederick Fung
Linda Stals
Quanling Deng
FaultTolerant Parallel Multigrid Method on Unstructured Adaptive Mesh
S145
S169
10.1137/23M1582904
https://epubs.siam.org/doi/abs/10.1137/23M1582904?af=R
© 2024 Society for Industrial and Applied Mathematics

Scalable Multiphysics Block Preconditioning for Low Mach Number Compressible Resistive MHD with Application to Magnetic Confinement Fusion
https://epubs.siam.org/doi/abs/10.1137/23M1582667?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. This study investigates multiphysics block preconditioners that are critical in devising scalable Newton–Krylov iterative solvers for longer timescale fully implicit fluid plasma models. The specific model of interest is the viscoresistive, low Mach number, compressible magnetohydrodynamics (MHD) model. This model describes the dynamics of conducting fluids in the presence of electromagnetic fields and can be used to study aspects of astrophysical phenomena, important science and technology applications, and basic plasma physics. The specific application of interest that motivates this study is the macroscopic simulation of longer timescale stability and disruptions of magnetic confinement fusion devices, specifically the ITER Tokamak. The computational solution of the governing balance equations for mass, momentum, heat transfer, and magnetic induction for resistive MHD systems can be extremely challenging. These difficulties arise from both the strong nonlinear, nonsymmetric coupling of fluid and electromagnetic phenomena as well as the significant range of time and length scales that the interactions of these physical mechanisms produce. To handle the range of time and spatial scales of interest, a fully implicit unstructured variational multiscale finite element formulation is employed. For the scalable solution of the Newton linearized systems, fully coupled block preconditioners are designed to leverage algebraic multigrid subsolves. Results are presented for the strong and weak scaling of the method as well as the robustness of these techniques for a large range of Lundquist numbers.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. This study investigates multiphysics block preconditioners that are critical in devising scalable Newton–Krylov iterative solvers for longer timescale fully implicit fluid plasma models. The specific model of interest is the viscoresistive, low Mach number, compressible magnetohydrodynamics (MHD) model. This model describes the dynamics of conducting fluids in the presence of electromagnetic fields and can be used to study aspects of astrophysical phenomena, important science and technology applications, and basic plasma physics. The specific application of interest that motivates this study is the macroscopic simulation of longer timescale stability and disruptions of magnetic confinement fusion devices, specifically the ITER Tokamak. The computational solution of the governing balance equations for mass, momentum, heat transfer, and magnetic induction for resistive MHD systems can be extremely challenging. These difficulties arise from both the strong nonlinear, nonsymmetric coupling of fluid and electromagnetic phenomena as well as the significant range of time and length scales that the interactions of these physical mechanisms produce. To handle the range of time and spatial scales of interest, a fully implicit unstructured variational multiscale finite element formulation is employed. For the scalable solution of the Newton linearized systems, fully coupled block preconditioners are designed to leverage algebraic multigrid subsolves. Results are presented for the strong and weak scaling of the method as well as the robustness of these techniques for a large range of Lundquist numbers. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Scalable Multiphysics Block Preconditioning for Low Mach Number Compressible Resistive MHD with Application to Magnetic Confinement Fusion
10.1137/23M1582667
SIAM Journal on Scientific Computing
20240705T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Peter Ohm
Jesus Bonilla
Edward Phillips
John N. Shadid
Michael Crockatt
Ray S. Tuminaro
Jonathan Hu
XianZhu Tang
Scalable Multiphysics Block Preconditioning for Low Mach Number Compressible Resistive MHD with Application to Magnetic Confinement Fusion
S170
S201
10.1137/23M1582667
https://epubs.siam.org/doi/abs/10.1137/23M1582667?af=R
© 2024 Society for Industrial and Applied Mathematics

Multigrid Methods Using Block Floating Point Arithmetic
https://epubs.siam.org/doi/abs/10.1137/23M1581819?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. Block floating point (BFP) arithmetic is currently seeing a resurgence in interest because it requires less power and less chip area and is less complicated to implement in hardware than standard floating point arithmetic. This paper explores the application of BFP to mixed and progressiveprecision multigrid methods, enabling the solution of linear elliptic partial differential equations (PDEs) in energy and hardwareefficient integer arithmetic. While most existing applications of BFP arithmetic tend to use small block sizes, the block size here is chosen to be maximal such that matrices and vectors share a single exponent for all entries. This is sometimes also referred to as a scaled fixed point format. We provide algorithms for BLASlike routines for BFP arithmetic that ensure exact vectorvector and matrixvector computations up to a specified precision. Using these algorithms, we study the asymptotic precision requirements for achieving discretizationerroraccuracy. We demonstrate that some computations can be performed using only 4bit integers, while the number of bits required to attain a certain target accuracy is similar to that of standard floating point arithmetic. Finally, we present a heuristic for full multigrid in BFP arithmetic based on saturation and truncation that still achieves discretizationerroraccuracy without the need for expensive normalization steps of intermediate results.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. Block floating point (BFP) arithmetic is currently seeing a resurgence in interest because it requires less power and less chip area and is less complicated to implement in hardware than standard floating point arithmetic. This paper explores the application of BFP to mixed and progressiveprecision multigrid methods, enabling the solution of linear elliptic partial differential equations (PDEs) in energy and hardwareefficient integer arithmetic. While most existing applications of BFP arithmetic tend to use small block sizes, the block size here is chosen to be maximal such that matrices and vectors share a single exponent for all entries. This is sometimes also referred to as a scaled fixed point format. We provide algorithms for BLASlike routines for BFP arithmetic that ensure exact vectorvector and matrixvector computations up to a specified precision. Using these algorithms, we study the asymptotic precision requirements for achieving discretizationerroraccuracy. We demonstrate that some computations can be performed using only 4bit integers, while the number of bits required to attain a certain target accuracy is similar to that of standard floating point arithmetic. Finally, we present a heuristic for full multigrid in BFP arithmetic based on saturation and truncation that still achieves discretizationerroraccuracy without the need for expensive normalization steps of intermediate results. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
Multigrid Methods Using Block Floating Point Arithmetic
10.1137/23M1581819
SIAM Journal on Scientific Computing
20240711T07:00:00Z
© 2024 Disney
Nils Kohl
Stephen F. McCormick
Rasmus Tamstorf
Multigrid Methods Using Block Floating Point Arithmetic
S202
S224
10.1137/23M1581819
https://epubs.siam.org/doi/abs/10.1137/23M1581819?af=R
© 2024 Disney

A FirstOrder Reduced Model for a Highly Oscillating Differential Equation with Application in Penning Traps
https://epubs.siam.org/doi/abs/10.1137/23M158351X?af=R
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. We derive a reduced firstorder model from a twoscale asymptotic expansion in a small parameter in order to approximate the solution of a stiff differential equation. The problem of interest is a multiscale Newton–Lorentz equation modeling the dynamics of a charged particle under the influence of a linear electric field and of a perturbed strong magnetic field. First, we show that in short times, the firstorder model provides a much better approximation than the zeroorder one, since it contains terms evolving at slow time scales. Then, thanks to the sourcefree property of the equations, we propose a volumepreserving method using a particular splitting technique to solve numerically the firstorder model. Finally, it turns out that the firstorder model does not systematically provide a satisfactory approximation in long times. To overcome this issue, we implement a recent strategy based on the Parareal algorithm, in which the firstorder approximation is used for the coarse solver. This approach allows one to perform efficient and accurate longtime simulations for any small parameter. Numerical results for two realistic Penning traps are provided to support these statements.
SIAM Journal on Scientific Computing, Ahead of Print. <br/> Abstract. We derive a reduced firstorder model from a twoscale asymptotic expansion in a small parameter in order to approximate the solution of a stiff differential equation. The problem of interest is a multiscale Newton–Lorentz equation modeling the dynamics of a charged particle under the influence of a linear electric field and of a perturbed strong magnetic field. First, we show that in short times, the firstorder model provides a much better approximation than the zeroorder one, since it contains terms evolving at slow time scales. Then, thanks to the sourcefree property of the equations, we propose a volumepreserving method using a particular splitting technique to solve numerically the firstorder model. Finally, it turns out that the firstorder model does not systematically provide a satisfactory approximation in long times. To overcome this issue, we implement a recent strategy based on the Parareal algorithm, in which the firstorder approximation is used for the coarse solver. This approach allows one to perform efficient and accurate longtime simulations for any small parameter. Numerical results for two realistic Penning traps are provided to support these statements. <p><img src="https://epubs.siam.org/na101/home/literatum/publisher/siam/journals/covergifs/sjoce3/cover.jpg" alttext="cover image"/></p>
A FirstOrder Reduced Model for a Highly Oscillating Differential Equation with Application in Penning Traps
10.1137/23M158351X
SIAM Journal on Scientific Computing
20240809T07:00:00Z
© 2024 Society for Industrial and Applied Mathematics
Sever A. Hirstoaga
A FirstOrder Reduced Model for a Highly Oscillating Differential Equation with Application in Penning Traps
S225
S245
10.1137/23M158351X
https://epubs.siam.org/doi/abs/10.1137/23M158351X?af=R
© 2024 Society for Industrial and Applied Mathematics