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

Designing effective reduced-order models (ROMs) for parametrized transport-dominated problems remains challenging because of the well-known Kolmogorov barrier. Autoencoder-based nonlinear ROMs have been developed to improve the compression ability for such systems. However, despite their stronger compression ability, autoencoder-based ROMs constructed in the Eulerian frame may fail to accurately predict future solutions, due to the poor coherence between historical and future solutions in the Eulerian frame. In contrast, we show that representing transport-dominated dynamics in the Lagrangian frame can lead to a significantly faster decay of the Kolmogorov n-width and improve coherence between historical and future solutions. Building on these insights, we develop two non-intrusive ROMs leveraging Lagrangian data: a Lagrangian autoencoder-based ROM and a Lagrangian parametric dynamic mode decomposition. Numerical experiments demonstrate that these Lagrangian ROMs achieve more accurate and stable future predictions than their Eulerian counterparts.

We present a numerical study of eigenvector deflation as a means of accelerating the WaveHoltz method for solving the Helmholtz equation. For energy-conserving (Dirichlet or Neumann) boundary conditions the WaveHoltz fixed-point iteration converges slowly at high frequency, requiring approximately 𝒪(ω^2d) iterations in d dimensions. We show that deflating the eigenvectors whose eigenvalues lie nearest the driving frequency substantially reduces iteration counts, and we examine two ways of incorporating the eigenvectors: direct eigenvector deflation (DEVD), in which the forcing and iterate are projected against the deflation set, and augmented-Krylov eigenvector deflation (AUKED) using deflated conjugate gradient (DCG), augmented GMRES (AGMRES), and augmented (recycled) BICGSTAB (ABICGSTAB). The required eigenpairs can be computed efficiently with the EigenWave approach, and we demonstrate, in two dimensions, that when the number of deflation vectors grows quadratically with ω the asymptotic convergence rate remains essentially constant. Because the eigenvectors on structured grids are naturally represented as matrices, we further apply SVD-based compression to reduce their storage. Numerical experiments on single curvilinear grids discretized with summation-by-parts operators, and on overset grids illustrate the robustness and efficiency of the approach, with the deflated solver breaking even against the undeflated solver after as few as two right-hand sides, when accounting for the cost of precomputing the eigenvectors.

We theoretically analyze the superconvergence of the upwind discontinuous Galerkin (DG) method for both the steady-state and time-dependent radiative transfer equation (RTE), and apply the Smooth-Increasing Accuracy-Conserving (SIAC) filters to enhance the accuracy order. Direct application of SIAC filters on low-dimensional macroscopic moments, often the quantities of practical interest, can effectively improve the approximation accuracy with marginal computational overhead. Using piecewise k-th order polynomials for the approximation and assuming constant cross sections, we prove (2k+2)-th order superconvergence for the steady-state problem at Radau points on each element and (2k+1/2)-th order superconvergence for the global L^2 and negative-order Sobolev norms for the time-dependent problem. Numerical experiments confirm the efficacy of the filtering, demonstrating post-filter convergence orders of 2k+2 for steady-state and 2k+1 for time-dependent problems. More significantly, the SIAC filter delivers substantial gains in computational efficiency. For a time-dependent problem, we observed an approximately 2.22 × accuracy improvement and a 19.94 × reduction in computational time. For the steady-state problems, the filter achieved a 4–9 × acceleration without any loss of accuracy.

The parametric radiative transfer equation (RTE) arises in multi-query applications, such as design optimization, inverse problems, and uncertainty quantification, which require solving the RTE multiple times for various parameters. Classical synthetic acceleration (SA) preconditioners are designed based on low-order approximations of a kinetic correction equation, e.g., its diffusion limit in diffusion synthetic acceleration (DSA). Despite their widespread success, these methods rely on empirical physical assumptions and do not leverage low-rank structures across parameters of the parametric problem. To address these limitations, our previous work introduced a reduced-order model (ROM) enhanced preconditioner called ROMSAD, which exploits low-rank structures across parameters and the original kinetic description of the correction equation. While ROMSAD improves overall efficiency compared with DSA, its efficiency reduces after the first iteration, because the construction of the underlying ROM ignores the preconditioner-dependence of the residual trajectory, leading to a mismatch between the offline and online residual trajectories. To overcome this issue, we introduce a trajectory-aware framework that iteratively constructs ROMs to eliminate the mismatch between offline and online residual trajectories. Numerical tests demonstrate superior efficiency over DSA, and substantial gains in both efficiency and robustness over ROMSAD. For a parametric lattice problem, trajectory-aware ROM preconditioners achieve rapid convergence within only 2-3 iterations online.

