Geometrically-exact time-integration mesh-free schemes for advection-diffusion problems derived from optimal transportation theory and their connection with particle methods

TL;DR

This paper introduces a novel mesh-free particle method for advection-diffusion problems based on optimal transportation theory, achieving geometrically exact updates and demonstrating robustness through 3D examples.

Contribution

The paper develops an innovative mesh-free scheme derived from optimal transportation principles, integrating entropy and Wasserstein distance for accurate advection-diffusion modeling.

Findings

Geometrically exact advection and volume preservation.

Robustness demonstrated through 3D numerical examples.

Effective mesh-free max-ent interpolation for density updates.

Abstract

We develop an Optimal Transportation Meshfree (OTM) particle method for advection-diffusion in which the concentration or density of the diffusive species is approximated by Dirac measures. We resort to an incremental variational principle for purposes of time discretization of the diffusive step. This principle characterizes the evolution of the density as a competition between the Wasserstein distance between two consecutive densities and entropy. Exploiting the structure of the Euler-Lagrange equations, we approximate the density as a collection of Diracs. The interpolation of the incremental transport map is effected through mesh-free max-ent interpolation. Remarkably, the resulting update is geometrically exact with respect to advection and volume. We present three-dimensional examples of application that illustrate the scope and robustness of the method.

Tables

Table 1. Table 1: Meshes to be used for the convergence analysis.

Mesh	Nodes	material-points	$h_{\min }$
Coarse	235	593	0.26
Medium	442	1,259	0.17
Fine	1,510	4,997	0.12

Equations

$$[left=\text{empheqlbrace}]$$

$$\kappa \nabla \rho \cdot n=0,$$

$$\rho (x,0)={\rho }_0(x),$$

$${\partial }_t\rho +\nabla \cdot (\rho v)=0,$$

$$\rho v=\rho u-\kappa \nabla \rho ,$$

$$\rho (\phi (x,t),t)=\frac{{\rho }_0(x)}{\det \nabla \phi (x,t)}$$

$${\partial }_t\phi (x,t)=v(\phi (x,t),t),$$

$${\int }_{\Omega }\eta (y)\rho (y,t)dy={\int }_{\Omega }\eta (\phi (x,t)){\rho }_0(x)dx,$$

$${\int }_{\Omega }\eta (y)d\mu (y,t)={\int }_{\Omega }\eta (\phi (x,t))d{\mu }_0(x),$$

$${\int }_{\Omega }{\partial }_t\phi ({\phi }^{-1}(y,t),t)\cdot \xi (y)dy={\int }_{\Omega }v(y,t)\cdot \xi (y)dy,$$

$$\xi \cdot n=0.$$

$$\int_{\Omega}\partial_{t}\varphi(\varphi^{-1}(y,t),t)\cdot\xi(y)\,dy=\int_{\Omega}\big{(}\rho(y,t)u(y,t)-\kappa\nabla\rho(y,t)\big{)}\cdot\xi(y)\,dy.$$

$$\int_{\Omega}\partial_{t}\varphi(\varphi^{-1}(y,t),t)\cdot\xi(y)\,dy=\int_{\Omega}\big{(}u(y,t)\cdot\xi(y)+\kappa\,{\rm div}\,\xi(y)\big{)}\rho(y,t)\,dy,$$

$$\int_{\Omega}\partial_{t}\varphi(\varphi^{-1}(y,t),t)\cdot\xi(y)\,dy=\int_{\Omega}\big{(}u(y,t)\cdot\xi(y)+\kappa\,{\rm div}\,\xi(y)\big{)}\,d\mu(y,t),$$

$${\partial }_t\rho +\nabla \cdot (\rho u)=0,$$

$${\partial }_t\chi (x,t)=u(\chi (x,t),t).$$

