Fluid Cohomology

The vorticity-streamfunction formulation for incompressible inviscid fluids is the basis for many fluid simulation methods in computer graphics, including vortex methods, streamfunction solvers, spectral methods, and Monte Carlo methods. We point out that current setups in the vorticity-streamfunction formulation are insufficient at simulating fluids on general non-simply-connected domains. This issue is critical in practice, as obstacles, periodic boundaries, and nonzero genus can all make the fluid domain multiply connected. These scenarios introduce nontrivial cohomology components to the flow in the form of harmonic fields. The dynamics of these harmonic fields have been previously overlooked. In this paper, we derive the missing equations of motion for the fluid cohomology components. We elucidate the physical laws associated with the new equations, and show their importance in reproducing physically correct behaviors of fluid flows on domains with general topology.

Fig. 1.With torus-shaped obstacles in the fluid domain, the cohomology of the fluid velocity becomes nontrivial.By incorporating the correct time-evolution of the harmonic parts, our method (middle) removes the unphysical behaviors of a general vorticity-streamfunction method (left).We compare our result against ground truth which is produced by a velocity-based method using pressure projection (right).

Our method Method with fixed harmonic part
Velocity-based method The vorticity-streamfunction formulation for incompressible inviscid fluids is the basis for many fluid simulation methods in computer graphics, including vortex methods, streamfunction solvers, spectral methods, and Monte Carlo methods.We point out that current setups in the vorticity-streamfunction formulation are insufficient at simulating fluids on general non-simplyconnected domains.This issue is critical in practice, as obstacles, periodic boundaries, and nonzero genus can all make the fluid domain multiply connected.These scenarios introduce nontrivial cohomology components to the flow in the form of harmonic fields.The dynamics of these harmonic fields have been previously overlooked.In this paper, we derive the missing equations of motion for the fluid cohomology components.We elucidate the physical laws associated with the new equations, and show their importance

INTRODUCTION
One of the main tasks in fluid simulation in computer graphics is to approximate the solutions to the incompressible Euler equations which govern the dynamics of the velocity fields of incompressible inviscid fluids.One approach is to instead simulate the equations in the vorticity-streamfunction formulation of fluid dynamics.This formulation takes the vorticity field as the primary variable, whereas the velocity field is reconstructed from the vorticity through a streamfunction [Bridson 2015, §14.2;Chorin and Marsden 1990, §1.2].Numerical methods based on this formulation have more direct controls over the preservation of local vorticity.The velocity fields derived from the streamfunction are automatically divergence Fig. 2. Fluid simulation on minimal surfaces plays an important role in texturing soap films.Costa's surface is a minimal surface with genus one and three boundaries, and therefore nontrivial cohomology.Our method is the first vorticity-streamfunction method to correctly simulate fluid on generic surfaces including those with complicated topology.
free.Researchers have adopted this formulation to simulate vortex dynamics [Gamito et al. 1995;Elcott et al. 2007;Zhang et al. 2015] and two-phase flows [Ando et al. 2015a].Other usage of the vorticity-streamfunction formulation in fluid simulation includes [De Witt et al. 2012] (a vorticity-based spectral method) and [Rioux-Lavoie et al. 2022] (a Monte Carlo method).Vorticity-streamfunction formulation is especially advantageous for fluids on 2D domains including curved surfaces, as the vorticity equation requires only a scalar advection, as opposed to a vector advection in velocity-based solvers [Yaeger et al. 1986;Elcott et al. 2007;Azencot et al. 2014].
The use of streamfunctions is also the foundation for divergencefree flow syntheses [Bridson et al. 2007;Sato et al. 2014Sato et al. , 2021] ] and interpolations [Chang et al. 2019[Chang et al. , 2022]].While vorticity-streamfunction formulation is well established, the equations therein are insufficient at describing fluid flows on general domains, in particular, those that are non-simply-connected.Such scenarios with nontrivial topology are common, as multiplyconnectedness can arise from the presence of obstacles (Fig. 1), periodic boundary conditions, or the domain being a surface1 with nonzero genus (Fig. 2).The issue with multiply-connected domains is that they support irrotational flows, known as harmonic fields, that cannot be represented by the vorticity variable.In topological terms, the harmonic fields represent the 1st de Rham cohomology components of the flow. 2 A rarely discussed fact is that harmonic parts have their own dynamics.The vorticity equation alone is not enough to describe an Euler fluid flow.One must consider a coupled system between the vorticity and the harmonic components.
Previous Approaches.Previous work fills in the harmonic components using various treatments during the recovery of the velocity from the vorticity variable.One approach is to reconstruct the velocity field from the vorticity data using the Biot-Savart integral (which works in the entire Euclidean space with no obstacles) followed by adding a pressure gradient [Lin 1941;Bridson 2015, §14.3.7;Method with fixed harmonic part

Our method
Velocity-based method   Golas et al. 2012;Brochu et al. 2012;Vines et al. 2013;Zhang and Bridson 2014;Xiong et al. 2021], or mirror-reflected vortices in the obstacles of specific shapes (e.g. a plane [Angelidis and Neyret 2005]) to prevent the velocity from penetrating the obstacles.These methods do ensure correct dynamics of the harmonic components implicitly.However, they work only when the fluid domain is a subset of a Euclidean space and when the circulation around every handle of the obstacle is zero.In particular, it does not apply to flows on surfaces, on periodic domains, or with nonzero circulations around obstacles.
A few methods skip the vorticity variable and directly reconstruct streamfunctions from the velocity data on a grid.This technique is intended to obtain the exact divergence-free flow interpolation facilitated by the streamfunction [Biswas et al. 2016;Chang et al. 2022].Unlike integrating from vorticity, the construction is by pathintegrating velocity, which will not lose the harmonic information of the velocity.However, as we show in our paper (Proposition 5), global streamfunctions exist only for special cases (e.g. when the multiply-connectedness arises from removing obstacles from a simply-connected domain). 3nother approach is to encapsulate the degrees of freedom of harmonic fields into inhomogeneous boundary conditions for the streamfunctions [Mizukami 1983].Determining the boundary conditions for the streamfunction on each "island" is a well-known problem in oceanography [Godfrey 1989;Pedlosky et al. 1997].This is commonly known as the island rule, which boils down to solving a global Dirichlet-to-Neumann map problem.Unfortunately, the treatment only resolves harmonic fields that arise from obstacles in 2D.The analogous formulations in 3D are more sophisticated as they involve keeping track of circulations along flowing loops.In general, the approach cannot account for harmonic fields arising from periodic boundary conditions or nonzero genus of a surface domain.
Lastly, a commonly adopted method accounting for harmonic fields is to extract the harmonic component of the initial velocity and add it to the streamfunction-represented velocity at each future time step [Elcott et al. 2007, §4.5;Azencot et al. 2014, Eq. (1); Ando et al. 2015a, §3].This method keeps the harmonic parts fixed, as if they were constant in time.If the initial harmonic components are all zero, then the treatment is equivalent to neglecting harmonic components on closed surfaces [Boatto and Koiller 2015, §3], or setting a (tangentially) zero boundary condition on the streamfunctions [Gamito et al. 1995;Rioux-Lavoie et al. 2022, §4.1].While most of these approaches address the nontrivial cohomology generated by both the obstacles and the genus, they ignore the dynamics of the harmonic parts.While it might seem that Kelvin's circulation theorem would imply the conservation of harmonic parts over time, this is in fact a misinterpretation of the theorem and a misgeneralization of the preservation of harmonic fields when there is no vorticity [Thomson 1868].The absence of changes in harmonic parts can lead to unrealistic fluid behaviors.For example, omitting the harmonic dynamics will lead to a fixed total flux over every cross-section of the domain, e.g. the cross-section of the tunnel in Fig. 1.Notably, a jet flow shot through the tunnel will be obstructed and disrupted by an artificial zero total flux condition persisting since the initial condition.Similar anomalies are pointed out in [Rioux-Lavoie et al. 2022, pp 12] and an errata [Ando et al. 2015b].
New Formulation.In this paper, we derive the missing dynamics of the harmonic parts for the vorticity-streamfunction formulation of the Euler equations.The system applies to all multiply-connected fluid domains.The dynamics are presented in a simple form that is easy to incorporate into previous vorticity-streamfunction methods.The dynamics can be summarized as a new physical law of Euler fluids: For each cohomology component associated with a cross-section of the domain, the difference between the fluid cross-sectional flux and the total linking number (or winding number in 2D) between the harmonic field streamlines and the fluid vortex lines (or vortex points in 2D) is a constant of motion (see Fig. 4).
In other words, the harmonic components, which are all directly related to the cross-sectional fluxes, evolve precisely at the rate at which a vortex line cuts through a harmonic streamline.Moreover, this seemingly intricate rate of topological changes boils down to a simple integral formula ∫  (u × w) • h  x in terms of the computationally available velocity u, vorticity w and harmonic basis h  , which makes incorporating the new dynamics light-weight and highly practical.This is the same as the time-evolution of coefficients of a vector spectral basis, which includes the harmonic basis, as proposed in the velocity-based model reduction work by [Liu et al. 2015].We point out that this equation for the harmonic basis is necessary to account for the missing dynamics of the harmonic components in the vorticity-streamfunction formulation, making the formulation suitable for fluid simulations on all general domains.We also describe our new conservation law from the perspective of Hamiltonian fluid mechanics as a Casimir invariant.Note that the only previously known nontrivial example of a Casimir invariant for 3D incompressible fluids is the helicity [Khesin et al. 2022], which has played important roles in meteorology, topological fluid mechanics, and plasma physics.The discovery of the new Casimir invariant will facilitate more exciting studies in these areas.
In practice, incorporating the new evolution equation of the harmonic components is quite simple.Express the fluid velocity field u = curl −1 w +    h  in terms of the vorticity field w and an orthonormal basis (h  )  for the harmonic fields (see Fig. 5).On top of a traditional vorticity equation solver that updates w, one only needs to add an extra step for updating the coefficients   = ∫  (u × w) • h  x (see Fig. 6).Through numerical examples, we show the importance of the dynamics of the harmonic components for reproducing physically correct behaviors of fluid flows on surfaces and on 3D domains with nontrivial topology (Fig. 3).

