Search

In recent years, there has been a renewed interest in preconditioning for multilevel Toeplitz systems, a research field that has been extensively explored over the past several decades. This work introduces novel preconditioning strategies using multilevel $\tau$ matrices for both symmetric and nonsymmetric multilevel Toeplitz systems. Our proposals constitute a general framework, as they are constructed solely based on the generating function of the multilevel Toeplitz coefficient matrix, when it can be defined. We begin with nonsymmetric systems, where we employ a symmetrization technique by permuting the coefficient matrix to produce a real symmetric multilevel Hankel structure. We propose a multilevel $\tau$ preconditioner tailored to the symmetrized system and prove that the eigenvalues of the preconditioned matrix sequence cluster at $\pm 1$, leading to rapid convergence when using the preconditioned minimal residual method. The high effectiveness of this approach is demonstrated through its application in solving space fractional diffusion equations. Next, for symmetric systems we introduce another multilevel $\tau$ preconditioner and show that the preconditioned conjugate gradient method can achieve an optimal convergence rate, namely a rate that is independent of the matrix size, when employed for a class of ill-conditioned multilevel Toeplitz systems. Numerical examples are provided to critically assess the effectiveness of our proposed preconditioners compared to several leading existing preconditioned solvers, highlighting their superior performance.

Variational models for image deblurring problems typically consist of a smooth term and a potentially non-smooth convex term. A common approach to solving these problems is using proximal gradient methods. To accelerate the convergence of these first-order iterative algorithms, strategies such as variable metric methods have been introduced in the literature. In this paper, we prove that, for image deblurring problems, the variable metric strategy can be reinterpreted as a right preconditioning method. Consequently, we explore an inexact left-preconditioned version of the same proximal gradient method. We prove the convergence of the new iteration to the minimum of a variational model where the norm of the data fidelity term depends on the preconditioner. The numerical results show that left and right preconditioning are comparable in terms of the number of iterations required to reach a prescribed tolerance, but left preconditioning needs much less CPU time, as it involves fewer evaluations of the preconditioner matrix compared to right preconditioning. The quality of the computed solutions with left and right preconditioning are comparable. Finally, we propose some non-stationary sequences of preconditioners that allow for fast and stable convergence to the solution of the variational problem with the classical $\ell^2$--norm on the fidelity term.

We introduce a two-level trust-region method (TLTR) for solving unconstrained nonlinear optimization problems. Our method uses a composite iteration step, which is based on two distinct search directions. The first search direction is obtained through minimization in the full/high-resolution space, ensuring global convergence to a critical point. The second search direction is obtained through minimization in the randomly generated subspace, which, in turn, allows for convergence acceleration. The efficiency of the proposed TLTR method is demonstrated through numerical experiments in the field of machine learning

Geological applications of phase-field methods for fracture are notably scarce. This work conducts a numerical examination of the applicability of standard phase-field models in reproducing jointing within sedimentary layers. We explore how the volumetric-deviatoric split alongside the AT1 and AT2 phase-field formulations have several advantages in simulating jointing, but also have intrinsic limitations that prevent a reliable quantitative analysis of rock fracture. The formulations qualitatively reproduce the process of joint saturation, along with the negative correlation between joint spacing and the height of the sedimentary layer. However, in quantitative comparison to alternative numerical methods and outcrop observations, the phase-field method overestimates joint spacings by a factor of 2 and induces unrealistic compressive fractures in the AT1 model, alongside premature shearing at layer interfaces for the AT2 model. The causes are identified to be intrinsic to the phase-field lengthscale and the unsuitable strength envelope arising from the Volumetric-Deviatoric split. Finally, our analysis elucidates on the phase-field lengthscale's distortion of the stress field around dilating fractures, causing the sedimentary layer to reach joint saturation prematurely, thereby stopping the nucleation of new fractures and leading to larger joint spacings than in natural examples. Decreasing the lengthscale results in gradual improvement but becomes prohibitively computationally expensive for very small lengthscales such that the limit of ``natural'' behavior was not reached in this study. Overall, our results motivate the development of constitutive phase-field models that are more suitable for geological applications and their benchmarking against geological observations., Comment: 24 pages

Simulation of the monodomain equation, crucial for modeling the heart's electrical activity, faces scalability limits when traditional numerical methods only parallelize in space. To optimize the use of large multi-processor computers by distributing the computational load more effectively, time parallelization is essential. We introduce a high-order parallel-in-time method addressing the substantial computational challenges posed by the stiff, multiscale, and nonlinear nature of cardiac dynamics. Our method combines the semi-implicit and exponential spectral deferred correction methods, yielding a hybrid method that is extended to parallel-in-time employing the PFASST framework. We thoroughly evaluate the stability, accuracy, and robustness of the proposed parallel-in-time method through extensive numerical experiments, using practical ionic models such as the ten-Tusscher-Panfilov. The results underscore the method's potential to significantly enhance real-time and high-fidelity simulations in biomedical research and clinical applications.

