Showing papers in "Mathematical Modelling and Numerical Analysis in 2015"

A Numerical Algorithm for L-2 Optimal Transport in 3D

Abstract: This paper introduces a numerical algorithm to compute the L2 optimal transport map between two measures µ and ν, where µ derives from a density ρ defined as a piecewise linear function (supported by a tetrahedral mesh), and where ν is a sum of Dirac masses. I first give an elementary presentation of some known results on optimal transport and then observe a relation with another problem (optimal sampling). This relation gives simple arguments to study the objective functions that characterize both problems. I then propose a practical algorithm to compute the optimal transport map between a piecewise linear density and a sum of Dirac masses in 3D. In this semi-discrete setting, Aurenhammer et.al [8th Symposium on Computational Geometry conf. proc., ACM (1992)] showed that the optimal transport map is determined by the weights of a power diagram. The optimal weights are computed by minimizing a convex objective function with a quasi-Newton method. To evaluate the value and gradient of this objective function, I propose an efficient and robust algorithm, that computes at each iteration the intersection between a power diagram and the tetrahedral mesh that defines the measure µ. The numerical algorithm is experimented and evaluated on several datasets, with up to hundred thousands tetrahedra and one million Dirac masses. Resume. Cet article decrit un algorithme numerique pour calculer l'application de transport optimal L2 entre deux mesures µ et ν, o` u µ derive d'une densite ρ lineaire par morceaux (supportee par un maillage tetraedrique), et o` u ν est une somme de masses de Dirac. Je donne tout d'abord une presentation elementaire de quelques resultats connus sur le transport optimal , et observe ensuite une relation avec un autreprobi eme (l'´ echantillonage optimal). Cette relation fournit des arguments simples pour etudier les fonctions objectifs caracterisant les deuxprobi emes. Je propose ensuite un algorithme pratique pour calculer le transport optimal entre une densite lineaire par morceaux et une somme de masses de Dirac en 3D. Dans ce cas semi-discret, Auren-hammer et.al [8th Symposium on Computational Geometry conf. proc., ACM (1992)] ont montre que l'application de transport optimal est determinee par les poids d'un diagramme de puissance. Les poids optimaux sont calcules en minimisant une fonction objectif convexe a l'aide d'une methode quasi-Newton. Pour evaluer cette fonction objectif et son gradient, je propose un algorithme efficace et robuste, qui calcule a chaque iteration l'intersection entre un diagramme de puissance et le maillage tetrahedrique qui definit la mesure µ. L'algorithme numerique est experimente et evalue sur plusieurs jeux de donnees , comportant jusqu'` a plusieurs centaines de milliers detetra edres et un million de masses de Dirac.

Multi-marginal optimal transport: theory and applications ∗

Abstract: Over the past five years, multi-marginal optimal transport, a generalization of the well known optimal transport problem of Monge and Kantorovich, has begun to attract considerable atten- tion, due in part to a wide variety of emerging applications. Here, we survey this problem, addressing fundamental theoretical questions including the uniqueness and structure of solutions. The answers to these questions uncover a surprising divergence from the classical two marginal setting, and reflect a delicate dependence on the cost function, which we then illustrate with a series of examples. We go on to describe some applications of the multi-marginal optimal transport problem, focusing primarily on matching in economics and density functional theory in physics.

Discrete least squares polynomial approximation with random evaluations − application to parametric and stochastic elliptic PDEs

Abstract: Motivated by the numerical treatment of parametric and stochastic PDEs, we analyze the least-squares method for polynomial approximation of multivariate functions based on random sampling according to a given probability measure. Recent work has shown that in the univariate case, the least-squares method is quasi-optimal in expectation in [A. Cohen, M A. Davenport and D. Leviatan. Found. Comput. Math. 13 (2013) 819–834] and in probability in [G. Migliorati, F. Nobile, E. von Schwerin, R. Tempone, Found. Comput. Math. 14 (2014) 419–456], under suitable conditions that relate the number of samples with respect to the dimension of the polynomial space. Here “quasi-optimal” means that the accuracy of the least-squares approximation is comparable with that of the best approximation in the given polynomial space. In this paper, we discuss the quasi-optimality of the polynomial least-squares method in arbitrary dimension. Our analysis applies to any arbitrary multivariate polynomial space (including tensor product, total degree or hyperbolic crosses), under the minimal requirement that its associated index set is downward closed. The optimality criterion only involves the relation between the number of samples and the dimension of the polynomial space, independently of the anisotropic shape and of the number of variables. We extend our results to the approximation of Hilbert space-valued functions in order to apply them to the approximation of parametric and stochastic elliptic PDEs. As a particular case, we discuss “inclusion type” elliptic PDE models, and derive an exponential convergence estimate for the least-squares method. Numerical results confirm our estimate, yet pointing out a gap between the condition necessary to achieve optimality in the theory, and the condition that in practice yields the optimal convergence rate.