Contributions. Highlights of this paper include
• Deriving the missing dynamics of the harmonic parts for incompressible inviscid fluids, and a simple and practical method for incorporating it to both 2D and 3D simulations (Section 4).• A new conservation law between harmonic streamlines and vortex lines in an Euler fluid (Theorem 2).• The first vorticity-streamfunction-based fluid solver on general surfaces that is consistent with the Euler equations.(Alg. 2 (Fig. 2)).• A new Hamiltonian formulation for incompressible Euler equations featuring new Casimir invariants (Section 6).

Additional Related Work
In addition to computer graphics and oceanography, we briefly review the fundamental research on the vorticity-streamfunction formulation in related areas.
1.1.1Computational Fluid Dynamics.The vorticity-streamfunction formulation is also extensively studied in the literature of computational fluid dynamics (CFD).There, the majority of attention is drawn to finding a compatible boundary condition in the presence of viscosity.Under viscosity, the boundary value of vorticity diffuses to the interior and can influence the resulting streamfunction from the Poisson solve.Therefore, the question becomes how to couple the boundary values of the vorticity and the streamfunction such that the resulting velocity satisfies both the no-through and no-slip conditions at the boundary.An initial pursuit for the boundary treatments by [Thom 1933] is followed by many variants [Taylor and Hood 1973;Orszag and Israeli 1974], revisions, and discussions [Quartapelle and Valz-Gris 1981;Anderson 1989;Quartapelle 1993;E and Liu 1996] throughout the remainder of the century.A comprehensive review on the highly debated topic of boundary conditions for the vorticity-streamfunction formulation can be found in [Rempfer 2006, §3.4].However, despite the large body of work searching for a proper vorticity-streamfunction formulation in CFD, few mentioned the effect of non-trivial topology.Different cohomology data of the flow can correspond to different sets of coupled vorticity-streamfunction boundary configurations, but these cohomological factors, which have their own dynamics, are irrotational and invisible to the vorticity variable.A few vorticity-streamfunction formulations that properly handle the cases of non-simply-connected domains are described in [Mizukami 1983;Tezduyar et al. 1988].However, these approaches only account for the non-simply-connectedness induced from obstacles rather than for general cases where cohomology components can arise.
1.1.2Geometric Fluid Dynamics.Geometric fluid dynamics is the mathematical discipline that discusses fluid dynamics in the geometric mechanics framework, including Hamiltonian systems [Salmon 1988;Morrison 1998] and geodesic equations on Lie groups [Arnold 1966].The mathematical foundation has inspired many structurepreserving fluid simulation algorithms in computer graphics [Elcott et al. 2007;Pavlov et al. 2011;Azencot et al. 2014;Liu et al. 2015;Chern et al. 2016;Yang et al. 2021;Nabizadeh et al. 2022].There, the Hamiltonian formulation of fluid dynamics often employs the vorticity-streamfunction formulation [Marsden and Weinstein 1983;Morrison 1998].Precisely, the phase space (with Poisson structure) is given by the space of vorticities, and the evaluation of Hamiltonian is defined with the aid of streamfunctions.Unfortunately, these descriptions would omit the cohomology components of the velocity.
A way to take cohomology components into account is to define the phase space as the quotient space Ω 1 / Ω 0 of velocity covectors modulo exact forms [Arnold and Khesin 1998, Theorem I.7.5;Oseledets 1989;Pavlov et al. 2011;Nabizadeh et al. 2022].However, this approach makes the phase space more abstract than the physically more understandable vorticity.To our knowledge, our work is the first to express the Hamiltonian formulation explicitly in terms of both the vorticity and the cross-sectional fluxes that parameterize the cohomology.Under this coordinate, we discover new Casimir invariants [Morrison 1998;Khesin and Chekanov 1989] that relate the cohomological flux and its linking number with the vortex lines.
A line of work that examines fluids' cohomology is the study of the commutant of the Lie group of volume-preserving diffeomorphisms [Arnold 1969;Banyaga 1978;Arnold and Khesin 1998, Definition I.7.12].They give an insightful characterization of the subgroup of flow maps generated by velocities that have no harmonic parts.This subgroup turns out to be the commutant of the group of all flow maps.However, the line of work did not comment on whether the physical fluid flow would stay or leave this commutant subgroup.The recent results on the helicity uniqueness conjecture [Khesin et al. 2022] were also restricted to this subgroup.Our work includes illustrative examples that show a fluid flow with initially zero harmonic components will later gain harmonic components.We show that fluid flows are generally not constrained in the commutant subgroup of the Lie group of volume-preserving flow maps.

BACKGROUND
In this section, we review the vorticity-streamfunction formulation for incompressible fluids.We present the theory of the paper in exterior calculus for its ability to unify 2D and 3D languages and to provide geometric intuitions.Readers who are less familiar with exterior calculus may find the following resources useful: (a) a comprehensive review of exterior calculus in Appendix A.1, (b) translations of exteior calculus operations to vector calculus in Table 1.Exterior calculus operations written in terms of vector calculus operations on a 2D domain .Here, ,  are scalar functions, a, b, v are vector fields, and J is the 90 • rotation operator.The map  :  ↩→  is the inclusion map of the boundary, whose normal vector is denoted by n.

Euler Equations
The incompressible inviscid fluid is governed by the following Euler equations.We assume the fluid density to be 1.On an -dimensional fluid domain  ( = 2 or 3), the time-dependent fluid velocity vector field u ∈ Γ( ) evolves under Here, n denotes the normal vector of the domain boundary , and the scalar function  ∈ Ω 0 () is the fluid pressure field.The covector (1-form) formulation [Nabizadeh et al. 2022] of (1) for the where  is the exterior derivative,  the codifferential, Lu the Lie derivative along u =  ♯ , and ★ the Hodge star.The Lagrangian pressure  L ∈ Ω 0 () is related to physical pressure by  L =  − 1 2 |u| 2 .The boundary condition is described using the pullback operator  * of the canonical inclusion map  :  ↩→ .Note that we call a differential form Dirichlet if it lies in ker(  * ), and co-Dirichlet if it lies in ker(  * ★). 4ll possible velocity fields form a linear subspace V 1 of all the co-closed (2b) and co-Dirichlet (2c) 1-forms.Definition 1 (Co-Dirichlet co-closed subspace).Let the subspace of co-Dirichlet co-closed -forms be denoted by The readers may verify that the co-Dirichlet co-closed subspace V  is  2 -perpendicular5 to the space im() of exact forms.Furthermore, together they span the space Ω  () of all k-forms.
As a result of this proposition, one can view the Lagrangian pressure term  L in (2a) as the  2 normal projection onto V 1 .Note that by Stokes theorem, the addition of an exact form  L does not affect the circulation along closed curves, i.e. ∮

