# How to simulate water in 2D physics with codea?

Looking for bestpractice advices on how to simulate water for interaction with contaos physical bodies.
Eg. Material point method for fluid simulation
Would like to implement the following with Lua:
This was originally implemented in java by http://grantkot.com/MPM/Liquid.html

Tried to port this script from processing to Lua, but failed dramatically.
Here the work from Golan Levin that I tried to port.
https://github.com/golanlevin/MPM-Fluid

The main porting work is done, but there are some caveats in the logic.
If one is interested to help i’ll try to post the code here (is there some way to share projects via this forum?)
Should be useful for all of us, so maybe we can put our brains together to solve this.

Edited main post

having that interacting with gravity would be pretty cool. I wonder, can you just use the physics api and lots of tiny bodies to get a good effect?

Sorry, I guess not… Too slow. It might need some kind of particles logic handling a hole grid of bodies but not sure how to achieve this within the physics engine.

Yeah, interacting with gravity should rock!

Hi @lttlfas,

I had trouble porting code from processing and I think it was because it was using a library. When the library was removed on processing I got the same result as with my own code in Codea.

Try posting your code - there may be something obvious, but check the library situation first. If you could point us to the processing code that may help too.

Bri_G

Hi @lttlfas,

Ignore my last comment found the processing code.

Bri_G

On your code the fps is?

@Bri_G
Didnt use any other library. But I had to fake the map() and millis() methods from processing.
Maybe thats the point do following functions make sense in Lua, are they equivalent to their respective processing methods?

``````function millis()
return ElapsedTime * 1000
end

function map(value,low1,high1,low2,high2)
if (value < low1) or (value > high1) then
return nil
end
local srcMax = high1 - low1
local dstMax = high2 - low2
local adjValue = value - low1
return (adjValue * dstMax / srcMax) + low2
end

``````

@Connorbot999
Difficult to say, the particles grid (table) is called out of bounds after half a second and throws an error

How would you check for framerate? Is that the right way to get it:

``````--setup
time_start=ElapsedTime
--draw
FPS = frames / ( ElapsedTime - time_start)
frames=frames+1
``````

Wow, will have to have a look on how to pretty print code on the forums, sorry
@Bri_G
Here’s the code for now, but i guess it’s useless in this form. Still looking for how to print the code in a more readable form.

–Here the code: init settings and grid

``````MPM = class()

function MPM:init()

self.settings = {}
self.settings[0]=2.5 --"Density" 2.5
self.settings[1]=0.5 --"Stiffness" 0.5
self.settings[2]=3.0--"Bulk Viscosity" 3.0
self.settings[3]=1 --"Elasticity" 1.0
self.settings[4]=1 --"Viscosity" 1.0
self.settings[5]=1 --"Yield Rate" 1.0
self.settings[6]=0.02 --"Gravity" //0.02 to 0.004
self.settings[7]=1 --"Smoothing" 1.0

self.particles = {}
self.particlesX=0
self.particlesY=0
self.nParticles=0

self.grid={} -- with [][]
self.activeNodeArray = {}
self.nActiveNodes=0

self.gsizeX=0
self.gsizeY=0

self.prevMousePressed = false
self.prevMouseX = 0
self.prevMouseY = 0
self.mode = -1
self.SC = 10 --scaling
self.timer=0
self.np =10
self.steps=0
self:setup()
end

function MPM:setup()
self.width=WIDTH
self.height=HEIGHT

self.np =15
self:initLiquid(1+math.floor(self.width/self.SC), 1+math.floor(self.height/self.SC), self.np,self.np)
self:initGridAndParticles()
end

function MPM:draw()
background(35, 45, 56, 255)
self:simulate()
self:render()
end

function MPM:initLiquid(ingsizeX, ingsizeY, inparticlesX, inparticlesY)
self.gsizeX     = ingsizeX
self.gsizeY     = ingsizeY
self.particlesX = inparticlesX
self.particlesY = inparticlesY