In this work, we introduce a novel Hermite method to handle Maxwell's equations for nonlinear dispersive media. The proposed method achieves high-order accuracy and is free of any nonlinear algebraic solver, requiring solving instead small local linear systems for which the dimension is independent of the order. The implementation of order adaptive algorithms is straightforward in this setting, making the resulting p-adaptive Hermite method appealing for the simulations of soliton-like wave propagation.

Nonlinear and nonaffine terms in parametric partial differential equations can potentially lead to a computational cost of a reduced-order model (ROM) that is comparable to the cost of the original full-order model (FOM). To address this, the Reduced Residual Reduced Over-Collocation method (R2-ROC) is developed as a hyper-reduction method within the framework of the reduced basis method in the collocation setting. R2-ROC greedily selects two sets of reduced collocation points based on the (generalized) empirical interpolation method for both solution snapshots and residuals, thereby avoiding the computational inefficiency. The vanilla R2-ROC method can face instability when applied to parametric fluid dynamic problems. To address this, an adaptive enrichment strategy has been proposed to stabilize the ROC method. However, this strategy can involve an excessive number of reduced collocation points, thereby negatively impacting online efficiency. To ensure both efficiency and accuracy, we propose an adaptive time partitioning and adaptive enrichment strategy-based ROC method (AAROC). The adaptive time partitioning dynamically captures the low-rank structure, necessitating fewer reduced collocation points being sampled in each time segment. Numerical experiments on the parametric viscous Burgers' equation and lid-driven cavity problems demonstrate the efficiency, enhanced stability, and accuracy of the proposed AAROC method.

We propose an inexact low-rank source iteration with diffusion synthetic acceleration (SI-DSA) for solving the multidimensional steady-state radiative transfer equation (RTE) in the second-order formulation. The angular flux is represented in either a low-rank matrix or hierarchical Tucker tensor (HTT) format, enabling substantial reductions in computational resources. Each SI step is solved using a preconditioned low-rank conjugate gradient (CG) method with a diffusion preconditioner. To further improve efficiency, we introduce an adaptive inexact strategy that dynamically relaxes the inner CG tolerance during early SI iterations. The method exploits the tensor-product structure of the discretized operators to perform all matrix-vector operations in low-rank form. Numerical experiments on 2D2V benchmark problems, including diffusion-dominated, transport-dominated, and multiscale problems, demonstrate that the proposed approach achieves errors on the order of 10^-4 to 10^-5 relative to full-rank reference solutions, while reducing the degrees of freedom by up to two orders of magnitude. In the diffusion-dominated case, the low-rank solver achieves speedups exceeding 90× over its full-rank counterpart and remains competitive in solving challenging transport-dominated and multiscale problems while providing substantial storage savings. To our knowledge, this work provides the first low-rank SI-DSA framework for multidimensional steady-state RTE.

Parametric radiative transfer equation (RTE) occurs in multi-query applications such as uncertainty quantification, inverse problems, and sensitivity analysis, which require solving RTE multiple times for a range of parameters. Consequently, efficient iterative solvers are highly desired. Classical Synthetic Acceleration (SA) preconditioners for RTE build on low order approximations to an ideal kinetic correction equation such as its diffusion limit in Diffusion Synthetic Acceleration (DSA). Their performance depends on the effectiveness of the underlying low order approximation. In addition, they do not leverage low rank structures with respect to the parameters of the parametric problem. To address these issues, we proposed a ROM-enhanced SA strategy, called ROMSAD, under the Source Iteration framework in Peng (2024). In this paper, we further extend the ROMSAD preconditioner to flexible general minimal residual method (FGMRES). The main new advancement is twofold. First, after identifying the ideal kinetic correction equation within the FGMRES framework, we reformulate it into an equivalent form, allowing us to develop an iterative procedure to construct a ROM for this ideal correction equation without directly solving it. Second, we introduce a greedy algorithm to build the underlying ROM for the ROMSAD preconditioner more efficiently. Our numerical examples demonstrate that FGMRES with the ROMSAD preconditioner (FGMRES-ROMSAD) is more efficient than GMRES with the right DSA preconditioner. Furthermore, when the underlying ROM in ROMSAD is not highly accurate, FGMRES-ROMSAD exhibits greater robustness compared to Source Iteration accelerated by ROMSAD.