Numerical methods for matching for teams and wasserstein barycenters

Abstract: Equilibrium multi-population matching (matching for teams) is a problem from mathematical economics which is related to multi-marginal optimal transport. A special but important case is the Wasserstein barycenter problem, which has applications in image processing and statistics. Two algorithms are presented: a linear programming algorithm and an efficient nonsmooth optimization algorithm, which applies in the case of the Wasserstein barycenters. The measures are approximated by discrete measures: convergence of the approximation is proved. Numerical results are presented which illustrate the efficiency of the algorithms.

Residual a posteriori error estimation for the Virtual Element Method for elliptic problems

Abstract: A posteriori error estimation and adaptivity are very useful in the context of the virtual element and mimetic discretization methods due to the flexibility of the meshes to which these numerical schemes can be applied. Nevertheless, developing error estimators for virtual and mimetic methods is not a straightforward task due to the lack of knowledge of the basis functions. In the new virtual element setting, we develop a residual based a posteriori error estimator for the Poisson problem with (piecewise) constant coefficients, that is proven to be reliable and efficient. We moreover show the numerical performance of the proposed estimator when it is combined with an adaptive strategy for the mesh refinement.

Reduced Basis Techniques For Nonlinear Conservation Laws

Abstract: In this paper we present a new reduced basis technique for parametrized nonlinear scalar conservation laws in presence of shocks. The essential ingredients are an efficient algorithm to approximate the shock curve, a procedure to detect the smooth components of the solution at the two sides of the shock, and a suitable interpolation strategy to reconstruct such smooth components during the online stage. The approach we propose is based on some theoretical properties of the solution to the problem. Some numerical examples prove the effectiveness of the proposed strategy.

On the Numerical Integration of Scalar Nonlocal Conservation Laws

Abstract: We study a rather general class of 1D nonlocal conservation laws from a numerical point of view. First, following [F. Betancourt, R. Burger, K.H. Karlsen and E.M. Tory, On nonlocal conservation laws modelling sedimentation. Nonlinearity 24 (2011) 855–885], we define an algorithm to numerically integrate them and prove its convergence. Then, we use this algorithm to investigate various analytical properties, obtaining evidence that usual properties of standard conservation laws fail in the nonlocal setting. Moreover, on the basis of our numerical integrations, we are led to conjecture the convergence of the nonlocal equation to the local ones, although no analytical results are, to our knowledge, available in this context.

A trace finite element method for a class of coupled bulk-interface transport problems ∗

Abstract: In this paper we study a system of advection-diffusion equations in a bulk domain coupled to an advection-diffusion equation on an embedded surface. Such systems of coupled partial differential equations arise in, for example, the modeling of transport and diffusion of surfactants in two-phase flows. The model considered here accounts for adsorption-desorption of the surfactants at a sharp interface between two fluids and their transport and diffusion in both fluid phases and along the interface. The paper gives a well-posedness analysis for the system of bulk-surface equations and introduces a finite element method for its numerical solution. The finite element method is unfitted, i.e. , the mesh is not aligned to the interface. The method is based on taking traces of a standard finite element space both on the bulk domains and the embedded surface. The numerical approach allows an implicit definition of the surface as the zero level of a level-set function. Optimal order error estimates are proved for the finite element method both in the bulk-surface energy norm and the L 2 -norm. The analysis is not restricted to linear finite elements and a piecewise planar reconstruction of the surface, but also covers the discretization with higher order elements and a higher order surface reconstruction.

