Victor M. Calo

Victor M. Calo, Y. Efendiev, Juan Galvis et al.

In this paper, we study the development of efficient multiscale methods for flows in heterogeneous media. Our approach uses the Generalized Multiscale Finite Element (GMsFEM) framework. The main idea of GMsFEM is to approximate the solution space locally using a few multiscale basis functions. This is typically achieved by selecting an appropriate snapshot space and a local spectral decomposition, e.g., the use of oversampled regions in order to achieve an efficient model reduction. However, the successful construction of snapshot spaces may be costly if too many local problems need to be solved in order to obtain these spaces. In this paper, we show that this efficiency can be achieved using a moderate quantity of local solutions (or snapshot vectors) with random boundary conditions on oversampled regions with zero forcing. Motivated by the randomized algorithm presented in [19], we consider a snapshot space which consists of harmonic extensions of random boundary conditions defined in a domain larger than the target region. Furthermore, we perform an eigenvalue decomposition in this small space. We study the application of randomized sampling for GMsFEM in conjunction with adaptivity, where local multiscale spaces are adaptively enriched. Convergence analysis is provided. We present representative numerical results to validate the method proposed.

Nathan Collier, Lisandro Dalcin, David Pardo et al.

In this paper we study how the use of a more continuous set of basis functions affects the cost of solving systems of linear equations resulting from a discretized Galerkin weak form. Specifically, we compare performance of linear solvers when discretizing using $C^0$ B-splines, which span traditional finite element spaces, and $C^{p-1}$ B-splines, which represent maximum continuity. We provide theoretical estimates for the increase in cost of the matrix-vector product as well as for the construction and application of black-box preconditioners. We accompany these estimates with numerical results and study their sensitivity to various grid parameters such as element size $h$ and polynomial order of approximation $p$. Finally, we present timing results for a range of preconditioning options for the Laplace problem. We conclude that the matrix-vector product operation is at most $\slfrac{33p^2}{8}$ times more expensive for the more continuous space, although for moderately low $p$, this number is significantly reduced. Moreover, if static condensation is not employed, this number further reduces to at most a value of 8, even for high $p$. Preconditioning options can be up to $p^3$ times more expensive to setup, although this difference significantly decreases for some popular preconditioners such as Incomplete LU factorization.

Victor Calo, Matteo Cicuttin, Quanling Deng et al.

We study the approximation of the spectrum of a second-order elliptic differential operator by the Hybrid High-Order (HHO) method. The HHO method is formulated using cell and face unknowns which are polynomials of some degree $k\geq0$. The key idea for the discrete eigenvalue problem is to introduce a discrete operator where the face unknowns have been eliminated. Using the abstract theory of spectral approximation of compact operators in Hilbert spaces, we prove that the eigenvalues converge as $h^{2t}$ and the eigenfunctions as $h^{t}$ in the $H^1$-seminorm, where $h$ is the mesh-size, $t\in [s,k+1]$ depends on the smoothness of the eigenfunctions, and $s>\frac12$ results from the elliptic regularity theory. The convergence rates for smooth eigenfunctions are thus $h^{2k+2}$ for the eigenvalues and $h^{k+1}$ for the eigenfunctions. Our theoretical findings, which improve recent error estimates for Hybridizable Discontinuous Galerkin (HDG) methods, are verified on various numerical examples including smooth and non-smooth eigenfunctions. Moreover, we observe numerically in one dimension for smooth eigenfunctions that the eigenvalues superconverge as $h^{2k+4}$ for a specific value of the stabilization parameter.

Michael Bartoň, Victor Manuel Calo

We introduce Gaussian quadrature rules for spline spaces that are frequently used in Galerkin discretizations to build mass and stiffness matrices. By definition, these spaces are of even degrees. The optimal quadrature rules we recently derived [5] act on spaces of the smallest odd degrees and, therefore, are still slightly sub-optimal. In this work, we derive optimal rules directly for even-degree spaces and therefore further improve our recent result. We use optimal quadrature rules for spaces over two elements as elementary building blocks and use recursively the homotopy continuation concept described in [6] to derive optimal rules for arbitrary admissible number of elements. We demonstrate the proposed methodology on relevant examples, where we derive optimal rules for various even-degree spline spaces. We also discuss convergence of our rules to their asymptotic counterparts, these are the analogues of the midpoint rule of Hughes et al. [16], that are exact and optimal for infinite domains.