Vorticity Equation
The vorticity 2-form is defined as  =  which measures local circulations.By taking the exterior derivative of (2a), applying  = 0, and using the commutativity between Lu and , we obtain the vorticity equation Eq. ( 4) is often regarded as a formulation "equivalent to (2a)."It is the basis of the class of fluid solvers known as vortex methods.
However, to fully establish the equivalence between (4) and (2a) (and to implement a vortex method), one must know how to reconstruct the velocity 1-form The reconstruction problem is posed as follows.
Unfortunately, the solution  may not be unique, depending on the fluid domain .If  ∈ ker() ∩ V  , readers can verify that  +  is also a solution to Problem 1.When  = 1, the subspace ker() ∩ V 1 is the space of all closed and co-closed 1-forms (corresponding to curl-free and divergence-free vector fields) satisfying the co-Dirichlet boundary condition (corresponding to the no-through boundary condition for vector fields).In general, the space ker() ∩ V  is the space of co-Dirichlet harmonic forms Here, the subscript "C" indicates the co-Dirichlet boundary condition  * ★ ℎ = 0.
2.2.1 Roadmap for Analyzing Problem 1.The standard method for reconstructing velocity from vorticity is through a streamfunction [Bridson 2015].We show that this streamfunction and its boundary condition emerge naturally in the context of solving a particular solution to Problem 1.We first introduce the concept of pseudoinverse of the exterior derivative in Section 2.3.Then we apply the idea to the problem of velocity reconstruction from vorticity in Section 2.4. im( )   argmin

Pseudoinverse of 𝑑
The pseudoinverse  +  ∈ Ω  () is given by the least-norm solution The general theory of pseudoinverses suggests that the pseudoinverse of a linear operator ( in our case) maps onto the orthogonal complement of the kernel of such a linear operator.Furthermore, the projection in (6) suggests that the pseudoinverse annihilates anything that is orthogonal to the image of the linear operator.Proposition 2. The space of -forms has an orthogonal decomposition The orthogonal projectors to im() ⊂ Ω +1 () and im( + ) ⊂ Ω  () are respectively Next, we characterize im( + ).Proposition 4. The image of  + is a subspace of V  .In particular: Putting together Propositions 1, 2, 4 and Eq. ( 5), we obtain Corollary 1.The space of -forms are orthogonally decomposed into:

Stream-Form
In the context of fluids, for each vorticity 2-form  ∈ im() ⊂ Ω 2 (), a particular solution to Problem 1 is  =  + .By Proposition 4, this velocity 1-form is given by the codifferential  of some  ∈ Ω 2 () that satisfies co-Dirichlet boundary condition.We call  the stream-form.
On a 2D domain, the stream-form  ∈ Ω 2 () is typically represented as  = ★ ψ by a scalar function ψ ∈ Ω 0 (), called the streamfunction.The co-Dirichlet boundary condition for  translates to the Dirichlet boundary condition ψ |  = 0 for ψ .For a streamfunctionrepresented velocity vector field, we have u = −J ∇ ψ where J is the counterclockwise 90 • rotation operator.
On a 3D domain, the stream-form  ∈ Ω 2 () is usually represented by a vector field  ∈ Γ( ) as  = ★ ♭ .The vector field  is called vector potential, stream vector field, or just streamfunction.The velocity is given by u = ∇ ×  and the boundary condition is that n ×  |  = 0, i.e. the stream vector field is normal to the boundary [Bridson et al. 2007].

Poisson Problem for Stream-Forms.
To concretely construct  ∈ Ω 2 () from  ∈ Ω 2 () one solves a Poisson problem.For details related to this Poisson problem and its boundary conditions, see Appendix B.
2.4.2Comments on 2D Streamfunctions.In many previous works involving 2D streamfunctions ψ ∈ Ω 0 (), the streamfunctions are allowed to have constant but nonzero boundary conditions.The velocity field  =  (★ ψ ) represented by this type of streamfunctions do carry the harmonic components.However, the existence of such streamfunctions only works for domains with special topology: Proposition 5.If a 2D domain  is the result of the removal of a few obstacles from a topological disk, then every  ∈ V 1 is coexact.
Proof.See Appendix E.5.□ In general, on a domain  that is a surface with a nonzero genus and with possibly a few obstacles removed, there are velocities that cannot be expressed by streamfunctions.

Summary
The vorticity data  ∈ im() ⊂ Ω 2 () is in one-to-one correspondence with a stream-form-represented velocity field  +  ∈ im( + ) ⊂ V 1 (Proposition 3).However, as shown in Corollary 1, the space im( + ) of stream-form-represented velocities is not the entirety of the space of all incompressible velocities V 1 .The gap is the space of co-Dirichlet harmonic forms H 1 C () that is isomorphic to the 1st de Rham cohomology  1 dR () = ker( ) /im( ), which is nontrivial whenever the domain is not simply-connected.In the context of Problem 1, H 1 C () is the space for non-uniqueness for the velocity reconstruction.
On a general non-simply-connected fluid domain, in order to pinpoint a velocity  ∈ V 1 in terms of vorticity, one must use both the vorticity data  and a co-Dirichlet harmonic form ℎ ∈ H 1 C (): While the evolution equation of the vorticity  is well-known, the evolution of the co-Dirichlet harmonic form ℎ has been overlooked in all vorticity-streamfunction-based methods.We investigate the time-evolution of this harmonic component of Euler equations in Section 3.

Summary in Vector Calculus
Here we reiterate the above background in vector calculus via Tables 1 and 2 to gain perspective in the vector counterparts of propositions about differential forms.
The fluid velocity that we wish to study is the space of vector fields that are divergence free and satisfying no-through boundary condition As documented in Tables 1 and 2, the vorticity  =  corresponds to the vorticity vector field w = ∇ × u in 3D or a vorticity scalar  = ∇ × u in 2D.Section 2.3 introduces the notion of  + , which can be seen as taking the inverse of the curl operator.This "curl −1 " is more precisely the pseudoinverse curl + of the non-invertible curl operator.By Proposition 4, the image of the curl + operator is given by velocities represented by streamfunctions (Section 2.4) subject to a specific boundary condition Known as the boundary-aware Helmholtz-Hodge decomposition, Corollary 1 asserts that the space Γ( ) of vector fields can be orthogonally decomposed into where  is the vector counterpart of (5) collecting harmonic vector fields satisfying the no-through boundary condition Proposition 3 further shows that there is a one-to-one correspondence im(curl + ) im(curl) between the space im(curl + ) and im(curl).This means that the vorticity vector field w in im(curl) (or scalar field  in 2D) is in one-to-one correspondence with a velocity field ) represented by a streamfunction field  (or a scalar function ψ in 2D) as characterized in ( 13).In particular, as demonstrated in ( 14), the space  of all incompressible velocities is larger than the space im(curl + ) im(curl) that can be captured by the vorticity data.The gap is the space  of harmonic vector fields with no-through boundary conditions.This gap  becomes nontrivial when the fluid domain is not simply-connected.Specifically,  is the kernel of the problem of reconstructing velocity from vorticity, which we discussed in Problem 1.
The decomposition ( 14) also reveals the difference between the velocity-based pressure projection and a vorticity-streamfunction solver in the context of advection-projection methods in fluid simulations.A pressure projection step removes im(grad) from Γ( ), keeping the information about the harmonic component  in  .This is consistent with the Euler equation.In contrast, a vorticitystreamfunction solver reconstructs the velocity using curl + from the vorticity data, effectively removing both im(grad) and  components from Γ( ).In particular, the vorticity-streamfunction solver leaves out the dynamics in the  component.As a result, different behaviors occurs between these two methods on non-simply-connected domains as demonstrated in Fig. 1.

THEORY
The goal of this section is to develop the full equations of motion for both , ℎ in (12).We first clarify the physical intuition to the harmonic component (Section 3.1).Next, we introduce the Lamb form (differential form counterpart of the Lamb vector), which is a central piece of the theory (Section 3.2).Finally, we derive the evolution equations for the harmonic part (Section 3.3), and explain the new physical law associated to it (Section 3.4).ACM Trans. Graph.,Vol. 42,No. 4