$${\rho }_{k+1}={{\chi }_{k\to k+1}}_{\#}{\rho }_k,$$

$${\chi }_{k\to k+1}(y)=\chi ({\chi }^{-1}(y,t_k),t_{k+1})$$

$${\partial }_t\rho =\kappa \Delta \rho ,$$

$$F({\phi }_{k\to k+1})=\frac{1}{2}\frac{d_W^2({\rho }_k,{\rho }_{k+1})}{t_{k+1}-t_k}+{\int }_{\Omega }\kappa {\rho }_{k+1}\log {\rho }_{k+1}dx,$$

$$\rho_{k+1}\circ\varphi_{k\to k+1}=\rho_{k}/\det\big{(}\nabla\varphi_{k\to k+1}\big{)}.$$

$$F({\phi }_{k\to k+1})\to \min !$$

$$\begin{array}{rl}⟨DF({\phi }_{k\to k+1}),{\xi }_{k+1}⟩& =\frac{1}{t_{k+1}-t_k}{\int }_{\Omega }⟨x-{\phi }_{k+1\to k}(x),{\xi }_{k+1}(x)⟩{\rho }_{k+1}(x)dx\\ & -{\int }_{\Omega }\kappa (\log {\rho }_{k+1}+1)\nabla \cdot ({\rho }_{k+1}{\xi }_{k+1})dx,\\ & =\frac{1}{t_{k+1}-t_k}{\int }_{\Omega }⟨x-{\phi }_{k+1\to k}(x),{\xi }_{k+1}(x)⟩{\rho }_{k+1}(x)dx\\ & +{\int }_{\Omega }\kappa \nabla {\rho }_{k+1}\cdot {\xi }_{k+1}dx,\end{array}$$

$${\phi }_{k+1\to k}={\phi }_{k\to k+1}^{-1},$$

$${\xi }_{n+1}\cdot n=0.$$

$${\rho }_{k+1}(x)\frac{x-{\phi }_{k+1\to k}(x)}{t_{k+1}-t_k}=-\kappa \nabla {\rho }_{k+1}(x),$$

$${\int }_{\Omega }{\rho }_{k+1}(y)\frac{y-{\phi }_{k+1\to k}(y)}{t_{k+1}-t_k}\cdot \xi (y)dy={\int }_{\Omega }\kappa {\rho }_{k+1}(y)\nabla \cdot \xi (y)dy.$$

$$\int\rho_{k}(x)\eta(x)\,dx=\int\rho_{k+1}(y)\eta\big{(}\varphi^{-1}_{k\to k+1}(y)\big{)}dy,$$

$$\rho_{h,k}(x)=\sum_{p=1}^{M}m_{p,k}\delta\big{(}x-x_{p,k}\big{)},$$

$$\sum _{p=1}^M{m}_{p,k}\eta (x_{p,k})=\sum _{p=1}^M{m}_{p,k+1}\eta (x_{p,k}),$$

$$m_{p,k}=m_{p,k+1}=m_p,$$

$$\rho_{h,k}(x)=\sum_{p=1}^{M}m_{p}\delta\big{(}x-x_{p,k}\big{)},$$

Full text

1 Introduction

Advection-diffusion problems have been traditionally regarded as defining time-evolutions in linear spaces. For instance, the classical analysis of the heat equation regards its solutions as elements of the Sobolev space $L^{2}(0,T;H^{1}_{0}(\Omega))$ with time derivatives in $L^{2}(0,T;H^{-1}(\Omega))$, [1]. The traditional formulation of numerical approximation schemes for advection-diffusion problems correspondingly works fully within this linear space structure and aims to approximate the governing PDEs by recourse to, e. g., finite-difference or finite-element discretization schemes. Such schemes are inevitably beset by a toxic combination of spurious oscillations, instabilities, poor accuracy, lack of monotonicity, numerical diffusion, leaking of mass and others (e. g., [2]). Despite massive and sustained efforts to overcome them, such deficiencies have shown remarkably stubborn staying power.

More recently, it has been recognized, starting with the seminal work of Jordan, Kinderlehrer and Otto [3, 4, 5], that the linear functional-analytical framework is not-at-all natural or convenient and that a more natural and productive framework is to regard advection-diffusion problems as transport problems defining time-evolutions in spaces of measures. From this new perspective, the essential difficulties experienced by linear-space-based approximation schemes for advection-diffusion problems maybe be regarded as the ’punishment’ for the ’crime’ of formulating the problems in an unnatural linear-space setting. Indeed, as we shall see, the overwhelming advantage of transport methods is that they are geometrically exact, i. e., they are exact with respect to advection and volume. This property effectively eliminates in one fell swoop the difficulties experienced by traditional linear-space methods in dealing with advection and volume, albeit at the expense of working in a somewhat more challenging non-linear space, or manifold, framework.

In a parallel development, particle methods have attained considerable acceptance in solid and fluid mechanics, e. g, the Smoothed Particle Hydrodynamics method (SPH) [6, 7], the Material Point Method [8], the Reproducing Kernel Particle method [9], the Corrective Smoothed Particle Method [10], the Modified Smoothed Particle Hydrodynamics method [11], the Optimal-Transportation Meshfree method (OTM) [12], and others. Particle methods have the appeal of reducing problems to the evolution of interacting discrete particles, in seeming analogy with molecular dynamics. There is also a clear connection between particle methods and transport of measures. Indeed, the particles maybe regarded as Dirac-delta approximations of otherwise continuous measures, or densities. However, continuum transport problems differ from true particle-dynamics problems in one important respect: in continuum transport problems, the velocity field that governs the instantaneous motion of the particles depends on local density gradients, not just particle positions. This differential structure introduces additional regularity requirements of the density interpolation, such as conformity, and is at the core of the difficulties that plague particle methods such as SPH. The differential structure of the mobility law of continuum transport problems in fact necessitates two types of representations: one for the measure itself, e. g., represented as a collection of Diracs; and a more regular representation for the transport velocity field, e. g., based on meshfree conforming interpolation. This double representational requirement was recognized in connection with the OTM method [12], but has remained obscure in much of the literature on particle methods.

The objective of the present work is to develop an OTM particle method for advection-diffusion in which the concentration or density of the diffusive species is approximated by Dirac measures. The new OTM formulation hybridizes elements of a Galerkin approximation with those of an updated Lagrangian approach in the context of optimal transport theory, to provide an alternative to traditional schemes formulated in linear spaces. We resort to the JKO variational principle [3, 4, 5] for purposes of time discretization of the diffusive step. This principle characterizes the evolution of the density as a competition between the Wasserstein distance, which penalizes departures from the initial conditions, and entropy, which tends to spread the density and make it uniform over the domain. Remarkably, the resulting update is geometrically exact with respect to advection and volume. We specifically regard the JKO incremental functional as a functional of an incremental transport map which rearranges the density over the time step. Finally, we proceed to discretize the JKO functional in space. This step requires two types of discretization, one for the density and another for the incremental transport map. Exploiting the structure of the Euler-Lagrange equations, which are linear in the density, we may treat said density as a measure and approximate it as a collection of Diracs. The interpolation of the transport map requires more regularity as, in particular, the transport map carries Jacobian information. We specifically resort to the max-ent interpolation scheme proposed by [13]. This interpolation scheme is meshfree and conforming and supplies converging approximations in general $W^{1,p}$ spaces [14].

The paper is organized as follows. In Sec. 2.1, we describe how to reformulate the classical diffusion equations as an optimal transportation problem, specifically how a variational structure naturally emerges from time discretization. In Sec. 2.3, we present the approximation used to discretize the density and the incremental transport map. In Sec. 3 we present numerical schemes that illustrate the scope and properties of the resulting OTM approach. Finally, in the Appendix we collect the main mathematical elements that underpin the present approach.

2 Formulation of the method of approximation

We consider the advection-diffusion initial-boundary-value problem

[TABLE]

where $\rho$ is the unknown density or concentration, $\rho_{0}$ is its initial value, $u$ is a given advection velocity field, $\kappa$ is a diffusion coefficient, $\Omega$ is a bounded domain in $\mathbb{R}^{d}$ and $n$ is the outward unit normal at the boundary.

As stated, problem (1) is defined for densities of sufficient regularity in an appropriate linear space (cf., e. g., [1]). In particular, equation (1) does not make sense for general measures such as Dirac masses and, therefore, it is not in a form suitable for the formulation of particle methods.

In order to eliminate this obstacle, we proceed to reformulate problem (1) as a problem of transport of measures. We begin by noting that (1a) can be equivalently reformulated as the pair of first-order partial-differential equations

[TABLE]

where $v$ is a velocity field that results from the combined effect of advection and diffusion. We identify (2a) as the continuity equation of Eulerian continuum mechanics and (2b) as a mobility law combining the effects of advection and diffusion. The corresponding Lagrangian formulation of the transport problem is

[TABLE]

where $\varphi:\Omega\times[0,T]\to\Omega$ is the transport map. Eq. (3a) simply states that $\rho(x,t)$ is the pushforward of $\rho_{0}(x)$ by the transport map $\varphi(\cdot,t)$, whereas eq. (3b) relates the transport map to the particle velocity.

As stated, the transport problem (2) still requires regularity of the density, but a further weakening of the equations permits extending them to general Radon measures. Thus, (3a) can be reformulated more generally in integral form as

[TABLE]

for all continuous test functions $\eta$ with compact support in $\Omega$. We note that the densities $\rho_{0}$ and $\rho$ now enter linearly into the integrals. Provided that $\varphi$ has sufficient regularity, (4) admits the extension to measures

[TABLE]

where, now, the distribution of the diffusive species is described by a time-dependent measure $\mu(\cdot,t)$ with initial value $\mu_{0}$. Evidently, the regular form (4) of the pushforward operation is recovered when $d\mu_{0}(x)=\rho_{0}(x)dx$ and $d\mu(y,t)=\rho(x,t)dy$ for absolutely continuous densities $\rho_{0}$ and $\rho$, respectively. In addition, expressing (3b) weakly, we obtain

[TABLE]

where $\xi$ are test functions satisfying the boundary condition

[TABLE]

Inserting the mobility law (2b) into (6), gives

[TABLE]

An integration by parts further gives

[TABLE]

where we have used (7). As in the case of the pushforward operation, this weak reformulation of the mobility law admits the extension to measures

[TABLE]

suitable to approximation by particle methods. Eqs. (5) and (10) may be thought as jointly defining an evolution for both the measure $\mu$ and the transport map $\varphi$. In particular, we note that the preceding extension to measures requires consideration of the transport map $\varphi$ as an additional unknown of the problem.

2.1 Time discretization

We begin by approximating problem (2) in time. To this end, let $t_{0}=0<t_{1}<$ $\cdots$ $<t_{N}=T$ be a discretization of the time interval $(0,T)$. The goal is now to determine corresponding discrete approximations $\rho_{0}$, $\rho_{1}$, $\dots$, $\rho_{N}$ of the densities and discrete approximations $\varphi_{0\to 1}$, $\varphi_{1\to 2}$, $\dots$, $\varphi_{N-1\to N}$ of the incremental transport maps.

We begin by noting that eq. (1a) may be regarded as an evolution equation in which the operator is the sum of two operators: advection and diffusion. This additive structure immediately suggests splitting the time integration into advective and diffusive fractional steps. The advective fractional step is governed by the advection equation

[TABLE]

which follows formally from (1a) by setting $\kappa=0$. Conveniently, the advection equation can be solved exactly. Thus, let $\chi:\Omega\times[0,T]\to\Omega$ be the flow corresponding to the given advection velocity field $u$, i. e.,

[TABLE]

Then, the incremental advection of the density $\rho_{k}$ at time $t_{k}$ is the pushforward (#) of $\rho_{k}$ through $\chi_{k\to k+1}$

[TABLE]

where

[TABLE]

is the incremental advection mapping.

Next, we turn to the diffusive fractional step, governed by the diffusive equation

[TABLE]

formally obtained from (1a) by setting $u=0$, and to the temporal discretization of the diffusion step. We characterize the incremental evolution of the density by means of the JKO functional

[TABLE]

where $d_{W}$ is the Wasserstein distance between measures and

[TABLE]

The incremental transport map $\varphi_{k\to k+1}$ then follows from the minimum problem

[TABLE]

We may formally verify that this formulation indeed supplies a time discretization of the mobility law (2b) in the limit of zero advection. Thus, taking variations we obtain (cf. Appendix, eq. (60), also [15])

[TABLE]

where we write

[TABLE]

and the variations satisfy the boundary condition

[TABLE]

Enforcing stationarity for all admissible variations yields

[TABLE]

which is indeed a time discretization of (2b).

We note that the density $\rho_{k+1}$ still enters (19) nonlinearly and differentiated. In order to facilitate particle approximations, we perform an integration by parts leading to the weak stationarity condition

[TABLE]

We verify that, indeed, the density $\rho_{k+1}$ now enters this weak stationarity condition linearly and undifferentiated.

Likewise, the incremental mass conservation relation (2a) must now be understood in a weak or distributional sense, i. e., as the requirement that

[TABLE]

for all test functions $\eta$. Again we verify that the density $\rho_{k+1}$ now enters this weak form of the transport equation linearly and undifferentiated.

2.2 Spatial discretization

Next, we turn to the question of spatial discretization of the weak semi-discrete Fick’s law (23) and the weak form of the transport equation (24). The structure of the latter reveals the need for two types of approximations: i) the discretization of the density $\rho_{k+1}$, and ii) the discretization of the incremental transport map $\varphi_{k+1\to k}$. We consider these two approximations in turn.

