Skip to content

Transport_fluid

Bruno Levy edited this page Sep 11, 2022 · 13 revisions

Simulation of incompressible fluids with Optimal Transport

Prerequisites

Before starting, if you do not have it already, you will need to install the Graphite software (and its WarpDrive plugin, that is included in the distribution):

  • if you are under Windows, install pre-compiled package (instructions here). The WarpDrive plugin is included.
  • if you are under Linux or MacOS, install Graphite from sources (instructions here). You will need to compile the WarpDrive plugin.

Programming with Graphite

In this tutorial, we will use the LUA programming language. LUA is a scripting language similar to Python, if you know Python already you will (nearly) find yourself at home. More information about LUA here.

Why LUA rather than Python ?

  • LUA is smaller/faster/more elegantly designed (in my own subjective opinion). A more objective measure: it has a super small footprint (in terms of source size, package size). It is directly included in Graphite (nothing to install).

  • If you really prefer Python (that is much more popular than LUA and that has many existing packages), there will be a Python version of this tutorial soon (Graphite also supports Python through the "gompy" plugin).

Now you can take a quick look at the Scripting Graphite tutorial to see how the editor works (see in particular automatic completion and help bubbles that can save you some time).

Step 1: a couple of particles

We will start by a super-simple simulatoin with a couple of particles, then add complexity, one element at a time.

First law of Newton

We first create a couple of objects and attributes, as follows:

N=10
scene_graph.clear()
Omega = scene_graph.create_object('OGF::MeshGrob','Omega')
Omega.I.Shapes.create_square()

points = scene_graph.create_object('OGF::MeshGrob','points')
scene_graph.current_object = 'points'

E = points.I.Editor
point    = E.find_attribute('vertices.point')
mass     = E.find_or_create_attribute('vertices.mass')
V        = E.find_or_create_attribute(
   {attribute_name='vertices.speed',dimension=2}
)
vertex_id = E.find_or_create_attribute('vertices.id')
points.shader.vertices_style='true; 0 1 0 1; 5'
  • the domain Omega of the simulation
  • the pointset points with the position (point), mass, speed vector V and id vertex_id used for visualization (display attribute so that each point as a different color).

Then we initialize the points, at random positions and with random initial speeds:

for v=0,N-1 do
   local x = 0.1 + 0.8*math.random()
   local y = 0.1 + 0.8*math.random()
   E.create_vertex(x,y)
   mass[v] = 1
   V[2*v]   = 3.0 * (math.random()-0.5)
   V[2*v+1] = 3.0 * (math.random()-0.5)
   vertex_id[v] = v
end

OK, so we can see our points, but it would be good to visualize the speed vectors as well. To do that, we create another mesh, and use its edges to visualize the speed vectors, as follows:

speeds_display=nil
function show_speeds()
   if speeds_display == nil then
      speeds_display = scene_graph.create_object('OGF::MeshGrob','speeds')
      E_speeds_display = speeds_display.I.Editor
      speeds_display_points = E_speeds_display.find_attribute('vertices.point')
   end
   speeds_display.clear()
   E_speeds_display.create_vertices(2*N)
   local off = 3*N
   local scale = 0.075
   for v=0,N-1 do
      local x = point[3*v]
      local y = point[3*v+1]
      local Vx = V[2*v]
      local Vy = V[2*v+1]
      speeds_display_points[3*v]   = x
      speeds_display_points[3*v+1] = y
      speeds_display_points[3*v  +off] = x + scale * Vx
      speeds_display_points[3*v+1+off] = y + scale * Vy
      E_speeds_display.create_edge(v, v+N)
   end
   speeds_display.visible=true
   speeds_display.update()
end

show_speeds()

Note remember that in Graphite, all meshes are in 3D, so even if we do a 2D simulation, point coordinates of vertex i are x = point[3*v],y=point[3*v+1] (instead of 2*v, 2*v+1 as one may have thought)!

OK, so now let us move our points ! For now we have no force, so all the points will move along a straight line (first law of Newton). We create a new function to do one step of simulation:

function Euler_step()
   local tau = 0.005 
   for v = 0,E.nb_vertices-1 do
      point[3*v]   = point[3*v]   + tau*V[2*v]
      point[3*v+1] = point[3*v+1] + tau*V[2*v+1]
   end
   points.update()
   show_speeds()
end

The function simply makes each point move along its speed vector. The parameter tau is the timestep, in seconds.