self.prevMousePressed = false
self.prevMouseX = 0
self.prevMouseY = 0
end

function MPM:initGridAndParticles()
self.grid = {} --new Node[gsizeX][gsizeY]
for i=0,self.gsizeX do
self.grid[i] = {}
for j=0,self.gsizeY do
self.grid[i][j] = Node()
end
end

self.activeNodeArray = {} --new Node[gsizeX * gsizeY];
self.nActiveNodes = 0

self.nParticles = self.particlesX * self.particlesY
self.particles = {} --new Particle [nParticles];
local np = 0
local fi = self.gsizeX / (self.particlesX+5)
for i=0,self.particlesX do
for j=0,self.particlesY do
local px = (i+1)*fi
local py = map(j, 0,self.particlesY-1, 2,self.gsizeY/3)
local p = Particle(px, py, 0.0, 0.0)
self.particles[np] = p
np = np + 1
end
end
end

function MPM:render()
noSmooth()
--stroke(0, 255, 239, 255)
--beginShape(POINTS)
fill(0, 211, 255, 255)
for ip=0,self.nParticles do
local p = self.particles[ip]
ellipse(self.SC*p.x, self.SC*p.y,10,10) --vertex()
end
--endShape()
end
``````

– simulate function part 1

``````function MPM:simulate ()

local mouseX=CurrentTouch.x
local mouseY=CurrentTouch.y
local drag = false
local mdx = 0.0
local mdy = 0.0
local mx = mouseX/self.SC
local my = mouseY/self.SC

local phi=0
local mousePressed=(not CurrentTouch.state == ENDED)

if mousePressed and self.prevMousePressed then
drag = true
mdx = (mouseX - self.prevMouseX)/self.SC
mdy = (mouseY - self.prevMouseY)/self.SC
end

self.prevMousePressed = mousePressed
self.prevMouseX = mouseX
self.prevMouseY = mouseY

for i=0,self.gsizeX do
for j=0,self.gsizeY do
self.grid[i][j]:clear()
end
end
self.nActiveNodes = 0

--Particles pass 1
local pcxTmp=0
local pcyTmp=0
for ip=0,self.nParticles do

local p = self.particles[ip]

pcxTmp = math.floor(p.x - 0.5)
pcxTmp = math.min(self.gsizeX - 4, math.max(0, pcxTmp))
local pcx = pcxTmp
p.cx = pcxTmp

pcyTmp = math.floor(p.y - 0.5)
pcyTmp = math.min(self.gsizeY - 4, math.max(0, pcyTmp))
local pcy = pcyTmp
p.cy = pcyTmp

local px = p.px
local py = p.py
local gx = p.gx
local gy = p.gy
local pu = p.u
local pv = p.v

local x = p.cx - p.x
px[0] = 0.5 * x * x + 1.5 * x + 1.125
gx[0] = x + 1.5
x = x + 1
px[1] = -x * x + 0.75
gx[1] = -2 * x
x = x + 1
px[2] = (0.5 * x * x - 1.5 * x) + 1.125
gx[2] = x - 1.5

local y = p.cy - p.y
py[0] = 0.5 * y * y + 1.5 * y + 1.125
gy[0] = y + 1.5
y = y + 1
py[1] = -y * y + 0.75
gy[1] = -2 * y
y = y + 1
py[2] = (0.5 * y * y - 1.5 * y) + 1.125
gy[2] = y - 1.5

for i = 0,3 do
local nrow = self.grid[pcx + i] -- potential for array index out of bounds here if simulation explodes.

local pxi = px[i]
local gxi = gx[i]

for j = 0,3 do
local n = nrow[pcy + j] -- potential for array index out of bounds here if simulation explodes.