A hybrid variational principle for the keller-segel system in r 2

Abstract: We construct weak global in time solutions to the classical Keller-Segel system cell movement by chemotaxis in two dimensions when the total mass is below the well-known critical value. Our construction takes advantage of the fact that the Keller-Segel system can be realized as a gradient flow in a suitable functional product space. This allows us to employ a hybrid variational principle which is a generalisation of the minimising implicit scheme for Wasserstein distances introduced by Jordan, Kinderlehrer and Otto (1998).

A convergent adaptive stochastic Galerkin finite element method with quasi-optimal spatial meshes

Abstract: We analyze a posteriori error estimation and adaptive refinement algorithms for stochastic Galerkin Finite Element methods for countably-parametric, elliptic boundary value problems. A residual error estimator which separates the effects of gpc-Galerkin discretization in parameter space and of the Finite Element discretization in physical space in energy norm is established. It is proved that the adaptive algorithm converges. To this end, a contraction property of its iterates is proved. It is shown that the sequences of triangulations which are produced by the algorithm in the FE discretization of the active gpc coefficients are asymptotically optimal. Numerical experiments illustrate the theoretical results.

Numerical analysis of augmented plane wave methods for full-potential electronic structure calculations

Abstract: This paper investigates the augmented plane wave methods which are widely used in full-potential electronic structure calculations. These methods introduce basis functions that describe different regions using different discretization schemes. We construct a nonconforming method based on this idea and present an a priori error analysis for both linear Schrodinger type equations and nonlinear Kohn−Sham equations. Some numerical experiments are presented to support our theory.

An analysis of hdg methods for convection-dominated diffusion problems

Abstract: We present the first a priori error analysis of the h -version of the hybridizable discontinuous Galkerin (HDG) methods applied to convection-dominated diffusion problems. We show that, when using polynomials of degree no greater than k , the L 2 -error of the scalar variable converges with order k + 1 / 2 on general conforming quasi-uniform simplicial meshes, just as for conventional DG methods. We also show that the method achieves the optimal L 2 -convergence order of k + 1 on special meshes. Moreover, we discuss a new way of implementing the HDG methods for which the spectral condition number of the global matrix is independent of the diffusion coefficient. Numerical experiments are presented which verify our theoretical results.

A two-phase shallow debris flow model with energy balance

Abstract: This paper proposes a thin layer depth-averaged two-phase model provided by a dissipative energy balance to describe avalanches of solid-fluid mixtures. This model is derived from a 3D two-phase model based on the equations proposed by Jackson [R. Jackson, The Dynamics of Fluidized Particles, 2000] which takes into account the force of buoyancy and the forces of interaction between the solid and fluid phases. Jackson's model is based on mass and momentum conservation within the two phases, i.e. two vector and two scalar equations. This system has five unknowns: the solid volume fraction, the solid and fluid pressures and the solid and fluid velocities, i.e. three scalars and two vectors. As a result, an additional equation is necessary to close the system. Surprisingly, this issue is inadequately accounted for in the models that have been developed on the basis of Jackson's work. In particular, Pitman and Le [E.B. Pitman, L. Le, Phil. Trans. R. Soc. A, 2005] replaced this closure simply by imposing an extra boundary condition. If the pressure is assumed to be hydrostatic, this condition can be considered as a closure condition. However, the corresponding model cannot account for a dissipative energy balance. We propose here a closure equation to complete Jackson's model, imposing incompressibility of the solid phase. We prove that the resulting whole 3D model is compatible with a dissipative energy balance. From this model, we deduce a 2D depth-averaged model and we also prove that the energy balance associated with this model is dissipative. Finally, we propose a numerical scheme to approximate the depth-averaged model. We present several numerical tests for the 1D case that are compared to the results of the model proposed by Pitman and Le.

Discretization of the 3D Monge-Ampere operator, between Wide Stencils and Power Diagrams

Abstract: We introduce a monotone (degenerate elliptic) discretization of the Monge-Ampere operator, on domains discretized on cartesian grids. The scheme is consistent provided the solution hessian condition number is uniformly bounded. Our approach enjoys the simplicity of the Wide Stencil method, but significantly improves its accuracy using ideas from discretizations of optimal transport based on power diagrams. We establish the global convergence of a damped Newton solver for the discrete system of equations. Numerical experiments, in three dimensions, illustrate the scheme efficiency.