Now we can create a new button to start our simulation:

Euler_dialog = {} 
Euler_dialog.visible = true
Euler_dialog.name = 'SimpleEuler' 
Euler_dialog.x = 100
Euler_dialog.y = 400
Euler_dialog.w = 150
Euler_dialog.h = 200
Euler_dialog.width = 400

function Euler_dialog.draw_window()
   if imgui.Button('TimeStep',-1,-1) then
      Euler_step()
   end
end

graphite_main_window.add_module(Euler_dialog)

Run the program (<F5>), each time you push the TimeStep button, each point moves along its speed vector, not super interesting for now. OK, to make things more interesting, let us add gravity.

Second law of Newton

Newton's second law states that F = ma, where F is the force applied to a particle, m its mass and a the acceleration. Put differently, it says that a = (1/m)*F, which means that when we move from t to t + tau, the speed V should be updated as V = V + tau*a, or V = V + (tau/m) * F.

Now which force do we have ? For now we have gravity, oriented towards the negative Y axis, or F = [0, -m*G] where G = 9.81. Hence we update our Euler_step() function as follows:

function Euler_step()
   local tau = 0.005 -- Timestep
   local G=9.81      -- Gravity in m/s^2 

   for v = 0,E.nb_vertices-1 do
      local Fx = 0.0
      local Fy = - mass[v] * G
      
      -- V += tau * a ; F = ma ==> V += tau * F / m
      V[2*v]   = V[2*v]   + tau * Fx / mass[v]
      V[2*v+1] = V[2*v+1] + tau * Fy / mass[v]
 
      -- position += tau * V
      point[3*v]   = point[3*v]   + tau*V[2*v]
      point[3*v+1] = point[3*v+1] + tau*V[2*v+1]
   end
   show_speeds()
end

Note here we used the simplest method (called "semi-explicit Euler") to approximate F=ma with discrete timesteps (what is called "time integration"). More sophisticated methods exist (leapfrog, Verlet), with the advantage of better conserving quantities supposed to be conserved (energy, momentum).

OK so now when you push the button TimeStep multiple times, you see your points moving along nice parabolas, as expected, but soon the points move outside the domain. Could we make them bounce on the domain boundaries ? Yes, it is simple, if the x (resp. y) coordinate of the point is outside [0,1], change the sign of Vx (resp. Vy). In addition, whenever a collision with the domain boundary occurs, you can multiply the speed vector by a damping coefficient:

function bounce_on_borders(v)
   local damp = 0.9
   
   if point[3*v] > 1.0 then
      point[3*v] = 1.0
      V[2*v]   = -damp * V[2*v]
      V[2*v+1] =  damp * V[2*v+1]
   end
      
   if point[3*v] < 0.0 then
      point[3*v] = 0.0
      V[2*v]   = -damp * V[2*v]
      V[2*v+1] =  damp * V[2*v+1]
   end

   if point[3*v+1] > 1.0 then
      point[3*v+1] = 1.0
      V[2*v]   =  damp * V[2*v]
      V[2*v+1] = -damp * V[2*v+1]
   end
      
   if point[3*v+1] < 0.0 then
      point[3*v+1] = 0.0
      V[2*v]   =  damp * V[2*v]	  
      V[2*v+1] = -damp * V[2*v+1]
   end
end

Then you call bounce_on_borders() for each particle at the beginning of Euler_step():

function Euler_step()
   local tau = 0.005 -- Timestep
   local G=9.81      -- Gravity in m/s^2 

   for v = 0,E.nb_vertices-1 do
      bounce_on_borders(v)
   end

   ... update forces, speeds and positions
end

The complete program is available here. You can also load it from the Examples... menu of the text editor (Examples.../WarpDrive/advanced/Euler_01.lua).

Note pushing the TimeStep button multiple times is certainly painful. The example program above lets you specify in the graphic interface how many timesteps should be computed each time the button is pushed.

Step 2: collisions

Let us now handle the collisions between the particles. To do that, we will modify the speeds V and store them in a new array V_new that we need to create (right after creating V):

  V        = E.find_or_create_attribute(
     {attribute_name='vertices.speed',dimension=2}
  )
  V_new    = E.find_or_create_attribute(
     {attribute_name='vertices.speed',dimension=2}
  )

Then we need a function to detect whether two particles have collided. We consider that two particles collide if they are nearer than a certain distance and if one of them is moving towards the other one:

radius2 = 0.001
function choc(v1,v2)
    local x1 = point[3*v1]
    local y1 = point[3*v1+1]
    local x2 = point[3*v2]
    local y2 = point[3*v2+1]
    local dx = x2-x1
    local dy = y2-y1
    local dist = dx*dx+dy*dy
    if dist < radius2 then
         local vx1 = V[2*v1]
         local vy1 = V[2*v1+1]
         local vx2 = V[2*v1]
         local vy2 = V[2*v1+1]
         if(
             dx*vx1+dy*vy1 > 0.0 or
             dx*vx2+dy*vy2 < 0.0
         ) then
	         return true
         end
    end
    return false
end

Then the function that makes two particles bounce if there was a choc. There are several possibilities that conserve momentum and kinetic energy. We do that as follows (and in addition, we multiply the speeds by a samping coefficient):

function bounce(v1,v2)
   -- Damping of speed when there is a choc
   local damp = 0.99

   if not choc(v1,v2) then
      return
   end
 
   local m1 = mass[v1]
   local m2 = mass[v2]

   local Vgx = (m1*V[2*v1]   + m2*V[2*v2]  ) / (m1+m2)
   local Vgy = (m1*V[2*v1+1] + m2*V[2*v2+1]) / (m1+m2)

   local Vgnorm = math.sqrt(Vgx*Vgx+Vgy*Vgy)

   local Nx=-Vgy/Vgnorm
   local Ny= Vgx/Vgnorm

   local dot = V[2*v1]*Nx + V[2*v1+1]*Ny
   V_new[2*v1  ] = damp * (V[2*v1]   - 2*dot*Nx)
   V_new[2*v1+1] = damp * (V[2*v1+1] - 2*dot*Ny)

   dot = V[2*v2]*Nx + V[2*v2+1]*Ny
   V_new[2*v2  ] = damp * (V[2*v2]   - 2*dot*Nx)
   V_new[2*v2+1] = damp * (V[2*v2+1] - 2*dot*Ny)

end

Then we simply add the following code at the beginning of the Euler_step() function:

   NL.blas.copy(V, V_new) 
   for v1=0,E.nb_vertices-1 do
     for v2=v1+1,E.nb_vertices-1 do
        bounce(v1,v2)
     end
   end
   NL.blas.copy(V_new,V)

The complete program is available here. You can also load it from the Examples... menu of the text editor (Examples.../WarpDrive/advanced/Euler_02.lua).

It is by default initialized with two particles and specific speeds so that we are sure collisions will occur (for testing). Try other values of N.

For N equal to 500 or larger, simulation starts to slow down a lot ! It is because we test all couples of particiles for collisions (O(N^2)). Can we do better ?

Step 3: many particles (faster collisions)

To avoid testing all N^2 particle pairs for collision, it is possible to decompose the domain into grid cells, and restrict the search for collision to each grid cell (plus adjacent cells if you do not want to miss collisions that occur exactly on cell edges).

In each grid cell, we need to store the list of points contained by the cell. We store the points of a grid cell in a "linked list", stored in arrays: the cell knows the index of the first point contained by the cell, and a vertex attribute next indicates the next point. The end of the list is indicated by -1:

search_grid = {}
search_grid.nU = math.floor(math.sqrt(N))*5
search_grid.nV = search_grid.nU
search_grid.first = gom.create('OGF::NL::Vector')
search_grid.first.resize({
   new_size=search_grid.nU * search_grid.nV,
   new_dim=1,
   element_meta_type=gom.resolve_meta_type('int')}
)
next = E.find_or_create_attribute({
    attribute_name='vertices.next',dimension=1,
    type=gom.resolve_meta_type('int')
})

Then we declare a couple of utility function. This one converts two grid indices U,V into a linear index (for search_grid.first[]):

function search_grid.linear_index(U,V)
   return V*search_grid.nU + U
end

Then one can get the index of the first point in grid cell (U,V):

function search_grid.first_by_UV(U,V)
   if U < 0 or V < 0 or
      U >= search_grid.nU or
      V >= search_grid.nV
   then
      return -1
   end
   local i = search_grid.linear_index(U,V)
   return search_grid.first[i]
end

This function gets the U and V indices of the grid cell that contains a point (x,y):

