A flux-interpolated advection scheme for fluid simulation

(1)

https://doi.org/10.1007/s00371-021-02187-2 O R I G I N A L A R T I C L E

A flux-interpolated advection scheme for fluid simulation

Naoyuki Hirasawa¹ ·Takashi Kanai¹·Ryoichi Ando²

Accepted: 31 May 2021 / Published online: 10 June 2021

Abstract

We propose a new advection scheme for fluid simulation that improves both conservation and numerical diffusion. Our work differs from previous works in that we re-formulate interpolation as cell-face flux of a vector field instantly constructed from a scalar field, rather than a per-point evaluation at back-traced positions. Our novel interpolation method enables excellent preservation of conservative quantities since the sum of flux exactly counteracts on cell faces, which eventually evaluates the boundary-flux of the whole domain. Our method can be implemented as a plug-and-play extension (or a temporary scratchpad) to the conventional semi-Lagrangian scheme; hence, our method naturally inherits all the benefits of semi-Lagrangian schemes and can be seamlessly integrated with existing fluid simulation pipelines together with other (black-boxed) solver components.

We conducted numerical experiments to verify the accuracy of our scheme (conservation and the improved numerical diffusion) and compared its qualitative results with state-of-the-art advection schemes that are in heavy use in the production environment such as MacCormack and the WENO6 interpolation.

Keywords Computer graphics·Fluid simulation·Conservative advection·Predictive diffusion correction·Flux interpolation

1 Introduction

Conservation laws are the fundamental principles for a broad range of physics phenomena. In Newtonian mechanics, total energy, mass and momentum preservation are some represen- tative examples. Conservation laws also play an important role for physically plausible animations [4]. For example, violation of the conservation laws causes numerous objec- tionable visual artifacts. These include simulation instability due to a persistent increase in total energy [8,12,22], ghost force introduced by the failure to preserve linear momentum [7,16], and disappearance of substances resulting from the loss of mass [32,62]. Interested readers should refer to the work of Söderström et al. [51] for the conservative version of the Navier–Stokes equations.

B

Naoyuki Hirasawa

hirasawa@g.ecc.u-tokyo.ac.jp Takashi Kanai

kanait@acm.org Ryoichi Ando rand@nii.ac.jp

1 Graduate School of Arts and Sciences, The University of Tokyo, Meguro-ku, Tokyo 153-8902, Japan

2 National Institute of Informatics, Chiyoda-ku, Tokyo 101-8430, Japan

In the production environment, semi-Lagrangian schemes and their extensions have been established as bread-and- butter tools for solving the advection term in the Navier–

Stokes equations (in conjunction with Chorin’s projection method [10]) [3,5,35,50] as they are easy to implement and unconditionally stable for large Courant–Friedrichs–

Lewy (CFL) numbers. However, as pointed out by several researchers, semi-Lagrangian schemes are not conservative by nature [31–33,55], and as such, cause perceivable fluctu- ations of the total mass and momentum.

In the past decade, previous works enforced conservation properties by applying the conservation form of partial differential equations (PDE) [13,38], the method of characteristics and its variants [18,42,45,59], corrective schemes [17,23,54, 65], energy-preserving integrators [12,34,37] or by simply increasing the order of accuracy [19,26,47]. Nevertheless, (reasonably) conservative advection methods have not yet been put into practical use for graphics although their importance is widely recognized by researchers [34]. We believe that this is due to the limited number of tools available and unforeseeable cost-effectiveness over off-the-shelf methods.

In this work, we first revisit a long-standing problem of non-conservative semi-Lagrangian methods and show that our highly mass conserving (not exactly as we discuss later) animations compare favorably to these schemes. Next, we

(2)

Fig. 1 From left to right: our method with trilinear interpolation (0.2% mass loss), semi-Lagrangian method with WENO6 interpolation (24.1%

mass loss), MacCormack method with WENO6 interpolation (700% mass gain). Grid resolution is 64×100×64 for velocity, 128×200×128 for density

introduce a novelflux-basedinterpolation scheme and exam- ine its numerical properties through several comparisons and examinations which include Taylor–Green vortex [58] and a 2D disk rotation. Finally, we run complex simulations on relatively large scales (e.g., up to 400×200×200) to validate the feasibility of our scheme.

1.1 Our contributions

The aim of this work is to develop a new physically consistent, practical, and robust method for handling advection such that a transported quantity is best preserved. To this end, we propose and develop the following novel schemes.

– A flux-based interpolation scheme for semi-Lagrangian advection that highly conserves the total mass (and momentum).

– A prediction-correction scheme to minimize numerical diffusion when the trilinear interpolation is employed.

– A local volume compensation scheme to counteract erro- neous compression or expansion of cells deformed by incompressible flow.

In addition to the above schemes, we provide self-contained (no dependencies on any external libraries) snippet C++

codes (MIT-licensed) to facilitate reproducibility for parts where the implementation may not be straightforward.

Note that currently our proposed scheme cannot be used to preserve total energy; the energy itself can also be inde- pendently advected via a conservative form [31], but it

usually does not exactly agree with the kinetic energy of the advected velocity (in a discrete sense) or the Kol- mogorov spectrum [43]. This discrepancy issue needs to be specially handled (e.g., a trapezoidal temporal integra- tor [37], advection-reflection solvers [39,63], an explicit correction [32]).

2 Previous work

Numerical algorithms for advection have a long history both in engineering and graphics. In engineering, a family of finite difference methods (FDM), which delivers both fast fixed- patterned memory access and required accuracy, is mainly used [14].

On the other hand, in graphics, semi-Lagrangian schemes are relatively widely adopted as they potentially allow for larger time steps while retaining stability [6,27]. Such advantages enable artists to perform multiple trial-and-error attempts within a limited amount of time, which is highly desirable for film production purposes. However, this also compromises accuracy due to excessive numerical diffusion or a lack of conservation. A number of extensions have therefore been proposed to address these shortcomings as reviewed below. For brevity, we only focus on works related to graphics or semi-Lagrangian schemes.

Conservative forms A number of conservative forms have been proposed for semi-Lagrangian schemes in computa- tional fluid dynamics (CFD) literature. For example, Iske and