A generalized model for optimal transport of images including dissipation and density modulation

Abstract: In this paper the optimal transport and the metamorphosis perspectives are combined. For a pair of given input images geodesic paths in the space of images are defined as minimizers of a resulting path energy. To this end, the underlying Riemannian metric measures the rate of transport cost and the rate of viscous dissipation. Furthermore, the model is capable to deal with strongly varying image contrast and explicitly allows for sources and sinks in the transport equations which are incorporated in the metric related to the metamorphosis approach by Trouve and Younes. In the non-viscous case with source term existence of geodesic paths is proven in the space of measures. The proposed model is explored on the range from merely optimal transport to strongly dissipative dynamics. For this model a robust and effective variational time discretization of geodesic paths is proposed. This requires to minimize a discrete path energy consisting of a sum of consecutive image matching functionals. These functionals are defined on corresponding pairs of intensity functions and on associated pairwise matching deformations. Existence of time discrete geodesics is demonstrated. Furthermore, a finite element implementation is proposed and applied to instructive test cases and to real images. In the non-viscous case this is compared to the algorithm proposed by Benamou and Brenier including a discretization of the source term. Finally, the model is generalized to define discrete weighted barycentres with applications to textures and objects.

Development and stability analysis of the inverse Lax−Wendroff boundary treatment for central compact schemes

Abstract: In this paper, we generalize the so-called inverse Lax−Wendroff boundary treatment [S. Tan and C.-W. Shu, J. Comput. Phys. 229 (2010) 8144–8166] for the inflow boundary of a linear hyperbolic problem discretized by the recently introduced central compact schemes [X. Liu, S. Zhang, H. Zhang and C.-W. Shu, J. Comput. Phys. 248 (2013) 235–256]. The outflow boundary is treated by the classical extrapolation and a stability analysis for the resulting scheme is provided. To ensure the stability of the considered schemes provided with the chosen boundaries, the G-K-S theory [B. Gustafsson, H.-O. Kreiss and A. Sundstrom, Math. Comput. 26 (1972) 649–686] is used, first in the semidiscrete case then in the fully discrete case with the third-order TVD Runge−Kutta time discretization. Afterwards, due to the high algebraic complexity of the G-K-S theory, the stability is analyzed by visualizing the eigenspectrum of the discretized operators. We show in this paper that the results obtained with these two different approaches are perfectly consistent. We also illustrate the high accuracy of the presented schemes on simple test problems.

An augmented mixed-primal finite element method for a coupled flow-transport problem

Abstract: In this paper we analyze the coupling of a scalar nonlinear convection-diffusion problem with the Stokes equations where the viscosity depends on the distribution of the solution to the transport problem. An augmented variational approach for the fluid flow coupled with a primal formulation for the transport model is proposed. The resulting Galerkin scheme yields an augmented mixed-primal finite element method employing Raviart−Thomas spaces of order k for the Cauchy stress, and continuous piecewise polynomials of degree ≤ k + 1 for the velocity and also for the scalar field. The classical Schauder and Brouwer fixed point theorems are utilized to establish existence of solution of the continuous and discrete formulations, respectively. In turn, suitable estimates arising from the connection between a regularity assumption and the Sobolev embedding and Rellich−Kondrachov compactness theorems, are also employed in the continuous analysis. Then, sufficiently small data allow us to prove uniqueness and to derive optimal a priori error estimates. Finally, we report a few numerical tests confirming the predicted rates of convergence, and illustrating the performance of a linearized method based on Newton−Raphson iterations; and we apply the proposed framework in the simulation of thermal convection and sedimentation-consolidation processes.

Energy Conservation and numerical stability for the reduced MHD models of the non-linear JOREK code

Abstract: In this paper we present a rigorous derivation of the reduced MHD models with and without parallel velocity that are implemented in the non-linear MHD code JOREK. The model we obtain contains some terms that have been neglected in the implementation but might be relevant in the non- linear phase. These are necessary to guarantee exact conservation with respect to the full MHD energy. For the second part of this work, we have replaced the linearized time stepping of JOREK by a non- linear solver based on the Inexact Newton method including adaptive time stepping. We demonstrate that this approach is more robust especially with respect to numerical errors in the saturation phase of an instability and allows to use larger time steps in the non-linear phase.