function search_grid.uv(x,y)
   local U = math.min(math.floor(x*search_grid.nU),search_grid.nU-1)
   local V = math.min(math.floor(x*search_grid.nV),search_grid.nV-1)
   return U,V   
end 

And this one gets the U and V indices of the grid ell that contains the point of index v: function search_grid.v_to_UV(v) local x = point[3v] local y = point[3v+1] local U,V = search_grid.uv(x,y) return U,V end

This function inserts the point of index v into grid cell (U,V):

function search_grid.insert(U,V,v)
   local i = search_grid.linear_index(U,V)
   next[v] = search_grid.first[i]
   search_grid.first[i] = v
end

Then one can insert all the points as follows:

function search_grid.init()
   for i=0,search_grid.nU*search_grid.nV-1 do
      search_grid.first[i] = -1
   end
   for v=0,N-1 do
      local x = point[3*v]
      local y = point[3*v+1]
      local U,V = search_grid.uv(x,y)
      search_grid.insert(U,V,v)
   end
end

Now it is easy to test all the intersections, by updating the function Euler_step() as follows:

   ...
   search_grid.init()
   NL.blas.copy(V, V_new)
   if(Euler_dialog.collisions) then
      local nb_tests = 0.0
      for v1=0,E.nb_vertices-1 do
          local U1,V1 = search_grid.v_to_UV(v1)
          for U=U1-1,U1+1 do
            for V=V1-1,V1+1 do
               local v2 = search_grid.first_by_UV(U,V)
               while v2 ~= -1 do
	              nb_tests = nb_tests+1
                  bounce(v1,v2)
                  v2 = next[v2]
               end
            end
          end
       end
    end
    NL.blas.copy(V_new,V)
    ...

The complete program is available here. You can also load it from the Examples... menu of the text editor (Examples.../WarpDrive/advanced/Euler_03.lua).

The new program is significantly faster than the previous one, and still runs at interactive speed with 1000 particles. What we can observe now corresponds to a gas. There are two different things we need to do to simulate a liquid:

  • liquids are incompressible, and we need to add new constraints for that
  • here we simulated each particle (or molecule) individually. If we want a realistic simulation, we need to simulate many many molecules, it is going to cost too much. So we need to "make a step backward" to have a more general view of the fluid, and simulate larger scale quantities.

This is where optimal transport enters the scene.

Step 4: simulating incompressible liquids

So what we have now is a simulation of a gas, with all the individual particles. To simulate a liquid, we need two things:

  • liquids are incompressible, and we need to add a new constraint for that;
  • we simulated each molecule individually. If we want a realistic simulation, we will need to simulate many many molecules, and it is going to cost too much. So we need to "make a step backward" to have a wider view of the fluid, and simulate quantities at a coarser scale.

We want to simulate what happens when one puts a heavy liquid on top of a lighter one, has here. Then both liquids want to exchange positions, but incompressibility makes things more difficult for them. First, at a microscopic scale, individual molecules exchange positions, which starts creating tiny vortices. These tiny vortices become larger and larger, at a faster and faster speed (it is called the Rayleigh-Taylor Instability).

In simulation, it will look like the small animation above.

Let us see now how it works. We are going to decompose the fluid into elements. Each element will represent a set of molecules of the fluid, which means that the elements will be able to change shape, but not to change volume: they are not allowed to exchange molecules. But we said we are not going to represent individual molecules, so how does it work ? In fact consider you represent many many molecules in each element. At the limit, there is an infinite number of molecules per element. So you can replace the individual molecules with integrated quantities over the cells.

So how do we describe the cells ? Each cell will depend on a set of parameters. These parameters are symbolized as big points on the figure above. For instance, if it was a Voronoi diagram, the parameters would be the side of the Voronoi diagram (but in fact it is slightly more complicated). So what we have to do now is just:

  • integrate Newton law F = ma with respect to space over the cells and transfer them to the parameters, which means take the limit when the number of molecules per cell tends to infinity;
  • integrate Newton law F = ma with respect to time, with discrete time-steps.

So we need to see which forces we have:

  • there is gravity, this one is easy, just like before;
  • there is also pressure, that can be seen either as a force or as a constraint. Let us elaborate on that.

Brenier's point of view on incompressible fluids

On the figure above, we are in an abstract space, where each point represents a configuration of the fluid. In other words, each point is a "fluid configuration", with all the coordinates of all the molecules of the fluid (the drawing is in 2D, but in fact it is a super-high-dimensional space !!!).