if not n.active then
n.active = true
self.activeNodeArray[self.nActiveNodes] = n
self.nActiveNodes = self.nActiveNodes + 1
end
phi   = pxi * py[j]
n.m  = n.m + phi
n.gx = n.gx +  gxi * py[j]
n.gy = n.gy +  pxi * gy[j]
n.u = n.u +  phi * pu
n.v = n.v +  phi * pv
end
end
end

for ni=0,self.nActiveNodes do
local n = self.activeNodeArray[ni]
if n and (n.m > 0) then
n.u = n.u / n.m
n.v = n.u / n.m
end
end
``````

– simulate function part 2

``````-- Particles pass 2
local densitySetting = self.settings[0]
local stiffness      = self.settings[1]
local bulkViscosity  = self.settings[2]
local elasticity     = self.settings[3]
local viscosity      = self.settings[4]
local yieldRate      = self.settings[5]

local T = millis()/1500
self.timer=T
for ip=0,self.nParticles do
local p = self.particles[ip]
local px = p.px
local py = p.py
local gx = p.gx
local gy = p.gy
local pcy = p.cy
local pcx = p.cx

local dudx = 0.0
local dudy = 0.0
local dvdx = 0.0
local dvdy = 0.0

local gxi=0
local pxi=0
local gxf=0
local gyf=0
for i = 0,3 do
local nrow = self.grid[pcx + i]
gxi = gx[i]
pxi = px[i]

local n0 = nrow[pcy]
gxf = gxi * py[0]
gyf = pxi * gy[0]
dudx = dudx +  n0.u * gxf
dudy = dudy +  n0.u * gyf
dvdx = dvdx +  n0.v * gxf
dvdy = dvdy +  n0.v * gyf

local n1 = nrow[pcy+1]
gxf = gxi * py[1]
gyf = pxi * gy[1]
dudx = dudx +  n1.u * gxf
dudy = dudy +  n1.u * gyf
dvdx = dvdx +  n1.v * gxf
dvdy = dvdy +  n1.v * gyf

local n2 = nrow[pcy+2]
gxf = gxi * py[2]
gyf = pxi * gy[2]
dudx = dudx +  n2.u * gxf
dudy = dudy +  n2.u * gyf
dvdx = dvdx +  n2.v * gxf
dvdy = dvdy +  n2.v * gyf
end

local w1  = dudy - dvdx
local wT0 = w1 * p.T01
local wT1 = 0.5 * w1 * (p.T00 - p.T11)
local D00 = dudx
local D01 = 0.5 * (dudy + dvdx)
local D11 = dvdy
local trace = 0.5 * (D00 + D11)
D00 = D00-trace
D11 = D11-trace

p.T00 = p.T00 + (-wT0 + D00) - yieldRate * p.T00
p.T01 = p.T01 + ( wT1 + D01) - yieldRate * p.T01
p.T11 = p.T11 + ( wT0 + D11) - yieldRate * p.T11
local norma = p.T00 * p.T00 + 2.0 * p.T01 * p.T01 + p.T11 * p.T11
if norma > 5.0 then
p.T00 = 0.0
p.T01 = 0.0
p.T11 = 0.0
end

local cx = math.abs(math.floor(p.x))
-- if cx>self.gsizeX then cx=5 end
local cy = math.abs(math.floor(p.y))
--  if cy>self.gsizeY then cy=5 end
--  cx=50
--   cy=50
-- test3=self.gsizeX.." / "..self.gsizeY.." / "..cx.." / "..cy

--test="uups "..cx.." - "..cy.." timer:"..self.timer.." steps:"..self.steps

local cxi = cx + 1
local cyi = cy + 1

local n00  = self.grid[cx][cy]
local n01  = self.grid[cx][cyi]
local n10  = self.grid[cxi][cy]
local n11  = self.grid[cxi][cyi]

local p00 = n00.m
local x00 = n00.gx
local y00 = n00.gy
local p01 = n01.m
local x01 = n01.gx
local y01 = n01.gy
local p10 = n10.m
local x10 = n10.gx
local y10 = n10.gy
local p11 = n11.m
local x11 = n11.gx
local y11 = n11.gy