Harmonic Part and Flux
In the fluid literature, the harmonic part of a velocity field is often loosely depicted as "flows around obstacles or holes." Some parameterize the harmonic components by the circulations on a set of 1st homology basis (loops around obstacles) [Marsden and Weinstein 1983, §4;Elcott et al. 2007, §4.6].While it is technically true that there is an isomorphism  1 ()  1 dR () H 1 C (), different choices of loops lead to different parameterizations of the harmonic part of the velocity, even if the loops are only re-chosen within the same homology class.This is because the mapping between the circulations and the strengths of harmonic components depends on an arbitrary choice of representative loops for  1 () where we measure the circulations.For a consistent parameterization for the harmonic component, one must fix this arbitrary choice of loops.This is in contrast to picturing flowing loops as in the setup of Kelvin's circulation theorem.In fact, Kelvin's theorem on circulation conservation along flowing loops implies nothing about the conservation of the harmonic components [Thomson 1868].
As detailed below, we clarify a precise mapping between the harmonic component and a physical quantity of the fluid that is independent of the choice of artificial test geometry.
The harmonic component ℎ of a velocity field  in ( 12) is directly related to the physical fluxes over crosssectional surfaces.(16) The definition of the Flux operator can be written in vector calculus as well.In 3D, Flux(u)() , or that  1 and  2 are homologous.Therefore, Flux()(•) of any velocity data  assigns a well-defined flux data on each homology class  −1 (, ) =  −1 (, ) / −1 (, ) of cross-sections.The specific surface representing the cross-section from the homology is not important.
The following proposition summarizes the above discussion and formally treats Flux as a map that sends velocity data to a set of flux data over the homology classes of cross-sections.Proposition 6.For each  ∈ V 1 , the linear functional Flux() :  −1 (, ) → R is well-defined over the relative ( − 1)-homology  −1 (, ) =  −1 (, ) / −1 (, ).As such we obtain a bilinear form Flux()([]) Flux()() overloading on the same name: By currying, we may view the Flux operator as a linear map that sends a velocity to the dual space of the relative ( − 1)-homology: Proof.See Appendix E.6.□ Now, we show that ( 18) gives an isomorphism between the space H 1 C () ⊂ V 1 of harmonic components and the space  −1 (, ) * of flux data on cross-sections.The first observation is that the flux data Flux() is independent of (and only of) the stream-form part im( + ) from vorticity.This property is particularly neat as im( + ) is the orthogonal complement of Moreover, every assignment of flux data is realizable by some co-Dirichlet harmonic form.
Diagrammatically, Propositions 7 and 8 can be summarized as a short exact sequence.The diagram shows that the space V 1 of velocities is split into two "coordinates": the stream-form part im( + ) and the harmonic part H 1 C ().The stream-form part is parameterized by the vorticity data (Proposition 3), and the harmonic part is parameterized by the cross-sectional fluxes In sum, each fluid state is described by two equally important pieces of information associated with concrete physical measurements.They are the vorticity field and the cross-sectional fluxes.With this endowment of physical meaning to fluid's cohomology, we can discuss its expected behavior by drawing on physical intuition.a much more challenging time passing through that surface since the total flux over the surface is maintained at zero.Fig. 7 shows a simple demonstration of this setup with the Kirchoff point vortex model computed using Biot-Savart integration.Kelvin's method of reflection is employed for handling circular obstacles.To obtain unphysical dynamics with conserved harmonic components, one adds an additional point vortex at the center of each obstacle so that the flux between the obstacles stays zero.These additional compensating vortices repel the vortex pair in the physical domain.
In general, in a fluid solver where the harmonic components are constant, the approaching vortices induce a repulsing momentum that deflects the vortices.Section 5.2 includes a quantitative study (Fig. 15) of this phenomenon for a general fluid solver.

Lamb 1-Form
To describe the full evolution equation in Section 3.3, we first introduce a relevant quantity called the Lamb form.Definition 4. For each velocity field  = u ♭ ∈ V 1 , define the Lamb 1-form as  − u , where  = . (21) The vector counterpart of the Lamb form is the Lamb vector l  ♯ given by in 3D, and l = − J u in 2D, where  is the scalar vorticity.
Geometrically, the Lamb 1-form, as a cloud of codimension-1 objects, is the resulting ribbon surfaces (or curve segments in 2D) from extruding vortex lines (or point vortices in 2D) along the (negative) velocity field for a unit of time.Since the vortices are just passively transported by the fluid flow as described by the Lie advection equation ( 4), one may think of  as the trailing, or the "motion blur, " of the moving vortices.(24) The lamb form in (23) corresponds to the convection term in the Euler equations.Using 3D vector calculus, we have l = u × w = u × (∇ × u).Recalling that the momentum equation of Euler equation is   u + u • ∇u = −∇, we can plug in the vector identity u Similarly, we can express (24) using 3D vector calculus and relate  to the advection and the stretching term of the vorticity equation in vector form.Since  is a 1-form, taking its exterior derivative corresponds to the curl of l.Expanding the curl we have Note that u and w are divergence-free.Therefore, (24) can be written as Rearranging terms, we obtain w  = (w • ∇)u.

Equations of Motion
We are ready to lay down the governing equations for fluids' harmonic components.C (), be a basis for the co-Dirichlet harmonic 1-forms.Note that the dual space of H 1 C () is the space of Dirichlet harmonic ( − 1)forms via the following dual pairing Determined uniquely by the basis ( 1 , . . .,   ) for With the above sets of bases, we express a velocity field To extract the vorticity component, take the exterior derivative  =  as usual; to extract the coefficients c, apply the dual pairing with the dual basis harmonic forms Theorem 1.  ,c evolves under the Euler equations (2) if and only if the vorticity  and the coefficients c evolve according to where  is the Lamb 1-form for  (,c) .One may take an orthonormal basis ( 1 , . . .,   ) for H 1 C () by applying a Gram-Schmidt process or an economic QR decomposition on any other basis.In that case,   = ★  .Let h  = (  ) ♯ be the vector counterparts of the harmonic forms.Then (30b) reads Proof.Eq. (30a) is the same as the vorticity equation ( 4) or (24).
To obtain (30b), apply   to (29) to get The last equality is due to

)). □
An earlier appearance of the time-evolution of the coefficients of any spectral basis under the Euler equations was in the velocitybased model reduction work by [Liu et al. 2015].Our equation ( 31) is a special case of equation ( 7) from [Liu et al. 2015], since a basis for the harmonic component is a subset of a full spectral basis.We show that (30b) (or (31) in vector calculus) is the essential element for making vorticity-based solvers applicable to general domains, including non-simply connected ones.In the following subsection, by elucidating the theory behind cross-sectional flux, harmonic streamlines, and vortex lines, we demonstrate that (30b) leads to a conservation law with concrete physical and geometric meaning.

Practical Note.
Remark 1.In practice, to account for the missing dynamics of the fluid's harmonic component, it is sufficient to "charge" the coefficients   's with the global aggregate sum of (l • h  ) across the domain, in conjunction with any vorticity-streamfunction solver.
Remark 2. Note that when the harmonic basis is orthonormal, the coefficients  1 , . . .,   do not directly represent the numeral values of the fluxes on domain's cross-sections as described in Section 3.1.In the following Section 3.4.1,we describe which alternative basis to take in order to let the coefficients represent the fluxes directly.
Remark 3 (Physical units in an orthonormal setup).Including physical units, each   of an orthonormal harmonic basis for ), and the dual basis is of type   ∈ Ω −1 (; R m  /2−1 ).That is, they are R m ∓  /2±1 -valued after being integrated over 1-chains and ( − 1)-chains respectively.The vector counterpart h  = (  ) ♯ has the unit of m −  /2 , and the coefficient   has the unit of m  /2+1 s −1 .

