SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1879-A1902, June 2024. Abstract. While extracting information from data with machine learning plays an increasingly important role, physical laws and other first principles continue to provide critical insights about systems and processes of interest in science and engineering. This work introduces a method that infers models from data with physical insights encoded in the form of structure and that minimizes the model order so that the training data are fitted well while redundant degrees of freedom without conditions and sufficient data to fix them are automatically eliminated. The models are formulated via solution matrices of specific instances of generalized Sylvester equations that enforce interpolation of the training data and relate the model order to the rank of the solution matrices. The proposed method numerically solves the Sylvester equations for minimal-rank solutions and so obtains models of low order. Numerical experiments demonstrate that the combination of structure preservation and rank minimization leads to accurate models with orders of magnitude fewer degrees of freedom than models of comparable prediction quality that are learned with structure preservation alone.
{"title":"Rank-Minimizing and Structured Model Inference","authors":"Pawan Goyal, Benjamin Peherstorfer, Peter Benner","doi":"10.1137/23m1554308","DOIUrl":"https://doi.org/10.1137/23m1554308","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1879-A1902, June 2024. <br/> Abstract. While extracting information from data with machine learning plays an increasingly important role, physical laws and other first principles continue to provide critical insights about systems and processes of interest in science and engineering. This work introduces a method that infers models from data with physical insights encoded in the form of structure and that minimizes the model order so that the training data are fitted well while redundant degrees of freedom without conditions and sufficient data to fix them are automatically eliminated. The models are formulated via solution matrices of specific instances of generalized Sylvester equations that enforce interpolation of the training data and relate the model order to the rank of the solution matrices. The proposed method numerically solves the Sylvester equations for minimal-rank solutions and so obtains models of low order. Numerical experiments demonstrate that the combination of structure preservation and rank minimization leads to accurate models with orders of magnitude fewer degrees of freedom than models of comparable prediction quality that are learned with structure preservation alone.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"18 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-05-29","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"141168810","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1850-A1878, June 2024. Abstract. For a nonlinear dynamical system that depends on parameters, this paper introduces a novel tensorial reduced-order model (TROM). The reduced model is projection-based, and for systems with no parameters involved, it resembles proper orthogonal decomposition (POD) combined with the discrete empirical interpolation method (DEIM). For parametric systems, TROM employs low-rank tensor approximations in place of truncated SVD, a key dimension-reduction technique in POD with DEIM. Three popular low-rank tensor compression formats are considered for this purpose: canonical polyadic, Tucker, and tensor train. The use of multilinear algebra tools allows the incorporation of information about the parameter dependence of the system into the reduced model and leads to a POD-DEIM type ROM that (i) is parameter-specific (localized) and predicts the system dynamics for out-of-training set (unseen) parameter values, (ii) mitigates the adverse effects of high parameter space dimension, (iii) has online computational costs that depend only on tensor compression ranks but not on the full-order model size, and (iv) achieves lower reduced space dimensions compared to the conventional POD-DEIM ROM. This paper explains the method, analyzes its prediction power, and assesses its performance for two specific parameter-dependent nonlinear dynamical systems.
{"title":"Tensorial Parametric Model Order Reduction of Nonlinear Dynamical Systems","authors":"Alexander V. Mamonov, Maxim A. Olshanskii","doi":"10.1137/23m1553789","DOIUrl":"https://doi.org/10.1137/23m1553789","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1850-A1878, June 2024. <br/> Abstract. For a nonlinear dynamical system that depends on parameters, this paper introduces a novel tensorial reduced-order model (TROM). The reduced model is projection-based, and for systems with no parameters involved, it resembles proper orthogonal decomposition (POD) combined with the discrete empirical interpolation method (DEIM). For parametric systems, TROM employs low-rank tensor approximations in place of truncated SVD, a key dimension-reduction technique in POD with DEIM. Three popular low-rank tensor compression formats are considered for this purpose: canonical polyadic, Tucker, and tensor train. The use of multilinear algebra tools allows the incorporation of information about the parameter dependence of the system into the reduced model and leads to a POD-DEIM type ROM that (i) is parameter-specific (localized) and predicts the system dynamics for out-of-training set (unseen) parameter values, (ii) mitigates the adverse effects of high parameter space dimension, (iii) has online computational costs that depend only on tensor compression ranks but not on the full-order model size, and (iv) achieves lower reduced space dimensions compared to the conventional POD-DEIM ROM. This paper explains the method, analyzes its prediction power, and assesses its performance for two specific parameter-dependent nonlinear dynamical systems.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"42 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-05-29","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"141187801","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Thomas Frachon, Peter Hansbo, Erik Nilsson, Sara Zahedi
SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1793-A1820, June 2024. Abstract. We study cut finite element discretizations of a Darcy interface problem based on the mixed finite element pairs [math], [math]. Here [math] is the space of discontinuous polynomial functions of degree less than or equal to [math] and [math] is the Raviart–Thomas space. We show that the standard ghost penalty stabilization, often added in the weak forms of cut finite element methods for stability and control of the condition number of the resulting linear system matrix, destroys the divergence-free property of the considered element pairs. Therefore, we propose new stabilization terms for the pressure and show that we recover the optimal approximation of the divergence without losing control of the condition number of the linear system matrix. We prove that with the new stabilization term the proposed cut finite element discretization results in pointwise divergence-free approximations of solenoidal velocity fields. We derive a priori error estimates for the proposed unfitted finite element discretization based on [math], [math]. In addition, by decomposing the computational mesh into macroelements and applying ghost penalty terms only on interior edges of macroelements, stabilization is applied very restrictively and active only where needed. Numerical experiments with element pairs [math], [math], and [math] (where [math] is the Brezzi–Douglas–Marini space) indicate that with the new method we have (1) optimal rates of convergence of the approximate velocity and pressure; (2) well-posed linear systems where the condition number of the system matrix scales as it does for fitted finite element discretizations; (3) optimal rates of convergence of the approximate divergence with pointwise divergence-free approximations of solenoidal velocity fields. All three properties hold independently of how the interface is positioned relative to the computational mesh. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/CutFEM/CutFEM-Library and in the supplementary materials (CutFEM-Library-master.zip [30.5MB]).
{"title":"A Divergence Preserving Cut Finite Element Method for Darcy Flow","authors":"Thomas Frachon, Peter Hansbo, Erik Nilsson, Sara Zahedi","doi":"10.1137/22m149702x","DOIUrl":"https://doi.org/10.1137/22m149702x","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1793-A1820, June 2024. <br/>Abstract. We study cut finite element discretizations of a Darcy interface problem based on the mixed finite element pairs [math], [math]. Here [math] is the space of discontinuous polynomial functions of degree less than or equal to [math] and [math] is the Raviart–Thomas space. We show that the standard ghost penalty stabilization, often added in the weak forms of cut finite element methods for stability and control of the condition number of the resulting linear system matrix, destroys the divergence-free property of the considered element pairs. Therefore, we propose new stabilization terms for the pressure and show that we recover the optimal approximation of the divergence without losing control of the condition number of the linear system matrix. We prove that with the new stabilization term the proposed cut finite element discretization results in pointwise divergence-free approximations of solenoidal velocity fields. We derive a priori error estimates for the proposed unfitted finite element discretization based on [math], [math]. In addition, by decomposing the computational mesh into macroelements and applying ghost penalty terms only on interior edges of macroelements, stabilization is applied very restrictively and active only where needed. Numerical experiments with element pairs [math], [math], and [math] (where [math] is the Brezzi–Douglas–Marini space) indicate that with the new method we have (1) optimal rates of convergence of the approximate velocity and pressure; (2) well-posed linear systems where the condition number of the system matrix scales as it does for fitted finite element discretizations; (3) optimal rates of convergence of the approximate divergence with pointwise divergence-free approximations of solenoidal velocity fields. All three properties hold independently of how the interface is positioned relative to the computational mesh. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/CutFEM/CutFEM-Library and in the supplementary materials (CutFEM-Library-master.zip [30.5MB]).","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"49 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-05-24","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"141150370","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Golo A. Wimmer, Ben S. Southworth, Thomas J. Gregory, Xian-Zhu Tang
SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1821-A1849, June 2024. Abstract. We present a novel solver technique for the anisotropic heat flux equation, aimed at the high level of anisotropy seen in magnetic confinement fusion plasmas. Such problems pose two major challenges: (i) discretization accuracy and (ii) efficient implicit linear solvers. We simultaneously address each of these challenges by constructing a new finite element discretization with excellent accuracy properties, tailored to a novel solver approach based on algebraic multigrid (AMG) methods designed for advective operators. We pose the problem in a mixed formulation, introducing the directional temperature gradient as an auxiliary variable. The temperature and auxiliary fields are discretized in a scalar discontinuous Galerkin space with upwinding principles used for discretizations of advection. We demonstrate the proposed discretization’s superior accuracy over other discretizations of anisotropic heat flux, achieving error [math] smaller for anisotropy ratio of [math], for closed field lines. The block matrix system is reordered and solved in an approach where the two advection operators are inverted using AMG solvers based on approximate ideal restriction, which is particularly efficient for upwind discontinuous Galerkin discretizations of advection. To ensure that the advection operators are nonsingular, in this paper we restrict ourselves to considering open (acyclic) magnetic field lines for the linear solvers. We demonstrate fast convergence of the proposed iterative solver in highly anisotropic regimes where other diffusion-based AMG methods fail.
{"title":"A Fast Algebraic Multigrid Solver and Accurate Discretization for Highly Anisotropic Heat Flux I: Open Field Lines","authors":"Golo A. Wimmer, Ben S. Southworth, Thomas J. Gregory, Xian-Zhu Tang","doi":"10.1137/23m155918x","DOIUrl":"https://doi.org/10.1137/23m155918x","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1821-A1849, June 2024. <br/> Abstract. We present a novel solver technique for the anisotropic heat flux equation, aimed at the high level of anisotropy seen in magnetic confinement fusion plasmas. Such problems pose two major challenges: (i) discretization accuracy and (ii) efficient implicit linear solvers. We simultaneously address each of these challenges by constructing a new finite element discretization with excellent accuracy properties, tailored to a novel solver approach based on algebraic multigrid (AMG) methods designed for advective operators. We pose the problem in a mixed formulation, introducing the directional temperature gradient as an auxiliary variable. The temperature and auxiliary fields are discretized in a scalar discontinuous Galerkin space with upwinding principles used for discretizations of advection. We demonstrate the proposed discretization’s superior accuracy over other discretizations of anisotropic heat flux, achieving error [math] smaller for anisotropy ratio of [math], for closed field lines. The block matrix system is reordered and solved in an approach where the two advection operators are inverted using AMG solvers based on approximate ideal restriction, which is particularly efficient for upwind discontinuous Galerkin discretizations of advection. To ensure that the advection operators are nonsingular, in this paper we restrict ourselves to considering open (acyclic) magnetic field lines for the linear solvers. We demonstrate fast convergence of the proposed iterative solver in highly anisotropic regimes where other diffusion-based AMG methods fail.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"38 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-05-24","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"141150418","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1747-A1763, June 2024. Abstract. A randomized algorithm for computing a compressed representation of a given rank-structured matrix [math] is presented. The algorithm interacts with [math] only through its action on vectors. Specifically, it draws two tall thin matrices [math] from a suitable distribution, and then reconstructs [math] from the information contained in the set [math]. For the specific case of a “Hierarchically Block Separable (HBS)” matrix (a.k.a. Hierarchically Semi-Separable matrix) of block rank [math], the number of samples [math] required satisfies [math], with [math] being representative. While a number of randomized algorithms for compressing rank-structured matrices have previously been published, the current algorithm appears to be the first that is both of truly linear complexity (no [math] factors in the complexity bound) and fully “black box” in the sense that no matrix entry evaluation is required. Further, all samples can be extracted in parallel, enabling the algorithm to work in a “streaming” or “single view” mode.
{"title":"Linear-Complexity Black-Box Randomized Compression of Rank-Structured Matrices","authors":"James Levitt, Per-Gunnar Martinsson","doi":"10.1137/22m1528574","DOIUrl":"https://doi.org/10.1137/22m1528574","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1747-A1763, June 2024. <br/> Abstract. A randomized algorithm for computing a compressed representation of a given rank-structured matrix [math] is presented. The algorithm interacts with [math] only through its action on vectors. Specifically, it draws two tall thin matrices [math] from a suitable distribution, and then reconstructs [math] from the information contained in the set [math]. For the specific case of a “Hierarchically Block Separable (HBS)” matrix (a.k.a. Hierarchically Semi-Separable matrix) of block rank [math], the number of samples [math] required satisfies [math], with [math] being representative. While a number of randomized algorithms for compressing rank-structured matrices have previously been published, the current algorithm appears to be the first that is both of truly linear complexity (no [math] factors in the complexity bound) and fully “black box” in the sense that no matrix entry evaluation is required. Further, all samples can be extracted in parallel, enabling the algorithm to work in a “streaming” or “single view” mode.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"26 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-05-17","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"141063181","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1714-A1746, June 2024. Abstract. Tensor wheel (TW) decomposition is an elegant compromise of the popular tensor ring decomposition and fully connected tensor network decomposition, and it has many applications. In this work, we investigate the computation of this decomposition. Three randomized algorithms based on random sampling or random projection are proposed. Specifically, by defining a new tensor product called the subwheel product, the structures of the coefficient matrices of the alternating least squares subproblems from the minimization problem of TW decomposition are first figured out. Then, using the structures and the properties of the subwheel product, a random sampling algorithm based on leverage sampling and two random projection algorithms respectively based on Kronecker subsampled randomized Fourier transform and TensorSketch are derived. These algorithms can implement the sampling and projection on TW factors and hence can avoid forming the full coefficient matrices of subproblems. We present the complexity analysis and numerical performance on synthetic data, real data, and image reconstruction for our algorithms. Experimental results show that, compared with the deterministic algorithm in the literature, they need much less computing time while achieving similar accuracy and reconstruction effect. We also apply the proposed algorithms to tensor completion and find that the sampling-based algorithm always has excellent performance and the projection-based algorithms behave well when the sampling rate is higher than 50%.
{"title":"Randomized Tensor Wheel Decomposition","authors":"Mengyu Wang, Yajie Yu, Hanyu Li","doi":"10.1137/23m1583934","DOIUrl":"https://doi.org/10.1137/23m1583934","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1714-A1746, June 2024. <br/> Abstract. Tensor wheel (TW) decomposition is an elegant compromise of the popular tensor ring decomposition and fully connected tensor network decomposition, and it has many applications. In this work, we investigate the computation of this decomposition. Three randomized algorithms based on random sampling or random projection are proposed. Specifically, by defining a new tensor product called the subwheel product, the structures of the coefficient matrices of the alternating least squares subproblems from the minimization problem of TW decomposition are first figured out. Then, using the structures and the properties of the subwheel product, a random sampling algorithm based on leverage sampling and two random projection algorithms respectively based on Kronecker subsampled randomized Fourier transform and TensorSketch are derived. These algorithms can implement the sampling and projection on TW factors and hence can avoid forming the full coefficient matrices of subproblems. We present the complexity analysis and numerical performance on synthetic data, real data, and image reconstruction for our algorithms. Experimental results show that, compared with the deterministic algorithm in the literature, they need much less computing time while achieving similar accuracy and reconstruction effect. We also apply the proposed algorithms to tensor completion and find that the sampling-based algorithm always has excellent performance and the projection-based algorithms behave well when the sampling rate is higher than 50%.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"25 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-05-15","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"141257897","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Ikrom Akramov, Sebastian Götschel, Michael Minion, Daniel Ruprecht, Robert Speck
SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1690-A1713, June 2024. Abstract. Spectral deferred corrections (SDC) are a class of iterative methods for the numerical solution of ordinary differential equations. SDC can be interpreted as a Picard iteration to solve a fully implicit collocation problem, preconditioned with a low-order method. It has been widely studied for first-order problems, using explicit, implicit, or implicit-explicit Euler and other low-order methods as preconditioner. For first-order problems, SDC achieves arbitrary order of accuracy and possesses good stability properties. While numerical results for SDC applied to the second-order Lorentz equations exist, no theoretical results are available for SDC applied to second-order problems. We present an analysis of the convergence and stability properties of SDC using velocity-Verlet as the base method for general second-order initial value problems. Our analysis proves that the order of convergence depends on whether the force in the system depends on the velocity. We also demonstrate that the SDC iteration is stable under certain conditions. Finally, we show that SDC can be computationally more efficient than a simple Picard iteration or a fourth-order Runge–Kutta–Nyström method. Reproducibility of computational results.This paper has been awarded the “SIAM Reproducibility Badge: code and data available,” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/Parallel-in-Time/pySDC/tree/master/pySDC/projects/Second_orderSDC.
{"title":"Spectral Deferred Correction Methods for Second-Order Problems","authors":"Ikrom Akramov, Sebastian Götschel, Michael Minion, Daniel Ruprecht, Robert Speck","doi":"10.1137/23m1592596","DOIUrl":"https://doi.org/10.1137/23m1592596","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1690-A1713, June 2024. <br/> Abstract. Spectral deferred corrections (SDC) are a class of iterative methods for the numerical solution of ordinary differential equations. SDC can be interpreted as a Picard iteration to solve a fully implicit collocation problem, preconditioned with a low-order method. It has been widely studied for first-order problems, using explicit, implicit, or implicit-explicit Euler and other low-order methods as preconditioner. For first-order problems, SDC achieves arbitrary order of accuracy and possesses good stability properties. While numerical results for SDC applied to the second-order Lorentz equations exist, no theoretical results are available for SDC applied to second-order problems. We present an analysis of the convergence and stability properties of SDC using velocity-Verlet as the base method for general second-order initial value problems. Our analysis proves that the order of convergence depends on whether the force in the system depends on the velocity. We also demonstrate that the SDC iteration is stable under certain conditions. Finally, we show that SDC can be computationally more efficient than a simple Picard iteration or a fourth-order Runge–Kutta–Nyström method. Reproducibility of computational results.This paper has been awarded the “SIAM Reproducibility Badge: code and data available,” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/Parallel-in-Time/pySDC/tree/master/pySDC/projects/Second_orderSDC.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"11 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-05-14","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"141063176","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Gabriel R. Barrenechea, Antonio Tadeu A. Gomes, Diego Paredes
SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1628-A1657, June 2024. Abstract. In this work we propose, analyze, and test a new multiscale finite element method called Multiscale Hybrid (MH) method. The method is built as a close relative to the Multiscale Hybrid Mixed (MHM) method, but with the fundamental difference that a novel definition of the Lagrange multiplier is introduced. The practical implication of this is that both the local problems to compute the basis functions, as well as the global problem, are elliptic, as opposed to the MHM method (and also other previous methods) where a mixed global problem is solved and constrained local problems are solved to compute the local basis functions. The error analysis of the method is based on a hybrid formulation, and a static condensation process is done at the discrete level, so the final global system only involves the Lagrange multipliers. We tested the performance of the method by means of numerical experiments for problems with multiscale coefficients, and we carried out comparisons with the MHM method in terms of performance, accuracy, and memory requirements.
{"title":"A Multiscale Hybrid Method","authors":"Gabriel R. Barrenechea, Antonio Tadeu A. Gomes, Diego Paredes","doi":"10.1137/22m1542556","DOIUrl":"https://doi.org/10.1137/22m1542556","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1628-A1657, June 2024. <br/> Abstract. In this work we propose, analyze, and test a new multiscale finite element method called Multiscale Hybrid (MH) method. The method is built as a close relative to the Multiscale Hybrid Mixed (MHM) method, but with the fundamental difference that a novel definition of the Lagrange multiplier is introduced. The practical implication of this is that both the local problems to compute the basis functions, as well as the global problem, are elliptic, as opposed to the MHM method (and also other previous methods) where a mixed global problem is solved and constrained local problems are solved to compute the local basis functions. The error analysis of the method is based on a hybrid formulation, and a static condensation process is done at the discrete level, so the final global system only involves the Lagrange multipliers. We tested the performance of the method by means of numerical experiments for problems with multiscale coefficients, and we carried out comparisons with the MHM method in terms of performance, accuracy, and memory requirements.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"153 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-05-13","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"140938967","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Nick Alger, Tucker Hartland, Noemi Petra, Omar Ghattas
SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1658-A1689, June 2024. Abstract. We present an efficient matrix-free point spread function (PSF) method for approximating operators that have locally supported nonnegative integral kernels. The PSF-based method computes impulse responses of the operator at scattered points and interpolates these impulse responses to approximate entries of the integral kernel. To compute impulse responses efficiently, we apply the operator to Dirac combs associated with batches of point sources, which are chosen by solving an ellipsoid packing problem. The ability to rapidly evaluate kernel entries allows us to construct a hierarchical matrix (H-matrix) approximation of the operator. Further matrix computations are then performed with fast H-matrix methods. This end-to-end procedure is illustrated on a blur problem. We demonstrate the PSF-based method’s effectiveness by using it to build preconditioners for the Hessian operator arising in two inverse problems governed by PDEs: inversion for the basal friction coefficient in an ice sheet flow problem and for the initial condition in an advective-diffusive transport problem. While for many ill-posed inverse problems the Hessian of the data misfit term exhibits a low-rank structure, and hence a low-rank approximation is suitable, for many problems of practical interest, the numerical rank of the Hessian is still large. The Hessian impulse responses, on the other hand, typically become more local as the numerical rank increases, which benefits the PSF-based method. Numerical results reveal that the preconditioner clusters the spectrum of the preconditioned Hessian near one, yielding roughly [math]–[math] reductions in the required number of PDE solves, as compared to classical regularization-based preconditioning and no preconditioning. We also present a comprehensive numerical study for the influence of various parameters (that control the shape of the impulse responses and the rank of the Hessian) on the effectiveness of the advection-diffusion Hessian approximation. The results show that the PSF-based method is able to form good approximations of high-rank Hessians using only a small number of operator applications.
{"title":"Point Spread Function Approximation of High-Rank Hessians with Locally Supported Nonnegative Integral Kernels","authors":"Nick Alger, Tucker Hartland, Noemi Petra, Omar Ghattas","doi":"10.1137/23m1584745","DOIUrl":"https://doi.org/10.1137/23m1584745","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1658-A1689, June 2024. <br/> Abstract. We present an efficient matrix-free point spread function (PSF) method for approximating operators that have locally supported nonnegative integral kernels. The PSF-based method computes impulse responses of the operator at scattered points and interpolates these impulse responses to approximate entries of the integral kernel. To compute impulse responses efficiently, we apply the operator to Dirac combs associated with batches of point sources, which are chosen by solving an ellipsoid packing problem. The ability to rapidly evaluate kernel entries allows us to construct a hierarchical matrix (H-matrix) approximation of the operator. Further matrix computations are then performed with fast H-matrix methods. This end-to-end procedure is illustrated on a blur problem. We demonstrate the PSF-based method’s effectiveness by using it to build preconditioners for the Hessian operator arising in two inverse problems governed by PDEs: inversion for the basal friction coefficient in an ice sheet flow problem and for the initial condition in an advective-diffusive transport problem. While for many ill-posed inverse problems the Hessian of the data misfit term exhibits a low-rank structure, and hence a low-rank approximation is suitable, for many problems of practical interest, the numerical rank of the Hessian is still large. The Hessian impulse responses, on the other hand, typically become more local as the numerical rank increases, which benefits the PSF-based method. Numerical results reveal that the preconditioner clusters the spectrum of the preconditioned Hessian near one, yielding roughly [math]–[math] reductions in the required number of PDE solves, as compared to classical regularization-based preconditioning and no preconditioning. We also present a comprehensive numerical study for the influence of various parameters (that control the shape of the impulse responses and the rank of the Hessian) on the effectiveness of the advection-diffusion Hessian approximation. The results show that the PSF-based method is able to form good approximations of high-rank Hessians using only a small number of operator applications.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"32 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-05-13","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"140939197","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Daniel Jodlbauer, Ulrich Langer, Thomas Wick, Walter Zulehner
SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1599-A1627, June 2024. Abstract. We consider the widely used continuous [math]-[math] quadrilateral or hexahedral Taylor–Hood elements for the finite element discretization of the Stokes and generalized Stokes systems in two and three spatial dimensions. For the fast solution of the corresponding symmetric, but indefinite system of finite element equations, we propose and analyze matrix-free monolithic geometric multigrid solvers that are based on appropriately scaled Chebyshev–Jacobi smoothers. The analysis is based on results by Schöberl and Zulehner (2003). We present and discuss several numerical results for typical benchmark problems.
{"title":"Matrix-Free Monolithic Multigrid Methods for Stokes and Generalized Stokes Problems","authors":"Daniel Jodlbauer, Ulrich Langer, Thomas Wick, Walter Zulehner","doi":"10.1137/22m1504184","DOIUrl":"https://doi.org/10.1137/22m1504184","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1599-A1627, June 2024. <br/> Abstract. We consider the widely used continuous [math]-[math] quadrilateral or hexahedral Taylor–Hood elements for the finite element discretization of the Stokes and generalized Stokes systems in two and three spatial dimensions. For the fast solution of the corresponding symmetric, but indefinite system of finite element equations, we propose and analyze matrix-free monolithic geometric multigrid solvers that are based on appropriately scaled Chebyshev–Jacobi smoothers. The analysis is based on results by Schöberl and Zulehner (2003). We present and discuss several numerical results for typical benchmark problems.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"21 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-05-09","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"140938884","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
扫码关注我们
求助内容：
应助结果提醒方式：
京公网安备 11010802042870号