local pdx = p10 - p00
local pdy = p01 - p00
local C20 =  3.0 * pdx - x10 - 2.0 * x00
local C02 =  3.0 * pdy - y01 - 2.0 * y00
local C30 = -2.0 * pdx + x10 + x00
local C03 = -2.0 * pdy + y01 + y00
local csum1 = p00 + y00 + C02 + C03
local csum2 = p00 + x00 + C20 + C30
local C21 =   3.0 * p11 - 2.0 * x01 - x11 - 3.0 * csum1  - C20
local C31 = (-2.0 * p11 +        x01 + x11 + 2.0 * csum1) - C30
local C12 =   3.0 * p11 - 2.0 * y10 - y11 - 3.0 * csum2  - C02
local C13 = (-2.0 * p11 +        y10 + y11 + 2.0 * csum2) - C03
local C11 = x01 - C13 - C12 - x00

local u  = p.x - cx*1.00 --float cx
local u2 = u * u
local u3 = u * u2
local v  = p.y - cy*1.00 -- float cy
local v2 = v * v
local v3 = v * v2
local density = p00 + x00 * u + y00 * v + C20 * u2 + C02 * v2 + C30 * u3 + C03 * v3 + C21 * u2 * v  + C31 * u3 * v + C12 * u  * v2 + C13 * u  * v3 + C11 * u  * v

local DS = densitySetting
--[[ spatially varying density field:
local DS = densitySetting * ( math.pow(  (p.y / self.gsizeY), 3.0) )
]]--

--if not DS then DS=2.5 end

local pressure = (stiffness / math.max(1.0, DS)) * (density - DS)

if pressure > 2.0 then
pressure = 2.0
end

local fx = 0.0
local fy = 0.0
if p.x < 3.0 then
fx = fx +  3.0 - p.x
elseif p.x > (self.gsizeX - 4) then
fx = fx +  (self.gsizeX - 4) - p.x
end

if p.y < 3 then
fy = fy +  3 - p.y
elseif p.y > (self.gsizeY - 4) then
fy = fy +  (self.gsizeY - 4) - p.y
end

trace = trace * self.settings[1]
local T00 = elasticity * p.T00 + viscosity * D00 + pressure + bulkViscosity * trace
local T01 = elasticity * p.T01 + viscosity * D01
local T11 = elasticity * p.T11 + viscosity * D11 + pressure + bulkViscosity * trace

local dx=nil
local dy=nil
for i=0,3 do
local nrow = self.grid[pcx + i]
local ppxi = px[i]
local pgxi = gx[i]

local n0 = nrow[pcy]
phi = ppxi * py[0]
dx  = pgxi * py[0]
dy  = ppxi * gy[0]
n0.ax = n0.ax +  fx * phi -(dx * T00 + dy * T01)
n0.ay = n0.ay + fy * phi -(dx * T01 + dy * T11)

local n1 = nrow[pcy + 1]
phi = ppxi * py[1]
dx  = pgxi * py[1]
dy  = ppxi * gy[1]
n1.ax = n1.ax + fx * phi -(dx * T00 + dy * T01)
n1.ay = n1.ay + fy * phi -(dx * T01 + dy * T11)

local n2 = nrow[pcy + 2]
phi = ppxi * py[2]
dx  = pgxi * py[2]
dy  = ppxi * gy[2]
n2.ax = n2.ax + fx * phi -(dx * T00 + dy * T01)
n2.ay = n2.ay + fy * phi -(dx * T01 + dy * T11)
end
end

for ni=0,self.nActiveNodes do
local n = self.activeNodeArray[ni]
if n and (n.m > 0.0) then
n.ax = n.ax/n.m
n.ay = n.ay/n.m
n.u = 0.0
n.v = 0.0
end
end
--test=self.settings[6]
``````