Flux Dynamics
In the remainder of the section, we expand on (30b) and elucidate its physical meaning by drawing a relation to cross-sectional fluxes.The purpose is to shed light to new physically and geometrically intuitive principles obeyed by incompressible fluids on multiply-connected domains.These principles may facilitate productions of qualitatively plausible fluid animations.Additional theoretic insights in terms of Hamiltonian mechanics are discussed in Section 6.
3.4.1When the   's Become the Fluxes.As described in Section 3.1, the cross-sectional fluxes directly reflect the harmonic component of a velocity field.To make this relation more explicit, we pick a basis ( 1 , . . .,   ) for H 1 C , different from an orthonormal one used in (31) and Section 3.3.1,so that the   's become the physical flux on cross-sections.
Let  1 , . . .,   ∈  −1 (, ) be a set of representative crosssectional surfaces that forms a basis for  −1 (, ).Construct closed curves  1 , . . .,   ∈  1 () such that their signed intersection products (denoted by Let    ∈ Ω −1 () and    ∈ Ω 1 () be the Dirac  forms concentrated on these curves   's and surfaces   's.Define The field   is the harmonic field "spread out from the concentrated current   ," as the result of a stream-form part removal (1 −  + ) from the current.The field   is the harmonic field "pumped out from the impulse at the cross-section   , " as the result of the pressure projection (1 −  + ) (exact form removal) of the impulse.Using vector calculus,   ∈ H 1 C () and   ∈ H −1 D () can be represented as harmonic vector fields Z  =   and X  = (★ −1   ) ♯ respectively.They are constructed as Note that ( − 1)-forms   's geometrically represent families of curves.These curves are the streamlines (integral curves) of the associated vector fields X  's.Therefore: Definition 5 (Harmonic stream).We call   in (33b) the harmonic stream(lines) associated to a cross-sectional surface   ∈  −1 (, ).Proposition 9.For   's and   's defined by (34), we have ∫    ∧  =    .Moreover, for each velocity  ∈ V 1 (u =  ♯ ∈  ), the coefficients   's in (28) are the fluxes through the cross-sectional surfaces: Remark 5. Once the basis  1 , . . .,   is constructed from a generator basis  1 , . . .,   , one can construct the dual basis  1 , . . .,   through a simple QR-decomposition followed by a small matrix inversion without the hassle of building  1 , . . .,   .The detail of this process is described in Section 4.1.2 Remark 6 (Physical units for the flux setup).Similar to Remark 3, we discuss the physical units for the harmonic bases (34) and the coefficients.Each harmonic stream   is of type   ∈ Ω −1 (; R m n−2 ), and its vector counterpart (★  ) ♯ has the unit of m −1 .The dual basis form   is of type   ∈ Ω 1 (; R m −n+2 ), whose vector counterpart (  ) ♯ has the unit of m −n+1 .The coefficients   has the unit of m n s −1 , which is the unit of a total flux.

Fluxes and Linkings between Vortices and Harmonic Streamlines.
Just as in the general theory (30b), the flux (35) through the -th cross-sectional surface   has a rate of change given by In the codimensional geometric picture of differential forms, recall Section 3.2 that the Lamb 1-form  is the collection of the "motion blur" ribbon surfaces trailing behind the flowing vortex lines, and   is the set of harmonic streamlines.The integrated wedge product ∫   ∧   is the total number of signed intersections between the motion blurs of vortex lines with the harmonic streamlines.In other words, ∫   ∧   is the rate at which vortex lines of  cut through the harmonic streamlines of   .Therefore, ∫   ∧   is the rate of change of the "linking number" between the flowing vortex lines of  and the static harmonic streamlines of   .We let this linking number be denoted by Link()(  ) (whose subtle mathematical definition for general manifolds will be discussed later).As such, ∫   ∧   =   Link()(  ).Therefore, ( 37) is stating about a balancing relation between two rates of changes, one about the flux through   , and the other about Link()(  ).We conclude this discovery in the following theorem.Theorem 2. In an Euler fluid, for each cross-sectional surface   which generates harmonic streamlines   , the difference between fluid's flux over   and the linking number between the vortex lines and harmonic streamlines is a constant of motion.(See Fig. 4.) 3.4.3Linking in 2D Domains.In 2D, the vorticity 2-form  is geometrically represented as a cloud of point vortices, instead of a cloud of vortex lines.In that case, Link()(  ) is understood as the total winding number of the harmonic streamlines about the point vortices.Such linking/winding numbers can be defined for the following special 2D domains.Suppose  has the topology of a disk with  obstacles removed.Then every cross-sectional curve   connects two of the ( + 1)boundary components.Moreover, the associated   ∈ H 1 D () is exact (Proposition 5): There exists harmonic functions   ∈ Ω 0 (), unique up to a constant, such that Geometrically, the harmonic streamlines of   are the level sets of the scalar potential   .Therefore, the winding number of these level lines around the vortex points of  admits an explicit formula: Corollary 2. For an Euler fluid on a disk with  obstacles, the difference between the fluid flux over   and the quantity (40) is a constant of motion.
Appendix D provides an analytical example for Corollary 2.
3.4.4Linking in Euclidean Domains.When  ⊂ R  , every Dirichlet closed forms, such as   , is exact   =   for some   ∈ Ω −2 () [Shonkwiler 2009;Zhao et al. 2019].Fixing any representative potential   , the linking number admits an explicit formula: In 2D,   corresponds to a scalar function   in a 2D Euclidean domain, and the harmonic vector field X  , which is the vector counterpart of   , satisfies X  = −J grad   .In a 3D Euclidean domain,   corresponds to a vector field a  with X  = curl a  .In terms of these vector calculus counterparts, we have Link 3.4.5Linking in General Domains.For general manifolds , defining the linking number is trickier.To make sense of Theorem 2 we define the linking number between  and   only through its variation with respect to .In that sense, Theorem 2 is understood as that this linking variation always balances out with the variation in the cross-sectional flux.For a detailed discussion on defining linking, see Appendix C.

IMPLEMENTATION
In the previous section, we derived the new dynamical equations for fluids' harmonic parts (30b) and the underlying physical law (Theorem 2).In Section 5, we use numerical examples to demonstrate the effect of restoring such harmonic dynamics in vorticitystreamfunction solvers.The results will show that (30b) is crucial for reproducing realistic fluid dynamics in both 2D and 3D fluid simulations.In this section, we describe the implementation details for incorporating (30b) into existing state-of-the-art vorticitystreamfunction solvers.Our algorithm is not limited to the particular choice of solver we integrate with, but could be applied to different advection schemes.
For 2D examples, we adopt the method of Functional Fluids [Azencot et al. 2014] with a 4th-order Runge-Kutta time integration.This method is an accurate vorticity-streamfunction solver that applies to general triangulated surfaces with boundaries.The algorithm is implemented in SideFX's Houdini 19.5, and the source code is available in the supplementary material.
For 3D examples, we employ Covector Fluids [Nabizadeh et al. 2022].The original Covector Fluids solver is a velocity-based advectionprojection method which provides base-line ground truth references.By replacing the pressure projection using a streamfunction solver [Ando et al. 2015a], the Covector Fluids solver becomes equivalent to a circulation-preserving vorticity-streamfunction method [Elcott et al. 2007] with a higher order advection bootstrapped by Backand-Forth Error Correction and Compensation (BFECC) [Kim et al. 2005;Selle et al. 2008].The algorithm is implemented in the C++ codebase provided by [Nabizadeh et al. 2022].
In each of these vorticity-streamfunction-based fluid simulators, we present algorithms for updating harmonic components.

2D Implementation on Triangle Meshes
We employ the method of lines for the 2D fluid simulations.That is, we spatially discretize the PDE (30), leaving a continuous-time ODE which is subsequently integrated using the 4th-order Runge-Kutta (RK4) method.The spatial discretization of the advective system on a triangle mesh follows [Azencot et al. 2014].
Let  = (P, E, F) be a triangle mesh.A discrete vector field is a piecewise constant vector u = (u f ∈  f ) f ∈F assigned on the triangles.Each vector field u is associated with a directional where  f is the area of triangle f ∈ F, and (grad  ) f is the gradient of the piecewise linear interpolated function  in face f ∈ F.8 An explicit formula for the gradient is given by where J f is the 90 • counterclockwise rotation within the tangent plane of the triangle f, and e f,p is the (unnormalized) edge vector opposite to the point p across the triangle f.
Next, we build the discrete ξ  = ★ −1   ∈ H 1 C () where   ∈ H −1 D () describes the harmonic stream described in Definition 5.For this process, we employ discrete exterior calculus (DEC) [Hirani 2003].For each   , build a discrete primal 1-form that represents the impulse    .That is, for a primal edge e ∈ E crossed by   we set (   ) e to be ±1.The sign is chosen so that d 1    = 0, where d 1 is the discrete exterior derivative operator (co-boundary operator) on 1-cochains.Each harmonic form ξ  is built by the pressure projection (1 −  + )   in the DEC sense.
Next, each ξ  is interpolated into a piecewise constant vector field ( X ∈  f ) f ∈F using Whitney interpolation [Bossavit 1998].Now, we  2 -orthonormalize X1 , . . ., X with respect to the inner product structure ⎷u, v⌄ = f ∈F ⟨u f , v f ⟩ f .For each  = 1, . . ., , pre-multiply the area factor X,f = Finally, we obtain an orthonormal basis (h 1 , . . ., h  ) for Alternatively, one can also use a randomized algorithm to obtain a set of orthonormal harmonic basis as explained in Section 4.2.1.

Harmonic Stream
Basis for Flux Coefficients.While an orthonormal basis (h 1 , . . ., h  ) for H 1 C () is convenient in computation, its associated coefficients do not have direct physical meaning.
Here, analogous to Section 3.4.1 we describe an alternative basis for H 1 C () so that  1 , . . .,   represent the total fluxes through  1 , . . .,   .Continuing (48), invert R and reassemble the orthogonal basis h1 , . . ., h into Define The bases (Z  ) and ( X ) are the discrete analogs of (  ) ♯ and (★ −1   ) ♯ defined in ( 34) respectively.Under these bases, one replaces h by Z in the velocity reconstruction (45) and Alg.1; and one replaces h by X in the c updates (46b).By doing so, the dynamics remains the same, but the coefficients  1 , . . .,   now represent the physical fluxes on the cross sections  1 , . . .,   .

3D Implementation on Staggered Grids
For 3D numerical examples, we integrate our algorithm to the codebase of Covector Fluids [Nabizadeh et al. 2022] for its equivalence to a circulation-preserving vorticity method.In particular, we replace the pressure projection step with a streamfunction solver applied to the vorticity.The discretization uses the standard MAC grid  = (V, E, F, C) to store the variables.Similar to many other grid-based streamfunction solvers [Ando et al. 2015a;Chang et al. 2022], we store the velocity on grid faces, and store vorticity and streamfunction on grid edges.We explain solving the streamfunction Poisson problem in details in Appendix B. So far, this is a classical 3D vorticity-streamfunction solver that does not have the dynamics of the harmonic components.
In our method, we incorporate (31).We store our harmonic basis as  vector fields on the MAC grid faces, similar to the velocity fields.To evaluate the right-hand side ∫  det(u, w, h  )  of ( 31), we interpolate the velocity ( f ) f ∈F , vorticity ( e ) e∈E , and harmonic basis (ℎ ,f ) f ∈F respectively into vector fields (u c ) c∈C , (w c ) c∈C , (h ,c ) c∈C that sit on the cell centers for easy local computations.With this discretization, we include the dynamics (31) by adding the following step in the main solver: where  c is the cell volume.The harmonic coefficients are also used to reconstruct velocity after the streamfunction solver step.We summarize our overall procedure in Alg. 3.  24) 11: 31) 13: 4.2.1 Generate Harmonic Basis in 3D.Methods for computing harmonic fields in a 3D domain are discussed in [Zhao et al. 2019].
In our case, we calculate an orthonormal basis for the harmonic components through a randomized algorithm.We generate  random vector fields where  is the dimension of H 1 C ().We store these vector fields as flux on the faces of the staggered grid.We subtract the exact and co-exact components of these vector fields by applying (1 −  + −  + ), i.e. by taking the difference between a standard pressure solve and a streamfunction solve.This procedure gives us  harmonic vector fields.These harmonic vector fields are almost surely linearly independent.We then perform a QR decomposition (e.g. a Gram-Schmidt or a Householder process) to obtain an orthonormal basis for the harmonic components.See Alg. 4.