Due to its optimal complexity, the multigrid (MG) method is one of the most popular approaches for solving large-scale linear systems arising from the discretization of partial differential equations. However, the parallel implementation of standard MG methods, which are inherently multiplicative, suffers from increasing communication complexity. In such cases, the additive variants of MG methods provide a good alternative due to their inherently parallel nature, although they exhibit slower convergence. This work combines the additive multigrid method with the multipreconditioned conjugate gradient (MPCG) method. In the proposed approach, the MPCG method employs the corrections from the different levels of the MG hierarchy as separate preconditioned search directions. In this approach, the MPCG method updates the current iterate by using the linear combination of the preconditioned search directions, where the optimal coefficients for the linear combination are computed by exploiting the energy norm minimization of the CG method. The idea behind our approach is to combine the $A$-conjugacy of the search directions of the MPCG method and the quasi $H_1$-orthogonality of the corrections from the MG hierarchy. In the numerical section, we study the performance of the proposed method compared to the standard additive and multiplicative MG methods used as preconditioners for the CG method.

The current study investigates the asymptotic spectral properties of a finite difference approximation of nonlocal Helmholtz equations with a Caputo fractional Laplacian and a variable coefficient wave number $\mu$, as it occurs when considering a wave propagation in complex media, characterized by nonlocal interactions and spatially varying wave speeds. More specifically, by using tools from Toeplitz and generalized locally Toeplitz theory, the present research delves into the spectral analysis of nonpreconditioned and preconditioned matrix-sequences. We report numerical evidences supporting the theoretical findings. Finally, open problems and potential extensions in various directions are presented and briefly discussed., Comment: 28 pages, 10 figures. arXiv admin note: text overlap with arXiv:2206.05171 by other authors

Fully explicit stabilized multirate (mRKC) methods are well-suited for the numerical solution of large multiscale systems of stiff ordinary differential equations thanks to their improved stability properties. To demonstrate their efficiency for the numerical solution of stiff, multiscale, nonlinear parabolic PDE's, we apply mRKC methods to the monodomain equation from cardiac electrophysiology. In doing so, we propose an improved version, specifically tailored to the monodomain model, which leads to the explicit exponential multirate stabilized (emRKC) method. Several numerical experiments are conducted to evaluate the efficiency of both mRKC and emRKC, while taking into account different finite element meshes (structured and unstructured) and realistic ionic models. The new emRKC method typically outperforms a standard implicit-explicit baseline method for cardiac electrophysiology. Code profiling and strong scalability results further demonstrate that emRKC is faster and inherently parallel without sacrificing accuracy.

Twenty years ago the anti-reflective numerical boundary conditions (BCs) were introduced in a context of signal processing and imaging, for increasing the quality of the reconstruction of a blurred signal/image contaminated by noise and for reducing the overall complexity to that of few fast sine transforms i.e. to $O(N\log N)$ real arithmetic operations, where $N$ is the number of pixels. Here for quality of reconstruction we mean a better accuracy and the elimination of boundary artifacts, called ringing effects. Now we propose numerical anti-reflective BCs in the context of nonlocal problems of fractional differential type: the goals are the same i.e. a smaller approximation error and the reduction of boundary artifacts. In the latter setting, we compare various types of numerical BCs, including the anti-symmetric ones considered in the case of fractional differential problems for modeling reasons. More in detail, given important similarities between anti-symmetric and anti-reflective BCs, we compare them from the perspective of computational efficiency, by considering nontruncated and truncated versions and also other standard numerical BCs such as reflective/Neumann. Several numerical tests, tables, and visualizations are provided and critically discussed. The conclusion is that the truncated numerical anti-reflective BCs perform better, both in terms of accuracy and low computational cost., Comment: 33 pages, 11 figures

We propose to train neural networks (NNs) using a novel variant of the ``Additively Preconditioned Trust-region Strategy'' (APTS). The proposed method is based on a parallelizable additive domain decomposition approach applied to the neural network's parameters. Built upon the TR framework, the APTS method ensures global convergence towards a minimizer. Moreover, it eliminates the need for computationally expensive hyper-parameter tuning, as the TR algorithm automatically determines the step size in each iteration. We demonstrate the capabilities, strengths, and limitations of the proposed APTS training method by performing a series of numerical experiments. The presented numerical study includes a comparison with widely used training methods such as SGD, Adam, LBFGS, and the standard TR method.

We propose to enhance the training of physics-informed neural networks (PINNs). To this aim, we introduce nonlinear additive and multiplicative preconditioning strategies for the widely used L-BFGS optimizer. The nonlinear preconditioners are constructed by utilizing the Schwarz domain-decomposition framework, where the parameters of the network are decomposed in a layer-wise manner. Through a series of numerical experiments, we demonstrate that both, additive and multiplicative preconditioners significantly improve the convergence of the standard L-BFGS optimizer, while providing more accurate solutions of the underlying partial differential equations. Moreover, the additive preconditioner is inherently parallel, thus giving rise to a novel approach to model parallelism., Comment: 23 pages, 7 figures

We propose a boundary element method for the accurate solution of the cell-by-cell bidomain model of electrophysiology. The cell-by-cell model, also called Extracellular-Membrane-Intracellular (EMI) model, is a system of reaction-diffusion equations describing the evolution of the electric potential within each domain: intra- and extra-cellular space and the cellular membrane. The system is parabolic but degenerate because the time derivative is only in the membrane domain. In this work, we adopt a boundary-integral formulation for removing the degeneracy in the system and recast it to a parabolic equation on the membrane. The formulation is also numerically advantageous since the number of degrees of freedom is sensibly reduced compared to the original model. Specifically, we prove that the boundary-element discretization of the EMI model is equivalent to a system of ordinary differential equations, and we consider a time discretization based on the multirate explicit stabilized Runge-Kutta method. We numerically show that our scheme convergences exponentially in space for the single-cell case. We finally provide several numerical experiments of biological interest., Comment: 20 pages, 9 figures, 1 table