2.2.1 Spatial discretization of the density.

As already noted, the density $\rho_{k+1}$ enters (23) and (24) linearly and undiffrentiated. Therefore, a natural and computationally convenient spatial discretization may be effected by considering mass densities of the particle type

[TABLE]

where $x_{p,k}$ represents the position at time $t_{k}$ of a material-point of mass $m_{p}$, \delta\big{(}x-x_{p,k}\big{)} is the Dirac-delta distribution centered at $x_{p,k}$, and $M$ is the number of material-points. For discrete mass distributions of the form (25), eq. (24) reduces to

[TABLE]

which must be satisfied for all test functions $\eta$. Hence,

[TABLE]

i. e., the material-points must have constant mass, and representation (25) reduces to

[TABLE]

with constant $\{m_{p},\ p=1,\dots,M\}$. Thus, the weak and transport reformulation of the problem results trivially in exact mass conservation, simply by keeping the mass of all particles constant.

2.2.2 Spatial discretization of the incremental transport map.

A full spatial discretization additionally requires the interpolation of the incremental transport map $\varphi_{k\to k+1}$. Since $\varphi_{k\to k+1}$ and its variations enter the governing equations (23) and (24) in differential form, the attendant interpolation must be conforming. To this end, we start by simplifying equation (23) by recourse to the change of variables