NUMERICAL EXAMPLES
In this section, we demonstrate the results from our incorporation of the dynamical harmonic components into vorticity-streamfunction fluid solvers for 2D surfaces (Alg.2) and 3D scenes (Alg.3).

2D Oscillating Fluxes with Closed Lamb Form
We design a setup so that the vortex dynamics (30a) vanish while a nontrivial evolution of the harmonic part (30b) is present.With this design, we single out the effect of the new equation (30b).We obtain the solution both numerically and analytically.In particular, the results demonstrate the necessity of (30b) for realistic fluid animation in contrast to previous methods.

Design
Rationale.How do we keep the vorticity equation steady while keeping the harmonic part dynamic?A classically known result in fluid mechanics directly following from ( 23) is that the flow is steady (   = 0) if and only if the Lamb form is the gradient of the Bernoulli pressure  =  B .Now, what if we relax the condition of  ∈ im() to just  ∈ ker()?The closedness of  will still ensure that the vorticity is steady (   = 0) using (30a).But such a non-exact closed , which contains harmonic components, can yield nontrivial dynamics in (30b), making the overall flow unsteady.
In 2D, a simple way to obtain a closed  is to set vorticity constant, say  = ★1.Then  = − u  = − ★ , which is closed since  = 0. Next, we need to make sure  is not exact.Observe that if the domain is the result of the removal of a few obstacles from a simplyconnected disk, every  = − ★  is exact (Proposition 5).Therefore, the simplest non-trivial example is a surface with nonzero genus.The surface must also have at least one boundary component since a closed surface must have zero total vorticity.

Analytic Solution.
For constant vorticity  =  ★ 1, the Lamb , where ψ = ★ is the streamfunction.Since  is closed, by (30a) the vorticity stays static.The only dynamics left is (30b), which now reads and similarly That is, ( 1 ,  2 ) ⊺ evolves in a simple harmonic oscillation which has an explicit solution The frequency  of the oscillation is proportional to the background vorticity  ∈ R( 1 /s) and a dimensionless number  depending on the geometry.Note that the construction of  relies only on the ★ on 1-forms.Hence  depends only on the conformal type of the surface [Soliman et al. 2021].

Numerical Results
. We apply our fluid solver (Alg.2) to the configuration of Section 5.1.2.We set  = 2 s −1 ,  1 (0) = 0.5 m 2 /s and  2 (0) = 0 m 2 /s.The solver is set with time discretization  = 0.1 s.The result is compared against the previous vorticity-streamfunction method on surfaces [Azencot et al. 2014], which is our algorithm without (30b).Fig. 10 shows the numerical values of  1 and  2 over time.In particular, they agree with the analytical solution  1 () = 0.5 cos(),  2 () = −0.5 sin() (cf.( 56)).Previous methods that keep ( 1 ,  2 ) = (0.5, 0) constant over time would lead to a vastly different flow both quantitatively and qualitatively.In Fig. 11, we advect colors to compare the resulting flow map of the previous method [Azencot et al. 2014] and ours.In [Azencot et al. 2014] the inertia is not correctly transferred, creating an unnatural laminar texture that is not turned by the background vorticity.The problem is solved in our method which incorporates the evolution of the harmonic components.

2D Vortex Pair Between Obstacles
Let the domain  be a flat hexagonal disk with two hexagonal holes removed, similar to an ocean basin with two islands.The dimension of its 1st relative homology is two, and therefore the dimension of the harmonic basis is two, as illustrated in Fig. 13.
We set a pair of vortices with opposite values as shown in Fig. 12 at  = 0 s.In principle, we expect the traveling vortex pair to pass through the space between the islands.However, using the method [Azencot et al. 2014] with a fixed harmonic part, the two vortices do not pass through.Instead, they turn around and go back in the opposite direction when they approach the islands.This unphysical behavior is due to a vanishing total flux along the curve connecting the two holes maintained by the fixed harmonic part.Using our method which includes (30b), the vortices pass through as expected.Fig. 12 shows the dynamics of the vortices over time for both methods.
In addition to inspecting the vortex motion, we also determine the correctness of the competing methods quantitatively.In Fig. 14 we measure the circulations around one of the islands over time for both methods.Since the flow is always tangential to the boundary, the boundary curve of the island that flows with the fluid will stay as the same boundary curve.In particular, Kelvin's circulation theorem or island rule [Godfrey 1989;Pedlosky et al. 1997] applies.The circulation along the island boundary should be conserved.As shown in Fig. 14, our method better conserves the circulation, while the method with a fixed harmonic part violates this conservation law.In Fig. 15, we consider a cross-section , and measure the corresponding Flux, Link, and Flux − Link.Using our method, Flux − Link is a conserved quantity.In contrast, using the method with a fixed harmonic part, the flux is always zero, and Flux − Link is not conserved.

3D Examples on Non-simply-connected Domains
Similar to the 2D examples, when one does not incorporate the dynamics of harmonic components into vorticity-streamfunction methods, unphysical behaviors occur in a domain with non-trivial topology.We set up three experiments by introducing different obstacles (fan, torus, and pillars) into a closed box domain, making the domain non-simply connected (see Fig. 16).
In Experiment 1 (Fan) we place a torus-shaped fan in the domain which raises the dimension of its 1st cohomology to one.Experiment 2 (Torus Rock) is staged with a torus-shaped rock submerged in an aquarium.Similar to the previous setup, the 1st cohomology's dimension is one.Experiment 3 (Pillars Rock) is also set in an aquarium, but the obstacles are two pillars.Since the pillars touch both the top and the bottom of the boxed domain, the dimension of the 1st cohomology is two.We illustrate the harmonic basis for Experiment 1 (Fan) and Experiment 3 (Pillars Rock) in Fig. 17.
We run our experiments on uniform staggered grids with a voxel resolution of 150 × 75 × 75 in a box size of 10 × 5 × 5 m 3 .The timestep duration is set to 1 /24 seconds with two substeps per frame.A vortex ring of radius 0.4 m is initialized to face in the positive -direction with strength 1.6 m 2 /s.As mentioned in Section 4, for 3D flows, we use Covector Fluids solver [Nabizadeh et al. 2022] as our velocity-based method reference (right).We replace the pressure  projection step with a streamfunction Poisson solver as the results which represent classical methods with fixed harmonic parts (left).
Our method includes the evolution of harmonic parts (30b) (middle).
Fig. 1, and Fig. 16 demonstrate the results of the three experiments detailed above.A significant amount of vorticity is trapped by the obstacles when one uses a traditional method that fixes harmonic parts (left).This behavior does not match the results in the velocity-based method reference (right).In particular, in the result of Experiment 1 (Fan) computed by the method fixing harmonic parts (left), only a small portion of the vortex passes through the fan while the majority of the vortex moves in the opposite direction to compensate for a vanishing total flux across the fan.Using our method (middle), the vortex passes through the obstacles with an overall dynamical behavior similar to the velocity-based method reference (see video).
In conclusion, the numerical examples demonstrate that a classical vorticity-streamfunction solver can generate incorrect flow patterns.Our method for including the dynamics of harmonic parts solves this problem.