In this work, we develop reduced order models (ROMs) to predict solutions to a multiscale kinetic transport equation with a diffusion limit under the parametric setting. When the underlying scattering effect is not sufficiently strong, the system governed by this equation exhibits transport-dominated behavior. Suffering from the Kolmogorov barrier for transport-dominant problems, classical linear ROMs may become inefficient in this regime. To address this issue, we first develop a piecewise linear ROM by introducing a novel goal-oriented adaptive time partitioning strategy. To avoid local over-refinement or under-refinement, we propose an adaptive coarsening and refinement strategy that remains robust with various initial empirical partitions. Additionally, for problems where a local linear approximation is not sufficiently efficient, we further develop a hybrid ROM, which combines autoencoder-based nonlinear ROMs and piecewise linear ROMs. Compared to previous autoencoder-based ROMs, this hybridized method reduces the offline autoencoder's training cost by only applying it to time intervals that are adaptively identified as the most challenging. Numerical experiments demonstrate that our proposed approaches successfully predict full-order solutions at unseen parameter values with both efficiency and accuracy. To the best of our knowledge, this is the first attempt to address the Kolmogorov barrier for multiscale kinetic transport problems with the coexistence of both transport- and diffusion-dominant behaviors.

Applications such as uncertainty quantification, shape optimization, and optical tomography, require solving the radiative transfer equation (RTE) many times for various parameters. Efficient solvers for RTE are highly desired. Source Iteration with Synthetic Acceleration (SISA) is a popular and successful iterative solver for RTE. Synthetic Acceleration (SA) acts as a preconditioning step to accelerate the convergence of Source Iteration (SI). After each source iteration, classical SA strategies introduce a correction to the macroscopic particle density by solving a low order approximation to a kinetic correction equation. For example, Diffusion Synthetic Acceleration (DSA) uses the diffusion limit. However, these strategies may become less effective when the underlying low order approximations are not accurate enough. Furthermore, they do not exploit low rank structures concerning the parameters of parametric problems. To address these issues, we propose enhancing SISA with data-driven ROMs for the parametric problem and the corresponding kinetic correction equation. First, the ROM for the parametric problem can be utilized to obtain an improved initial guess. Second, the ROM for the kinetic correction equation can be utilized to design a low rank approximation to it. Unlike the diffusion limit, this ROM-based approximation builds on the kinetic description of the correction equation and leverages low rank structures concerning the parameters. We further introduce a novel SA strategy called ROMSAD. ROMSAD initially adopts our ROM-based approximation to exploit its greater efficiency in the early stage of SISA, and then automatically switches to DSA to leverage its robustness in the later stage. Additionally, we propose an approach to construct the ROM for the kinetic correction equation without directly solving it. Through a series of numerical tests, we compare the proposed methods with SI-DSA and DSA preconditioned Krylov solver. Particularly, for a multiscale parametric pin-cell problem, ROMSAD achieves approximately 10 times the acceleration compared to SI-DSA and 4 times acceleration compared to DSA preconditioned GMRES.

The Kolmogorov n-width of the solution manifolds of transport-dominated problems can decay slowly. As a result, it can be challenging to design efficient and accurate reduced order models (ROMs) for such problems. To address this issue, we propose a new learning-based projection method to construct nonlinear adaptive ROMs for transport problems. The construction follows the offline-online decomposition. In the offline stage, we train a neural network to construct adaptive reduced basis dependent on time and model parameters. In the online stage, we project the solution to the learned reduced manifold. Inheriting the merits from both deep learning and the projection method, the proposed method is more efficient than the conventional linear projection-based methods, and may reduce the generalization error of a solely learning-based ROM. Unlike some learning-based projection methods, the proposed method does not need to take derivatives of the neural network in the online stage.

We develop a universally applicable embedded boundary finite difference method, which results in a symmetric positive definite linear system and does not suffer from small cell stiffness. Our discretization is efficient for the wave, heat and Poisson equation with Dirichlet boundary conditions. When the system needs to be inverted we can use the conjugate gradient method, accelerated by algebraic multigrid techniques. A series of numerical tests for the wave, heat and Poisson equation and applications to shape optimization problems verify the accuracy, stability, and efficiency of our method. Our fast computational techniques can be extended to moving boundary problems (e.g. Stefan problem), to the Navier–Stokes equations, and to the Grad-Shafranov equations for which problems are posed on domains with complex geometry and fast simulations are of great interest.

Quantum computing has received significant amounts of interest from many different research communities over the last few years. Although there are many introductory texts that focus on the algorithmic parts of quantum computing, there is a dearth of publications that describe the modeling, calibration and operation of current quantum computing devices. One aim of this report is to fill that void by providing a case study that walks through the entire procedure from the characterization and optimal control of a qudit device at Lawrence Livermore National Laboratory (LLNL) to the validation of the results. A goal of the report is to provide an introduction for students and researchers, especially computational mathematicians, who are interested in but new to quantum computing. Both experimental and mathematical aspects of this procedure are discussed. We present a description of the LLNL QuDIT testbed, the mathematical models that are used to describe it, and the numerical methods that are used to to design optimal controls. We also present experimental and computational methods that can be used to characterize a quantum device. Finally, an experimental validation of an optimized control pulse is presented, which relies on the accuracy of the characterization and the optimal control methodologies.