[TABLE]

with the result

[TABLE]

We consider general linear interpolation schemes of the form

[TABLE]

where $a$ indexes a nodal point set within $\Omega$, $\{N_{a,k},\ a=1,\dots,N\}$ are the corresponding first-order consistent nodal shape functions at time $t_{k}$, and $\{x_{a,k+1},\ a=1,\dots,N\}\equiv x_{k+1}$ is the array of unknown nodal coordinates at time $t_{k+1}$. Consistency here means, specifically, that the shape functions satisfy the identities

[TABLE]

In addition, we interpolate the weight functions as

[TABLE]

and

[TABLE]

Inserting these representations into (30) defines the fully-discrete problem

[TABLE]

However

[TABLE]

or

[TABLE]

whereupon (35) reduces to

[TABLE]

Enforcing this equation for all nodal weights $\{\xi_{a},\ a=1,\dots,N\}$ gives the system of equations

[TABLE]

Finally, inserting the discrete density (28) into (39) gives the particularly simple semi-discrete evolution equation

[TABLE]

with always positive definite mass matrix

[TABLE]

abbreviated $M_{k}\equiv\{\{M_{k,ab},\ a=1,\dots,N\},\ b=1,\dots,N\}$, and nodal flux array

[TABLE]

abbreviated $f_{k}\equiv\{f_{k,a},\ a=1,\dots,N\}$. We note that the mass matrices $M_{k}$ are symmetric.