– simulate function part 3

``````  -- Particles pass 3
local gravity =self.settings[6]
for ip=0,self.nParticles do
local p = self.particles[ip]

local px = p.px -- with []
local py = p.py -- with []
local pcy = p.cy
local pcx = p.cx
for i=0,3 do
local nrow = self.grid[pcx + i]
local ppxi = px[i]

local n0 = nrow[pcy]
phi = ppxi * py[0]
p.u = p.u +  phi * n0.ax
p.v = p.v +  phi * n0.ay

local n1 = nrow[pcy + 1]
phi = ppxi * py[1]
p.u = p.u + phi * n1.ax
p.v = p.v + phi * n1.ay

local n2 = nrow[pcy + 2]
phi = ppxi * py[2]
p.u = p.u + phi * n2.ax
p.v = p.v + phi * n2.ay
end

p.v = p.v + gravity
if drag then
local vx = abs(p.x - mx)
local vy = abs(p.y - my)
if (vx < 10) and (vy < 10) then
local weight = (1.0 - vx / 10) * (1.0 - vy / 10)
p.u = p.u + weight * (mdx - p.u)
p.v = p.v + weight * (mdy - p.v)
end
end

local xf = p.x + p.u
local yf = p.y + p.v
if xf < 2.0 then
p.u = p.u + (2.0 - xf) + math.random()*0.01
elseif (xf > (self.gsizeX - 3.0)) then
p.u = p.u + (self.gsizeX - 3.0) - xf - math.random()*0.01
end
if (yf < 2.0) then
p.v = p.v + (2.0 - yf) + math.random()*0.01
elseif (yf > (self.gsizeY - 3.0)) then
p.v = p.v + (self.gsizeY - 3.0) - yf - math.random()*0.01
end

local pu = p.u
local pv = p.v
for i=0,3 do
local nrow = self.grid[pcx + i]
local ppxi = px[i]

local n0 = nrow[pcy]
phi = ppxi * py[0]
n0.u = n0.u + phi * pu
n0.v = n0.v + phi * pv

local n1 = nrow[pcy + 1]
phi = ppxi * py[1]
n1.u = n1.u + phi * pu
n1.v = n1.v + phi * pv

local n2 = nrow[pcy + 2]
phi = ppxi * py[2]
n2.u = n2.u + phi * pu
n2.v = n2.v + phi * pv
end
end

for ni=0,self.nActiveNodes do
local n = self.activeNodeArray[ni]
if n and (n.m > 0.0) then
n.u = n.u/n.m
n.v = n.v/n.m
end
end

--Particles pass 4
local gu=nil
local gv=nil
local smoothing = 1.00--self.settings[7]
for ip=0,self.nParticles do
local p = self.particles[ip]

gu = 0.0
gv = 0.0

local px = p.px
local py = p.py
local pcy = p.cy
local pcx = p.cx
for i=0,3 do
local nrow = self.grid[pcx + i]
local ppxi = px[i]

local n0 = nrow[pcy]
phi = ppxi * py[0]
gu = gu + phi * n0.u
gv = gv + phi * n0.v

local n1 = nrow[pcy + 1]
phi = ppxi * py[1]
gu = gu + phi * n1.u
gv = gv + phi * n1.v

local n2 = nrow[pcy + 2]
phi = ppxi * py[2]
gu = gu + phi * n2.u
gv = gv + phi * n2.v
end

p.gu = gu
p.gv = gv
p.x = p.x + gu
p.y = p.y + gv

p.u = p.u + smoothing * (gu - p.u)
p.v = p.v + smoothing * (gv - p.v)
self.steps = self.steps + 1
--    test2=p.x.." - "..p.y.." steps:"..self.steps