In the higher levels of superconducting transmon devices, and more generally charge sensitive devices, $T_2^*$ measurements made in the presence of low-frequency time-correlated $1/f$ charge noise and quasiparticle-induced parity flips can give an underestimation of the total dephasing time. The charge variations manifest as beating patterns observed in the overlay of several Ramsey fringe curves, and are reproduced with a phenomenological Ramsey curve model which accounts for the charge variations. $T_2^*$ dephasing times which more accurately represent the total dephasing time are obtained. The phenomenological model is compared with a Lindblad master equation model. Both models are found to be in agreement with one another and the experimental data. Finally, the phenomenological formulation enables a simple method in which the power spectral density (PSD) for the low-frequency noise can be inferred from the overlay of several Ramsey curves.

Motivated by the noisy and fluctuating behavior of current quantum computing devices, this paper presents a data-driven characterization approach for estimating transition frequencies and decay times in a Lindbladian dynamical model of a superconducting quantum device. The data includes parity events in the transition frequency between the first and second excited states. A simple but effective mathematical model, based upon averaging solutions of two Lindbladian models, is demonstrated to accurately capture the experimental observations. A deterministic point estimate of the device parameters is first performed to minimize the misfit between data and Lindbladian simulations. These estimates are used to make an informed choice of prior distributions for the subsequent Bayesian inference. An additive Gaussian noise model is developed for the likelihood function, which includes two hyper-parameters to capture the noise structure of the data. The outcome of the Bayesian inference are posterior probability distributions of the transition frequencies, which for example can be utilized to design risk neutral optimal control pulses. The applicability of our approach is demonstrated on experimental data from the Quantum Device and Integration Testbed (QuDIT) at Lawrence Livermore National Laboratory, using a tantalum-based superconducting transmon device.

By defining spaces and differential operators on mesh surfaces, we extend the total generalized variation (TGV) model on the 2-dimensional space to triangulated surfaces. Based on the new definition of TGV model on triangulated surfaces, we introduce total generalized variation restoration optimization problem for data (image pixel/surface normal) restoration over triangulated surfaces. The optimization problem is solved effectively by augmented Lagrangian method (ALM). Closed form solutions for subproblems of the ALM method are obtained. Convergence analysis of the ALM algorithm is presented. Through series of numerical experiments, we show that the TGV method can alleviate the staircase effect and recover more structures and details. As a result, the TGV model outperforms several existing models visually and quantitatively. The robustness of the TGV method is also confirmed numerically.

A novel approach to computing time-harmonic solutions of Maxwell's equations by time-domain simulations is presented.The method, EM-WaveHoltz, results in a positive definite system of equations which makes it amenable to iterative solution with the conjugate gradient method or with GMRES.Theoretical results guaranteeing the convergence of the method away from resonances are presented.Numerical examples illustrating the properties of EM-WaveHoltz are given.

Linear kinetic transport equations play a critical role in optical tomography, radiative transfer and neutron transport. The fundamental difficulty hampering their efficient and accurate numerical resolution lies in the high dimensionality of the physical and velocity/angular variables and the fact that the problem is multiscale in nature. Leveraging the existence of a hidden low-rank structure hinted by the diffusive limit, in this work, we design and test the angular-space reduced order model for the linear radiative transfer equation, the first such effort based on the celebrated reduced basis method (RBM). Our method is built upon a high-fidelity solver employing the discrete ordinates method in the angular space, an asymptotic preserving upwind discontinuous Galerkin method for the physical space, and an efficient synthetic accelerated source iteration for the resulting linear system. Addressing the challenge of the parameter values (or angular directions) being coupled through an integration operator, the first novel ingredient of our method is an iterative procedure where the macroscopic density is constructed from the RBM snapshots, treated explicitly and allowing a transport sweep, and then updated afterwards. A greedy algorithm can then proceed to adaptively select the representative samples in the angular space and form a surrogate solution space. The second novelty is a least-squares density reconstruction strategy, at each of the relevant physical locations, enabling the robust and accurate integration over an arbitrarily unstructured set of angular samples toward the macroscopic density. Numerical experiments indicate that our method is effective for computational cost reduction in a variety of regimes.

We present EM-WaveHoltz, a new method for computing frequency domain solutions to Maxwell's equations in linear complex media. EM-WaveHoltz is based on a fixed point iteration which is applied to the initial data of a time domain solution. The fixed point iteration converges to a single time-harmonic component and can be formulated as a symmetric positive definite linear system of equations. The method is agnostic to the material model and to the time-domain solver used.