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.

Approximated forms of the RII and RIII redistribution matrices are frequently applied to simplify the numerical solution of the radiative transfer problem for polarized radiation, taking partial frequency redistribution (PRD) effects into account. A widely used approximation for RIII is to consider its expression under the assumption of complete frequency redistribution (CRD) in the observer frame (RIII CRD). The adequacy of this approximation for modeling the intensity profiles has been firmly established. By contrast, its suitability for modeling scattering polarization signals has only been analyzed in a few studies, considering simplified settings. In this work, we aim at quantitatively assessing the impact and the range of validity of the RIII CRD approximation in the modeling of scattering polarization. Methods. We first present an analytic comparison between RIII and RIII CRD. We then compare the results of radiative transfer calculations, out of local thermodynamic equilibrium, performed with RIII and RIII CRD in realistic 1D atmospheric models. We focus on the chromospheric Ca i line at 4227 A and on the photospheric Sr i line at 4607 A.

The efficient construction of anatomical models is one of the major challenges of patient-specific in-silico models of the human heart. Current methods frequently rely on linear statistical models, allowing no advanced topological changes, or requiring medical image segmentation followed by a meshing pipeline, which strongly depends on image resolution, quality, and modality. These approaches are therefore limited in their transferability to other imaging domains. In this work, the cardiac shape is reconstructed by means of three-dimensional deep signed distance functions with Lipschitz regularity. For this purpose, the shapes of cardiac MRI reconstructions are learned to model the spatial relation of multiple chambers. We demonstrate that this approach is also capable of reconstructing anatomical models from partial data, such as point clouds from a single ventricle, or modalities different from the trained MRI, such as the electroanatomical mapping (EAM).

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.

Images conveniently capture the result of physical processes, representing rich source of information for data driven medicine, engineering, and science. The modeling of an image as a graph allows the application of graph-based algorithms for content analysis. Amongst these, one of the most used is the Dijkstra Single Source Shortest Path algorithm (DSSSP), which computes the path with minimal cost from one starting node to all the other nodes of the graph. However, the results of DSSSP remains unknown for nodes until they are explored. Moreover, DSSSP execution is associated to frequent jumps between distant locations in the graph, which results in non-optimal memory access, reduced parallelization, and finally increased execution time. Therefore, we propose AnyDijkstra, an iterative implementation of the Dijkstra SSSP algorithm optimized for images, that retains anytime properties while accessing memory following a cache-friendly scheme and maximizing parallelization., Comment: 7 pages, 4 figures

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.

Computational models of atrial fibrillation have successfully been used to predict optimal ablation sites. A critical step to assess the effect of an ablation pattern is to pace the model from different, potentially random, locations to determine whether arrhythmias can be induced in the atria. In this work, we propose to use multi-fidelity Gaussian process classification on Riemannian manifolds to efficiently determine the regions in the atria where arrhythmias are inducible. We build a probabilistic classifier that operates directly on the atrial surface. We take advantage of lower resolution models to explore the atrial surface and combine seamlessly with high-resolution models to identify regions of inducibility. When trained with 40 samples, our multi-fidelity classifier shows a balanced accuracy that is 10% higher than a nearest neighbor classifier used as a baseline atrial fibrillation model, and 9% higher in presence of atrial fibrillation with ablations. We hope that this new technique will allow faster and more precise clinical applications of computational models for atrial fibrillation., Comment: 22 pages, 7 figures

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.