HAMILTONIAN FORMULATION
One of the highlights in the theory of Hamiltonian Fluid Mechanics is interpreting the vorticity equation ( 4) or (30a) as a reduced infinite-dimensional Hamiltonian system in classical mechanics.In this final discussion section, we extend the Hamiltonian description to include our new equation (30b) of the time-evolution of the harmonic coefficients.Readers can find introductions to the most common Hamiltonian formulation for fluid dynamics in [Salmon 1988] using elementary continuum mechanics, [Morrison 1998] using non-canonical transformations, or [Arnold 1966;Marsden and Weinstein 1983;Arnold and Khesin 1998] using group theory.
Background.The setup for a Hamiltonian formulation involves several steps.One first establishes a set of variables that describe the state of the physical system.The space of all possible states is called the phase space M. One then describes a symplectic structure.A symplectic structure is a non-degenerate closed 2-form  ∈ Ω 2 (M) that encodes the interrelation among the variables such as positionmomentum paring.When the phase space has a symplectic structure, it is called a symplectic manifold.In many weaker cases including the case of fluids, the phase space M is merely foliated into many symplectic submanifolds called symplectic leaves.Such a phase space is called a non-canonical space or a Poisson manifold [Weinstein Note that using the method with a fixed harmonic part, the vorticity appears to be trapped by the obstacles, especially in Experiment 1 (Fan).Our method resolves this issue and produces visually similar results to ground truth (velocity-based advection-projection scheme).The above setup is typically easy to describe for particle systems and their continuum limits, as one may define the positions and momenta as the obvious ones for each particle.The non-trivial aspect about Hamiltonian Fluid Mechanics is after the so-called reduction by sorting out the particle-relabeling symmetry [Salmon 1988].After the reduction, the position and momentum of each individual particle are no-longer parts of the coordinates of the phase space.Instead, what naturally emerges is that the phase space for incompressible fluids is the quotient space M = Ω 1 ()/ Ω 0 () of velocity 1-forms modulo an exact form ( of Lagrangian pressure) [Arnold and Khesin 1998, Theorem I.7.5].In fact, this phase space is a Poisson manifold with symplectic leaves given by the orbits of Lie transportations of velocity 1-forms under any volume-preserving flow maps (i.e. the pullback of velocity 1-forms under the inverse of the flow maps).In particular, Kelvin's circulation theorem can be thought of as a Casimir-typed conservation law.A more concrete example of Casimir in a 3D incompressible fluid is the helicity ∫   ∧  that measures the self-linking number of vortex lines.
To make the quotient space M = Ω 1 ()/ Ω 0 () less abstract, a common practice is to describe M = im( 1 ) ⊂ Ω 2 () as the space of vorticity fields (by assuming there is no cohomology component) [Marsden and Weinstein 1983, §4].Such a treatment is attractive, since vorticity has a concrete physical meaning, and the Hamiltonian flow ( = total kinetic energy) on the vorticity variable directly yields the familiar vorticity equation ( 4) or (30a).The drawback is of course it does not describe fluids on non-simply-connected domains.
Our Modification.Similar to our Theorem 1, the fluid state is described by both vorticity  ∈ im( 1 ) ⊂ Ω 2 () and coefficients of harmonic components c = (  )  =1 ∈ R  .Importantly, as described in ( 19) and Section 3.4.1, the coefficients c have the concrete physical meaning of flux through cross-sections  1 , . . .,   of the domain.Recall that each cross section   is associated with a harmonic field   defined in (33b).Our phase space M is given by the coordinate of vorticity data and the flux data: This phase space is a Poisson manifold.We describe the Poisson structure by defining its symplectic foliation as follows.A tangent vector ( ω, c)| (,c) ∈  (,c) M at state (, c) ∈ M is tangent to the symplectic leaf if it takes the following form for some divergence-free vector field  ∈ Γ( ).One can check that the distribution of tangent subspaces (58) is integrable10 and hence form a foliation.We define the symplectic form  on the symplectic leaf: For every two tangent vectors ( ω, c) = (−  , ∫  (−  ) ∧ ), ( , c) = (−  , ∫  (−  ) ∧  ) tangent to the leaf (where ,  ∈ Γ( ) are arbitrary divergence-free vector fields), This symplectic form is the same as the one defined in the literature [Marsden and Weinstein 1983] but now lifted to a larger space (57).Note that the definition (58) ensures that Flux()(  ) − Link()(  ) is always invariant under any variation along a symplectic leaf.We elaborate on the variation of Link in Appendix C. In particular, Flux() (  ) − Link()(  ) is a Casimir for each  = 1, . . .,  (cf.Theorem 2).Finally, one may verify that our full equation ( 30) is the Hamiltonian flow of the total kinetic energy In summary, the above equations complete the vorticity-based Hamiltonian description of incompressible fluid dynamics in domains with general topology.This mathematical framework may serve as a foundation for future research in mathematical fluid mechanics.

CONCLUSION
In this paper, we tackle a long-existing problem in the vorticitystreamfunction formulation of incompressible Euler fluids on nonsimply-connected domains.We demonstrate that the dynamics for the harmonic (cohomology) components of incompressible Euler fluids can be described by (30b) (Theorem 1).This have been overlooked by previous vorticity-streamfunction solvers.We give a simple and practical algorithm in Section 4 that easily incorporates such dynamics into previous vorticity-streamfunction solvers.
Our numerical examples (Section 5) demonstrate that it is necessary to include the additional equations in a vortex solver to avoid unphysical artifacts in fluid animations.
The proposed algorithm for 2D surface simulation (Alg.2) is particularly significant because while vorticity-based methods (involving only scalar advections) are much more straightforward to compute on a 2D surface domain compared to a velocity-based solver (involving covariant vector advection), previous vorticity-based surface fluid solvers do not include the dynamics of the harmonic components.Our Alg. 2 is the first vortex solver that is consistent with the Euler equations on surfaces with arbitrary topology.
We also find a new physical conservation law associated with the new evolution of harmonic components (Theorem 2).We describe our new equation as a Hamiltonian system on the state space coordinated by the vorticity and the cross-sectional fluxes (Section 6).Compared to previous editions of Hamiltonian formulation for incompressible Euler flows, the new framework does not omit the cohomology and each variable is a physically meaningful measurement.The mathematical investigations have also led us to discover new Casimir invariants as well as interesting analytic flows, e.g. the example presented in Section 5.1.It would be exciting to explore if this new conserved quantity that we discovered could serve as the foundation for new algorithms in the future.
This paper leaves a few open questions that are beyond our current investigation.We only work on inviscid Euler fluids on a fixed oriented Riemannian manifold.In particular, we do not consider moving obstacles or that the domain is an evolving surface.Expanding our analysis to moving domains can yield new perspectives to the studies of underwater swimmers and solid-vortex interactions [Vankerschaver et al. 2009;Weißmann 2014].Lastly, what is the effect of viscosity?What if the domain has changing topology?What if the domain is non-orientable?We expect exciting research answering these and related questions in the near future.

A PRELIMINARIES A.1 Differential Forms
Let  be an -dimensional Riemannian manifold.The space of vector fields is denoted by Γ( ), and the space of differential -forms, 0 ≤  ≤ , is denoted by Ω  ().
A differential -form  ∈ Ω  () is a formal object to-be-integrated over an oriented -dimensional test surface .This evaluation method is denoted by ∫  .For example, in 3D, a 3-form describes a density field or a measure that awaits being integrated over a volumetric region.A 2-form describes a flux that is to be evaluated over an oriented surface.A 1-form is to be line-integrated into circulations, and 0-forms are synonyms for scalar functions that are to be evaluated over points.
A.1.1 Forms are Distributions of Codimensional Geometries.Differential -forms should not be confused with vector fields in their visual representations.Geometrically, a vector field v is an assignment of an infinitesimal "arrow" v  ∈    at every point  ∈ , whose directions and magnitudes depict an instantaneous flow velocity within the domain .A -form  ∈ Ω  (), on the other hand, is a distribution of ( − )-dimensional (codimension-) oriented planes over .For example, in 3D, a 3-form is illustrated as a point cloud, a 2-form is a line segment cloud, a 1-form is a plane field, and a 0-form is a superposition of sublevel sets of the corresponding scalar function.The orientations and densities of the codimension- plane cloud are given so that the integration ∫   over a test -surface  is the total signed intersection between  and the codimensional- plane cloud.A 1-form can be converted into a vector, and vice versa, by the sharp ♯ : Ω 1 () → Γ( ) and flat ♭ = ♯ −1 : Γ( ) → Ω 1 () operators.For a 1-form  ∈ Ω 1 (), the vector field  ♯ is pointwise an arrow whose direction is orthogonal to the hyperplane of  at the point, and whose magnitude is set as | ♯ | = | |.An ( − 1)form  (typically representing a flux) can also be converted into a vector by (★) ♯ .Both 1-forms and ( − 1)-forms are often identified as vector fields;  ♯ is a vector field orthogonal to the hyperplane field  ∈ Ω 1 (), and (★) ♯ is a vector tangential to the line field  ∈ Ω −1 ().See also Tables 1 and 2.
A.1.3Exterior Derivatives are Boundary Operators.The exterior derivative operator  : Ω  () → Ω +1 () takes the boundary of Remark 7. The condition (71d) in the Coulomb gauge ensures the uniqueness of (71).In practice, the condition (71d) is sometimes neither mentioned nor enforced.This is fine, as the remaining system of (71) only has a small finite dimensional kernel of H +1  (), which does not jeopardize most linear solvers and does not influence the result of  in (70).
On a triangle mesh, we set both  and ψ on vertices following [Azencot et al. 2014].Another possible approach is to let the values of  and ψ sit on edge centers, and represent ψ using the Crouzeix-Raviart elements [Poelke and Polthier 2016].

B.2 Solving Streamfunctions on a 3D MAC Grid
Here, we describe an implementation of (71) in 3D under the Marker-And-Cell (MAC) discretization scheme [Harlow and Welch 1965;Bridson 2015].The streamfunction solve is a simplified special case of [Ando et al. 2015a] without a varying density.In the MAC scheme, the velocity field is given as fluxes assigned to the faces on a regular grid.The vorticity and the streamfunction are defined on edges.That is, the MAC scheme is equivalent to a Discrete Exterior Calculus (DEC) scheme for our exterior calculus formulation up to a Hodge dual.For notation distinction, let η = ★ ∈ Ω 2 () be the velocity flux defined on faces, ω = ★ ∈ Ω 1 () be the vorticity 1-form defined on edges, and ψ = ★ be the stream 1-form on edges.
Let the grid (V, E, F, C) be organized by the vertex set V, edge set E, face set F and cell set C. We assume every edge has length ℎ.The adjacency and orientation between the cubical complex are summarized in the discrete exterior derivative matrices (a.k.a.coboundary operators) forming a (co)chain complex Note that each matrix d  has entries either −1, 0, 1.The division by ℎ in a finite difference derivative will be denoted separately, keeping d  's purely topological.
To assign obstacles, for each cell in C, label "fluid (interior  )" or "solid (boundary )." Next, label fluid and solid on V, E, F by taking the closure of the solid cells.That is, first initialize each element of V, E, F as fluid, and then label every boundary face of a solid cell a "solid" face; subsequently, label every boundary edge of every solid face as solid, and finally label every boundary vertex of every boundary edge as solid.After consistently labeling of fluid ( ) and solid () elements and similarly for the other d  's.One checks that the interior (fluid) blocks form a cochain complex: (77) Take "curl" 1 ℎ d ⊺ 1,  on both sides of the equation: We label ω ∈ R |E  | with subscript  to acknowledge that ω = Here we provide a well-defined linking number (Theorem 2) between two fields on a general manifold .This linking number between smooth fields is also called the cross-helicity.For general domain, define the linking Link()() between a closed -form  and a closed Dirichlet ( −  + 1)-form  by its variation with respect to local transportations ω = L  while fixing : Now, we verify this variational gradient is indeed a gradient by checking that the second derivatives commute.Take two vector fields  1 ,  2 with Lie commutivity [ 1 ,  2 ] = 0 (i.e. two coordinate directions on the diffeomorphism group Diff ()).The mixed second derivative is (87) In particular, we find that which is indeed constant over time.
The splittings (8) imply that ,  + are isomorphisms between the subspaces im( + ) and im().Moreover, when being restricted to im( + ), im(), they are the inverse of each other by construction.These properties imply (9) and Proposition 3. □

E.5 Proof of Proposition 5
To show the coexactness of , or equivalently the exactness of η ★, it suffices to check that ∮  η = 0 for all closed curve  ⊂ .If  is the complement of a few obstacles in a simply-connected domain, then every  is homologous to a boundary curve along which ∮  η = 0 using the no-through condition  * ★ .□

Fig. 3 .
Fig. 3. Side-view of smoke passing through a torus-shaped obstacle.Also see the top row of Fig. 1 (see video 00:31).
Weißmann and Pinkall 2012; Ishida et al. 2022], a vortex sheet over the obstacle surface [Park and Kim 2005; Weißmann and Pinkall 2010; Fig. 4. The change of linking between harmonic streamlines (blue) and fluid vortex lines (green) in a fluid domain exterior to the obstacle (torus).

Fig. 6 .
Fig.6.The rate of change of the harmonic coefficients (  1/, . . .,  / ) is computed by projecting the Lamb vector onto the harmonic basis.The Lamb vector is the cross-product of velocity and vorticity fields.
Fig.7.A vortex pair approaching two circular obstacles in the plane.The physically correct behavior for the vortices is to pass through the space between the two obstacles (top row).An artificial fluid dynamics that preserves the harmonic component of the flow repels the vortices from passing between the obstacles (bottom row).The simulation is computed using Kirchoff's point vortex dynamics with Kelvin's method of reflection for handling circular obstacles.

Fig. 8 .
Fig. 8. Left: vortices on a surface with genus = 2 simulated using our method.Right: dye advected by the vortices.
Fig. 9.A surface with genus = 1 and one boundary component.

Fig. 10 .FluidFig. 11 .
Fig.10.Evolution of the coefficients of the harmonic components for the setup of Section 5.1.4,simulated by our method described in Section 4.1.

Fig. 12 .
Fig. 12.The time-evolution (from left to right) of a vortex pair moving through a fluid domain with two small hexagonal obstacles, simulated by [Azencot et al. 2014] (top) and our method (bottom) (see video 01:00).
Fig.16.Visualizations of vorticities from Experiment 1 (Fan), Experiment 2 (Torus Rock), and Experiment 3 (Pillars Rock).Note that using the method with a fixed harmonic part, the vorticity appears to be trapped by the obstacles, especially in Experiment 1 (Fan).Our method resolves this issue and produces visually similar results to ground truth (velocity-based advection-projection scheme).

Fig. 17 .
Fig. 17.Harmonic bases for Experiment 1 (Fan) in the left figure, and Experiment 3 (Pillars Rock) in the middle and right figures.Note that our domain is inside a closed box, but we only include the bottom side of the box to better visualize the harmonic fields inside.
Conversions.The Riemannian metric over the manifold defines a notion of distance and orthogonality.The magnitude | | of a -form  is the throughput of the codimension cloud across an orthogonal cross section of unit -area.The Hodge star converts a -form to a ( − )-form, ★: Ω  () → Ω +1 (), so that the codimension cloud of ★ is pointwise the orthogonal complement of the codimension cloud of , and that they share the same magnitude |★ | = | |.The parity is chosen so that ★ ★  = (−1)  (− )  for a -form . ★  ♯ 76) which is known as the relative cochain complex.Now, the no-through boundary condition  * η = 0 and the Dirichlet condition  * ψ = 0 translates to the fact that η and ψ only live on F  and E  respectively, and vanish over F  , E  .The discrete relation between the streamfunction ψ ∈ R |E  | and η ∈ R |F  | is given by η = 1 ℎ d 1,  ψ .