We propose a nonlinear additive Schwarz method for solving nonlinear optimization problems with bound constraints. Our method is used as a "right-preconditioner" for solving the first-order optimality system arising within the sequential quadratic programming (SQP) framework using Newton's method. The algorithmic scalability of this preconditioner is enhanced by incorporating a solution-dependent coarse space, which takes into account the restricted constraints from the fine level. By means of numerical examples, we demonstrate that the proposed preconditioned Newton methods outperform standard active-set methods considered in the literature.

It is known that the solution of a conservative steady-state two-sided fractional diffusion problem can exhibit singularities near the boundaries. As consequence of this, and due to the conservative nature of the problem, we adopt a finite volume elements discretization approach over a generic non-uniform mesh. We focus on grids mapped by a smooth function which consist in a combination of a graded mesh near the singularity and a uniform mesh where the solution is smooth. Such a choice gives rise to Toeplitz-like discretization matrices and thus allows a low computational cost of the matrix-vector product and a detailed spectral analysis. The obtained spectral information is used to develop an ad-hoc parameter free multigrid preconditioner for GMRES, which is numerically shown to yield good convergence results in presence of graded meshes mapped by power functions that accumulate points near the singularity. The approximation order of the considered graded meshes is numerically compared with the one of a certain composite mesh given in literature that still leads to Toeplitz-like linear systems and is then still well-suited for our multigrid method. Several numerical tests confirm that power graded meshes result in lower approximation errors than composite ones and that our solver has a wide range of applicability.

One of the state-of-the-art strategies for predicting crack propagation, nucleation, and interaction is the phase-field approach. Despite its reliability and robustness, the phase-field approach suffers from burdensome computational cost, caused by the non-convexity of the underlying energy functional and a large number of unknowns required to resolve the damage gradients. In this work, we propose to solve such nonlinear systems in a monolithic manner using the Schwarz preconditioned inexact Newton's (SPIN) method. The proposed SPIN method leverages the field split approach and minimizes the energy functional separately with respect to displacement and the phase-field, in an additive and multiplicative manner. In contrast to the standard alternate minimization, the result of this decoupled minimization process is used to construct a preconditioner for a coupled linear system, arising at each Newton's iteration. The overall performance and the convergence properties of the proposed additive and multiplicative SPIN methods are investigated by means of several numerical examples. Comparison with widely-used alternate minimization is also performed and we show a reduction in the execution time up to a factor of 50. Moreover, we also demonstrate that this reduction grows even further with increasing problem size and larger loading increments.

Context. Numerical solutions to transfer problems of polarized radiation in solar and stellar atmospheres commonly rely on stationary iterative methods, which often perform poorly when applied to large problems. In recent times, stationary iterative methods have been replaced by state-of-the-art preconditioned Krylov iterative methods for many applications. However, a general description and a convergence analysis of Krylov methods in the polarized radiative transfer context are still lacking. Aims. We describe the practical application of preconditioned Krylov methods to linear transfer problems of polarized radiation, possibly in a matrix-free context. The main aim is to clarify the advantages and drawbacks of various Krylov accelerators with respect to stationary iterative methods. Methods. We report the convergence rate and the run time of various Krylov-accelerated techniques combined with different formal solvers when applied to a 1D benchmark transfer problem of polarized radiation. In particular, we analyze the GMRES, BICGSTAB, and CGS Krylov methods, preconditioned with Jacobi, or (S)SOR. Results. Krylov methods accelerate the convergence, reduce the run time, and improve the robustness of standard stationary iterative methods. Jacobi-preconditioned Krylov methods outperform SOR-preconditioned stationary iterations in all respects. In particular, the Jacobi-GMRES method offers the best overall performance for the problem setting in use. Conclusions. Krylov methods can be more challenging to implement than stationary iterative methods. However, an algebraic formulation of the radiative transfer problem allows one to apply and study Krylov acceleration strategies with little effort. Furthermore, many available numerical libraries implement matrix-free Krylov routines, enabling an almost effortless transition to Krylov methods.

In this paper, we investigate the combination of multigrid methods and neural networks, starting from a Finite Element discretization of an elliptic PDE. Multigrid methods use interpolation operators to transfer information between different levels of approximation. These operators are crucial for fast convergence of multigrid, but they are generally unknown. We propose Deep Neural Network models for learning interpolation operators and we build a multilevel hierarchy based on the output of the network. We investigate the accuracy of the interpolation operator predicted by the Neural Network, testing it with different network architectures. This Neural Network approach for the construction of grid operators can then be extended for an automatic definition of multilevel solvers, allowing a portable solution in scientific computing, Comment: To appear in Springer Journal: "The 26th International Domain Decomposition Conference (DD26)"

We propose a globally convergent multilevel training method for deep residual networks (ResNets). The devised method can be seen as a novel variant of the recursive multilevel trust-region (RMTR) method, which operates in hybrid (stochastic-deterministic) settings by adaptively adjusting mini-batch sizes during the training. The multilevel hierarchy and the transfer operators are constructed by exploiting a dynamical system's viewpoint, which interprets forward propagation through the ResNet as a forward Euler discretization of an initial value problem. In contrast to traditional training approaches, our novel RMTR method also incorporates curvature information on all levels of the multilevel hierarchy by means of the limited-memory SR1 method. The overall performance and the convergence properties of our multilevel training method are numerically investigated using examples from the field of classification and regression.