A. F. Sarmiento, A. M. A. Cortes, D. A. Garcia et al.

We describe the development of a high-performance solution framework for isogeometric discrete differential forms based on B-splines: PetIGA-MF. Built on top of PetIGA, PetIGA-MF is a general multi-field discretization tool. To test the capabilities of our implementation, we solve different viscous flow problems such as Darcy, Stokes, Brinkman, and Navier-Stokes equations. Several convergence benchmarks based on manufactured solutions are presented assuring optimal convergence rates of the approximations, showing the accuracy and robustness of our solver.

Victor Calo, Quanling Deng, Vladimir Puzyrev

We develop and analyze quadrature blending schemes that minimize the dispersion error of isogeometric analysis up to polynomial order seven with maximum continuity in the span ($C^{p-1}$). The schemes yield two extra orders of convergence (superconvergence) on the eigenvalue errors, while the eigenfunction errors are of optimal convergence order. Both dispersion and spectrum analysis are unified in the form of a Taylor expansion for eigenvalue errors. As a consequence, the schemes increase the accuracy and robustness of isogeometric analysis for wave propagation as well as the differential eigenvalue problems. We analyze the methods' robustness and efficacy and utilize numerical examples to verify our analysis of the performance of the proposed schemes.

Victor M. Calo, Nathaniel O. Collier, Antti H. Niemi

We analyze the discontinuous Petrov-Galerkin (DPG) method with optimal test functions when applied to solve the Reissner-Mindlin model of plate bending. We prove that the hybrid variational formulation underlying the DPG method is well-posed (stable) with a thickness-dependent constant in a norm encompassing the $L_2$-norms of the bending moment, the shear force, the transverse deflection and the rotation vector. We then construct a numerical solution scheme based on quadrilateral scalar and vector finite elements of degree $p$. We show that for affine meshes the discretization inherits the stability of the continuous formulation provided that the optimal test functions are approximated by polynomials of degree $p+3$. We prove a theoretical error estimate in terms of the mesh size $h$ and polynomial degree $p$ and demonstrate numerical convergence on affine as well as non-affine mesh sequences.

Antti H. Niemi, Nathaniel O. Collier, Victor M. Calo

We investigate the application of the discontinuous Petrov-Galerkin (DPG) finite element framework to stationary convection-diffusion problems. In particular, we demonstrate how the quasi-optimal test space norm can be utilized to improve the robustness of the DPG method with respect to vanishing diffusion. We numerically compare coarse-mesh accuracy of the approximation when using the quasi-optimal norm, the standard norm, and the weighted norm. Our results show that the quasi-optimal norm leads to more accurate results on three benchmark problems in two spatial dimensions. We address the problems associated to the resolution of the optimal test functions with respect to the quasi-optimal norm by studying their convergence numerically. In order to facilitate understanding of the method, we also include a detailed explanation of the methodology from the algorithmic point of view.

Michael Bartoň, Vladimir Puzyrev, Quanling Deng et al.

Calabro et al. (2017) changed the paradigm of the mass and stiffness computation from the traditional element-wise assembly to a row-wise concept, showing that the latter one offers integration that may be orders of magnitude faster. Considering a B-spline basis function as a non-negative measure, each mass matrix row is integrated by its own quadrature rule with respect to that measure. Each rule is easy to compute as it leads to a linear system of equations, however, the quadrature rules are of the Newton-Cotes type, that is, they require a number of quadrature points that is equal to the dimension of the spline space. In this work, we propose weighted quadrature rules of Gaussian type which require the minimum number of quadrature points while guaranteeing exactness of integration with respect to the weight function. The weighted Gaussian rules arise as solutions of non-linear systems of equations. We derive rules for the mass and stiffness matrices for uniform $C^1$ quadratic and $C^2$ cubic isogeometric discretizations. Our rules further reduce the number of quadrature points by a factor of $(\frac{p+1}{2p+1})^d$ when compared to Calabro et al. (2017), $p$ being the polynomial degree and $d$ the dimension of the problem, and consequently reduce the computational cost of the mass and stiffness matrix assembly by a similar factor.