A Finite Element Method with Singularity Reconstruction for Fractional Boundary Value Problems

Abstract: We consider a two-point boundary value problem involving a Riemann−Liouville fractional derivative of order α ∈ (1, 2) in the leading term on the unit interval (0, 1). The standard Galerkin finite element method can only give a low-order convergence even if the source term is very smooth due to the presence of the singularity term x α−1 in the solution representation. In order to enhance the convergence, we develop a simple singularity reconstruction strategy by splitting the solution into a singular part and a regular part, where the former captures explicitly the singularity. We derive a new variational formulation for the regular part, and show that the Galerkin approximation of the regular part can achieve a better convergence order in the L 2 (0, 1), H α/2 (0, 1) and L ∞ (0, 1)-norms than the standard Galerkin approach, with a convergence rate for the recovered singularity strength identical with the L 2 (0, 1) error estimate. The reconstruction approach is very flexible in handling explicit singularity, and it is further extended to the case of a Neumann type boundary condition on the left end point, which involves a strong singularity x α−2 . Extensive numerical results confirm the theoretical study and efficiency of the proposed approach.

Motivations, ideas and applications of ramified optimal transportation

Abstract: In this survey article, the author summarizes the motivations, key ideas and main applications of ramified optimal transportation that the author has studied in recent years.

Optimal transport with Coulomb cost. Approximation and duality