Internal interfaces in a domain could exist as a material defect or they can appear due to propagations of cracks. Discretization of such geometries and solution of the contact problem on the internal interfaces can be computationally challenging. We employ an unfitted Finite Element (FE) framework for the discretization of the domains and develop a tailored, globally convergent, and efficient multigrid method for solving contact problems on the internal interfaces. In the unfitted FE methods, structured background meshes are used and only the underlying finite element space has to be modified to incorporate the discontinuities. The non-penetration conditions on the embedded interfaces of the domains are discretized using the method of Lagrange multipliers. We reformulate the arising variational inequality problem as a quadratic minimization problem with linear inequality constraints. Our multigrid method can solve such problems by employing a tailored multilevel hierarchy of the FE spaces and a novel approach for tackling the discretized non-penetration conditions. We employ pseudo-$L^2$ projection-based transfer operators to construct a hierarchy of nested FE spaces from the hierarchy of non-nested meshes. The essential component of our multigrid method is a technique that decouples the linear constraints using an orthogonal transformation of the basis. The decoupled constraints are handled by a modified variant of the projected Gauss-Seidel method, which we employ as a smoother in the multigrid method. These components of the multigrid method allow us to enforce linear constraints locally and ensure the global convergence of our method. We will demonstrate the robustness, efficiency, and level independent convergence property of the proposed method for Signorini's problem and two-body contact problems.

The dual formulation for linear elasticity, in contrast to the primal formulation, is not affected by locking, as it is based on the stresses as main unknowns. Thus it is quite attractive for nearly incompressible and incompressible materials. Discretization with mixed finite elements will lead to -- possibly large -- linear saddle point systems with a particular structure. Whereas efficient multigrid methods exist for solving problems in mixed plane elasticity, to the knowledge of the authors, no multigrid methods are readily available for the general dual formulation. Two are the main challenges in constructing a multigrid method for the dual formulation for linear elasticity. First, in the incompressible limit, the matrix block related to the stress is semi-positive definite. Second, the stress belongs to $\textbf{H}_{\text{div}}$ and standard smoothers, working for $\textbf{H}^1$ regular problems, cannot be applied. We present a novel patch-based smoother for the dual formulation for linear elasticity. We discuss different types of local boundary conditions for the patch subproblems. Based on our patch-smoother, we build a multigrid method for the solution of the resulting saddle point problem and investigate its efficiency and robustness. Numerical experiments show that Robin conditions best fit the multigrid framework, leading eventually to multigrid performance., Comment: 27 pages, 88 figures

We present a novel approach which aims at high-performance uncertainty quantification for cardiac electrophysiology simulations. Employing the monodomain equation to model the transmembrane potential inside the cardiac cells, we evaluate the effect of spatially correlated perturbations of the heart fibers on the statistics of the resulting quantities of interest. Our methodology relies on a close integration of multilevel quadrature methods, parallel iterative solvers and space-time finite element discretizations, allowing for a fully parallelized framework in space, time and stochastics. Extensive numerical studies are presented to evaluate convergence rates and to compare the performance of classical Monte Carlo methods such as standard Monte Carlo (MC) and quasi-Monte Carlo (QMC), as well as multilevel strategies, i.e. multilevel Monte Carlo (MLMC) and multilevel quasi-Monte Carlo (MLQMC) on hierarchies of nested meshes. Finally, we employ a recently suggested variant of the multilevel approach for non-nested meshes to deal with a realistic heart geometry.

The parallel solution of large scale non-linear programming problems, which arise for example from the discretization of non-linear partial differential equations, is a highly demanding task. Here, a novel solution strategy is presented, which is inherently parallel and globally convergent. Each global non-linear iteration step consists of asynchronous solutions of local non-linear programming problems followed by a global recombination step. The recombination step, which is the solution of a quadratic programming problem, is designed in a way such that it ensures global convergence. As it turns out, the new strategy can be considered as a globalized additively preconditioned inexact Newton (ASPIN) method. However, in our approach the influence of ASPIN's non-linear preconditioner on the gradient is controlled in order to ensure a sufficient decrease condition. Two different control strategies are described and analyzed. Convergence to first-order critical points of our non-linear solution strategy is shown under standard trust-region assumptions. The strategy is investigated along difficult minimization problems arising from non-linear elasticity in 3D solved on a massively parallel computer with several thousand cores.

We introduce a novel variant of the recursive multilevel trust-region (RMTR) method, called MASTR. The method is designed for solving non-convex bound-constrained minimization problems, which arise from the finite element discretization of partial differential equations. MASTR utilizes an active-set strategy based on the truncated basis approach in order to preserve the variable bounds defined on the finest level by the coarser levels. Usage of this approach allows for fast convergence of the MASTR method, especially once the exact active-set is detected. The efficiency of the method is demonstrated by means of numerical examples.