2.3 Update algorithm

Eq. (40) can be solved for the updated nodal coordinates, with the result

[TABLE]

This scheme amounts to transporting the nodes ballistically over each time step with nodal velocities

[TABLE]

The resulting update for the diffusive fractional step is summarized in the following algorithm.

This forward solution has the usual structure of explicit time-integration and updated-Lagrangian schemes. However, the scheme differs fundamentally from the classical schemes for advection-diffusion equations, which have a linear-space structure. Thus, the present scheme relies on the transport of nodes and mass points to evolve the system. In particular, all the finite kinematics of the motion, including advection and the mass density and volume updates, are geometrically exact.

As in common practice with explicit schemes, the linear solver implied in (43) can be rendered trivial by mass lumping. Thus, for instance, row-sum lumping yields the diagonal mass matrix

[TABLE]

Although this additional approximation strictly deviates from the minimum principle (16) and, therefore, constitutes a variational crime, it is nevertheless adopted for convenience in all numerical examples.

We note that we are left with considerable latitude in the choice of shape functions. In the examples that follow, we choose to use meshfree max-ent shape functions [13] computed from the convected nodal coordinates. The continuous reconstruction of the shape functions has the effect of automatically reconnecting the material-points and the nodal set, with no need of remapping the information carried by the material-points. This property of the method is particularly convenient for nonlinear and history-dependent extensions of the method, in which the local material state often includes additional internal variable information.

From the computational point of view, a challenging aspect of methods based on optimal transportation theory is the periodic redefinition of the nodal and material point neighbor lists, which requires the use of efficient search algorithms [12]. For large-scale problems, the reconstruction of neighbor-list adds overhead to the calculations. We mitigate this overhead by performing neighbor-list updates only when the distortion of a local neighborhood reaches a certain prespecified tolerance. The shape function update is performed at every step.

3 Examples of Application

We present selected numerical simulations that showcase the range and robustness of the approach in examples involving pure diffusion, pure advection, and combined advection-diffusion. In all calculations, the initial nodal-point set is obtained through an auxiliary Delaunay triangulation of the support of the initial density, with material points then placed at the barycenters of the corresponding tetrahedra.

The critical step for the advection step is determined by a Courant condition bases on particle velocity and mean-free path between particles. For the small particle velocities in the examples under consideration, diffusion controls the critical time step. The time-step requirement is estimated with the following relation

[TABLE]

where $\Delta x$ is the minimum distance between two nodes belonging to the same material-point neighborhood. It is also possible to advance the advection and diffusive steps with different time steps using subcycling, but such extensions of the method are beyond the scope of the paper. The time step is recomputed frequently, in particular after the redefinition of the connectivity table. In all examples, diffusion and possibly advection take place concurrently within confined volumes. Neumann boundary conditions are enforced by controlling the position of the nodes. Thus, nodes exiting the domain are repositioned on the boundary. No constraints are applied to the material-points.

3.1 Pure diffusion in a spherical volume

We consider mass diffusion in a spherical volume from an initial configuration in the form of a sphere of radius $R_{0}=1$ with uniform density $\rho_{0}=1$. The mass is discretized into particle sets of three different sizes ranging from $593$ to $4,997$ material points, cf. Table 1. Material-points are subsequently free to diffuse within a sphere of radius $R_{1}=7$.