In this drawing, the "blob" represents the set of fluid configurations that correspond to an incompressible motion (and here I mean the boundary of the blob, because this set of incompressible fluid configurations is of lower dimension than the set of all possible configurations). Hence a "valid" (that is, incompressible) fluid configuration is always on the boundary of the blob.

In this drawing, the "velocity" arrow represents all the velocities of all the molecules of the fluid (it is also of super high dimension). Since we want the fluid to remain incompressible, this means that the point that represents the fluid configuration is supposed to remain on the curve. Thus, at any time, the velocity should be tangent to the curve. Put differently, the fluid configuration follows a geodesic on the blob.

But we are going to move through discrete timesteps, which means that when moving between t and t + delta t, the fluid with "take off" a bit from the surface of the blob, so the idea is to reproject it onto the blob at each timesteps. How can we do that ? We need to find the orthogonal projection of any fluid configuration onto the set of incompressible configurations. Brenier's polar factorization theorem pdf relates this projection with optimal transport (this is where optimal transport enters the scene).

The Gallouet-Merigot numerical scheme

Projecting the fluid onto the blob at each timestep would be too "brutal", so the idea in Gallouet and Merigot's numerical scheme is to attach a "spring" between the fluid configuration (that may be outside the blob) and its projection onto the blob. Then the fluid configuration will oscillate around the blob (see figure above, still in abstract space).

What was said was rather abstract, how does it work in practice ? In practice, we start with a set of points (left figure above). Then, as shown in the center figure, we compute the unique Laguerre diagram that has prescribed areas for the cells, using semi-discrete optimal transport see this tutorial. In practice, the "projection onto the blob" corresponds to the centroids of the Laguerre cells. They are shown as white dots on the right figure. But instead of directly moving the points to the centroids of the Laguerre cell, we attach a "spring". It can be written as follows in Lua:

function Euler_step()
   local OT = points.query_interface('OGF::MeshGrobTransport')
 
   -- Timestep
   local tau = 0.001    

   -- Stiffness of the 'spring' pressure force that pulls the 
   -- points towards the centroids
   local epsilon = 0.004

   -- Gravity on earth in m/s^2
   local G = 9.81

   local inveps2 = 1.0/(epsilon*epsilon)

   -- Compute the centroids of the unique Laguerre diagram defined
   -- from the points that has constant areas.
   OT.compute_optimal_Laguerre_cells_centroids(
      {Omega=Omega,centroids=centroid,weights=weight,mode='EULER_2D'}
   )

   -- Update forces, speeds and positions (Explicit Euler scheme, super simple !)
   for v = 0,E.nb_vertices-1 do
      -- Compute forces: F = spring_force(point, centroid) - m G Z
      local Fx = inveps2 * (centroid[2*v]   - point[3*v]  )
      local Fy = inveps2 * (centroid[2*v+1] - point[3*v+1]) - mass[v] * G
      -- V += tau * a ; F = ma ==> V += tau * F / m
      V[2*v]   = V[2*v]   + tau * Fx / mass[v]
      V[2*v+1] = V[2*v+1] + tau * Fy / mass[v]
      -- position += tau * V
      point[3*v]   = point[3*v]   + tau*V[2*v]
      point[3*v+1] = point[3*v+1] + tau*V[2*v+1]
   end
end

The complete program is available here. You can also load it from the Examples... menu of the text editor (Examples.../WarpDrive/advanced/TaylorRaileigh_2d_in_lua.lua).

References

  • [1] Polar Factorizartion and Monotone Rearrangement of Vector-Valued Functions, Yann Brenier, Comm. on Pure and Applied Math. PDF.
  • [2] A Lagrangian scheme for the incompressible Euler equation using optimal transport, Thomas Gallouet and Quentin Merigot, Foundations of Computational Math. arXiv
  • [3] Partial Optimal Transport for a Constant-Volume Lagrangian Mesh with Free Boundaries, Bruno Lévy, J. of Computational Physics arXiv
  • [4] Power particles: an incompressible fluid solver based on power diagrams, F. de Goes, C. Wallez, J. Huang, D. Pavlov and M. Desbrun, ACM SIGGRAPH ACM DL PDF
  • [5] The power particle-in-cell method, Z. Qu, M. Li, F. de Goes and C. Jiang, ACM SIGGRAPH ACM DL