Context. The numerical modeling of the generation and transfer of polarized radiation is a key task in solar and stellar physics research and has led to a relevant class of discrete problems that can be reframed as linear systems. In order to solve such problems, it is common to rely on efficient stationary iterative methods. However, the convergence properties of these methods are problem-dependent, and a rigorous investigation of their convergence conditions, when applied to transfer problems of polarized radiation, is still lacking. Aims. After summarizing the most widely employed iterative methods used in the numerical transfer of polarized radiation, this article aims to clarify how the convergence of these methods depends on different design elements, such as the choice of the formal solver, the discretization of the problem, or the use of damping factors. The main goal is to highlight advantages and disadvantages of the different iterative methods in terms of stability and rate of convergence. Methods. We first introduce an algebraic formulation of the radiative transfer problem. This formulation allows us to explicitly assemble the iteration matrices arising from different stationary iterative methods, compute their spectral radii and derive their convergence rates, and test the impact of different discretization settings, problem parameters, and damping factors. Conclusions. The general methodology used in this article, based on a fully algebraic formulation of linear transfer problems of polarized radiation, provides useful estimates of the convergence rates of various iterative schemes. Additionally, it can lead to novel solution approaches as well as analyses for a wider range of settings, including the unpolarized case.

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 train deep residual networks with a stochastic variant of the nonlinear multigrid method MG/OPT. To build the multilevel hierarchy, we use the dynamical systems viewpoint specific to residual networks. We report significant speed-ups and additional robustness for training MNIST on deep residual networks. Our numerical experiments also indicate that multilevel training can be used as a pruning technique, as many of the auxiliary networks have accuracies comparable to the original network.

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.

Electroanatomical maps are a key tool in the diagnosis and treatment of atrial fibrillation. Current approaches focus on the activation times recorded. However, more information can be extracted from the available data. The fibers in cardiac tissue conduct the electrical wave faster, and their direction could be inferred from activation times. In this work, we employ a recently developed approach, called physics informed neural networks, to learn the fiber orientations from electroanatomical maps, taking into account the physics of the electrical wave propagation. In particular, we train the neural network to weakly satisfy the anisotropic eikonal equation and to predict the measured activation times. We use a local basis for the anisotropic conductivity tensor, which encodes the fiber orientation. The methodology is tested both in a synthetic example and for patient data. Our approach shows good agreement in both cases, with an RMSE of 2.2ms on the in-silico data and outperforming a state of the art method on the patient data. The results show a first step towards learning the fiber orientations from electroanatomical maps with physics-informed neural networks., Comment: 10 pages, 3 figures

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

Electroanatomical mapping, a keystone diagnostic tool in cardiac electrophysiology studies, can provide high-density maps of the local electric properties of the tissue. It is therefore tempting to use such data to better individualize current patient-specific models of the heart through a data assimilation procedure and to extract potentially insightful information such as conduction properties. Parameter identification for state-of-the-art cardiac models is however a challenging task. In this work, we introduce a novel inverse problem for inferring the anisotropic structure of the conductivity tensor, that is fiber orientation and conduction velocity along and across fibers, of an eikonal model for cardiac activation. The proposed method, named PIEMAP, performed robustly with synthetic data and showed promising results with clinical data. These results suggest that PIEMAP could be a useful supplement in future clinical workflows of personalized therapies., Comment: 12 pages, 4 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 propose a novel training method based on nonlinear multilevel minimization techniques, commonly used for solving discretized large scale partial differential equations. Our multilevel training method constructs a multilevel hierarchy by reducing the number of samples. The training of the original model is then enhanced by internally training surrogate models constructed with fewer samples. We construct the surrogate models using first-order consistency approach. This gives rise to surrogate models, whose gradients are stochastic estimators of the full gradient, but with reduced variance compared to standard stochastic gradient estimators. We illustrate the convergence behavior of the proposed multilevel method to machine learning applications based on logistic regression. A comparison with subsampled Newton's and variance reduction methods demonstrate the efficiency of our multilevel method.

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.

Complex geometries as common in industrial applications consist of multiple patches, if spline based parametrizations are used. The requirements for the generation of analysis-suitable models are increasing dramatically since isogeometric analysis is directly based on the spline parametrization and nowadays used for the calculation of higher-order partial differential equations. The computational, or more general, the engineering analysis necessitates suitable coupling techniques between the different patches. Mortar methods have been successfully applied for coupling of patches and for contact mechanics in recent years to resolve the arising issues within the interface. We present here current achievements in the design of mortar technologies in isogeometric analysis within the Priority Program SPP 1748, Reliable Simulation Techniques in Solid Mechanics. Development of Non-standard Discretisation Methods, Mechanical and Mathematical Analysis.