For the fine mesh, Fig. 1(a) shows the initial configuration of material-points (red) and nodes (green), whereas Figs. 1(b)-1(c) show the final configuration attained by the material points and nodes, respectively. We note that the mass distribution evolves towards a uniform density in which the material points are uniformly distributed and carry ostensibly the same volume, Fig. 1(b). For the three particle-set sizes, Fig. 2(a) displays the time history of the average density of the system and the volume occupied by the material-points. The expected monotonic trend towards uniform density is also clear from these figures. In terms of the minimum principle (16), this trend may be understood as the result of a competition between entropy, which favors a uniform distribution, and the Wasserstein distance, which penalizes the movement of the particles. Fig. 2(b) shows the time history of the radius of the support of the particles, which experiences a steady expansion from its initial value until such time as the particles collide with the boundary of the container. Note that freezing the nodal points at the boundary of the sphere enforces condition (1b). These results are nearly independent of the number of particles over the range considered, which demonstrates the fast convergence of the method and the high accuracy of relatively coarse approximations.

3.2 Advection-diffusion in a hollow cylinder

As a second example, we consider the problem of advection-diffusion in a circular annulus with square cross section and periodic boundary conditions, Fig. 3a. The external radius of the channel is $0.5$ and the internal radius $0.25$. The initial mass density is uniform within a spherical region spanning the cross section of the channel. The density is discretized into 723 nodes and 2,215 material-points. The advection velocity field corresponds to a rigid-body motion at a constant angular velocity of 4. We consider the case of pure advection, $\kappa=0$, and the case of combined advection-diffusion with diffusivity $\kappa=0.001$.

Successive snapshots of the density evolution in the pure-advection problem spanning the short-term part of the solution are shown in Fig. 3. Remarkably, the density function is advected exactly by the algorithm with no numerical diffusion or spurious noise of any type. This type of behavior is, of course, a manifestation of the geometrically exact character of the algorithm and is in sharp contrast to linear-space schemes, which experience difficulty in dealing with advection.

Fig. 4 collects successive snapshots of the advection-diffusion solution. In this case, a certain amount of diffusion tending to spread out the mass distribution is superimposed on the overall circular motion imparted by the advective flow. Eventually, the particles fill the entire channel and the mass density attains a uniform value. The robustness and stability with which the algorithm transports the mass particles and accounts for complex particle-particle and particle-boundary interactions is noteworthy.

Fig. 5 shows the time history of the average material-point density and of the volume occupied by the material-point set for the pure advection and the advection-diffusion problems. As expected, the average density and occupied volume remain exactly constant for pure-advection, again a manifestation of the geometrically exact nature of the method. In the advection-diffusion problem, the average density steadily decreases and, correspondingly, the occupied volume increases with time. The smoothness of the average quantities in time is noteworthy.

A number of considerations of spatial and time resolution need to be carefully born in mind. Firstly, max-ent meshless shape functions such as formulated in [13] are defined through a convex programming problem that is sure to yield solutions when the domain is convex. If, contrariwise, the boundary has concave regions, the nodal point set must resolve the minimum radius of curvature $R$, or feature size, for the max-ent shape functions to be well defined. Thus, we must have

[TABLE]

which places spatial resolution requirements in terms of the geometry of the domain. In addition, as already noted the explicit nature of present formulation places restrictions on the time step, cf. eq. 46. Combining (46) and (47) we additionally have

[TABLE]

which places temporal resolution requirements in terms of the geometry of the domain. The consequences of insufficient resolution are illustrated in Fig. 6, corresponding to the problem of advection diffusion in a circular channel of square cross-section discussed in the foregoing with a larger diffusivity of $\kappa=0.01$, which tightens the resolution requirements. Thus, whereas the material and nodal point sets provide sufficient resolution for a diffusivity $\kappa=0.001$, they fail to do so for the larger diffusivity, with the result that the material and nodal point sets spill over the concave part of the boundary and exit the domain of analysis.

3.3 Advection and diffusion in a bucket

The long-term behavior of the method is illustrated by means of the problem of advection-diffusion in a circular cylinder, Fig. 7. The radius of the cylinder is $0.5$ and its length is $0.25$. The initial mass density is uniform within a spherical region spanning the length of the cylinder. The density is discretized into 723 nodes and 2,215 material-points. The advection velocity field corresponds to a rigid-body motion at a constant angular velocity of 4. We consider three diffusivities of values $\kappa=0.005$, $\kappa=0.0075$, and $\kappa=0.01$.

Fig. 7 shows successive snapshots of the material-point distribution for the case of diffusivity $\kappa=0.0075$. As may be observed from the figure, the material and nodal point set spreads from its initial location and tends towards a configuration of uniform density. The long-term stages of the solution are thus characterized by the fine mixing of the material points. Again, the robust manner in which the method effects the requisite mixing of material points is noteworthy.