In this work, we propose a multigrid preconditioner for Jacobian-free Newton-Krylov (JFNK) methods. Our multigrid method does not require knowledge of the Jacobian at any level of the multigrid hierarchy. As it is common in standard multigrid methods, the proposed method also relies on three building blocks: transfer operators, smoothers, and a coarse level solver. In addition to the restriction and prolongation operator, we also use a projection operator to transfer the current Newton iterate to a coarser level. The three-level Chebyshev semi-iterative method is employed as a smoother, as it has good smoothing properties and does not require the representation of the Jacobian matrix. We replace the direct solver on the coarsest level with a matrix-free Krylov subspace method, thus giving rise to a truly Jacobian-free multigrid preconditioner. We will discuss all building blocks of our multigrid preconditioner in detail and demonstrate the robustness and the efficiency of the proposed method using several numerical examples.

The identification of the initial ventricular activation sequence is a critical step for the correct personalization of patient-specific cardiac models. In healthy conditions, the Purkinje network is the main source of the electrical activation, but under pathological conditions the so-called earliest activation sites (EASs) are possibly sparser and more localized. Yet, their number, location and timing may not be easily inferred from remote recordings, such as the epicardial activation or the 12-lead electrocardiogram (ECG), due to the underlying complexity of the model. In this work, we introduce GEASI (Geodesic-based Earliest Activation Sites Identification) as a novel approach to simultaneously identify all EASs. To this end, we start from the anisotropic eikonal equation modeling cardiac electrical activation and exploit its Hamilton--Jacobi formulation to minimize a given objective function, e.g. the quadratic mismatch to given activation measurements. This versatile approach can be extended to estimate the number of activation sites by means of the topological gradient, or fitting a given ECG. We conducted various experiments in 2D and 3D for in-silico models and an in-vivo intracardiac recording collected from a patient undergoing cardiac resynchronization therapy. The results demonstrate the clinical applicability of GEASI for potential future personalized models and clinical intervention., Comment: 38 pages, 17 figures

In this work, we consider the modeling of inclusions in the material using an unfitted finite element method. In the unfitted methods, structured background meshes are used and only the underlying finite element space is modified to incorporate the discontinuities, such as inclusions. Hence, the unfitted methods provide a more flexible framework for modeling the materials with multiple inclusions. We employ the method of Lagrange multipliers for enforcing the interface conditions between the inclusions and matrix, this gives rise to the linear system of equations of saddle point type. We utilize the Uzawa method for solving the saddle point system and propose preconditioning strategies for primal and dual systems. For the dual systems, we review and compare the preconditioning strategies that are developed for FETI and SIMPLE methods. While for the primal system, we employ a tailored multigrid method specifically developed for the unfitted meshes. Lastly, the comparison between the proposed preconditioners is made through several numerical experiments.

