Abstract: We present two new Lagrangian methods for hydrodynamics, in a systematic comparison with moving-mesh, smoothed particle hydrodynamics (SPH), and stationary (non-moving) grid methods. The new methods are designed to simultaneously capture advantages of both SPH and grid-based/adaptive mesh refinement (AMR) schemes. They are based on a kernel discretization of the volume coupled to a high-order matrix gradient estimator and a Riemann solver acting over the volume ‘overlap’. We implement and test a parallel, second-order version of the method with self-gravity and cosmological integration, in the code gizmo:1 this maintains exact mass, energy and momentum conservation; exhibits superior angular momentum conservation compared to all other methods we study; does not require ‘artificial diffusion’ terms; and allows the fluid elements to move with the flow, so resolution is automatically adaptive. We consider a large suite of test problems, and find that on all problems the new methods appear competitive with moving-mesh schemes, with some advantages (particularly in angular momentum conservation), at the cost of enhanced noise. The new methods have many advantages versus SPH: proper convergence, good capturing of fluid-mixing instabilities, dramatically reduced ‘particle noise’ and numerical viscosity, more accurate sub-sonic flow evolution, and sharp shock-capturing. Advantages versus non-moving meshes include: automatic adaptivity, dramatically reduced advection errors and numerical overmixing, velocity-independent errors, accurate coupling to gravity, good angular momentum conservation and elimination of ‘grid alignment’ effects. We can, for example, follow hundreds of orbits of gaseous discs, while AMR and SPH methods break down in a few orbits. However, fixed meshes minimize ‘grid noise’. These differences are important for a range of astrophysical problems.

Abstract: The aim of this paper is to provide the key theoretical and numerical aspects of the NEMOH software. NEMOH is a numerical solver for computations of first order hydrodynamic coefficients in frequency domain. It is based on linear free surface potential flow theory. It has been developed in Ecole Centrale de Nantes for 30 years. It has been released in open source in January 2014. This paper starts with recalling the assumptions and the resulting boundary value problem for the free surface flow around a body with an arbitrary condition for the normal velocity on the body surface. Then, it is recalled how the 3D flow problem can be transformed into the 2D problem of a source distribution on the body surface using Green's second identity and the appropriate Green function. Discretization of the mathematical problem using the constant panel method leads to a linear matrix problem whose coefficients are the influence coefficients. Depending on the distance between the source and field point, the influence coefficients are calculated analytically or approximately. Irregular frequencies are removed by adding constraints on the interior flow problem. Hydrodynamic coefficients for usual diffraction and radiation problems can be obtained by specifying to the solver the appropriate body conditions. The flexibility in specifying the body conditions allows dealing with complex structures, e.g deformable wave energy converters. Eventually, it is explained how far field coefficients and free surface elevation can be obtained from the solver outputs.

Abstract: This work develops new minimization and saddle point principles for the coupled problem of Darcy–Biot-type fluid transport in porous media at fracture. It shows that the quasi-static problem of elastically deforming, fluid-saturated porous media is related to a minimization principle for the evolution problem. This two-field principle determines the rate of deformation and the fluid mass flux vector. It provides a canonically compact model structure, where the stress equilibrium and the inverse Darcy's law appear as the Euler equations of a variational statement. A Legendre transformation of the dissipation potential relates the minimization principle to a characteristic three field saddle point principle, whose Euler equations determine the evolutions of deformation and fluid content as well as Darcy's law. A further geometric assumption results in modified variational principles for a simplified theory, where the fluid content is linked to the volumetric deformation. The existence of these variational principles underlines inherent symmetries of Darcy–Biot theories of porous media. This can be exploited in the numerical implementation by the construction of time- and space-discrete variational principles, which fully determine the update problems of typical time stepping schemes. Here, the proposed minimization principle for the coupled problem is advantageous with regard to a new unconstrained stable finite element design, while space discretizations of the saddle point principles are constrained by the LBB condition. The variational principles developed provide the most fundamental approach to the discretization of nonlinear fluid–structure interactions, showing symmetric systems in algebraic update procedures. They also provide an excellent starting point for extensions towards more complex problems. This is demonstrated by developing a minimization principle for a phase field description of fracture in fluid-saturated porous media. It is designed for an incorporation of alternative crack driving forces, such as a convenient criterion in terms of the effective stress. The proposed setting provides a modeling framework for the analysis of complex problems such as hydraulic fracture. This is demonstrated by a spectrum of model simulations.

Abstract: We study accelerated mirror descent dynamics in continuous and discrete time. Combining the original continuous-time motivation of mirror descent with a recent ODE interpretation of Nesterov's accelerated method, we propose a family of continuous-time descent dynamics for convex functions with Lipschitz gradients, such that the solution trajectories converge to the optimum at a O(1/t2) rate. We then show that a large family of first-order accelerated methods can be obtained as a discretization of the ODE, and these methods converge at a O(1/k2) rate. This connection between accelerated mirror descent and the ODE provides an intuitive approach to the design and analysis of accelerated first-order algorithms.