Fig. 8 shows the time history of the average material-point density and of the volume occupied by the material-point set for the pure advection and the advection-diffusion problems. As expected, the the average density steadily decreases and, correspondingly, the occupied volume increases, with time. However, two distinct regimes are clearly discernible in the evolution of the system. The first stage corresponds to the free spreading of the initial distribution of mass, followed by a second stage corresponding to the mixing of the material points. Interestingly, mixing is observed to accelerate the diffusion process during the mixing phase of the solution.

4 Concluding remarks

We have developed an Optimal Transportation Meshfree (OTM) particle method for advection-diffusion problems regarded as optimal transport of measures. Thus, in sharp contrast to traditional methods of approximation, which regard advection-diffusion as an evolution of a density function taking place is a suitable linear space, here mass distribution is viewed as a measure evolving in a non-linear space, or manifold. The method follows after the OTM formulation of [12] for fluid and solid flows and represents a scalar version of that method. As in the vector OTM method, we resort to the incremental variational principle of Jordan, Kinderlehrer and Otto, [3, 4, 5], for purposes of time discretization of the diffusive step. This principle characterizes the evolution of the density as a competition between the Wasserstein distance between two consecutive densities and entropy. By virtue of the structure of the Euler-Lagrange equations, sufficiently weakened so as to be linear and undifferentiated in the density, we may approximate the density as a collection of Dirac masses. The interpolation of the incremental transport map, which needs to be conforming, is effected through mesh-free max-ent interpolation [13]. Remarkably, the resulting update is geometrically exact with respect to advection and volume, which overcomes the chronic problems of spurious diffusion and oscillations that afflict linear-space methods. We present three-dimensional examples of application that illustrate the scope and robustness of the method.

The addition of distributed sources and sinks is straightforward and it simply requires allowing the mass of each particle to change with time. Non-zero Neumann boundary conditions can also be implemented simply by injecting particles or allowing particles to escape through the boundary. Finally, Dirichlet boundary conditions can be enforce by placing fixed particles consistent with a prescribed density on the boundary.

The accuracy and robustness of the OTM method have been assessed by means of selected advection-diffusion test problems. The numerical tests exhibit clearly the geometrically exact nature of the method with respect to advection and volume. Thus, advection is handled directly by means of a pushforward of the mass density by the advection flow. In addition, the volume of the material points is updated by means of the Jacobian of the incremental transport map. Both operations are geometrically exact, in contrast to the conventional Eulerian treatment of advection and volume through upwinding and the solution of a Poisson problem. The numerical solutions generally exhibit a remarkable robustness and accuracy, even for coarse approximations, enabling, for instance, the simulation of mixing of phases and the propagation of density fronts.

We have noted that the method, as presented here, is subject to spatial and temporal resolution requirements, in particular as regards the need to resolve concave regions of the boundary. However, it should be carefully noted that said resolution requirements are directly tied to the particular implementation adopted here and are not inherent to the optimal transport treatment of advection-diffusion problems. Thus, the convex-programming nature of max-ent interpolation, and its relation to convex domains, may be relaxed by considering signed shape functions [16]. In addition, the explicit character of the implementation presented here can be generalized by considering incremental minimum principles, other than JKO (16), in the spirit of variational integrators [17, 18, 19, 20]. This extension enables, in particular, the formulation of implicit and asynchronous time-integration schemes.

Finally, we note the connection between transport problems and gradient flows Wasserstein spaces (e. g., [21, 22]). This connection suggests the applicability of the present approach to broad classes of problems in mechanics and physics, including flow in porous media, Oswald ripening, the Cahn-Hilliard equation, and others. These generalizations and extensions define worthwhile directions for further study but are beyond the scope of this paper.

\ack

LF and MO gratefully acknowledge the support of the U.S. National Science Foundation through the Partnership for International Research and Education (PIRE) on Science at the Triple Point Between Mathematics, Mechanics and Materials Science, Award Number 0967140.

Appendix A Elements of optimal transportation theory

For completeness, in this appendix we collect basic elements of optimal transportation theory that underpin the Optimal Transportation Meshfree (OTM) method presented in the foregoing. A complete account of the theory may be found, e. g., in [15, 23].