Victor M. Calo, Nathan Collier, Matthias Gehre et al.

We study the numerical recovery of Manning's roughness coefficient for the diffusive wave approximation of the shallow water equation. We describe a conjugate gradient method for the numerical inversion. Numerical results for one-dimensional model are presented to illustrate the feasibility of the approach. Also we provide a proof of the differentiability of the weak form with respect to the coefficient as well as the continuity and boundedness of the linearized operator under reasonable assumptions using the maximal parabolic regularity theory.

Quanling Deng, Pouria Behnoudfar, Victor M. Calo

The generalized-$α$ method encompasses a wide range of time integrators. The method possesses high-frequency dissipation while minimizing unwanted low-frequency dissipation and the numerical dissipation can be controlled by the user. The method is unconditionally stable and is of second-order accuracy in time. We extend the second-order generalized-$α$ method to third-order in time while the numerical dissipation can be controlled in a similar fashion. We establish that the third-order method is unconditionally stable. We discuss a possible path to the generalization to higher order schemes. All these high-order schemes can be easily implemented into programs that already contain the second-order generalized-$α$ method.

Pouria Behnoudfar, Victor M. Calo, Quanling Deng et al.

We present a variationally separable splitting technique for the generalized-$α$ method for solving parabolic partial differential equations. We develop a technique for a tensor-product mesh which results in a solver with a linear cost with respect to the total number of degrees of freedom in the system for multi-dimensional problems. We consider finite elements and isogeometric analysis for the spatial discretization. The overall method maintains user-controlled high-frequency dissipation while minimizing unwanted low-frequency dissipation. The method has second-order accuracy in time and optimal rates ($h^{p+1}$ in $L^2$ norm and $h^p$ in $L^2$ norm of $\nabla u$) in space. We present the spectrum analysis on the amplification matrix to establish that the method is unconditionally stable. Various numerical examples illustrate the performance of the overall methodology and show the optimal approximation accuracy.

Quanling Deng, Vladimir Puzyrev, Victor Calo

We approximate the spectra of a class of $2n$-order differential operators using isogeometric analysis in mixed formulations. This class includes a wide range of differential operators such as those arising in elliptic, biharmonic, Cahn-Hilliard, Swift-Hohenberg, and phase-field crystal equations. The spectra of the differential operators are approximated by solving differential eigenvalue problems in mixed formulations, which require auxiliary parameters. The mixed isogeometric formulation when applying classical quadrature rules leads to an eigenvalue error convergence of order $2p$ where $p$ is the order of the underlying B-spline space. We improve this order to be $2p+2$ by applying optimally-blended quadrature rules developed in \cite{20,52} and this order is an optimum in the view of dispersion error. We also compare these results with the mixed finite elements and show numerically that mixed isogeometric analysis leads to significantly better spectral approximations.

Quanling Deng, Vladimir Puzyrev, Victor Calo

We study the spectral approximation of a second-order elliptic differential eigenvalue problem that arises from structural vibration problems using isogeometric analysis. In this paper, we generalize recent work in this direction. We present optimally blended quadrature rules for the isogeometric spectral approximation of a diffusion-reaction operator with both Dirichlet and Neumann boundary conditions. The blended rules improve the accuracy and the robustness of the isogeometric approximation. In particular, the optimal blending rules minimize the dispersion error and lead to two extra orders of super-convergence in the eigenvalue error. Various numerical examples (including the Schr$\ddot{\text{o}}$dinger operator for quantum mechanics) in one and three spatial dimensions demonstrate the performance of the blended rules.

Judit Muñoz-Matute, David Pardo, Victor M. Calo et al.