(3)

Käser [21] proposed a conservative semi-Lagrangian scheme for Voronoi cells, where conservation is achieved by subdividing advected cells with overlapping background cells. A similar scheme can also be used to achieve our goal, but in practice, the generation of cut-cells on deformed hexahedral meshes in 3D space requires dedicated implementation [57]

and can be computationally expensive. Meanwhile, Qiu and Shu [41] instead focused on the flux of a derivative of a quantity and integrated it along with the moving position to compute the temporal evolution of a variable within a time step. Such an approach conserves mass by construction, but requires an accurate integration of flux along a long trajec- tory path for large time steps. Note that our proposed scheme also uses flux, but we emphasize that we use the flux of aspa- tially integrated variableof interest (not a spatial derivative of a quantity itself). As a result, our scheme is still considered aninterpolation-based semi-Lagrangian scheme, which does not requireanyflux integration over back-trace trajectories.

Therefore, our method robustly handles cases with large time steps.

Method of characteristics The method of characteristics (MC) is a promising approach for reducing numerical diffusion and realizing long-term conservation of total energy, mass and momentum. It was first introduced by Stam [53]

in graphics and subsequently extended to long-term advection methods by several researchers [18,59]. The method has been further investigated [45] and improved for higher visual quality and accuracy [42]. The above MC schemes deliver excellent conservation over a long time, but conservation is slightly lost in map initialization. In this sense, MC schemes are not exactly conservative. On the other hand, our proposed scheme provides better conservation without such MC tracking (except for cases where artificial manipulations are made such as velocity extrapolation to the solid).

High-order methods Higher-order schemes are promising approaches for increasing both the spatial and temporal accuracies of advection schemes. In graphics, the MacCor- mack method [47], BFECC [26], QUICK [36], CIP [28] and WENO [49] methods have been investigated. However, as also mentioned by Qu et al. [42], these higher-order schemes nevertheless accumulate error on conservative properties and eventually those that must be conserved can noticeably fade- out in the long run. This issue is demonstrated throughout the paper.

Lagrangian methods Lagrangian methods (e.g., particles, meshes and vortex sheets) are often used to handle physical quantities (e.g., liquid body [20,61], velocity [23,67], vorticity [48] and deformation gradient [24]). Unlike Eulerian methods, particle-based methods deliver excellent conservation since particles do not dissipate; however, an explicit

viscosity term may be included for stability [20]. It should be noted that dedicated particle re-sampling is required for uni- formly distributed particles [1,2,30,52,56]. Despite the fact that particles may be used to simulate smoke [17,44,46,52, 64,66], Eulerian methods are often preferred for smoke since grids are suitable for smoothly representing varying density and velocity fields.

Other methods For level-set based approaches, the volume of liquid body can be lost during advection [6]. Such loss may be fixed by a volume correction scheme [25] or switching to the volume of fluid (VOF) method [9,38]. Several works attempt to recover vorticity that is lost during advection [15, 29,40,65]. Alternatively, the vorticity may be fully preserved via the form of explicit filaments [60] or circulation [13]. In the work of Elcott et al. [13], the preservation of vorticity is realized similarly to our work in the sense that the per-cell circulation counteracts on cell boundaries when integrating circulations over the whole domain.

3 Method overview

Advection is written as the following time-dependent hyper- bolic equation

∂φ

∂t + ∇ ·(φu)=0, (1) where t,φ,u and∇ denote time, quantity to be advected, velocity, and spatial derivative ∇ = [∂/∂x∂/∂y∂/∂z], respectively. For incompressible flow (that is∇ ·u=0), (1) can be alternatively written as

∂φ

∂t +(u· ∇) φ=0, (2) or more simply,

Dφ

Dt =0, (3)