p
84) where we have used  = 0 and the Dirichlet closedness of  for ∫   (•) ∧  = 0 in the second to last equality.Therefore, for each fixed , Link()() is well-defined locally in the space of 's Lie-transported by diffeomorphisms.We demonstrate an analytical example of Corollary 2 by considering point vortices with strengths  1 , . . .,   located at p 1 , . . ., p  in the plane exterior to a circular obstacle   of radius  centered at the origin.Note that p  's move over time while   's remain constant.Using the method of images, the vortex dynamics can be represented without the obstacle by introducing mirrored vortices reflected in the disk.These mirrored vortices have strengths   = −  and positions p  =  2 |p  | 2 p  .The domain R 2 \   admits a harmonic vector field given by the clockwise 90 • rotated gradient of  (x) = 1 2 log |x|.The corresponding vortex-harmonic streamline linking number is Link =  =1    (p  ) =  =1   2 log |p  |. (85) The flux over any cross-section extending from the obstacle boundary to infinity can be computed by the difference of the streamfunction at the endpoints of the cross-section.In this case, the flux is the negative streamfunction value − (p) = −  =1 (   (p, p  ) +    (p, p  )) at any p ∈   , where  (x, y) = − 1 2 log |x − y| is the R 2 -Green's function.Using the geometric identity that the ratio |p−p  | /|p−p  | is the constant |p  | /, we obtain Flux = − (p) =  =1   2 log |p − p  | −   2 log |p − p  | show that ker(Flux) = im( + ).For notation convenience, we write η = ★ for each  ∈ V 1 to absorb the Hodge star.Let us also call Ω  D () = { ∈ Ω  () |  *  = 0} the space of Dirichletcondition-satisfying -forms.Note that the Dirichlet de Rham com- ∫  η = ∫   ψ = ∫  ψ = ∫   * ψ = 0. Now letus show the other direction ("⇒") using linear algebraic techniques.The statement we want to show is that if η ∈ Ω −1 D () satyisfies ∫  η = 0 for all cycles  ∈  −1 (, ), then η is exact in the Dirichlet de Rham complex.To see this, inspect the following diagram in which  and  are adjoint of each other  −1 (, )

Table 2 .
Same as Table 1 but on a 3D domain.