In electrocardiography, the "classic" inverse problem is the reconstruction of electric potentials at a surface enclosing the heart from remote recordings at the body surface and an accurate description of the anatomy. The latter being affected by noise and obtained with limited resolution due to clinical constraints, a possibly large uncertainty may be perpetuated in the inverse reconstruction. The purpose of this work is to study the effect of shape uncertainty on the forward and the inverse problem of electrocardiography. To this aim, the problem is first recast into a boundary integral formulation and then discretised with a collocation method to achieve high convergence rates and a fast time to solution. The shape uncertainty of the domain is represented by a random deformation field defined on a reference configuration. We propose a periodic-in-time covariance kernel for the random field and approximate the Karhunen-Lo\`eve expansion using low-rank techniques for fast sampling. The space-time uncertainty in the expected potential and its variance is evaluated with an anisotropic sparse quadrature approach and validated by a quasi-Monte Carlo method. We present several numerical experiments on a simplified but physiologically grounded 2-dimensional geometry to illustrate the validity of the approach. The tested parametric dimension ranged from 100 up to 600. For the forward problem the sparse quadrature is very effective. In the inverse problem, the sparse quadrature and the quasi-Monte Carlo method perform as expected, except for the total variation regularisation, where convergence is limited by lack of regularity. We finally investigate an $H^{1/2}$ regularisation, which naturally stems from the boundary integral formulation, and compare it to more classical approaches., Comment: 23 pages, 6 figures, 1 table

Non-linear phase field models are increasingly used for the simulation of fracture propagation models. The numerical simulation of fracture networks of realistic size requires the efficient parallel solution of large coupled non-linear systems. Although in principle efficient iterative multi-level methods for these types of problems are available, they are not widely used in practice due to the complexity of their parallel implementation. Here, we present Utopia, which is an open-source C++ library for parallel non-linear multilevel solution strategies. Utopia provides the advantages of high-level programming interfaces while at the same time a framework to access low-level data-structures without breaking code encapsulation. Complex numerical procedures can be expressed with few lines of code, and evaluated by different implementations, libraries, or computing hardware. In this paper, we investigate the parallel performance of our implementation of the recursive multilevel trust-region (RMTR) method based on the Utopia library. RMTR is a globally convergent multilevel solution strategy designed to solve non-convex constrained minimization problems. In particular, we solve pressure-induced phase-field fracture propagation in large and complex fracture networks. Solving such problems is deemed challenging even for a few fractures, however, here we are considering networks of realistic size with up to 1000 fractures., Comment: CCF Trans. HPC (2021)

We consider two parallel-in-time approaches applied to a (reaction) diffusion problem, possibly non-linear. In particular, we consider PFASST (Parallel Full Approximation Scheme in Space and Time) and space-time multilevel strategies. For both approaches, we start from an integral formulation of the continuous time-dependent problem. Then, a collocation form for PFASST and a discontinuous Galerkin discretization in time for the space-time multigrid are employed, resulting in the same discrete solution at the time nodes. Strong and weak scaling of both multilevel strategies is compared for varying order of the temporal discretization. Moreover, we investigate the respective convergence behavior for non-linear problems and highlight quantitative differences.

We present a tailored multigrid method for linear problems stemming from a Nitsche-based extended finite element method (XFEM). Our multigrid method is robust with respect to highly varying coefficients and the number of interfaces in a domain. It shows level independent convergence rates when applied to different variants of Nitsche's method. Generally, multigrid methods require a hierarchy of finite element (FE) spaces which can be created geometrically using a hierarchy of nested meshes. However, in the XFEM framework, standard multigrid methods might demonstrate poor convergence properties if the hierarchy of FE spaces employed is not nested. We design a prolongation operator for the multigrid method in such a way that it can accommodate the discontinuities across the interfaces in the XFEM framework and recursively induces a nested FE space hierarchy. The prolongation operator is constructed using so-called pseudo-$L^2$-projections; as common, the adjoint of the prolongation operator is employed as the restriction operator. The stabilization parameter in Nitsche's method plays an important role in imposing interface conditions and also affects the condition number of the linear systems. We discuss the requirements on the stabilization parameter to ensure coercivity and review selected strategies from the literature which are used to implicitly estimate the stabilization parameter. Eventually, we compare the impact of different variations of Nitsche's method on discretization errors and condition number of the linear systems. We demonstrate the robustness of our multigrid method with respect to varying coefficients and the number of interfaces and compare it with other preconditioners., Comment: 21 pages, 6 figures

We present a novel approach aimed at high-performance uncertainty quantification for time-dependent problems governed by partial differential equations. In particular, we consider input uncertainties described by a Karhunen-Loeeve expansion and compute statistics of high-dimensional quantities-of-interest, such as the cardiac activation potential. Our methodology relies on a close integration of multilevel Monte Carlo methods, parallel iterative solvers, and a space-time discretization. This combination allows for space-time adaptivity, time-changing domains, and to take advantage of past samples to initialize the space-time solution. The resulting sequence of problems is distributed using a multilevel parallelization strategy, allocating batches of samples having different sizes to a different number of processors. We assess the performance of the proposed framework by showing in detail its application to the solution of nonlinear equations arising from cardiac electrophysiology. Specifically, we study the effect of spatially-correlated perturbations of the heart fibers conductivities on the mean and variance of the resulting activation map. As shown by the experiments, the theoretical rates of convergence of multilevel Monte Carlo are achieved. Moreover, the total computational work for a prescribed accuracy is reduced by an order of magnitude with respect to standard Monte Carlo methods.

The simulation of crack initiation and propagation in an elastic material is difficult, as crack paths with complex topologies have to be resolved. Phase-field approach allows to simulate crack behavior by circumventing the need to explicitly model crack paths. However, the underlying mathematical model gives rise to a non-convex constrained minimization problem. In this work, we propose a recursive multilevel trust region (RMTR) method to efficiently solve such a minimization problem. The RMTR method combines the global convergence property of the trust region method and the optimality of the multilevel method. The solution process is accelerated by employing level dependent objective functions, minimization of which provides correction to the original/fine-level problem. In the context of the phase-field fracture approach, it is challenging to design efficient level dependent objective functions as the underlying mathematical model relies on the mesh dependent parameters. We introduce level dependent objective functions that combine fine level description of the crack path with the coarse level discretization. The overall performance and the convergence properties of the proposed RMTR method are investigated by means of several numerical examples in three dimensions.

Multifidelity Monte Carlo methods rely on a hierarchy of possibly less accurate but statistically correlated simplified or reduced models, in order to accelerate the estimation of statistics of high-fidelity models without compromising the accuracy of the estimates. This approach has recently gained widespread attention in uncertainty quantification. This is partly due to the availability of optimal strategies for the estimation of the expectation of scalar quantities-of-interest. In practice, the optimal strategy for the expectation is also used for the estimation of variance and sensitivity indices. However, a general strategy is still lacking for vector-valued problems, nonlinearly statistically-dependent models, and estimators for which a closed-form expression of the error is unavailable. The focus of the present work is to generalize the standard multifidelity estimators to the above cases. The proposed generalized estimators lead to an optimization problem that can be solved analytically and whose coefficients can be estimated numerically with few runs of the high- and low-fidelity models. We analyze the performance of the proposed approach on a selected number of experiments, with a particular focus on cardiac electrophysiology, where a hierarchy of physics-based low-fidelity models is readily available., Comment: Manuscript submitted to Journal of Computational Physics

This work describes a domain embedding technique between two non-matching meshes used for generating realizations of spatially correlated random fields with applications to large-scale sampling-based uncertainty quantification. The goal is to apply the multilevel Monte Carlo (MLMC) method for the quantification of output uncertainties of PDEs with random input coefficients on general, unstructured computational domains. We propose a highly scalable, hierarchical sampling method to generate realizations of a Gaussian random field on a given unstructured mesh by solving a reaction-diffusion PDE with a stochastic right-hand side. The stochastic PDE is discretized using the mixed finite element method on an embedded domain with a structured mesh, and then the solution is projected onto the unstructured mesh. This work describes implementation details on how to efficiently transfer data from the structured and unstructured meshes at coarse levels, assuming this can be done efficiently on the finest level. We investigate the efficiency and parallel scalability of the technique for the scalable generation of Gaussian random fields in three dimensions. An application of the MLMC method is presented for quantifying uncertainties of subsurface flow problems. We demonstrate the scalability of the sampling method with non-matching mesh embedding, coupled with a parallel forward model problem solver, for large-scale 3D MLMC simulations with up to $1.9\cdot 10^9$ unknowns., Comment: To appear in Numerical Linear Algebra with Applications

In this paper we combine the Parareal parallel-in-time method together with spatial parallelization and investigate this space-time parallel scheme by means of solving the three-dimensional incompressible Navier-Stokes equations. Parallelization of time stepping provides a new direction of parallelization and allows to employ additional cores to further speed up simulations after spatial parallelization has saturated. We report on numerical experiments performed on a Cray XE6, simulating a driven cavity flow with and without obstacles. Distributed memory parallelization is used in both space and time, featuring up to 2,048 cores in total. It is confirmed that the space-time-parallel method can provide speedup beyond the saturation of the spatial parallelization.

We present a new class of iterative schemes for solving initial value problems (IVP) based on discontinuous Galerkin (DG) methods. Starting from the weak DG formulation of an IVP, we derive a new iterative method based on a preconditioned Picard iteration. Using this approach, we can systematically construct explicit, implicit and semi-implicit schemes with arbitrary order of accuracy. We also show that the same schemes can be constructed by solving a series of correction equations based on the DG weak formulation. The accuracy of the schemes is proven to be $\min\{2p+1, K+1\}$ with $p$ the degree of the DG polynomial basis and $K$ the number of iterations. The stability is explored numerically; we show that the implicit schemes are $A$-stable at least for $0 \leq p \leq 9$. Furthermore, we combine the methods with a multilevel strategy to accelerate their convergence speed. The new multilevel scheme is intended to provide a flexible framework for high order space-time discretizations and to be coupled with space-time multigrid techniques for solving partial differential equations (PDEs). We present numerical examples for ODEs and PDEs to analyze the performance of the new methods. Moreover, the newly proposed class of methods, due to its structure, is also a competitive and promising candidate for parallel in time algorithms such as Parareal, PFASST, multigrid in time, etc., Comment: Submitted to Math. of Comp

Multiscale models allow for the treatment of complex phenomena involving different scales, such as remodeling and growth of tissues, muscular activation, and cardiac electrophysiology. Numerous numerical approaches have been developed to simulate multiscale problems. However, compared to the well-established methods for classical problems, many questions have yet to be answered. Here, we give an overview of existing models and methods, with particular emphasis on mechanical and bio-mechanical applications. Moreover, we discuss state-of-the-art techniques for multilevel and multifidelity uncertainty quantification. In particular, we focus on the similarities that can be found across multiscale models, discretizations, solvers, and statistical methods for uncertainty quantification. Similarly to the current trend of removing the segregation between discretizations and solution methods in scientific computing, we anticipate that the future of multiscale simulation will provide a closer interaction with also the models and the statistical methods. This will yield better strategies for transferring the information across different scales and for a more seamless transition in selecting and adapting the level of details in the models. Finally, we note that machine learning and Bayesian techniques have shown a promising capability to capture complex model dependencies and enrich the results with statistical information; therefore, they can complement traditional physics-based and numerical analysis approaches.

The applicability of the Parareal parallel-in-time integration scheme for the solution of a linear, two-dimensional hyperbolic acoustic-advection system, which is often used as a test case for integration schemes for numerical weather prediction (NWP), is addressed. Parallel-in-time schemes are a possible way to increase, on the algorithmic level, the amount of parallelism, a requirement arising from the rapidly growing number of CPUs in high performance computer systems. A recently introduced modification of the "parallel implicit time-integration algorithm" could successfully solve hyperbolic problems arising in structural dynamics. It has later been cast into the framework of Parareal. The present paper adapts this modified Parareal and employs it for the solution of a hyperbolic flow problem, where the initial value problem solved in parallel arises from the spatial discretization of a partial differential equation by a finite difference method. It is demonstrated that the modified Parareal is stable and can produce reasonably accurate solutions while allowing for a noticeable reduction of the time-to-solution. The implementation relies on integration schemes already widely used in NWP (RK-3, partially split forward Euler, forward-backward). It is demonstrated that using an explicit partially split scheme for the coarse integrator allows to avoid the use of an implicit scheme while still achieving speedup.

In this technical report we study the convergence of Parareal for 2D incompressible flow around a cylinder for different viscosities. Two methods are used as fine integrator: backward Euler and a fractional step method. It is found that Parareal converges better for the implicit Euler, likely because it under-resolves the fine-scale dynamics as a result of numerical diffusion., Comment: 16 pages, 7 figures

In-silico investigation of skin permeation is an important but also computationally demanding problem. To resolve all scales involved in full detail will not only require exascale computing capacities but also suitable parallel algorithms. This article investigates the applicability of the time-parallel Parareal algorithm to a brick and mortar setup, a precursory problem to skin permeation. The C++ library Lib4PrM implementing Parareal is combined with the UG4 simulation framework, which provides the spatial discretization and parallelization. The combination's performance is studied with respect to convergence and speedup. It is confirmed that anisotropies in the domain and jumps in diffusion coefficients only have a minor impact on Parareal's convergence. The influence of load imbalances in time due to differences in number of iterations required by the spatial solver as well as spatio-temporal weak scaling is discussed., Comment: 11 pages, 8 figures

In view of the rapid rise of the number of cores in modern supercomputers, time-parallel methods that introduce concurrency along the temporal axis are becoming increasingly popular. For the solution of time-dependent partial differential equations, these methods can add another direction for concurrency on top of spatial parallelization. The paper presents an implementation of the time-parallel Parareal method in a C++ domain specific language for stencil computations (STELLA). STELLA provides both an OpenMP and a CUDA backend for a shared memory parallelization, using the CPU or GPU inside a node for the spatial stencils. Here, we intertwine this node-wise spatial parallelism with the time-parallel Parareal. This is done by adding an MPI-based implementation of Parareal, which allows us to parallelize in time across nodes. The performance of Parareal with both backends is analyzed in terms of speedup, parallel efficiency and energy-to-solution for an advection-diffusion problem with a time-dependent diffusion coefficient.

Spectral deferred correction (SDC) methods are an attractive approach to iteratively computing collocation solutions to an ODE by performing so-called sweeps with a low-order time stepping method. SDC allows to easily construct high order split methods where e.g. stiff terms of the ODE are treated implicitly. This requires the solution to full accuracy of multiple linear systems of equations during each sweep, e.g. with a multigrid method. In this paper, we present an inexact variant of SDC, where each solve of a linear system is replaced by a single multigrid V-cycle and thus significantly reduces the cost for each sweep. For the investigated examples, this strategy results only in a small increase of the number of required sweeps and we demonstrate that "inexact spectral deferred corrections" can provide a dramatic reduction of the overall number of multigrid V-cycles required to complete an SDC time step.

For the time-parallel Parareal method, there exists both numerical and analytical proof that it converges very well for diffusive problems like the heat equation. Many applications, however, do not lead to simple homogeneous diffusive scenarios but feature strongly inhomogeneous and possibly anisotropic coefficients. The paper presents results from a numerical study of how space- and time-dependent coefficients in a diffusion setup affect Parareal's convergence behaviour. It is shown that, for the presented examples, non-constant diffusion coefficients have only marginal influence on how fast Parareal converges. Furthermore, an example is shown that illustrates how for linear problems the maximum singular value of the Parareal iteration matrix can be used to estimate convergence rates.

The paper presents first a linear stability analysis for the time-parallel Parareal method, using an IMEX Euler as coarse and a Runge-Kutta-3 method as fine propagator, confirming that dominant imaginary eigenvalues negatively affect Parareal's convergence. This suggests that when Parareal is applied to the nonlinear Navier-Stokes equations, problems for small viscosities could arise. Numerical results for a driven cavity benchmark are presented, confirming that Parareal's convergence can indeed deteriorate as viscosity decreases and the flow becomes increasingly dominated by convection. The effect is found to strongly depend on the spatial resolution.

An accurate assessment for fatigue damage as a function of activation and deactivation cycles is vital for the design of many engineering parts. In this paper we extend the probabilistic and local approach to this problem proposed in [1,2] and [3] to the case of non-constant temperature fields and thermomechanical loading. The method has been implemented as a finite element postprocessor and applied to an example case of a gas-turbine blade which is made of a conventionally cast nickel base superalloy., Comment: 8 pages, 3 figures

The paper presents a combination of the time-parallel "parallel full approximation scheme in space and time" (PFASST) with a parallel multigrid method (PMG) in space, resulting in a mesh-based solver for the three-dimensional heat equation with a uniquely high degree of efficient concurrency. Parallel scaling tests are reported on the Cray XE6 machine "Monte Rosa" on up to 16,384 cores and on the IBM Blue Gene/Q system "JUQUEEN" on up to 65,536 cores. The efficacy of the combined spatial- and temporal parallelization is shown by demonstrating that using PFASST in addition to PMG significantly extends the strong-scaling limit. Implications of using spatial coarsening strategies in PFASST's multi-level hierarchy in large-scale parallel simulations are discussed., Comment: 10 pages

The spectral deferred correction (SDC) method is an iterative scheme for computing a higher-order collocation solution to an ODE by performing a series of correction sweeps using a low-order timestepping method. This paper examines a variation of SDC for the temporal integration of PDEs called multi-level spectral deferred corrections (MLSDC), where sweeps are performed on a hierarchy of levels and an FAS correction term, as in nonlinear multigrid methods, couples solutions on different levels. Three different strategies to reduce the computational cost of correction sweeps on the coarser levels are examined: reducing the degrees of freedom, reducing the order of the spatial discretization, and reducing the accuracy when solving linear systems arising in implicit temporal integration. Several numerical examples demonstrate the effect of multi-level coarsening on the convergence and cost of SDC integration. In particular, MLSDC can provide significant savings in compute time compared to SDC for a three-dimensional problem., Comment: Appears in BIT Numerical Mathematics, 2014

An accurate risk assessment for fatigue damage is of vital importance for the design and service of today's turbomachinery components. We present an approach for quantifying the probability of crack initiation due to surface driven low-cycle fatigue (LCF). This approach is based on the theory of failure-time processes and takes inhomogeneous stress fields and size effects into account. The method has been implemented as a finite-element postprocessor which uses quadrature formulae of higher order. Results of applying this new approach to an example case of a gas-turbine compressor disk are discussed., Comment: 12 pages, 6 figures