Let $\Omega\subseteq\mathbb{R}^{d}$ an open convex set and let $\mathscr{P}(\Omega)$, the space of Borel probability measures on $\Omega$. If $\mu\in\mathscr{P}(\Omega)$ and $\varphi:\Omega\rightarrow\Omega$ is a Borel map, $\varphi_{\#}\mu\in\mathscr{P}(\Omega)$ denotes the push-forward of $\mu$ through $\varphi$, defined by

[TABLE]

Let $\pi^{i},\,\,i=1\dots n$ be the canonical projection operator from a product space $\Omega\times\dots\times\Omega$ onto $\Omega$,

[TABLE]

Given $\mu_{1}\in\mathscr{P}(\Omega)$ and $\mu_{2}\in\mathscr{P}(\Omega)$ the class $\Gamma(\mu_{1},\mu_{2})$ of transport plans or couplings between $\mu_{1}$ and $\mu_{2}$ is defined as

[TABLE]

We denote by $\mathscr{P}_{2}(\Omega)$ the space of Borel probability measures with finite second moment: $\mu\in\mathscr{P}(\Omega)$ belongs to $\mathscr{P}_{2}(\Omega)$ iff:

[TABLE]

For every pair of measures $\mu,\nu\in\mathscr{P}_{2}(\Omega)$ we consider the transport cost

[TABLE]

By the direct method of calculus of variations, the minimum problem (53) admits at least a solution. The minimum value $d_{W}(\mu,\nu)$ defines a distance between the measures $\mu,\nu\in\mathscr{P}_{2}(\Omega)$ and the metric space $(\mathscr{P}_{2}(\Omega),d_{W})$ is referred to as the ($L^{2}$-) Wasserstein space on $\Omega$. We denote by $\mathscr{P}_{2}^{a}(\Omega)$ the subset of $\mathscr{P}_{2}(\Omega)$ formed by the absolutely continuous measures with respect to the Lebesgue measure,

[TABLE]

The following result establishes the existence and uniqueness of optimal transport plans induced by maps, or optimal transport maps, got the case in which the initial measure $\mu$ belongs to $\mathscr{P}_{2}^{a}(\Omega)$.

Theorem A.1** (Existence and uniqueness of optimal transport maps [24, 25] ).**

For any $\mu\in\mathscr{P}_{2}^{a}(\Omega)$ and $\nu\in\mathscr{P}_{2}(\Omega)$, the Kantorovich optimal transport problem (53) has a unique solution $\gamma$, which is concentrated on the graph of a transport map, namely, the unique minimizer of Monge’s optimal transport problem

[TABLE]

The map $\varphi$ is cyclically monotone and there exists a convex function $\phi:\Omega\rightarrow\mathbb{R}$ such that $\varphi(x)=\nabla\phi(x)$ for $\mu-a.e.\,x\in\Omega$.

We also recall that the manifold $\mu\in\mathscr{P}_{2}(\Omega)$ can be endowed with a canonical differential structure. In particular, the tangent space to $\mathscr{P}_{2}(\Omega)$ at the point $\mu$ is defined as

[TABLE]

where the elements of $\textnormal{Tan}_{\mu}\mathscr{P}_{2}(\Omega)$ may be regarded as velocity fields. Indeed, for any absolutely continuous curve $t\mapsto\mu_{t}\in\mathscr{P}_{2}(\Omega)$ with $\mu_{0}=\mu$, there exists a unique $\xi\in\textnormal{Tan}_{\mu}\mathscr{P}_{2}(\Omega)$ such that

[TABLE]

In particular, if $\mu_{t}\in\mathscr{P}_{2}^{a}(\Omega)$, and $\rho_{t}$ is its density with $\rho_{0}=\rho$, then

[TABLE]

for all test functions $\eta\in\mathcal{D}(\Omega)$. Thus, changes in densities in $\mathscr{P}_{2}^{a}(\Omega)$ are to be regarded as mass fluxes induced by velocity fields.

We also recall the following differentiation formula. Let $b>a$, $\rho_{a}$, $\rho_{b}\in\mathscr{P}_{2}^{a}(\Omega)$ and let $\rho_{t}\in\mathscr{P}_{2}^{a}(\Omega)$, $t\in[a,b]$, be such $\rho(\cdot,a)=\rho_{a}$, $\rho(\cdot,b)=\rho_{b}$ and

[TABLE]

for some $C^{1}$, globally bounded, velocity field $\xi$. Then [15],

[TABLE]

where $\varphi_{\rho_{b}}$ is the unique optimal transport map between $\rho_{a}$ and $\rho_{b}$.

We conclude by showing that the first variation in $\mathscr{P}_{2}(\Omega)$ of the entropy

[TABLE]

gives (19) and (22). Thus, let $\xi\in C^{\infty}_{c}({\Omega})$ be a smooth vector field with compact support. Define the variation $y=\psi(x,t)=\psi_{t}(x)$ as

[TABLE]

Let $\rho_{t}$ be the pushforward measure of $\rho_{k+1}$ by $\psi_{t}$, i. e.,

[TABLE]

for all $\eta\in C^{\infty}_{0}(\Omega)$. Since $\psi_{t}$ is regular, (63) is equivalent to the relation between densities

[TABLE]

In particular,

[TABLE]

and $\rho_{t}$ is in $\mathscr{P}_{2}(\Omega)$. In addition, we have

[TABLE]

and

[TABLE]

Write

[TABLE]

Then, we have

[TABLE]

where we have used that $\psi(x,0)=\text{id}$. Taking the limit $t\rightarrow 0$, we find

[TABLE]

as advertised.