Variational space-time formulations for Partial Differential Equations have been of great interest in the last decades. While it is known that implicit time marching schemes have variational structure, the Galerkin formulation of explicit methods in time remains elusive. In this work, we prove that the explicit Runge-Kutta methods can be expressed as discontinuous Petrov-Galerkin methods both in space and time. We build trial and test spaces for the linear diffusion equation that lead to one, two, and general stage explicit Runge-Kutta methods. This approach enables us to design explicit time-domain (goal-oriented) adaptive algorithms

Michael Barton, Victor Calo

We introduce optimal quadrature rules for spline spaces that are frequently used in Galerkin discretizations to build mass and stiffness matrices. Using the homotopy continuation concept [6] that transforms optimal quadrature rules from source spaces to target spaces, we derive optimal rules for splines generated above finite domains. Starting with the classical Gaussian quadrature for polynomials, which is an optimal rule for a discontinuous odd-degree space, we derive rules for target spaces of higher continuity. We further show how the homotopy methodology handles cases where the source and target rules require different numbers of optimal quadrature points. We demonstrate it by deriving optimal rules for various odd-degree spline spaces, particularly with non-uniform knot sequences and non-uniform multiplicities. We also discuss convergence of our rules to their asymptotic counterparts, that is, the analogues of the midpoint rule of Hughes et al. [34], that are exact and optimal for infinite domains. For spaces of low continuities, we numerically show that the derived rules quickly converge to their asymptotic counterparts as the weights and nodes of a few boundary elements differ from the asymptotic values.

Mehdi Ghommem, Michael Presho, Victor M. Calo et al.

In this paper, we combine concepts of the generalized multiscale finite element method and mode decomposition methods to construct a robust local-global approach for model reduction of flows in high-contrast porous media. This is achieved by implementing proper orthogonal decomposition (POD) and dynamic mode decomposition (DMD) techniques on a coarse grid. The resulting reduced-order approach enables a significant reduction in the flow problem size while accurately capturing the behavior of fully resolved solutions. We consider a variety of high-contrast coefficients and present the corresponding numerical results to illustrate the effectiveness of the proposed technique. This paper is a continuation of the first part where we examine the applicability of POD and DMD to derive simplified and reliable representations of flows in high-contrast porous media. In the current paper, we discuss how these global model reduction approaches can be combined with local techniques to speed-up the simulations. The speed-up is due to inexpensive, while sufficiently accurate, computations of global snapshots.

Michael Bartoň, Rachid Ait-Haddou, Victor Manuel Calo

We provide explicit expressions for quadrature rules on the space of $C^1$ quintic splines with uniform knot sequences over finite domains. The quadrature nodes and weights are derived via an explicit recursion that avoids an intervention of any numerical solver and the rule is optimal, that is, it requires the minimal number of nodes. For each of $n$ subintervals, generically, only two nodes are required which reduces the evaluation cost by $2/3$ when compared to the classical Gaussian quadrature for polynomials. Numerical experiments show fast convergence, as $n$ grows, to the "two-third" quadrature rule of Hughes et al. for infinite domains.

Leonardo A. Poveda, Sebastian Huepo, Victor M. Calo et al.

We study linear elasticity problems with high contrast in the coefficients using asymptotic limits recently introduced. We derive an asymptotic expansion to solve heterogeneous elasticity problems in terms of the contrast in the coefficients. We study the convergence of the expansion in the $H^1$ norm.

Maciej Paszynski, Victor Calo, David Pardo

This paper describes a direct solver algorithm for a sequence of finite element meshes that are h-refined towards one or several point singularities. For such a sequence of grids, the solver delivers linear computational cost O(N) in terms of CPU time and memory with respect to the number of unknowns N. The linear computational cost is achieved by utilizing the recursive structure provided by the sequence of h-adaptive grids with a special construction of the elimination tree that allows for reutilization of previously computed partial LU factorizations over the entire unrefined part of the computational mesh. The reutilization technique reduces the computational cost of the entire sequence of h-refined grids from O(N^2) down to O(N). Theoretical estimates are illustrated with numerical results on two- and three-dimensional model problems exhibiting one or several point singularities.

Nathan Collier, David Pardo, Maciej Paszynski et al.