where D/(Dt is the material derivative: D/(Dt ≡ ∂/∂t + u·∇. Using the method of characteristics (and the first-order approximation of the temporal integration), (3) may be solved by

φ(x,t+Δt)=φ(x−Δtu,t). (4) In practice,φ(x−Δtu,t)is evaluated by interpolation. This is known as the semi-Lagrangian scheme [53]. For a continuous cases, (3) is conservative

Ωall

φ(x,t+Δt)dV =

Ωall

φ(x,t)dV, (5)

(4)

whereΩalldenotes the whole simulation domain and dV = d xd yd z. Unfortunately, for discrete cases, (4) is usually not conservative

i

φ(xi,t+Δt)=

i

φ(xi−Δtu,t), (6)

wherei,xi denote a cell index and set of all grid cell positions, respectively. As discussed by Lentine et al. [32], this issue arises from the lack of partition of unity of (transposed) interpolation. For example, suppose that interpolation (4) is expressed as

φ(xi,t+Δt)=

j

ωi jφ(xj,t) (7)

whereωi j denotes an interpolation coefficient. Notice that

jωi j =1 for anyi but

iωi j =1 is generallynottrue for any j. As a consequence, not all theφ(xi,t−Δt)is carried to the next time stept, leading to the loss ofφ. Such an issue can be fixed by normalizingωi j (together with partial forward advection of the missingφ) [32] using a fundamen- tally different approach explained in the next subsection.

3.1 Flux-based interpolation

Letφ(x)be a continuous scalar field ofφ. Subsequently, we defineΦi as a corresponding scalar value stored on a grid pointigiven by

Φi =

Ωi

φ(x)dV, (8)

whereΩidenotes a local support cell domain. In actual simulations, (8) may be scaled by a large constant number (e.g., 1/ΔV) to avoid extremely small numbers. In this setting, the material derivative ofΦi is given as

DΦi

Dt = D Dt

Ωi

φ(x)dV. (9)

Next, suppose that a vector variableΨ is provided such that φ(x)= ∇ ·Ψ(x). (10) Note that at this momentΨ(x)can be any vector satisfy- ing (10). Substituting (10) to (9) yields

D (Dt

Ωi

φ(x)dV = D Dt

Ωi

∇ ·Ψ(x)dV (11)

= D Dt

∂Ωi

n·Ψ(x)dA, (12)

where∂Ωidenotes the surface boundary of a cell (essentially faces) and dA infinitesimal face area. Approximating (12) using a quadrature rule gives

D Dt

∂Ωi

n·Ψ(x)dA= D Dt

k

wkni(xk)·Ψ(xk)ΔA,

(13) where wk, ni(xk), ΔA denote quadrature weight, surface normal with respect to celliat positionxk, and cell face area, respectively. Applying the semi-Lagrangian scheme to (13) gives

Φi(t+Δt)=

k

wkni(xk−uΔt,t)·Ψ(xk−uΔt,t)ΔA,

(14) whereni(xk−uΔt,t)denotes surface normal on the back- tracked face. Note that for two adjacent cellsi,jand position xk on a face sharing the two cells, it holds that ni(xk) =

−nj(xk). Therefore, (13) counteracts on cell boundaries regardless the accuracy of the quadrature rule when computing the sum ofΦi, leading to

i

Φi(t+Δt)

=

i

k

wkni(xk−uΔt,t)·Ψ(xk−uΔt,t)ΔA (15)

=

k∈∂Ωall

wkn^∗(xk−uΔt,t)·Ψ(xk−uΔt,t)ΔA (16)

=

i

Φi(t), (17)

where k ∈ ∂Ωall and n^∗ in (16) denote all the sets of quadrature points on the domain boundary and the associated (axis-aligned) normal vector, respectively. Step from (16) to (17) is clear if we assume thatu=0 (the no-slip condition) on the (rectangular) simulation domain boundary. For the case ofn·u=0 (the free-slip condition), (17) is true if (16) is accurately computed on the domain boundary. Therefore, our method exactly conservesφprovided that a correspond- ingΨ is found and the flux integration is exact on the domain boundary. In this paper, we particularly choose

Ψ(x,y,z)= _y

0

[0 φ(x,y,z) 0]^Tdy. (18) For ease of subsequent calculations, we defineψ(x)and re- defineφ(x)as

ψ(x,y,z)= _y

0 φ(x,y,z)dy, (19)

(5)

φ(x)= ∂

∂yψ(x,y,z). (20)

Note that with this definition,

Ω∇·ΨdV =

∂ΩψnydA whereny denotes the ycomponent of the normaln. Equa- tion (19) can be efficiently computed on grids. An example would be

ψ_i_,_j₊¹

2,k =ψ_i_,_j₋¹

2,k+φi,j,kΔx, (21)

ψ_i_,−¹

2,k =0, (22)

where index j±¹₂denotes a cell face on a xz-plane.

3.2 Algorithm overview

Algorithm 1Our Flux-Interpolated Advection 1: Convertφtoψusing (21)

2: Back-trace faces

3:fori=1,2, . . .in cellsdo 4: Φi=

Flux ofψusing (14) 5:end for

Our algorithm is outlined in Algorithm1. We first begin with convertingφtoψ(Line 1). This may be done by simply accumulatingΦi along theyaxis as in (21). Next, we back- trace all the cell faces (Line 2), by moving all the grid points backward in time. Finally, we loop over the cells and compute the flux of ψi on the six surrounding faces (Line 4). The remainder of this paper primarily focuses on how to complete each step in Algorithm1.

4 Computing the flux

The most important (precision sensitive) part of our scheme is how to adequately evaluate the flux ofψon faces as in (13).

Although the face flux eventually counteracts in the bulk of the domain, a single face-centered quadrature rule would fail since it produces severe ghost popping artifacts as shown in Fig.3.

To illustrate the source of this problem, we assume that the simulation domain is 2D. Suppose thatψ(x,y)is given as a function depending on only thexcoordinates:ψ(x,y)= f(x). In this scenario, the analytical solution to∂ψ(x,y)/∂y is zero everywhere. However, when a single-point quadrature rule is applied, the quadrature points on different faces do not lie on the same vertical line (Fig.2), which can cause ghost flux.

A feasible workaround is to subdivide the integration region on a face into two (or more) parts. Take for example an edge (face in 2D) connecting two vertices(x1,y1),

Fig. 2 Single-pointed estimation of flux. Suppose that an underlying scalar fieldψis given asψ(x,y)= f(x). In this specific setting, the fluxψ(x,y)n_y should be net zero over the surfaces of deformed cell.

However, when a single-pointed (area weighted) evaluation is used, this will not be generally satisfied

Fig. 3 Artifacts induced by a single-pointed evaluation (Fig. 2) per face. This issue can be minimized by the piece-wise flux computation followed by our face subdivision as we discuss in Sect.4.2

(x2,y2)and an integer a (typically, an integer grid coordinate). Ifa−1 <x1<aanda <x2<a+1, we compute the flux as

_x₂

x1

ψnydx= _a

x1

ψnydx+ _x₂

a

ψnydx. (23)

Assuming thatψis a piece-wise function, (23) can be evaluated at high precision using a multi-point quadrature rule.

For three dimensions, this involves a face subdivision with both yz and yx planes, which we will discuss in the next subsection.

4.1 Change of variables

Directly subdividing faces in three dimensions is not straightforward. We facilitate this task by changing variables.

Suppose a mapping

X:(ξ, η)→(x,y,z), (24)

where(ξ, η)denotes a unit 2D coordinate 0 ≤ξ ≤1,0 ≤ η ≤ 1 and(x,y,z)3D coordinate on a face. Therefore,X translates a unit coordinate to the corresponding 3D coordinate on a face. For convenience, we call this mapping as (x,y,z)=X(ξ, η)orx =X(ξ). For a face with four vertices (face of a regular cell), this is given by

(6)

X(ξ, η)=

⎡

⎢⎢

⎣

(1−ξ)(1−η) ξ(1−η) (1−ξηξ)η

⎤

⎥⎥

⎦

T⎡

⎢⎢

⎣ v1

v2

v3

v4

⎤

⎥⎥

⎦, (25)

where v1,v2,v3 and v4 are four face vertices enumerated counter-clockwise. With this setting, flux on a face is re- interpreted as

∂Ωface

ψnydA

= ₁

0

₁

0 ψ (X(ξ, η))ny(X(ξ, η))S(ξ, η)dξdη, (26) whereS(ξ, η)denotes a relative area scaling factor.

ny(X(ξ, η))S(ξ, η)is given by

ny(X(ξ, η))S(ξ, η)=

⎡

⎣0 1 0

⎤

⎦· ∂

∂ξX(ξ, η)

× ∂

∂ηX(ξ, η)

. (27)

The cross-product term in (27) is explicitly given by ∂

∂ξX(ξ, η)

× ∂

∂ηX(ξ, η)

=

⎡

⎢⎢

⎣

(1−ξ)(1−η) ξ(1−η)

ξη (1−ξ)η

⎤

⎥⎥

⎦

T⎡

⎢⎢

⎣

(v2−v1)×(v4−v1) (v2−v1)×(v3−v2) (v3−v4)×(v4−v1) (v3−v4)×(v3−v2)

⎤

⎥⎥

⎦. (28)

By substituting (28) to (27) and (27) to (26), then to get ψ (X(ξ, η)) by interpolation, we can conveniently inter- pret (26) as an integration in a unit(ξ, η)coordinate using only four vertex variables. Consequently, our flux computation could be simplified from 3D to a simpler 2D problem.

4.2 Face subdivision

Our next task is to subdivide faces in(ξ, η)space (extend (23) to two dimensions). Conceptually, this is trivially done by recursive 2D contouring. The core idea is that any non-self- intersected polygon can be subdivided into multiple parts using the devised Marching Squares algorithm. We outline our subdivision procedures in Algorithm2. Details of each step are discussed in the subsections.

4.2.1 Finding integer numbers

The subdivision of a face starts with a unit square polygon and the target coordinatex. For each vertex in the polygon,

Algorithm 2Recursive 2D Contouring 1:Subdivide(unit square,x)

2:functionSubdivide(polygon,i∈(x,z)) 3: x_i^min,x_i^max←Find Min/Max Coordinate §4.2.1 4: a_i^min← x^min_i

5: ifa_i^min=x_i^minthen 6: a_i^min←a^min_i +1 7: end if

8: ifa_i^min≤ x_i^maxthen

9: P←Contour(pol ygon,a_i^min) §4.2.2 10: foreachpol y_kinPdo

11: Subdivide(pol y_k,i) 12: end for

13: else if i=x then 14: Subdivide(pol ygon,z) 15: else

16: set polygon as fully subdivided 17: end if

18:end function

we first find a set of integer numbers with which we need to split (Line 3). We do this by first interpolating a pair of 3D coordinatesx,zon four vertices. Next, for each coordinate we find the minimal x^min and maximalx^max. Finally, the minimal integer is given by

a^min_x = x^min. (29)

Ifx^minis an exact integer, we seta^min_x =x^min+1 (Line 6).

4.2.2 Recursive contouring

Once the minimal a_x^min is found, we first check ifa_x^min ≤ x^max. In this case, we subdivide the polygon into multiple (in most cases, two) sub-polygons by devised contouring. If this is not the case, we switch the target coordinate fromx toz(Line 13) and recursively attempt the same subdivision procedure (Line 14). If the target coordinate is already set as z, we set the polygon as fully subdivided (Line 16).

Our contouring scheme is as follows. Letq1. . .qn ∈R²be sequential polygon vertices in(ξ, η)space ordered counter- clockwise withq1=qn. For each index pairkandk+1, we check if the signs ofx(qk)−a^min_x andx(qk+1)−a^min_x are different wherex(qk)denotes thexcoordinate interpolated at qk. The interpolation may be done by interpolatingv1. . .v4

(four vertices of the face) at its(ξ, η)coordinate. The different signs mean that there is a crossing iso-pointcbetween kandk+1. The crossing pointcon the edge can be found using the two signed distance values, similar to the Marching Squares algorithm. Finally, we connect pointscwith edges e1 = (c1,c2) ,e2 = (c3,c4) . . .en = (c2n−1,c2n). Since edgesendo not intersect one another, we can subdivide this polygon with these edges without complication. Although our strategy can support multiple cuts for onea^min_x ; except for some rare cases (e.g., face inversion), only one edge is

(7)

Algorithm 32D Contouring of A Polygon 1:functionContour(poly,a_i^min)

2: P← {},C← {}

3: fork=1,2, . . .in poly verticesdo 4: D₁←x_i(q_k)−a_i^min

5: D₂←x_i(q_k₊₁)−a_i^min 6: ifD₁D₂<0then 7: θ1←D₂/(D₂−D₁) 8: θ2← −D₁/(D₂−D₁) 9: c←θ1q_k+θ2q_k+1 10: addctoC 11: end if 12: end for

13: for(c_2n−1,c_2n)∈Cdo 14: e_n←create edge(c_2n−1,c_2n)

15: {poly,sub-poly}←cut poly withe_ninto two 16: add sub-poly toP

17: end for 18: add poly toP 19: returnP 20:end function

likely to be generated. We outline our contouring algorithm in Algorithm3.

4.2.3 Integrating the flux

After the above face subdivision, we proceed to compute the integration in (26). To this end, we use a single-pointed quadrature rule because it is easy to implement and fast to evaluate. However, a single quadrature point is exactly accurate only for either piece-wise constant or linear func- tions. In this paper, we use the trilinear interpolant (as we discuss later), which may introduces some error. If such an error is not tolerable, one may alternatively analytically inte- grate (26), which would be slightly more expensive than the single-pointed quadrature.

5 Extending to the trilinear interpolant

In the above exposition, we assume that analyticalφ(x)such thatΦi =

Ωiφ(x)dV is available at any timing. How- ever, in numerical simulation, we only have discrete variables φ1, φ2. . . φn(stored on grids), andφ(x)is defined as a function of f(x, φ1, φ2. . . φn). In most cases, this is given as interpolation:φ(x)=

iαi(x)φi whereαi(x)denotes the shape function with respect to interpolation. In this case,φ is not self-producible

Φi−

Ωi

φ(x)dV =Φi−

Ωi

j

αj(x)φj dV

=i, (30)

Fig. 4 Examples of subdivided polygons on a unit square. The flux of ψis computed per subdivided polygons in theξ, ηcoordinate

where i denotes an error induced by interpolation, which is generally not zero except for piece-wise constant interpolation. Therefore, the aforementioned algorithm works only for the piece-wise constant approximation ofφ(x); other- wise, our method introduces excessive numerical diffusion and loss of mass (momentum) (i.e.,

iφi =

φ(x)dV).

Note that our piece-wise constant approximation can still capture subtle sub-cell changes since single face flux computation may contain contributions from multiple overlapping cells. Therefore, with respect to the interpolation accuracy, our piece-wise approximation is conceptually similar to the trilinearly interpolated semi-Lagrangian method. However, in many cases, high-order interpolants are preferable because they produce superior visual fidelity [15,19]. In this paper, we extended our method to support the trilinear interpolant, which we find less diffusive than trilinarly (or even WENO6) interpolated semi-Lagrangian schemes. Higher-order interpolants are left for future work.

5.1 Predictive diffusion correction

Except where noted, we assume thatφ(x)is trilinear polynomial. To enforce self-reproducibility, we pre-computei

every time step using (30), and re-defineφ(x)as φ(x)=

i

αi(x)φi

+(x), (31)

where (x) denotes a piece-wise constant interpolatedi. A straightforward way to computei would be to define a functionφ^∗(x) =

iαi(x)φi and apply our flux scheme withu = 0 to get an intermediateφ_i^∗. Finally,i is given asi =φi −φ_i^∗. This can be done efficiently since all the face subdivisions and splits are pre-computable. Since (31) is conservative (i.e.,

iφi =

φ(xdV), the newly advected

iφi(t+Δt)will be conservative again.

5.2 Vertical integration

For a piece-wise approximated φ(x), ψ(x) is computed by (21). However, for trilinear polynomials, (21) is no longer valid. For this reason, we use a slightly different strategy for

(8)

trilinear interpolantsφ(x). First, note that for any arbitrary interpolant,

_R

−R

αi(x,y,z)dy=1, (32)

for a reasonably large constant R. Such R is known as a support radius for the interpolant kernel. Next, the objective function that we wish to solve is

ψ(x)= _y

0 φ(x,y,z)dy

= _y

0

i

αi(x,y,z)φi+(x,y,z)dy. (33)

Using Fubini’s theorem, (33) is re-written as _y

0

i

αi(x,y,z)φi +(x,y,z)dy= (34)

i

φi

y 0

αi(x,y,z)dy+ y

0

(x,y,z)dy. (35)

The second term of (35) can be computed by (21) since (x,y,z)is piece-wise constant. The first term in (35) needs to be treated carefully. If both yi > R and yi < y−R, _y

0 αi(x,y,z)dy=1. If not, _y

0

αi(x,y,z)dy= _y₂

y₁

αi(x,y,z)dy, (36) y1=max(0,yi −R), y2=min(y,yi+R). (37) Equation (36) can be calculated either symbolically or numerically using a quadrature rule according to the prac- titioners’ preference. Whenyis large, most of (36) are equal to one, and only a few of (36) needs dedicated calculation.

Therefore, the first term in (35) can be simplified as

i

φi

_y

0

αi(x,y,z)dy

=ψ_j∗+¹2 +

j^∗+W+2 i=j^∗−W+1

φi

_y

y_j∗αi(x,y,z)dy, (38)

j^∗= y/Δx, yj^∗=(j^∗+1

2)Δx, W =R/Δx (39) where W is half a stencil width for the interpolant and ψj^∗+1/2 is vertically accumulatedφi using (21). Note that in (38) we used an implicit rule thatφi =0 wheni <1 ori exceeds the maximal index in theycoordinate.

6 Cell incompressibility correction

For incompressible flow, cell deformation due to the advection should not change its local volume. Unfortunately, for discrete simulation where cell faces cannot be distorted (faces remain planar), the cell volume can change as a numerical error. As a result, our advection scheme can yield an erro- neous sink or source.

We correct this error by a cell incompressibility correction as follows. After the advection, we compute volume errorσi

per celli. Next, we solve the following linear system

j∈N

p_i^∗−p^∗_j

=σi, (40)

where p^∗and jdenotes artificial pressure and six neighbor cellsN =(i±1,j,k), (i,j±1,k), (i,j,k±1), respectively (we setp^∗_j = p_i^∗when jis outside the domain). Finally, we move substances across faces as

φi^corrected=

j∈N

p^∗_i −p^∗_j 1

2

φi+φj

. (41)

Solving (40) accurately is costly. Fortunately, we found that when an iterative solver (e.g., a pre-conditioned conjugate gradient) is used, the L² relative residual of around e.g., 10⁻¹is sufficient. This can be achieved by one or a few itera- tions using an algebraic multigrid pre-conditioned BiCGStab solver provided by AMGCL [11]. If possible, the same pre- conditioner could be re-used for the pressure solver later.

7 Results

We implemented our method by replacing the advection rou- tine in a conventional fluid solver; as such, the majority parts of solver components are inherited. In doing so, we chose not to extrapolate density towards solids as this changes the total mass but instead introduced an artificial potential force to push the density out of solids when the density is erro- neously carried to inside the solids.

Regarding the performance, we emphasize that our current implementation is designed as a concept of proof; and is not yet optimized for any reasonably fast runtime. For example, we have generated all the quadrature points in the first pass and used them in the second pass to interpolate variables. As a result, our implementation is significantly memory intensive.

We chose to do this because in this way we can provide more detailed timings for individual components of our method.

However, we believe that this bottleneck is easily alle- viated by the on-the-fly computations of face subdivision, quadrature point assembly, and the flux integrationperface in a massively parallel fashion with a single pass. We also note

(9)

Table 1 Timings break down for Fig.1

timing(ms) Computation ofψusing (19)

(piecewise-constant)

25.036 Computation ofψusing (19) (trilinear) 343.338 2D contouring of A Polygon

(piecewise-constant)

1433.398 2D contouring of A Polygon (trilinear) 2529.101 Computation ofΦi(14) (piecewise-constant) 133.900 Computation ofΦi(14) (trilinear) 284.656 Cell incompressibility correction (40)

and (41)

1019.864

Table 2 Mass changes in the last step. In the legend, Flux corresponds to our flux-interpolated method, MC MacCormack, SL semi-Lagrangian, MOC the method of characteristics [59] with 20 steps for the initialization interval, Lentine et al. [32]

Advection scheme Figure1 Figure5

Flux + Tri 0.998 0.942

Lentine + Tri 0.999 0.999

MC + Tri 25.251 2.249

MC + WENO6 8.007 1.958

SL + Tri 2.031 0.172

SL + WENO6 0.759 0.080

MOC N/A 0.0735

Tri indicates trilinear interpolation and WENO6 likewise

that since our method can be seen as an (extended) Marching Squares method on all the cell faces, it could be efficiently implemented on the (multiple) GPUs without complexity.

For this matter, in the next subsections will primarily focus on the qualitative results of our method. Timings breakdown

is provided (Table1) as a reference to indicate which task is more comparably expensive. The overview of the accuracy of mass conservation is listed in Table2.

7.1 Smoke rising past a sphere

Figure1shows a one-way coupled smoke simulation using three different methods. When a semi-Lagrangian method is used, the total mass tends to lose despite an accurate WENO6 scheme is used for interpolation. For this simulation, the total mass was lost 24% in the last video frame. MacCormack method on the other hand improves spatiotemporal accuracy of the semi-Lagrangian method, we found that the mass tends to persistently increase. Eventually, the mass gain was approximately 8 times. Our method, on the other hand, preserves the mass within the tolerance of 0.3%.

7.2 Two smoky spheres colliding

Figure 5shows two smoky spheres colliding at the center using five different methods. For the case of the semi- Lagrangian method, we observed that 92% of mass is lost in the course of the simulation. The method of characteristics better preserves mass, but due to the initialization at every 20 steps, the mass is eventually 92.6% lost. MacCormack similarly to Fig.1gain mass, and the mass was amplified to 1.95 times in this simulation. Lentine et al. [32] exactly preserves mass but introduces excessive numerical diffusion. It may be possible to extend the method of Lentine et al. [32] to use WENO6 interpolation instead, but the algorithms for tracking all the interpolation coefficients for the WENO6 interpolation might not be straightforward or efficient. Our method on the other hand preserves mass within the tolerance of 6%

retaining the sharp edges of smokes.

Fig. 5 Two spheres collision. From left to right: semi-Lagrangian method with WENO6 interpolation (92% mass loss), method of characteristics with WENO6 (92.6% mass loss), MacCormack method with

WENO6 (95.8% mass gain), the method of Lentine et al. [32] (no mass change), our method with trilinear interpolation (6% mass loss). Grid resolution is 120×120×192 for velocity, 240×240×384 for density

(10)

Fig. 6 2D disk rotation. Left: the method of Lentine et al. [32]. Right:

our method. Notice that our method retains sharper edges

Fig. 7 Taylor–Green Vortex experiment. Left: without our cell-wise volume correction. Right: our volume correction. Notice that a uni- formly distributed density fields remain (unnoticeably) the same for a long period of time

7.3 2D disk rotation

Figure6shows a 2D disk rotation to clarify the difference of numerical diffusion using the method of Lentine et al. [32]

and our scheme. Mainly due to our predictive diffusion correction scheme, our method better retains sharp swirly lines of density fields.

7.4 Taylor–Green vortex test

Figure7illustrates the importance of our cell-wise volume correction scheme using a Taylor–Green Vortex config- uration. Without the volume correction scheme, the cell deformation introduces sink and source and they accumulate over time. When our correction scheme is employed, the error is nearly invisible.

7.5 Additional results and discussions

We also have conducted momentum advection experiments, which are discussed in the supplemental materials.

Fig. 8 Total mass change of Fig.1

Fig. 9 Total mass change of Fig.5

8 Conclusions

In this paper, we presented a novel flux-interpolated advection scheme for fluid simulation. Our novel interpolation scheme provides both high conservation of advected quantities and improved numerical diffusion without resorting to high-order interpolation schemes such as WENO6. Our method can be implemented as an extended face-wise March- ing Squares algorithms per cell-face, and seamlessly inte- grates to the existing simulation pipelines. In the future work, we are interested in extending our method to volume-of-fluid method for tracking liquid interfaces while preserving the total volume.

Supplementary Information The online version contains supplementary material available athttps://doi.org/10.1007/s00371-021-02187- 2.

Acknowledgements This study is supported by JSPS Grant-in-Aid for Young Scientists (18K18060).

References

1. Adams, B., Pauly, M., Keiser, R., Guibas, L.J.: Adaptively sampled particle fluids. ACM Trans. Graph.26(3) (2007)

2. Ando, R., Thürey, N., Tsuruno, R.: Preserving fluid sheets with adaptively sampled anisotropic particles. IEEE Trans. Visual Com- put. Graphics18(8), 1202–1214 (2012)

3. Autodesk: Maya (2020)

(11)

4. Bargteil, A.W., Shinar, T.: An introduction to physics-based animation. In: ACM SIGGRAPH 2019 Courses, SIGGRAPH ’19.

Association for Computing Machinery, New York, NY, USA (2019)

5. Blender: Blender—a 3D modelling and rendering package. Blender Foundation, Stichting Blender Foundation, Amsterdam (2020) 6. Bridson, R.: Fluid Simulation for Computer Graphics, 2nd edn. AK

Peters/CRC Press, Boca Raton (2015)

7. Brochu, T.: Fluid animation with explicit surface meshes and boundary-only dynamics. Master’s thesis, The University Of British Columbia (2006)

8. Brochu, T., Batty, C., Bridson, R.: Matching fluid simulation ele- ments to surface geometry and topology. ACM Trans. Graph.29(4), 1–9 (2010)

9. Chentanez, N., Müller, M.: Mass-conserving eulerian liquid simulation. In: Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation, SCA ’12, pp. 245–254. Euro- graphics Association, Goslar, DEU (2012)

10. Chorin, A.J.: Numerical solution of the Navier–Stokes equations.

Math. Comput.22(104), 745–762 (1968)

11. Demidov, D.: Amgcl - a c++ library for efficient solution of large sparse linear systems. Softw. Impacts6,100037 (2020)

12. Dinev, D., Liu, T., Kavan, L.: Stabilizing integrators for real-time physics. ACM Trans. Graph.37(1), 1–19 (2018)

13. Elcott, S., Tong, Y., Kanso, E., Schröder, P., Desbrun, M.: Stable, circulation-preserving, simplicial fluids. ACM Trans. Graph.26(1), 4–es (2007)

14. Ewing, R.E., Wang, H.: A summary of numerical methods for time- dependent advection-dominated partial differential equations. J.

Comput. Appl. Math.128(1–2), 423–445 (2001)

15. Fedkiw, R., Stam, J., Jensen, H.W.: Visual simulation of smoke. In:

Proceedings of the 28th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH ’01, p. 15–22. Association for Computing Machinery, New York, NY, USA (2001)

16. Fei, Y.R., Maia, H.T., Batty, C., Zheng, C., Grinspun, E.: A multi- scale model for simulating liquid-hair interactions. ACM Trans.

Graph.36(4), 1–17 (2017)

17. Fu, C., Guo, Q., Gast, T., Jiang, C., Teran, J.: A polynomial particle- in-cell method. ACM Trans. Graph.36(6), 1–12 (2017)

18. Hachisuka, T.: Combined Lagrangian–Eulerian approach for accurate advection. In: ACM SIGGRAPH 2005 Posters, SIGGRAPH

’05, pp. 114–es. Association for Computing Machinery, New York, NY, USA (2005)

19. Heo, N., Ko, H.S.: Detail-preserving fully-Eulerian interface tracking framework. In: ACM SIGGRAPH Asia 2010 Papers, SIG- GRAPH ASIA ’10. Association for Computing Machinery, New York, NY, USA (2010)

20. Ihmsen, M., Orthmann, J., Solenthaler, B., Kolb, A., Teschner, M.:

SPH Fluids in computer graphics. In: S. Lefebvre, M. Spagnuolo (eds.) Eurographics 2014—State of the Art Reports, pp. 21–42.

The Eurographics Association (2014)

21. Iske, A., Käser, M.: Conservative semi-Lagrangian advection on adaptive unstructured meshes. Numer. Methods Partial Differ. Equ.

20(3), 388–411 (2004)

22. Jeschke, S., Wojtan, C.: Water wave animation via wavefront parameter interpolation. ACM Trans. Graph.34(3), 1–14 (2015) 23. Jiang, C., Schroeder, C., Selle, A., Teran, J., Stomakhin, A.: The

affine particle-in-cell method. ACM Trans. Graph.34(4), 1–10 (2015)

24. Jiang, C., Schroeder, C., Teran, J., Stomakhin, A., Selle, A.: The material point method for simulating continuum materials. In:

ACM SIGGRAPH 2016 Courses, SIGGRAPH ’16. Association for Computing Machinery, New York, NY, USA (2016)

25. Kim, B., Liu, Y., Llamas, I., Jiao, X., Rossignac, J.: Simulation of bubbles in foam with the volume control method. ACM Trans.

Graph.26(3), 98–es (2007)

26. Kim, B., Liu, Y., Llamas, I., Rossignac, J.: Flowfixer: Using bfecc for fluid simulation. In: Proceedings of the First Eurographics Con- ference on Natural Phenomena, NPH’05, pp. 51–56. Eurographics Association, Goslar, DEU (2005)

27. Kim, D.: Fluid Engine Development. AK Peters/CRC Press, Boca Raton (2016)

28. Kim, D., Song, O.Y., Ko, H.S.: A semi-Lagrangian cip fluid solver without dimensional splitting. Comput. Graphics Forum 27(2), 467–475 (2008)

29. Kim, T., Thürey, N., James, D., Gross, M.: Wavelet turbulence for fluid simulation. In: ACM SIGGRAPH 2008 Papers, SIGGRAPH

’08. Association for Computing Machinery, New York, NY, USA (2008)

30. Kugelstadt, T., Longva, A., Thuerey, N., Bender, J.: Implicit density projection for volume conserving liquids. IEEE Transactions on Visualization and Computer Graphics (2019)

31. Kwatra, N., Grétarsson, J.T., Fedkiw, R.: Practical animation of compressible flow for shock waves and related phenomena. In:

Proceedings of the 2010 ACM SIGGRAPH/Eurographics Sympo- sium on Computer Animation, SCA ’10. Eurographics Association, Goslar, DEU (2010)

32. Lentine, M., Aanjaneya, M., Fedkiw, R.: Mass and momentum conservation for fluid simulation. In: Proceedings of the 2011 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, SCA ’11, pp. 91–100. Association for Computing Machinery, New York, NY, USA (2011)

33. Li, G.S., Tricoche, X., Hansen, C.: Physically-based dye advection for flow visualization. Comput. Graphics Forum27(3), 727–734 (2008)

34. Li, W., Chen, Y., Desbrun, M., Zheng, C., Liu, X.: Fast and scal- able turbulent flow simulation with two-way coupling. ACM Trans.

Graph.39(4), 47:1–47:20 (2020) 35. MAXON: Cinema 4d (2020)

36. Molemaker, J., Cohen, J.M., Patel, S., Noh, J.: Low viscosity flow simulations for animation. In: Proceedings of the 2008 ACM SIG- GRAPH/Eurographics Symposium on Computer Animation, SCA

’08, pp. 9–18. Eurographics Association, Goslar, DEU (2008) 37. Mullen, P., Crane, K., Pavlov, D., Tong, Y., Desbrun, M.: Energy-

preserving integrators for fluid animation. ACM Trans. Graph.

28(3), 1–8 (2009)

38. Mullen, P., McKenzie, A., Tong, Y., Desbrun, M.: A variational approach to eulerian geometry processing. ACM Trans. Graph.

26(3), 66–es (2007)

39. Narain, R., Zehnder, J., Thomaszewski, B.: A second-order advection-reflection solver. Proc. ACM Comput. Graph. Interact.

Tech.2(2), 1–14 (2019)

40. Pfaff, T., Thuerey, N., Gross, M.: Lagrangian vortex sheets for animating fluids. ACM Trans. Graph.31(4), 1–8 (2012)

41. Qiu, J.M., Shu, C.W.: Conservative high order semi-Lagrangian finite difference weno methods for advection in incompressible flow. J. Comput. Phys.230(4), 863–889 (2011)

42. Qu, Z., Zhang, X., Gao, M., Jiang, C., Chen, B.: Efficient and conservative fluids using bidirectional mapping. ACM Trans. Graph.

38(4), 1–12 (2019)

43. Rasmussen, N., Nguyen, D.Q., Geiger, W., Fedkiw, R.: Smoke simulation for large scale phenomena. ACM Trans. Graph.22(3), 703–707 (2003)

44. Ren, B., Yan, X., Yang, T., Li, C.F., Lin, M.C., Hu, S.M.: Fast sph simulation for gaseous fluids. Vis. Comput.32(4), 523–534 (2016) 45. Sato, T., Batty, C., Igarashi, T., Ando, R.: Spatially adaptive long- term semi-Lagrangian method for accurate velocity advection.

Comput. Vis. Media4, 220–230 (2018)

46. Schechter, H., Bridson, R.: Evolving sub-grid turbulence for smoke animation. In: Proceedings of the 2008 ACM SIG- GRAPH/Eurographics Symposium on Computer Animation, SCA

’08. Eurographics Association, Goslar, DEU (2008)

(12)

47. Selle, A., Fedkiw, R., Kim, B., Liu, Y., Rossignac, J.: An unconditionally stable Maccormack method. J. Sci. Comput.35, 350–371 (2008)

48. Selle, A., Rasmussen, N., Fedkiw, R.: A vortex particle method for smoke, water and explosions. ACM Trans. Graph.24(3), 910–914 (2005)

49. Shu, C.W.: High Order ENO and WENO Schemes for Computa- tional Fluid Dynamics, pp. 439–582. Springer, Berlin (1999) 50. SideFX: Houdini (2020)

51. Söderström, A., Karlsson, M., Museth, K.: A pml-based nonreflec- tive boundary for free surface fluid animation. ACM Trans. Graph.

29(5), 1–17 (2010)

52. Solenthaler, B., Gross, M.: Two-scale particle simulation. ACM Trans. Graph.30(4), 1–8 (2011)

53. Stam, J.: Stable fluids. In: Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH ’99, pp. 121–128. ACM Press/Addison-Wesley Pub- lishing Co., USA (1999)

54. Su, J., Sheth, R., Fedkiw, R.: Energy conservation for the simulation of deformable bodies. IEEE Trans. Visual Comput. Graphics19(2), 189–200 (2013)

55. Sussman, M., Ohta, M.: A stable and efficient method for treating surface tension in incompressible two-phase flow. SIAM J. Sci.

Comput.31(4), 2447–2471 (2009)

56. Takahashi, T., Lin, M.C.: A geometrically consistent viscous fluid solver with two-way fluid-solid coupling. Comput. Graph. Forum 38(2), 49–58 (2019)

57. Tao, M., Batty, C., Fiume, E., Levin, D.I.W.: Mandoline: robust cut- cell generation for arbitrary triangle meshes. ACM Trans. Graph.

38(6), 1–17 (2019)

58. Taylor, G.I., Green, A.E.: Mechanism of the production of small eddies from large onesmechanism of the production of small eddies from large ones. Proc. R. Soc. Lond. Ser. A Math. Phys. Sci.

158(895), 499–521 (1937)

59. Tessendorf, J., Pelfrey, B.: The characteristic map for fast and efficient vfx fluid simulations. In: Proc. Computer Graphics Inter- national Workshop on VFX, Computer Animation, and Stereo Movies. Ottawa, Canada (2011)

60. Weißmann, S., Pinkall, U.: Filament-based smoke with vortex shedding and variational reconnection. ACM Trans. Graph.29(4), 1–12 (2010)

61. Wojtan, C., Müller-Fischer, M., Brochu, T.: Liquid simulation with mesh-based surface tracking. In: ACM SIGGRAPH 2011 Courses, SIGGRAPH ’11. Association for Computing Machinery, New York, NY, USA (2011)

62. Wojtan, C., Turk, G.: Fast viscoelastic behavior with thin features.

ACM Trans. Graph.27(3), 1–8 (2008)

63. Zehnder, J., Narain, R., Thomaszewski, B.: An advection-reflection solver for detail-preserving fluid simulation. ACM Trans. Graph.

37(4), 1–8 (2018)

64. Zhang, M., Liu, S., Sun, H., Si, W., Qian, Y.: Hybrid vortex model for efficiently simulating turbulent smoke. In: Proceedings of the 13th ACM SIGGRAPH International Conference on Virtual- Reality Continuum and Its Applications in Industry, VRCAI ’14, pp. 71–79. Association for Computing Machinery, New York, NY, USA (2014)

65. Zhang, X., Bridson, R., Greif, C.: Restoring the missing vorticity in advection-projection fluid solvers. ACM Trans. Graph.34(4), 1–8 (2015)

66. Zhu, B., Yang, X., Fan, Y.: Creating and Preserving Vortical Details in SPH Fluid. Comput. Graph. Forum29(7), 2207–2214 (2010) 67. Zhu, Y., Bridson, R.: Animating sand as a fluid. ACM Trans. Graph.

24(3), 965–972 (2005)

Publisher’s Note Springer Nature remains neutral with regard to juris- dictional claims in published maps and institutional affiliations.

Naoyuki Hirasawa is a Master’s student at the Graduate School of Arts and Sciences, the University of Tokyo, Japan. He received BS degree from the Department of Physics, Faculty of Science and Technology, Keio University, Japan.

His research interests include physics-based animations and fluid simulations for computer graphics.

Takashi Kanai is an associate professor in the Graduate School of Arts and Sciences, the Univer- sity of Tokyo, Japan. His research interests include geometry processing and physics-based animation in computer graphics. He received his doctor degree of engineering from the University of Tokyo in 1998. He is a member of ACM, ACM SIGGRAPH, IPSJ (Infor- mation Processing Society of Japan), JSPE (Japan Society for Precision Engineering).

Ryoichi Ando is an Assistant Professor at National Institute of Informatics, Japan. Prior to that, he has been a postdoctoral resea- rcher at IST Austria. His research interests include physics-based animations for computer graphics with a strong emphasis on fluid. Many of my research animations can be seen athttps://ryichando.graphics