{'fraction_non_alphanumeric': 0.04140038414527676, 'fraction_numerical': 0.011838658983761132, 'mean_word_length': 4.556514989812749, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Constrained particle-mesh projections in a hybridized discontinuous Galerkin framework with applications to advection-dominated flows

By combining concepts from particle-in-cell (PIC) and hybridized discontinuous Galerkin (HDG) methods, we present a particle-mesh scheme which allows for diffusion-free advection, satisfies mass and momentum conservation principles in a local sense, and allows the extension to high-order spatial accuracy. To achieve this, we propose a novel particle-mesh projection operator required for the exchange of information between the particles and the mesh. Key is to cast these projections as a PDE-constrained $\ell^2$-optimization problem to allow the advective field naturally located on Lagrangian particles to be expressed as a mesh quantity. By expressing the control variable in terms of single-valued functions at cell interfaces, this optimization problem seamlessly fits in a HDG framework. Owing to this framework, the resulting scheme can be implemented efficiently via static condensation. The performance of the scheme is demonstrated by means of various numerical examples for the linear advection-diffusion equation and the incompressible Navier-Stokes equations. The results show that optimal spatial accuracy can be achieved, and given the particular time-stepping strategy, second-order time accuracy is confirmed. The robustness of the scheme is illustrated by considering benchmarks for advection of discontinuous fields and the Taylor-Green vortex instability in the high Reynolds number regime.


Introduction
Flow and transport phenomena in which advection plays a dominant role, and that can be modeled as incompressible, arise in many engineering applications, such as the advection and mixing of pollutants, turbulent flows, and multiphase flows. Accurately simulating such problems requires a discretization of the underlying mass and momentum conservation laws having minimal artificial dissipation relative to the physical damping, as well as imposing the incompressibility constraint. Whereas Lagrangian particlebased methods offer a diffusion-free treatment of advective phenomena, pure, Lagrangian particle motion consistent with the incompressibility constraint is difficult, if not impossible to achieve. On the other hand, Eulerian mesh-based methods are amenable to imposition of incompressibility and possess global or local (i.e. cellwise) conservation properties. However, in the advection dominated regime, the stabilization of mesh-based methods is non-trivial and typically introduces artificial diffusion. By combining concepts from Lagrangian particle-based and Eulerian mesh-based methods, this paper formulates a novel particle-mesh scheme that reconciles a diffusion-free discretization of the advective part of the problem with an efficient discretization of the constitutive part, while rigorously satisfying local conservation principles. Furthermore, the method extends to arbitrarily high-order spatial accuracy on structured or unstructured meshes.
Tracing back to the 1960s, the particle-in-cell (PIC) method [1] was the first to consider Lagrangian particles in the advective part of a problem, while conveniently using a mesh to account for the dynamic interaction between the particles. The coupling between the particles and the mesh is established by means of two auxiliary projection steps: a particle-mesh projection for reconstructing field data from the scattered particles, and a mesh-particle projection to update the properties of the Lagrangian particles from the solution on the mesh. Conservation in PIC and related methods, such as the material-point method (MPM) [2], is generally examined in the context of global conservation [3,4,5,6]. The salient principle is that the particle discretization of the solution and the mesh-based discretization of the solution should contain the same total mass, momentum and energy, as these are two representations of the same material. Indeed, exact conservation of these quantities is possible when using basic, low-order methods in which particles are used as moving point masses or volumes [5]. However, this does not yet insure that the particle fields are locally consistent with the mesh description of the flow. For example, the mass or momentum represented by the particles within a discretization mesh cell does not typically match the mass and momentum associated with the cell. Extension to higher-order accuracy is possible by considering the particles as moving sampling points of the material at the mesh [7,8,9,10]. While being consistent, this compromises however exact conservation [8,10]. This trade-off between conservation and accuracy is also observed in other particle methods such as SPH [11,12].
In this work, as in related previous work [10], we frame our particle-mesh scheme as an operator splitting method and use particles as moving sampling points for reconstructing fields at the mesh. By extending the concepts explored in [10], this paper formulates particle-mesh projection operators which provide high-order accuracy, while simultaneously satisfying rigorous global and local conservation principles on the mesh.
Not surprisingly, a minimal requirement for obtaining the latter is to employ a mesh-based framework exhibiting local conservation properties. This requirement is met by many meshed-based discretizations, and in particular is met by discontinuous Galerkin (DG) methods (see, e.g., [13,14] and references) and hybridized discontinuous Galerkin (HDG) methods (see, e.g., [15,16,17,18,19,20]). Typical to (H)DG is to pose the problem in terms of cellwise balances augmented with properly defined numerical fluxes at cell interfaces. These numerical fluxes are usually expressed in terms of the element field variables, and require upwinding (e.g., [15,17,18,19,20]) or other carefully designed formulations (e.g., [16,21]) in order to obtain stable solutions in advection dominated regimes. What distinguishes HDG from DG is that the former employs single-valued interface functions which are defined on cell facets only. Typically, the interface flux in HDG is partly formulated in terms of these interface fields, thereby providing the conditions to enforce inter-element flux continuity in a weak sense. Essential to our novel particle-mesh projection is the reconstruction of HDG interface fields from the advection of Lagrangian particles. In this way, we inherit the diffusion-free properties of a particle-based method while retaining the local conservation properties of a HDG method. The key idea in formulating the particle-mesh projection is that the particle motion has to satisfy a mesh-based advection operator.
This advection operator is used as a constraint to augment the local 2 -projections used in [10], resulting in a PDE-constrained minimization problem for the particle-mesh projections. By expressing the control variable in terms of single-valued interface functions, this approach seamlessly fits in a HDG framework.
The HDG framework provides an additional benefit; it enables an efficient implementation procedure via static condensation. We will present our particle-mesh scheme for the advection-diffusion equation and the incompressible Navier-Stokes equations, where the inclusion of the diffusion terms and/or the incompressibility constraint is formulated as an operator splitting scheme. Considering the simpler case of advection-diffusion allows analysis of our scheme, before using our approach for more demanding problems, including the challenging case of vortex instability in high Reynolds number flows. Other potential applications, beyond the scope of this paper, include the sharp representation, using particles, of fluid-fluid interfaces in multiphase flows.
The remainder of this work is structured as follows. In Section 2 the problems of interest are formulated and a spatial-temporal operator splitting technique is introduced. Formulating the PDE-constrained particlemesh projection, and proving consistency and conservation constitute the main part of Section 3, leading to semi-discrete formulations for the advection-diffusion equation and the incompressible Navier-Stokes equations. Section 4 presents the corresponding fully-discrete formulations and discusses various algorithmic aspects of the resulting numerical schemes. The performance of the method is assessed for various benchmark tests in Section 5. The tests demonstrate, for an appropriate choice of basis functions, optimal second and third-order convergence rates in space, and, given the particular time stepping strategy, second-order accuracy is obtained in time. Conclusions and an outlook for future research are presented in Section 6.

Problem formulation and definitions
We define the scalar-valued linear advection-diffusion equation and the incompressible Navier-Stokes equations on a domain Ω ⊂ R d (with d = 2, 3), having Lipschitz continuous boundary Γ = ∂Ω. The boundary Γ is partitioned into Dirichlet and Neumann boundaries Γ D and Γ N , respectively, satisfying The outward pointing unit vector normal to Γ is denoted by n. The time interval of interest is I = t 0 , t N , where t 0 and t N are the start and end time of the simulation, respectively.

Advection-diffusion equation
On the space-time domain Ω × I, the scalar-valued linear advection-diffusion equation is defined as a system of two first-order equations as follows: given a solenoidal velocity field a : Ω × I → R d , an initial condition φ 0 : Ω → R, diffusivity κ, and boundary conditions h : Γ N × I → R and g : Γ D × I → R, find the scalar quantity φ : Ω × I → R such that The parameter γ in Eq. (1c) is equal to one at inflow boundaries (where a · n < 0) and equal to zero on outflow Neumann boundaries (where a · n ≥ 0), with h specifying the total flux on inflow Neumann boundaries and the diffusive flux on outflow Neumann boundaries.

Incompressible Navier-Stokes equations
The incompressible Navier-Stokes equations on the space-time domain Ω × I are stated as follows: given the kinematic viscosity ν, the forcing term f : Ω × I → R d , the boundary conditions g : Γ D × I → R d and h : Γ N × I → R d , and a solenoidal initial condition u 0 : Ω → R d , find the velocity field u : Ω × I → R d and pressure field p : Ω × I → R such that where σ is the momentum flux tensor, I is the identity tensor, and ∇ s u = 1 2 ∇(u)+ 1 2 ∇(u) is the symmetric velocity-gradient tensor. Similar to its use in the advection-diffusion equation, the parameter γ is equal to one on inflow parts of Γ N (that is, where u · n < 0), imposing in this way the total momentum flux. On the outflow parts of Γ N (where u · n ≥ 0) γ equals zero, thus prescribing the diffusive part of the momentum flux only.

Operator splitting
As in [10], the particle-mesh method is conceived as an operator splitting procedure. To this end, let the time interval of interest I be partitioned using a sequence of N + 1 discrete time levels {t 0 , t 1 , . . . , t N −1 , t N } which for n = 0, N − 1 define the half-open subintervals I n = (t n , t n+1 ] such that n I n = I, while I := {I n } defines the ordered sequence of subintervals. Furthermore, let the total fluxes (given by Eqs. (1b) and (2b)) be decomposed into advective parts σ a and diffusive parts σ d .
A spatial-temporal operator splitting procedure for the advection-diffusion problem, Eq. (1), now involves a scalar field ψ : Ω × I n → R satisfying an advection problem, and a scalar field φ : Ω × I n → R satisfying a diffusion problem, to be applied sequentially for every I n ∈ I. In the advection stage the Dirichlet boundary condition can only be prescribed on inflow Dirichlet boundaries, denoted with Γ − D . The flux prescribed on Γ N is split into an advective flux h a and a diffusive flux h d , such that h = h a + h d . Note that the advective flux cannot be specified at outflow Neumann boundaries, which is automatically taken care of by virtue of Eq. (3c). Furthermore, P L and P E are projection operators, which are introduced in order to couple the fields ψ and φ, with these fields being naturally defined on the particles and the mesh, respectively, in PIC methods.
More precisely, the projection operator P L provides the initial condition at t n to advance the Lagrangian advection problem to t n+1 , and the projection operator P E provides the initial condition at t n to advance the Eulerian diffusion problem to t n+1 . To remain consistent with Eq. (1) it is required that P E • P L equals the identity operator.
In a similar fashion, an operator splitting procedure for the incompressible Navier-Stokes equations (Eq. (2)) is formulated by finding a vector-valued field v : Ω × I n → R d satisfying an advection problem, and a velocity field u : Ω × I n → R d satisfying an incompressible Stokes problem, with notation similar to Eqs. (3)(4). The advective field U in Eq. (5b) is not yet specified, other than to require this field to be a consistent approximation to u which is piecewise constant on I.
We aim to formulate the projection operators P L (in Eqs. (3e) and (5e)) and P E (in Eqs. (4e) and (6f)) such that exact conservation is guaranteed, in the sense that the projected field satisfies Eq. (3a) or Eq. (5a) in an integral sense over each cell.

Definitions
We next introduce some notation needed for the spatial discretizations of the Eulerian and Lagrangian subproblems, respectively, and the definition of the associated projections P E and P L . The fixed Eulerian is denoted by F.

Function spaces
The following scalar finite element spaces are defined: in which P (K) and P (F ) denote the spaces spanned by Lagrange polynomials on K and F , respectively,

Particle definitions
The Lagrangian particle configuration in the domain Ω at a fixed time instant t is defined as follows in which x p denotes the spatial coordinates of particle p and N p is the number of particles.
Furthermore, the index set of all particles and the index set of particles hosted by cell K, at a fixed time instant t are defined as, respectively Finally, Lagrangian scalar and vector fields on the particles are defined as, respectively, where ψ p and v p (t) denote the corresponding scalar and vector quantities associated with particle p.
Importantly, subscripts p and h are consistently used to distinguish between Lagrangian particle data and Eulerian mesh fields, respectively.

Semi-discrete formulations
We provide now a collection of semi-discrete forms that together constitute a particle-mesh operator splitting method for the advection-diffusion equation (Eqs. (3)(4)) and for the incompressible Navier-Stokes equations (Eqs. (5)(6)). The splitting procedure involves the following sequence of steps (referring to the relevant equations for the advection-diffusion operator splitting for the sake of illustration): 1. Lagrangian discretization of the advection problem, in order to solve Eqs. (3a -3d) at the particle level; 2. particle-mesh projection, in order to project the scattered particle data onto a (scalar) field on the Eulerian mesh using the operator P E : Ψ t → W h (see Eq. (4e)); 3. Eulerian discretization of the diffusion equation, in order to solve Eqs. (4a -4d) on the mesh; 4. mesh-particle projection, in order to update the particle properties from the scalar field on the Eulerian mesh using the operator P L : W h → Ψ t (see Eq. (3e)).
The resulting semi-discrete components for the advection-diffusion problem and the incompressible Navier-Stokes problem are presented in Sections 3.1 and 3.2, respectively. Since the formulation of conservative projection operators is our main objective we primarily focus on steps 2 and 4, referring to Maljaars et al. [10] for a discussion of steps 1 and 3.

Advection-diffusion equation
3.1.1. Lagrangian discretization of the advection problem In a Lagrangian, particle-based frame of reference, solving the advection problem Eqs. (3a-3d) on the time interval I n is straightforward by decomposing the problem into two ordinary differential equations for the particle scalar quantity and the particle position, given bẏ whereψ p (t) andẋ p (t) are the total derivatives at time t ∈ I n of the scalar quantity and the position of particle p, respectively. Furthermore, a(x, t n ) is a prescribed solenoidal velocity field at time t n . An important observation is that ψ p stays constant throughout the particle advection stage by virtue of Eq. (17a).

PDE-constrained particle-mesh interaction
In [10] the projection operator P E , transferring the scattered particle data to piecewise continuous fields at the mesh, is formulated in terms of local least-squares projections obtained from the minimization problem With W h a discontinuous function space, this approach allows for an efficient cellwise implementation, and gives accurate results provided that the particle locations satisfy unisolvency with respect to W h . However, conservation of linear quantities is not guaranteed after projection. We now extend the variational framework of that approach by imposing constraints as to ensure conservation throughout the particle-mesh projection.
Recalling that the Lagrangian particles are used for solving an advection operator, it is argued that the projected particle motion should also satisfy an advection operator from the perspective of the Eulerian mesh. Therefore, we constrain the particle-mesh projection by a hyperbolic conservation law.
To this end, we augment the functional in Eq. (18) with terms multiplying Eq. (3) with a Lagrange multiplier λ h ∈ T h . Integration by parts leaves an unknown flux on interior cell facets which is formulated in terms of a control variableψ h ∈W h,g , while the Neumann boundary condition is substituted on exterior facets. For a given particle field ψ p ∈ Ψ t , an advective velocity field a : Ω × I n → R d , the initial condition ψ n h ∈ W h , and an advective Neumann boundary condition h a : Γ N × I n → R, the minimization problem then involves finding the stationary points of the Lagrangian functional for every t ∈ I n . The collection of terms containing λ h constitutes a weak form of the advection subproblem, Eq. (3). Furthermore, the unknown facet-based fieldψ h ∈W h,g determines the interface flux, the implications of which will become clear later on. The additional term containing β > 0 penalizes the jumps between ψ h andψ h on cell interfaces thereby avoiding the problem of becoming singular in cases with vanishing normal velocity a · n.
Equating the variations of Eq. (19) with respect to the three unknowns ψ h , λ h ,ψ h ∈ W h , T h ,W h,g to zero, results in the following system of variational equations. At time t ∈ I n , variation with respect to the scalar field ψ h gives the co-state equation Variation with respect to the Lagrange multiplier λ h gives the state equation, Finally, variation with respect to the control variableψ h leads to the optimality condition, Note that the optimality system Eq. (20) is naturally formulated using a HDG approach, since we require the optimal solution for the state variable ψ h to be in the discontinuous function space W h , and the optimal control to be provided by the single-valued, facet-based control variableψ h ∈W h,g . After an appropriate discretization of the time derivatives in Eqs. (20a and 20b) -to be discussed in Section 4 -a field ψ h ∈ W h can be reconstructed from the particle data ψ p ∈ Ψ t by solving the optimality system Eq. (20). This reconstructed field is used to provide the initial condition for the subsequent diffusion subproblem.

Eulerian discretization of the diffusion equation
We next present the semi-discrete formulation of the diffusion problem Eq. (4). To this end, we will use the HDG method presented by Labeur & Wells [15] to seek solutions φ h ∈ W h . This choice allows for trivial projections between fields ψ h and φ h , which seamlessly fits in the approach used for solving the optimality system Eq. (20). Referring to [15] for further details, the HDG discretization results in a set of local and global problems which are respectively stated as: at time t ∈ I n , given the initial condition φ n h ∈ W h , the diffusive Neumann boundary condition h d : and whereσ d,h is a diffusive flux at cell facets defined bŷ in which α is a dimensionless parameter as is typical to interior penalty methods [22]. Solving Eq. (21) yieldsφ h and φ h , where the latter is to be used for updating the particle properties ψ p in the subsequent mesh-particle projection step.

Mesh-particle projection
The mesh-particle projection P L : W h → Ψ t is based on the following minimization problem min ψp(t) where we emphasize that the objective functional J is also at the basis of the particle-mesh projection Eq. (19). Carrying out the minimization yields the particularly simple result The mesh-particle projection is not restricted to the mapping of the scalar field φ h ∈ W h only, but can be applied to project arbitrary fields in W h -e.g. the temporal increments of φ h , see Section 4.1.4 -onto the particles. As such, it provides the tool for updating the particle quantities, which completes the semi-discrete sequence of steps comprising the particle-mesh operator splitting of the advection-diffusion equation.

Incompressible Navier-Stokes equations
In the same vein, we consider a particle-mesh operator splitting for the incompressible Navier-Stokes equations. Details of the semi-discrete operator splitting procedure, which involves a non-linear advection step followed by an unsteady Stokes step, can be found in [10]. The main difference between the approach in the aforementioned reference and the method presented here, is the PDE-constrained particle-mesh projection performed in step 2. For completeness, we also briefly recapitulate the remaining model components.
To set the stage, some additional notation is introduced. The momentum field at the Eulerian mesh, to be reconstructed from the particle field, is denoted by v h ∈ W h . Throughout the advection stage we will make use of advective fields

Lagrangian discretization of the advection problem
The particle advection step for the incompressible Navier-Stokes equations proceeds along similar lines as for the linear advection-diffusion equation (Eq. (17)). Particle specific momenta and positions are updated by integratingv over the time interval I n . In this equation x p (t) denotes the particle position at time t ∈ I n . We still leave the advective field U h ∈ W h unspecified, but recall that it is a consistent approximation of u h which is piecewise constant in I. Similar to the scalar-valued case, v p stays constant throughout the advection stage by virtue of Eq. (24).

PDE-constrained particle-mesh interaction
The formulation of the vector-valued PDE-constrained projection operator P E : V t → W h of step 2 proceeds along similar lines as for the linear advection-diffusion equation. That is, the functional constituting the least-squares minimization problem (as in Eq. (18)) is augmented with Lagrange multiplier terms that enforce weak satisfaction of the momentum advection equation by the projected field v h . Integration by parts of the constraint introduces a facet-based variablev h controlling the flux across cell interfaces. The resulting optimality system can be stated as follows: at time t ∈ I n , given the initial condition v n h ∈ W h , the particle field v p (t) ∈ V t , the advective field U h ,Ū h ∈ W h ,W h and the advective Neumann boundary It is worth mentioning the similarity of Eq. (26) with Eq. (20) for the linear scalar advection problem, owing to the use of the advective fields U h andŪ h which linearizes the momentum fluxes, thereby simplifying the optimality system considerably.
The PDE-constrained projection step leads to a field v h ∈ W h to be used in the subsequent unsteady Stokes problem.

Unsteady Stokes equations
The semi-discrete Stokes problem is formulated in terms of a velocity field u h ∈ W h which enables a trivial projection between the fields u h and v h . This practically implies that we will resort to the class of (H)DG methods for the spatial discretization of the unsteady Stokes equations. More specifically, the HDG method from Rhebergen & Wells [23,20] is employed. In a particle-mesh setting, the restriction on the choice of function spaces advocated by these authors has some distinct advantages compared to the general formulation of Labeur & Wells in earlier work [18]. This will be further discussed in Section 3.3.3.
Following [20], the semi-discrete form of the unsteady incompressible Stokes problem on the time interval I n is stated as follows: at time t ∈ I n , given the initial condition u n h ∈ W h , a forcing term f : Ω → R d , the diffusive Neumann boundary condition h d : Γ N → R d and the viscosity ν, find (u h ,ū h , p h ,p h ) ∈ W h ,W h,g , Q h ,Q h satisfying the local and global momentum balances, respectively, and, simultaneously, the local and global mass balances, respectively, in which the diffusive interface fluxσ d,h in Eq. (27a) and Eq. (27b) is defined aŝ with I the identity tensor, and α a parameter as is typical to interior penalty methods [22].
Solving Eq. (27) for (u h ,ū h , p h ,p h ) yields the divergence-free velocity field u h ∈ W h which is to be used in the subsequent mesh-particle projection step.

Mesh-particle projection
Next, the vector-valued counterpart of Eq. (22) is used for the mesh-particle projection P L : Thus, the specific momentum v p associated with particle p is updated via This mesh-particle projection concludes the sequence of steps for a particle-mesh operator splitting of the incompressible Navier-Stokes equations.

Properties of the semi-discrete formulation
To better illustrate the properties of the PDE-constrained particle-mesh projection, we now prove consistency and conservation of the above presented scheme. In addition, we prove consistency and conservation of the resulting particle-mesh operator splitting schemes for the advection-diffusion equation and the incompressible Navier-Stokes equations.

Consistency
Consistency of the particle-mesh operator splitting scheme entails three different aspects, (i) the operator splitting has to be consistent with the unsplit governing equations (Eqs. (1) or (2)), (ii) the Eulerian part (step 3) has to be consistent with the diffusion subproblem (Eqs. (4) or (6)), and (iii), the constraint imposed weakly in the projection operator P E (step 2) has to be consistent with the advection subproblem (Eqs. (3) or (5)).
To start with (i), splitting of the advection-diffusion equation (Eq. (1)) and the incompressible Navier-Stokes equations (Eq. (2)) into a kinematic part (advection problem) and a constitutive part (diffusion equation or Stokes problem) has been the subject of numerous studies, and is known to be consistent up to a time step dependent splitting error which vanishes in the continuous time limit (see, e.g., [24,25,26] among many others). However, specific to a particle-mesh framework, we also have to require consistency of the projection operators P E (particle -mesh) and P L (mesh -particle). This implies that subsequent application of the projection operators P L and P E must recover the initial scalar-valued mesh field ψ h exactly. Similarly, provided that the particle field ψ p is initialized by means of the mesh-particle projection P L , subsequent application of the projection operators P E and P L must recover the initial particle field ψ p exactly. Mathematically, these compatibility conditions can be formulated as, respectively which implies mutual consistency of steps 2 and 4 in the absence of steps 1 and 3. To verify these conditions for the presented formulations for P E and P L , consider the scalar-valued projections. For a repeated backand-forth mapping between particles and mesh, hence omitting step 1 (i.e. a = 0) and step 3 (i.e. κ = 0), the constraint will be inactive in Eq. (20), i.e. λ h = 0 everywhere, andψ h and ψ h coincide in a weak sense at cell facets. As a result, the only term remaining in the particle-mesh projection stems from the 2 -term in the Lagrangian functional L given by Eq. (19). Hence, the functional L underpinning the particle-mesh projection is in this case similar to the objective function J in Eq. (18). By recalling that a trivial projection between ψ h ∈ W h and φ h ∈ W h can be established, it follows that the mesh-particle projection (Eq. (22)) is based on the same objective function, and only differs in terms of the minimizer. Owing to this symmetry Eq. (29) is satisfied, which shows consistency of the projections in the semi-discrete setting.
Concerning (ii), consistency of the HDG method used in the diffusion step was proven in [15]. For a consistency proof of the HDG discretization of the Stokes equations, reference is made to [18,23].
It remains to prove consistency of the PDE-constrained projection (iii). To this end, consider a sufficiently smooth scalar field ψ. Substitution into Eq. (20b) gives, after integration by parts, which demonstrates consistency with the strong form of the advection problem Eqs. (3a-3b) and the Neumann boundary condition Eq. (3c), with the enforcement ofψ = ψ on interior cell facets and the Dirichlet boundary Γ D .

Conservation
We will prove that the scheme is mass and momentum conservative from the perspective of the Eulerian mesh. Global and local conservation of the HDG method for the Eulerian diffusion step was demonstrated in [15]. Global and local conservation of the HDG method for the unsteady Stokes step was demonstrated in [18,23,20]. It therefore remains to prove conservation of the PDE-constrained particle-mesh projection.
Setting δλ h = 1 in Eq. (20b), and rearranging, leads to For a point-wise divergence free vector field a the boundary integral on the union of interior cell facets vanishes, due toψ h being single-valued on facets F ∈ F. The right-hand side therefore equals the total ingoing advective flux at the exterior boundary Γ , thereby proving global mass conservation.
For local mass conservation, setting δλ h = 1 on cell K and δλ h = 0 on Ω \ K gives, after rearrangement, The right-hand side of Eq. In this work particles carry point evaluations of the underlying field with no notion of mass or volume.
As such, a notion of conservation is not well-defined at the particle level. Therefore, the current approach can only pursue a notion of conservation at the mesh level.

Function spaces
To are H(div)-conforming and pointwise divergence-free [20]. In the scope of particle-mesh methods, this has a distinct advantage in that particles can be advected through velocity fields being pointwise-divergence free within a cell and having a continuous normal velocity component across cell facets. In [10], it was shown that satisfying these two criteria is paramount in maintaining a uniform particle distribution over time. We will further exploit this feature in Section 4.2.1.

Fully-discrete formulations
The final step preceding a computer implementation involves discretization of the respective semi-discrete model components in time. The particular time stepping strategy presented in this section largely follows the approach used in [10], albeit special care is required to render the particle-mesh projection compatible with the mesh-particle projection. This section concludes by discussing some algorithmic aspects. In particular, we evaluate the choice for the Lagrange multiplier space and present an efficient solution procedure for the PDE-constrained particle-mesh interaction by using static condensation.

PDE-constrained particle-mesh projection
A fully-discrete PDE-constrained particle-mesh projection is formulated with the objective to optimize the scalar field ψ h at time level n+1 given the particle field ψ p at time t n+1 . Employing the θ-method, where 1/2 ≤ θ ≤ 1, the constraint in the Lagrangian functional L is evaluated at time t n+θ := (1 − θ)t n + θt n+1 using linear interpolation between discrete time levels. The scalar field ψ h at time level n + θ is then approximated by in which ψ * ,n h ∈ W h is an initial field given by where ∆t n = t n+1 − t n is the time step size, θ L is an additional time stepping parameter (1/2 ≤ θ L ≤ 1, but possibly different from θ), and the incrementsφ m h (with m = n − 1, n) are defined bẏ with φ m h and φ * ,m h being fields in W h related to the Eulerian diffusion step, to be formulated shortly. Using these increments in Eq. (34), insures compatibility of the fully-discrete projection operator P L with P E which is required for consistency, as discussed in Section 3.3.1 and to be elaborated further in Section 4.3.2.
The time derivative of the scalar field ψ h at time level n + θ is now given by since variations δψ * ,n h ∈ W h vanish. Given these approximations, the fully-discrete co-state equation reads: given the particle field ψ n p ∈ Ψ t , the post-advection particle positions x n+1 p ∈ X t , and the intermediate field Correspondingly, the fully-discrete counterpart of the state equation Eq. (20b) reads: Finally, the fully-discrete optimality condition becomes In these equations, The Lagrange multiplier λ h and the control variableψ h are conveniently taken at time level n + 1, which is allowed since these variables are fully-implicit, not requiring differentiation in time. and in which the fully-discrete diffusive fluxσ n+1 d,h is given bŷ Solving Eq. (39) for φ n+1 h ,φ n+1 h provides the ingredients necessary to update the particle field ψ p in the subsequent mesh-particle projection step.

Mesh-particle projection
Finally, a fully-discrete mesh-particle projection is formulated by mapping the increments of the mesh related field (φ h ) to the particles using where 1/2 ≤ θ L ≤ 1, andφ n h ∈ W h is defined according to Eq. (35). We emphasize the consistency of the formulation for the mesh-particle update (Eq. (40)) with the definition of ψ * ,n h ∈ W h in Eq. (34), which is required to respect the compatibility condition, Eq. (29), in the fully-discrete setting. Furthermore, it readily follows that in the advective limit (i.e. for κ = 0) it holds that ψ n+1 p = ψ n p , sinceφ n h =φ n+1 h = 0.

A fully-discrete formulation for the incompressible Navier-Stokes equations
We now briefly present the fully-discrete formulations for the four model components constituting the particle-mesh operator splitting method for the incompressible Navier-Stokes equations. The formulations are quite similar to those for the advection-diffusion problem but, in addition, involve the incompressibility constraint and the discretization of the advective fields U h ∈ W h andŪ h ∈W h,g .

Particle advection
Once the advective velocity field U h ∈ W h has been chosen, the particle advection for the incompressible Navier-Stokes equations proceeds along similar lines as for the linear advection-diffusion equation. To this end, a judicious choice for the advective field U h (x, t) in Eq. (25) for t ∈ I n , where I n = (t n , t n+1 ], is to use the mesh related velocity field u h ∈ W h frozen at the old time level t n , i.e.
The H(div)-conformity of the velocity field u h ∈ W h following the solution of the Stokes problem Eq. (27) guarantees that the advective field U h is pointwise divergence-free (see

PDE-constrained particle-mesh projection
The fully-discrete vector-valued PDE-constrained particle-mesh projection is formulated similarly to its in which, analogously to Eq.
where the acceleration fieldu n ∈ W h (i.e., the vector-valued counterpart of Eq. (35) which de facto linearizes the PDE-constrained particle-mesh projection problem.
Given these definitions, we now present a fully-discrete implementation of the optimality system Eq. (26).
The time-discrete counterpart of the co-state equation (Eq. (26a)) reads as follows: given the vector-valued particle field v n p ∈ V t , the particle positions x n+1 p ∈ X t , and the field v * ,n Simultaneously the fully-discrete state equation (Eq. (26b)) is to be satisfied, given the advective Neumann boundary condition h a , together with the fully-discrete optimality condition (Eq. (26c)), which reads, gives the reconstructed field v n+1 h , which will serve as an input to the Stokes solver (to be discussed next). We emphasize that the fully-discrete state equation (45b) does not involve upwinding of the (interface) fluxes.

Stokes problem
The unsteady Stokes equations are discretized in time by using the backward Euler method, similarly to the method used in [10]. The fully-discrete problem then reads as follows: given the intermediate field , and (local and global mass balances), are satisfied. In these equations, the diffusive flux σ n+1 d,h is given by and the diffusive flux at cell facets,σ n+1 d,h , is defined aŝ For a detailed analysis of the unsteady Stokes solver, we refer to [23,20]. Here we note that the solver is inf − sup stable, conserves momentum globally and locally, is energy stable, produces H(div)-conforming velocity fields u n+1 h ∈ W h , and last but not least, is amenable to a static condensation procedure.

Mesh-particle projection
For the incompressible Navier-Stokes equations, the mesh-particle projection is the vector-valued counterpart of Eq. (40), i.e. the momentum field on particles is updated according to v n+1 requiring differentiation in time. Hence, the major advantage of choosing l = 0 is that it avoids time stepping dependencies in the PDE-constrained particle-mesh projection, thus rendering the particle-mesh projection independent of θ.
Exploiting this evident advantage, we will restrict the discussion in the sequel to the choice l = 0. In fact, it is only for this particular choice of the Lagrange multiplier space that second-order accuracy in time can be expected for the given time integration scheme. Furthermore, independent of the choice for l, we expect optimal spatial convergence rates of the scheme of order k + 1, provided a sufficiently accurate particle advection scheme is used.

Accuracy of the time stepping scheme
For an overview, let us briefly summarize the essential ingredients required for obtaining second-order accuracy in time. In [10] -with the particle-mesh projection being formulated in terms of local 2 -projectionsit was shown that overall second-order accuracy in time can be expected when combining a backward Euler method for the Eulerian step, Eq. (39) or Eq. (46), with a second-order mesh-particle update as in Eq. (40) or Eq. (49), using θ L = 1/2. Note that the latter requires saving an earlier particle field. Backward Euler suffices on the Eulerian step since φ n+1 h and u n+1 h are advanced over one time step from the second-order fields φ * ,n h and u * ,n h , respectively. Thus, the local, second-order accuracy of backward Euler pertains. To extend this approach to the PDE-constrained particle-mesh projection pursued here, we have to use l = 0 as argued above. In addition, we have to ensure compatibility between the mesh-particle update, Eq. (40) or Eq. (49), and the PDE-constrained particle-mesh interaction in the fully-discrete setting, Eq. (38) or Eq. (45). As further detailed in Appendix A, the latter requirement leads to definitions for ψ * ,n h or v * ,n h via Eq. (34) or Eq. (43), respectively.

Numerical implementation
Similarly to the fully-discrete diffusion equation and the fully-discrete Stokes equations, the fully-discrete PDE-constrained particle-mesh projections can be implemented efficiently by eliminating the unknowns local

Advection-diffusion: Gaussian hump
The accuracy of the presented method is assessed by considering a (rotating) Gaussian pulse in the diffusive limit, for moderate diffusion, and in the advective limit. The domain of interest is the circular disk Ω := {(x, y) | x 2 + y 2 ≤ 0.5} and the velocity field is either set to a = 0 (diffusive limit), or given by The corresponding analytical solution for a rotating Gaussian pulse is given by in which x c = (x c , y c ) is the position vector of the center, σ is the initial standard deviation, and κ is a constant diffusivity. Furthermore, · 2 denotes the square of the Euclidean norm andx(x, t) is a coordinate in the rotating frame of reference, given by (x cos(πt) + y sin(πt), −x sin(πt) + y cos(πt)) . The initial condition φ(x, 0) is deduced from Eq. (51), with the initial standard deviation σ set to 0.1 and the Gaussian pulse initially being centered at (x c , y c ) = (−0.15, 0). The disk-shaped domain Ω is triangulated using a sequence of mesh refinements, and particles are randomly seeded in Ω such that each cell contains on average 30 particles, initially.
Four different cases (listed in Table 1) are considered for three values of the diffusivity, κ = 0.01, κ = 0.001, and κ = 0 (advective limit). The Dirichlet boundary condition for the diffusion step is deduced from the analytical solution Eq. (51), and we emphasize that in the advective limit the particle specific mass ψ p stays constant by virtue of Eq. (40). The same meshes, initial conditions, and particle distributions are used in all four cases.  Results obtained after a full revolution are presented in Table 2 for the different test cases. For Cases 1-3, at least second-order convergence is obtained. More precisely, for the largest value of the diffusivity (i.e. κ = 0.01), the convergence rate tends to second-order, whereas near-optimal convergence is obtained for moderate diffusion, with diffusivity κ = 0.001, resulting in third-order convergence for Case 1 and Case 3, and second-order convergence for Case 2. In the advective limit, optimal convergence rates are obtained for Case 2 and Case 3.  Mesh

Advection: rigid body rotation
In order to qualitatively assess the scheme's ability to preserve point and line singularities, we next consider the advection test proposed in [28]. This test comprises the rigid body rotation of a pointy cone (initially centered at (x, y) = (−0.3, 0)), a slotted disk (initially centered at (x, y) =(0, -0.3)), and a Gaussian hump (initially centered at (x, y) = (0.15, 0.15)) on the the circular disk Ω := {(x, y) |x 2 + y 2 ≤ 0.5}. The velocity field is given by Eq. (50). By virtue of Eq. (40) it follows that in the advective limit the particle specific masses need not be updated, so that any discontinuities at the particle level persist. A mesh with moderate spatial resolution (Mesh 1, containing 16189 cells), and a fine resolution mesh (Mesh 2, containing 64561 cells) are considered. Approximately 30 particles are assigned per cell initially, resulting in a total number of slightly over 5 × 10 5 and 2 × 10 6 particles in total for Mesh 1 and Mesh 2, respectively. Time step sizes are chosen such to keep the CFL-number approximately the same on both meshes, see Table 3.
The results obtained for the two different configurations (see Table 3) are assessed visually after a half and a full rotation in Fig. 1. Since the particle values are not updated, the initial discontinuities are maintained at the particle level. As such, the test provides insight into how well the mesh-based solution follows the particle field. The shapes of the pointy cone and the Gaussian hump are well-preserved (without numerical damping), both for a half rotation as well as for a full rotation. Although the shape is well-preserved for the slotted-disk, localized overshoot is observed near the discontinuities for a half rotation. Rather than being a dispersion artifact, the over-and undershoot should be interpreted as a resolution issue with the mesh resolution being too coarse to capture the sharp discontinuity at the particle level monotonically. This is clearly illustrated by Fig. 1c and Fig. 1f, showing that the initial condition is accurately recovered for a full rotation at t = 2, thus indicating that the mesh-based solution is able to follow the moving particle field.
We also investigate the mass conservation errors for the two different configurations. To this end, a measure for the relative global mass conservation error at time T is defined as in which φ h (x, 0) and φ h (x, T ) are the mesh related fields at time 0 and time T , respectively.
The local mass conservation error is investigated via the L 2 -norm of the time-discrete counterpart of the local conservation statement, Eq. (32). For the problem under consideration, this local mass conservation error norm at time level n + 1 is given by The mass conservation errors as defined by Eqs. (52) and (53) after a half rotation (T = 1) and a full rotation (T = 2) are tabulated in Table 4. The results confirm global and local mass conservation of the PDE-constrained projection to machine precision, irrespective of the resolution of the Eulerian mesh.   accordingly. Similar to the preceding rigid body rotation benchmark, the particle specific masses need not be updated, so that any discontinuities at the particle level persist without artificial diffusion. Hence, this test is well-suited for assessing the behavior of the scheme in the presence of steep gradients. Furthermore, in the particle-mesh setting of this study, the test case can be used to assess the global conservation statement in the presence of inflow and outflow boundary conditions. The fields φ h at t = 2.0 are plotted in Fig. 2. For a characteristic direction of 45 • , the discontinuity is captured exactly at the mesh level. For the other characteristic directions, an overshoot is observed near the discontinuity. However, this overshoot remains strictly localized to one mesh cell upstream of the discontinuity. This behavior can be expected since no attempts are presently made to preserve monotonicity at the mesh, while the discontinuity at the particle level is inherently maintained without any diffusion. We leave the introduction of limiters as a fruitful area for future research, and refer to the work of Bochev and coworkers on bound-preserving remaps as an interesting starting point [29].
The global mass conservation property of the scheme is verified by virtue of the time-discrete equivalent of Eq. (31). The global mass conservation error, denoted by ∆φΩ , is the residual after subtracting the right-hand side from the left-hand side in this equation. The values for ∆φΩ thus obtained at t = 2 are tabulated in Table 5 for the four characteristic angles, confirming that global mass conservation is satisfied to machine precision.

Navier-Stokes equations: Poiseuille flow
Plane Poiseuille flow is considered to assess the behavior of the particle-mesh operator splitting scheme

Convergence study
A coarse, regular triangular mesh is constructed using 8 × 4 × 2 = 64 cells. The mesh size is halved three times, giving a sequence of mesh refinements with the finest mesh containing 4096 cells. On average, 30 particles are initially assigned per cell. This number is kept constant upon mesh refinement. Throughout, we will use a θ L -value of 1/2 for the particle updating, and the time step size corresponds to a CFL-number of approximately 0.25. Results at the dimensionless time instant t * = tU/2d = 100 are presented in Table 6 for polynomial orders (k, l) = (1, 0) and (k, l) = (2, 0), respectively. The observed convergence is near-optimal (i.e. order k + 1 in the velocity and order k in the pressure).

Assessing the particle resolution
Evidently, it is desirable from an efficiency perspective to keep the number of particles as low as possible without compromising accuracy. We therefore investigate the influence of the particle resolution on the accuracy by considering the Poiseuille flow benchmark, using l = 0 and k = 1, 2, 3, 4, combined with a variable particle resolution. We restrict the discussion to the mesh containing 256 cells and a time step of ∆t = 0.1 and only vary the particle resolution for this configuration so that the number of particles per cell (denoted withS K 0 ) is in the range 2-50, initially. In order to have full control over the initial particle configuration, particles are placed on a regular lattice. As expected, many of the low particle resolution tests fail prematurely (i.e. before reaching the end time t * = tU/2d = 100). In Fig. 3, results are visualized by plotting the L 2 -error of the velocity field against the (average) number of particles per cell for the different polynomial orders. For convenience, we assign an error-value of 1 to all low particle resolution runs failing prematurely.
Evident from the figure is the sharp transition in the error levels. That is, for low particle resolutions the computations are prone to fail prematurely, whereas from a certain particle resolution threshold onward, accurate results are obtained with error values being independent of the particle resolution. Interestingly, this threshold particle resolution depends on the polynomial order of the basis functions and is approximately equal toS K 0 = 8, 8, 10, 15 for k = 1, 2, 3, 4, respectively. By recognizing that the particle-mesh projection is based on the local least-squares minimization problem (Eq. (18)), we note that unisolvency of the particle locations with respect to the discontinuous function space W h is a necessary condition for the particle-mesh projection to be accurate. Empirically, we thus observe this condition to be satisfied if the minimal number of particles in a cell is at least equal to the number of local basis functions, except for linear basis functions (k = 1) and to a lesser extent for quadratic basis functions (k = 2), which require a higher number of particles per cell.

Navier-Stokes equations: Taylor-Green flow
The Taylor This problem features a periodic sequence of decaying vortices. Provided the Reynolds number is sufficiently small, closed analytical expressions for the velocity and the pressure are given by, respectively, u(x, t) = U exp −2νπ 2 t (− cos(k x x) sin(k y y) , sin(k x x) cos(k y y)) , p(x, t) = 1 4 exp −4νπ 2 t (cos(2k x x) + cos(2k y y)) , in which U is the initial velocity amplitude, ν is the kinematic viscosity, and k x = 2π/L x and k y = 2π/L y are wave numbers in the x-and y-direction, respectively, with L x and L y being the associated wave lengths.
Since only an initial condition needs to be specified, this test is particularly suited to study the accuracy of the time stepping scheme. Furthermore, due to the absence of body forces and boundary tractions the test is useful to investigate momentum conservation. To this end, a mesh-related measure for the global momentum conservation error is defined as

Convergence study
To investigate the convergence properties of the scheme, we consider a time interval of interest I = (0, 2]. The initial peak velocity U in Eq. (54) is set to 1, the wave lengths L x = L y = 2, and the kinematic viscosity is either ν = 2 · 10 −2 or ν = 2 · 10 −3 , resulting in Reynolds numbers of Re = U L/ν = 100 and Re = 1000, respectively. The chosen time step corresponds to a CFL-number of approximately 0.4. Different function space combinations are considered, see Table 7. For comparison, Case 3 uses the (non-conservative) 2 particle-mesh interaction from [10].  Velocity and pressure errors are tabulated in Table 8 for different model runs. In all these runs, the average number of particles per cell is 28. For Case 1, consistent second-order convergence is observed for the velocity and first-order convergence is obtained for the pressure. Given the function spaces (piecewise linear velocity, and piecewise constant pressure), no better convergence rates would be expected.
For Case 2, we obtain approximately second-order convergence in the velocity and the pressure, both for the Re = 100 test case and the Re = 1000 test case. This indicates that the time stepping error becomes dominant over the spatial error, where the former is expected to converge with second-order. As expected, momentum is conserved globally up to machine precision for the PDE-constrained particle-mesh interaction (Case 1, Case 2), whereas this is clearly not so for the unconstrained projection method (Case 3), using local 2 -projections. Noteworthy to mention is that the errors in the velocity and pressure, respectively, are almost identical for Case 2 and Case 3.

High Reynolds number test
Finally, we show that the developed scheme is robust for high Reynolds numbers. To this end, we investigate an 8 × 8 array of Taylor-Green vortices (i.e. using L x = L y = 0.5). The initial velocity amplitude is U = 1 and the kinematic viscosity is chosen such that the Reynolds number is Re = U L/ν = 2.5 · 10 5 . A 200 × 200 × 2 regular triangular mesh is used, the polynomial order are set to k = 2 and l = 0, and 2.25 · 10 6 particles are seeded in the domain resulting in approximately 28 particles per cell, initially. Furthermore, the time step size is ∆t = 5 · 10 −3 , resulting in a CFL-number of approximately 1.  The vorticity field at various time instants is plotted in Fig. 4. For this high Reynolds number, the regular initial lattice of vortices collapses into a chaotic, turbulent motion. Two features characteristic of 2D turbulence are qualitatively observed in Fig. 4. First, comparing the vorticity fields at t = 7.5 and t = 20, the emergence of large-scale coherent structures is clearly observed for longer simulation runtimes.
These large-scale vortices are only weakly dissipative and hence persist for many eddy-turnover times [31].
Second, the viscous dissipation takes place in the thin vortex filaments, clearly observed in Fig. 4. Somewhat remarkable are the small-scale, high-vorticity spots, best observed for the longer simulation runtimes (e.g. for t = 15 or t = 20). As reported in literature [32], such small-scale vortices can be generated by Kelvin-Helmholtz type instabilities of the vortex filaments for sufficiently high Reynolds numbers.
The energy density spectra for different time instants are plotted in Fig. 5. At time instant 0, the energy density is concentrated at one wave-number. At later time instances, an energy cascade develops with a slope of approximately k −3 covering almost two decades. This is in good correspondence with the k −3 slope for the direct enstrophy cascade as predicted by Kraichnan [33]. A slight steepening of the spectrum is observed at later time instances. Corroboration for this behavior is found in the literature, where the coherent vortices are reported to destroy scale invariance and can produce steeper spectra than predicted by the k −3 direct enstrophy cascade [31]. The increase in the energy density contained in the low wave numbers at the expense of the energy density contained in the higher wave number range, further illustrates the formation of large-scale coherent structures.

Conclusions and outlook
This paper presented a particle-mesh scheme which allows for diffusion-free advection, satisfies mass and momentum conservation properties in a local (i.e. cellwise) sense, and allows the extension to high-order spatial accuracy. Central in obtaining this, is to formulate the particle-mesh projections in terms of a PDE-constrained minimization problem. The key idea in formulating the constraint is that from a meshperspective the particle motion must satisfy an advection operator. By expressing the control variable in terms of single-valued functions at cell interfaces, it was shown that the HDG method naturally provides the necessary ingredients for formulating the optimality system. Consistency and conservation of the constrained interaction were proven. An analysis of the resulting discrete optimality system, revealed that a specific choice for the Lagrange multiplier field renders the PDE-constrained optimization procedure independent of the time integration method, thus resulting in a particularly attractive and robust scheme.
The PDE-constrained particle-mesh interaction is embedded in a Lagrangian-Eulerian operator splitting approach for the advection-diffusion and the incompressible Navier-Stokes equations. Consistency requirements specific to particle-mesh methods were formulated, and a particular fully-discrete formulation was formulated to obey these requirements.
A range of numerical experiments unveils the potential of the presented approach. Standard convergence rates in space and second-order convergence rates in time were confirmed for a range of examples for the linear advection-diffusion problem and the incompressible Navier-Stokes problem. The method turns out to be promising, albeit monotonicity is not guaranteed. A high Reynolds number test case gives further evidence of the robustness of the scheme.
The presented method opens many avenues for further investigation. Future work can include a more rigorous mathematical analysis of the method, finding alternative (high-order time-accurate) operator splitting approaches and/or time stepping strategies, and the reconstruction of monotonicity preserving background fields from scattered particle data by, e.g., exploiting the subgrid information available at the particle level.
Future work will also investigate the applicability of the presented projection method for (mass conservative) interface tracking in, e.g., multiphase flows.