The multi-frontal direct solver is the state-of-the-art algorithm for the direct solution of sparse linear systems. This paper provides computational complexity and memory usage estimates for the application of the multi-frontal direct solver algorithm on linear systems resulting from B-spline-based isogeometric finite elements, where the mesh is a structured grid. Specifically we provide the estimates for systems resulting from $C^{p-1}$ polynomial B-spline spaces and compare them to those obtained using $C^0$ spaces.

Victor M. Calo, Quanling Deng, Sergio Rojas et al.

Most variational forms of isogeometric analysis use highly-continuous basis functions for both trial and test spaces. For a partial differential equation with a smooth solution, isogeometric analysis with highly-continuous basis functions for trial space results in excellent discrete approximations of the solution. However, we observe that high continuity for test spaces is not necessary. In this work, we present a framework which uses highly-continuous B-splines for the trial spaces and basis functions with minimal regularity and possibly lower order polynomials for the test spaces. To realize this goal, we adopt the residual minimization methodology. We pose the problem in a mixed formulation, which results in a system governing both the solution and a Riesz representation of the residual. We present various variational formulations which are variationally-stable and verify their equivalence numerically via numerical tests.

Quanling Deng, Victor Calo

We study the generalized finite element methods (GFEMs) for the second-order elliptic eigenvalue problem with an interface in 1D. The linear stable generalized finite element methods (SGFEM) were recently developed for the elliptic source problem with interfaces. We first generalize SGFEM to arbitrary order elements and establish the optimal error convergence of the approximate solutions for the elliptic source problem with an interface. We then apply the abstract theory of spectral approximation of compact operators to establish the error estimation for the eigenvalue problem with an interface. The error estimations on eigenpairs strongly depend on the estimation of the discrete solution operator for the source problem. We verify our theoretical findings in various numerical examples including both source and eigenvalue problems.

Leonardo A. Poveda, Sebastian Huepo, Victor M. Calo et al.

We solve elliptic systems of equations posed on highly heterogeneous materials. Examples of this class of problems are composite structures and geological processes. We focus on a model problem which is a second-order elliptic equation with discontinuous coefficients. These coefficients represent the conductivity of a composite material. We assume a background with low conductivity that contains inclusions with different thermal properties. Under this scenario we design a multiscale finite element method to efficiently approximate solutions. The method is based on an asymptotic expansions of the solution in terms of the ratio between the conductivities. The resulting method constructs (locally) finite element basis functions (one for each inclusion). These bases that generate the multiscale finite element space where the approximation of the solution is computed. Numerical experiments show the good performance of the proposed methodology.

Jesus Gonzalez-Sieiro, David Pardo, Vincenzo Nava et al.

We propose a method for reducing the spatial discretization error of coarse computational fluid dynamics (CFD) problems by enhancing the quality of low-resolution simulations using deep learning. We feed the model with fine-grid data after projecting it to the coarse-grid discretization. We substitute the default differencing scheme for the convection term by a feed-forward neural network that interpolates velocities from cell centers to face values to produce velocities that approximate the down-sampled fine-grid data well. The deep learning framework incorporates the open-source CFD code OpenFOAM, resulting in an end-to-end differentiable model. We automatically differentiate the CFD physics using a discrete adjoint code version. We present a fast communication method between TensorFlow (Python) and OpenFOAM (c++) that accelerates the training process. We applied the model to the flow past a square cylinder problem, reducing the error from 120% to 25% in the velocity for simulations inside the training distribution compared to the traditional solver using an x8 coarser mesh. For simulations outside the training distribution, the error reduction in the velocities was about 50%. The training is affordable in terms of time and data samples since the architecture exploits the local features of the physics.

M. Shahriari, D. Pardo, J. A. Rivera et al.