Abstract: We revisit the duality theorem for multimarginal optimal transportation problems. In par- ticular, we focus on the Coulomb cost. We use a discrete approximation to prove equality of the extremal values and some careful estimates of the approximating sequence to prove existence of maximizers for the dual problem (Kantorovich's potentials). Finally we observe that the same strategy can be applied to a more general class of costs and that a classical results on the topic cannot be applied here.

Multi-physics optimal transportation and image interpolation

Abstract: Optimal transportation theory is a powerful tool to deal with image interpolation. This was first investigated by [Benamou and Brenier, Numer. Math. 84 (2000) 375–393.] where an algorithm based on the minimization of a kinetic energy under a conservation of mass constraint was devised. By structure, this algorithm does not preserve image regions along the optimal interpolation path, and it is actually not very difficult to exhibit test cases where the algorithm produces a path of images where high density regions split at the beginning before merging back at its end. However, in some applications to image interpolation this behaviour is not physically realistic. Hence, this paper aims at studying how some physics can be added to the optimal transportation theory, how to construct algorithms to compute solutions to the corresponding optimization problems and how to apply the proposed methods to image interpolation.

A Nitsche finite element method for dynamic contact : 1. Semi-discrete problem analysis and time-marching schemes

Abstract: This paper presents a new approximation of elastodynamic frictionless contact problems based both on the finite element method and on an adaptation of Nitsche's method which was initially designed for Dirichlet's condition. A main interesting characteristic is that this approximation produces well-posed space semi-discretizations contrary to standard finite element discretizations. This paper is then mainly devoted to present an analysis of the semi-discrete problem in terms of consistency, well-posedness and energy conservation, and also to study the well-posedness of some time-marching schemes (theta-scheme, Newmark and a new hybrid scheme). The stability properties of the schemes and the corresponding numerical experiments can be found in a second paper.

Energy consistent discontinuous Galerkin methods for a quasi-incompressible diffuse two phase flow model

Abstract: We design consistent discontinuous Galerkin finite element schemes for the approximation of a quasi-incompressible two phase flow model of Allen–Cahn/Cahn–Hilliard/Navier–Stokes–Korteweg type which allows for phase transitions. We show that the scheme is mass conservative and monotonically energy dissipative. In this case the dissipation is isolated to discrete equivalents of those effects already causing dissipation on the continuous level, that is, there is no artificial numerical dissipation added into the scheme. In this sense the methods are consistent with the energy dissipation of the continuous PDE system.

Eulerian models and algorithms for unbalanced optimal transport

Abstract: Benamou and Brenier formulation of Monge transportation problem (Numer. Math. 84:375-393, 2000) has proven to be of great interest in image processing to compute warpings and distances between pair of images (SIAM J. Math. Analysis, 35:61-97, 2003). One requirement for the algorithm to work is to interpolate densities of same mass. In most applications to image interpolation, this is a serious limitation. Existing approaches to overcome this caveat are reviewed, and discussed. Due to the mix between transport and $L^2$ interpolation, these models can produce instantaneous motion at finite range. In this paper we propose new methods, parameter-free, for interpolating unbalanced densities. One of our motivations is the application to interpolation of growing tumor images.

A priori error estimates to smooth solutions of the third order Runge–Kutta discontinuous Galerkin method for symmetrizable systems of conservation laws

Abstract: In this paper we present an a priori error estimate of the Runge–Kutta discontinuous Galerkin method for solving symmetrizable conservation laws, where the time is discretized with the third order explicit total variation diminishing Runge–Kutta method and the finite element space is made up of piecewise polynomials of degree k ≥ 2. Quasi-optimal error estimate is obtained by energy techniques, for the so-called generalized E-fluxes under the standard temporal-spatial CFL condition τ ≤ γh , where h is the element length and τ is time step, and γ is a positive constant independent of h and τ . Optimal estimates are also considered when the upwind numerical flux is used.

Stabilized mixed approximation of axisymmetric Brinkman flows

Abstract: This paper is devoted to the numerical analysis of an augmented finite element approximation of the axisymmetric Brinkman equations. Stabilization of the variational formulation is achieved by adding suitable Galerkin least-squares terms, allowing us to transform the original problem into a formulation better suited for performing its stability analysis. The sought quantities (here velocity, vorticity, and pressure) are approximated by Raviart−Thomas elements of arbitrary order k ≥ 0, piecewise continuous polynomials of degree k + 1, and piecewise polynomials of degree k , respectively. The well-posedness of the resulting continuous and discrete variational problems is rigorously derived by virtue of the classical Babuska–Brezzi theory. We further establish a priori error estimates in the natural norms, and we provide a few numerical tests illustrating the behavior of the proposed augmented scheme and confirming our theoretical findings regarding optimal convergence of the approximate solutions.

Analysis of a high-order space and time discontinuous Galerkin method for elastodynamic equations. Application to 3D wave propagation

Abstract: In this paper, we introduce a high-order discontinuous Galerkin method, based on centered fluxes and a family of high-order leap-frog time schemes, for the solution of the 3D elastodynamic equations written in velocity-stress formulation. We prove that this explicit scheme is stable under a CFL type condition obtained from a discrete energy which is preserved in domains with free surface or decreasing in domains with absorbing boundary conditions. Moreover, we study the convergence of the method for both the semi-discrete and the fully discrete schemes, and we illustrate the convergence results by the propagation of an eigenmode. We also propose a series of absorbing conditions which allow improving the convergence of the global scheme. Finally, several numerical applications of wave propagation, using a 3D solver, help illustrating the various properties of the method.

Finite Element Decomposition and Minimal Extension for Flow Equations

Abstract: In the simulation of flows, the correct treatment of the pressure variable is the key to stable time-integration schemes. This paper contributes a new approach based on the theory of differential- algebraic equations. Motivated by the index reduction technique of minimal extension, a remodelling of the flow equations is proposed. It is shown how this reformulation can be realized for standard finite elements via a decomposition of the discrete spaces and that it ensures stable and accurate approxi- mations. The presented decomposition preserves sparsity and does not call on variable transformations which might change the meaning of the variables. Since the method is eventually an index reduction, high index effects leading to instabilities are eliminated.

Computing quantities of interest for random domains with second order shape sensitivity analysis

Abstract: We consider random perturbations of a given domain. The characteristic amplitude of these perturbations is assumed to be small. We are interested in quantities of interest which depend on the random domain through a boundary value problem. We derive asymptotic expansions of the first moments of the distribution of this output function. A simple and efficient method is proposed to compute the coefficients of these expansions provided that the random perturbation admits a low- rank spectral representation. By numerical experiments, we compare our expansions with Monte-Carlo simulations. Mathematics Subject Classification. 60G35, 65N75, 65N99.