Abstract: Although the simulation of quantum chemistry is one of the most anticipated applications of quantum computing, the scaling of known upper bounds on the complexity of these algorithms is daunting. Prior work has bounded errors due to discretization of the time evolution (known as ``Trotterization'') in terms of the norm of the error operator and analyzed scaling with respect to the number of spin orbitals. However, we find that these error bounds can be loose by up to 16 orders of magnitude for some molecules. Furthermore, numerical results for small systems fail to reveal any clear correlation between ground-state error and number of spin orbitals. We instead argue that chemical properties, such as the maximum nuclear charge in a molecule and the filling fraction of orbitals, can be decisive for determining the cost of a quantum simulation. Our analysis motivates several strategies to use classical processing to further reduce the required Trotter step size and estimate the necessary number of steps, without requiring additional quantum resources. Finally, we demonstrate improved methods for state preparation techniques which are asymptotically superior to proposals in the simulation literature.

Abstract: The aim of this paper is to develop a virtual element method for the two-dimensional Steklov eigenvalue problem. We propose a discretization by means of the virtual elements presented in [L. Beirao da Veiga et al., Basic principles of virtual element methods, Math. Models Methods Appl. Sci.23 (2013) 199–214]. Under standard assumptions on the computational domain, we establish that the resulting scheme provides a correct approximation of the spectrum and prove optimal-order error estimates for the eigenfunctions and a double order for the eigenvalues. We also prove higher-order error estimates for the computation of the eigensolutions on the boundary, which in some Steklov problems (computing sloshing modes, for instance) provides the quantity of main interest (the free surface of the liquid). Finally, we report some numerical tests supporting the theoretical results.