Deep learning (DL) is a numerical method that approximates functions. Recently, its use has become attractive for the simulation and inversion of multiple problems in computational mechanics, including the inversion of borehole logging measurements for oil and gas applications. In this context, DL methods exhibit two key attractive features: a) once trained, they enable to solve an inverse problem in a fraction of a second, which is convenient for borehole geosteering operations as well as in other real-time inversion applications. b) DL methods exhibit a superior capability for approximating highly-complex functions across different areas of knowledge. Nevertheless, as it occurs with most numerical methods, DL also relies on expert design decisions that are problem specific to achieve reliable and robust results. Herein, we investigate two key aspects of deep neural networks (DNNs) when applied to the inversion of borehole resistivity measurements: error control and adequate selection of the loss function. As we illustrate via theoretical considerations and extensive numerical experiments, these interrelated aspects are critical to recover accurate inversion results.

Victor M. Calo, Albert Romkes, Eirik Valseth

We introduce an automatic variationally stable analysis (AVS) for finite element (FE) computations of scalar-valued convection-diffusion equations with non-constant and highly oscillatory coefficients. In the spirit of least squares FE methods, the AVS-FE method recasts the governing second order partial differential equation (PDE) into a system of first-order PDEs. However, in the subsequent derivation of the equivalent weak formulation, a Petrov-Galerkin technique is applied by using different regularities for the trial and test function spaces. We use standard FE approximation spaces for the trial spaces, which are C0, and broken Hilbert spaces for the test functions. Thus, we seek to compute pointwise continuous solutions for both the primal variable and its flux (as in least squares FE methods), while the test functions are piecewise discontinuous. To ensure the numerical stability of the subsequent FE discretizations, we apply the philosophy of the discontinuous Petrov-Galerkin (DPG) method by Demkowicz and Gopalakrishnan, by invoking test functions that lead to unconditionally stable numerical systems (if the kernel of the underlying differential operator is trivial). In the AVS-FE method, the discontinuous test functions are ascertained per the DPG approach from local, decoupled, and well-posed variational problems, which lead to best approximation properties in terms of the energy norm. We present various 2D numerical verifications, including convection-diffusion problems with highly oscillatory coefficients and extremely high Peclet numbers, up to a billion. These show the unconditional stability without the need for any upwind schemes nor any other artificial numerical stabilization. The results are not highly diffused for convection-dominated problems ...

Lisandro Dalcin, Nathan Collier, Philippe Vignal et al.

We present PetIGA, a code framework to approximate the solution of partial differential equations using isogeometric analysis. PetIGA can be used to assemble matrices and vectors which come from a Galerkin weak form, discretized with Non-Uniform Rational B-spline basis functions. We base our framework on PETSc, a high-performance library for the scalable solution of partial differential equations, which simplifies the development of large-scale scientific codes, provides a rich environment for prototyping, and separates parallelism from algorithm choice. We describe the implementation of PetIGA, and exemplify its use by solving a model nonlinear problem. To illustrate the robustness and flexibility of PetIGA, we solve some challenging nonlinear partial differential equations that include problems in both solid and fluid mechanics. We show strong scaling results on up to 4096 cores, which confirm the suitability of PetIGA for large scale simulations.

Michael Bartoň, Victor Manuel Calo

We introduce a new concept for generating optimal quadrature rules for splines. Given a target spline space where we aim to generate an optimal quadrature rule, we build an associated source space with known optimal quadrature and transfer the rule from the source space to the target one, preserving the number of quadrature points and therefore optimality. The quadrature nodes and weights are, considered as a higher-dimensional point, a zero of a particular system of polynomial equations. As the space is continuously deformed by modifying the source knot vector, the quadrature rule gets updated using polynomial homotopy continuation. For example, starting with $C^1$ cubic splines with uniform knot sequences, we demonstrate the methodology by deriving the optimal rules for uniform $C^2$ cubic spline spaces where the rule was only conjectured heretofore. We validate our algorithm by showing that the resulting quadrature rule is independent of the path chosen between the target and the source knot vectors as well as the source rule chosen.

Rachid Ait-Haddou, Michael Bartoň, Victor Manuel Calo

We provide explicit expressions for quadrature rules on the space of $C^1$ cubic splines with non-uniform, symmetrically stretched knot sequences. The quadrature nodes and weights are derived via an explicit recursion that avoids an intervention of any numerical solver and the rule is optimal, that is, it requires minimal number of nodes. Numerical experiments validating the theoretical results and the error estimates of the quadrature rules are also presented.