--       if p.x<0 then p.x=0 end
--  if p.x>WIDTH/self.SC-4 then p.x=WIDTH/self.SC-5 end
--   if p.y<0 then p.y=0 end
--   if p.y>HEIGHT/self.SC-4 then p.y=HEIGHT/self.SC-4 end
end

end
``````

– The particle and node classes

``````Particle = class()

function Particle:init(inx, iny, inu, inv)
if not inx then inx=0 end
if not iny then iny=0 end
if not inu then inu=0 end
if not inv then inv=0 end

self.x = inx
self.y = iny
self.u = inu
self.v = inv

self.cx = 0
self.cy = 0

self.gu=0
self.gv=0
self.T00=0
self.T01=0
self.T11=0

self.px={0,0,0}
self.py={0,0,0}
self.gx={0,0,0}
self.gy={0,0,0}
end

Node = class()

function Node:init()
self.m  = 0
self.d  = 0
self.gx = 0
self.gy = 0
self.u  = 0
self.v  = 0
self.ax = 0
self.ay = 0
self.active = false
end

function Node:clear()
self.m = 0
self.d = 0
self.gx = 0
self.gy = 0
self.u = 0
self.v = 0
self.ax = 0
self.ay = 0
self.active = false
end

``````

–and finally the main class

``````function setup()
displayMode(FULLSCREEN)
liquid=MPM()
end

function draw()
liquid:draw()
end
``````

@connorbot, thanks, the 3 tilde kept code formatting

Is not running the error is

``````error: error: [string "MPM = class()..."]:86: attempt to call global 'map' (a nil value)

Pausing playback
``````

@connorbot
You will have to include the map() and millis() functions from above

Notice, even then, this is a broken implementation. However, changing the settings makes the simulation crash more or less faster

@ruilov
Here is a quick and dirty test with lots of tiny bodies (there is a weird ?bug? though after few seconds, breaking the world bounds)
With 100 bodies it’s ok, with 200 it becomes slow, with 300 it fails

``````
function setup()
maxitems=100
drops={}
supportedOrientations(LANDSCAPE_LEFT)
displayMode(FULLSCREEN)
local x=WIDTH/2
local y=HEIGHT/2
for i=0,maxitems do
end
boundingbox()
end

function boundingbox()
ground = physics.body(POLYGON, vec2(0,1), vec2(0,0), vec2(WIDTH,0), vec2(WIDTH,1))
ground.type = STATIC
groundleft = physics.body(POLYGON, vec2(1,0), vec2(0,0), vec2(0,HEIGHT), vec2(1,HEIGHT))
groundleft.type = STATIC
groundright = physics.body(POLYGON, vec2(WIDTH - 1,0), vec2(WIDTH,0), vec2(WIDTH,HEIGHT), vec2(WIDTH - 1,HEIGHT))
groundright.type = STATIC
ceiling = physics.body(POLYGON, vec2(0,HEIGHT-1), vec2(0,HEIGHT), vec2(WIDTH,HEIGHT), vec2(WIDTH,HEIGHT-1))
ceiling.type = STATIC
end

local r=5
local circle = physics.body(CIRCLE, r)
circle.gravityScale = 10+math.random()
circle.interpolate = true
circle.x = x
circle.y = y
circle.restitution = 0.7
circle.sleepingAllowed = false
return circle
end

function draw()
local g=vec3(Gravity.x,Gravity.y,Gravity.z)
physics.gravity(g)
background(40, 40, 50)
noStroke()
for i,body in ipairs(drops) do
fill(0, 200, 255, 126)
ellipse(body.x,body.y,5)
end
end

``````

@Bri_G
Many thanks for the bug fix, updated the above test code.
Maybe the garbage collector is the issue in the MPM class, I’ll have a look later

Hi @littfass,

Just tried out your small tester, floor falls away after several cycles. Resolved that by removing local definition on your floor, left, right and ceiling. By that I mean delete the preceding word ‘local’ on each line. Then it runs fine.

Hope that helps.

Bri_G