Abstract: This work concerns linearization methods for efficiently solving the Richards` equation,a degenerate elliptic-parabolic equation which models flow in saturated/unsaturated porous media.The discretization of Richards` equation is based on backward Euler in time and Galerkin finite el-ements in space. The most valuable linearization schemes for Richards` equation, i.e. the Newtonmethod, the Picard method, the Picard/Newton method and theLscheme are presented and theirperformance is comparatively studied. The convergence, the computational time and the conditionnumbers for the underlying linear systems are recorded. The convergence of theLscheme is theo-retically proved and the convergence of the other methods is discussed. A new scheme is proposed,theLscheme/Newton method which is more robust and quadratically convergent. The linearizationmethods are tested on illustrative numerical examples.

Abstract: 1. Introduction.- Part I: Analytical and Numerical Foundations.- 2. Analytical Background.- 3. FEM for Linear Problems.- 4. Concepts for Discretized Problems.- Part II: Approximation of Classical Formulations.- 5. The Obstacle Problem.- 6. The Allen-Cahn Equation.- 7. Harmonic Maps.- 8. Bending Problems.- Part III: Methods for Extended Formulations.- 9. Nonconvexity and Microstructure.- 10. Free Discontinuities.- 11. Elastoplasticity.- Auxiliary Routines.- Frequently Used Notation.- Index.

Abstract: The fractional Fokker--Planck equation is an important physical model for simulating anomalous diffusions with external forces. Because of the nonlocal property of the fractional derivative an interesting problem is to explore high accuracy numerical methods for fractional differential equations. In this paper, a space-time spectral method is presented for the numerical solution of the time fractional Fokker--Planck initial-boundary value problem. The proposed method employs the Jacobi polynomials for the temporal discretization and Fourier-like basis functions for the spatial discretization. Due to the diagonalizable trait of the Fourier-like basis functions, this leads to a reduced representation of the inner product in the Galerkin analysis. We prove that the time fractional Fokker--Planck equation attains the same approximation order as the time fractional diffusion equation developed in [X. Li and C. Xu, SIAM J. Numer. Anal., 47 (2009), pp. 2108--2131] by using the present method. That indicates an e...

Abstract: Introduction.- Basic Theory of Partial Differential Equations and Their Discretization.- Theory of PDE-constrained Optimization.- Numerical Optimization Methods.- Box-constrained Problems.- Nonsmooth PDE-constrained Optimization.

Abstract: Many problems from mass transport can be reformulated as variational problems under a prescribed divergence constraint (static problems) or subject to a time-dependent continuity equation, which again can be formulated as a divergence constraint but in time and space. The variational class of mean field games, introduced by Lasry and Lions, may also be interpreted as a generalization of the time-dependent optimal transport problem. Following Benamou and Brenier, we show that augmented Lagrangian methods are well suited to treat such convex but non-smooth problems. They include in particular Monge historic optimal transport problem. A finite-element discretization and implementation of the method are used to provide numerical simulations and a convergence study.

Abstract: . We present a new method and implementation (Instaseis) to store global Green's functions in a database which allows for near-instantaneous (on the order of milliseconds) extraction of arbitrary seismograms. Using the axisymmetric spectral element method (AxiSEM), the generation of these databases, based on reciprocity of the Green's functions, is very efficient and is approximately half as expensive as a single AxiSEM forward run. Thus, this enables the computation of full databases at half the cost of the computation of seismograms for a single source in the previous scheme and allows to compute databases at the highest frequencies globally observed. By storing the basis coefficients of the numerical scheme (Lagrange polynomials), the Green's functions are 4th order accurate in space and the spatial discretization respects discontinuities in the velocity model exactly. High-order temporal interpolation using Lanczos resampling allows to retrieve seismograms at any sampling rate. AxiSEM is easily adaptable to arbitrary spherically symmetric models of Earth as well as other planets. In this paper, we present the basic rationale and details of the method as well as benchmarks and illustrate a variety of applications. The code is open source and available with extensive documentation at www.instaseis.net .

Abstract: This article details a novel numerical scheme to approximate gradient flows for optimal transport (i.e. Wasserstein) metrics. These flows have proved useful to tackle theoretically and numerically non-linear diffusion equations that model for instance porous media or crowd evolutions. These gradient flows define a suitable notion of weak solutions for these evolutions and they can be approximated in a stable way using discrete flows. These discrete flows are implicit Euler time stepping according to the Wasserstein metric. A bottleneck of these approaches is the high computational load induced by the resolution of each step. Indeed, this corresponds to the resolution of a convex optimization problem involving a Wasserstein distance to the previous iterate. Following several recent works on the approximation of Wasserstein distances, we consider a discrete flow induced by an entropic regularization of the transportation coupling. This entropic regularization allows one to trade the initial Wasserstein fidelity term for a Kulback-Leibler divergence, which is easier to deal with numerically. We show how KL proximal schemes, and in particular Dykstra's algorithm, can be used to compute each step of the regularized flow. The resulting algorithm is both fast, parallelizable and versatile, because it only requires multiplications by a Gibbs kernel. On Euclidean domains discretized on an uniform grid, this corresponds to a linear filtering (for instance a Gaussian filtering when $c$ is the squared Euclidean distance) which can be computed in nearly linear time. On more general domains, such as (possibly non-convex) shapes or on manifolds discretized by a triangular mesh, following a recently proposed numerical scheme for optimal transport, this Gibbs kernel multiplication is approximated by a short-time heat diffusion.

Abstract: We have investigated the use of the adaptive high-order finite-element method (FEM) for geoelectromagnetic modeling. Because high-order FEM is challenging from the numerical and computational points of view, most published finite-element studies in geoelectromagnetics use the lowest order formulation. Solution of the resulting large system of linear equations poses the main practical challenge. We have developed a fully parallel and distributed robust and scalable linear solver based on the optimal block-diagonal and auxiliary space preconditioners. The solver was found to be efficient for high finite element orders, unstructured and nonconforming locally refined meshes, a wide range of frequencies, large conductivity contrasts, and number of degrees of freedom (DoFs). Furthermore, the presented linear solver is in essence algebraic; i.e., it acts on the matrix-vector level and thus requires no information about the discretization, boundary conditions, or physical source used, making it readily ef...

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.

Abstract: Hydrogels are capable of coupled mass transport and large deformation in response to external stimuli. In this paper, a nonlinear, transient finite element formulation is presented for initial boundary value problems associated with swelling and deformation of hydrogels, based on a nonlinear continuum theory that is consistent with classical theory of linear poroelasticity. A mixed finite element method is implemented with implicit time integration. The incompressible or nearly incompressible behavior at the initial stage imposes a constraint to the finite element discretization in order to satisfy the Ladyzhenskaya–Babuska–Brezzi (LBB) condition for stability of the mixed method, similar to linear poroelasticity as well as incompressible elasticity and Stokes flow; failure to choose an appropriate discretization would result in locking and numerical oscillations in transient analysis. To demonstrate the numerical method, two problems of practical interests are considered: constrained swelling and flat-punch indentation of hydrogel layers. Constrained swelling may lead to instantaneous surface instability for a soft hydrogel in a good solvent, which can be regulated by assuming a stiff surface layer. Indentation relaxation of hydrogels is simulated beyond the linear regime under plane strain conditions, in comparison with two elastic limits for the instantaneous and equilibrium states. The effects of Poisson’s ratio and loading rate are discussed. It is concluded that the present finite element method is robust and can be extended to study other transient phenomena in hydrogels.

Abstract: In this paper we present a quasi-static formulation of a phase-field model for a pressurized crack in a poroelastic medium. The mathematical model represents a linear elasticity system with a fading Gassman tensor as the crack grows, that is coupled with a variational inequality for the phase-field variable containing an entropy inequality. We introduce a novel incremental approximation that decouples displacement and phase-field problems. We establish convergence to a solution of the quasi-static problem, including Rice's condition, when the time discretization step goes to zero. Numerical experiments confirm the robustness and efficiency of this approach for multidimensional test cases.