depot-fluid.py

#!/usr/bin/env python
import pyfefp
from pyfefp.legendre import *
from pyfefp import  scontact2Interface
import scontact2
import numpy as np
from pyfefp import mesh
import os
import time
import random

def genInitialPosition(filename, r, ox, oy, lx, ly, rhop) :
    p = scontact2.ParticleProblem()
    scontact2Interface.MeshLoader(p, "mesh.msh", ["Top", "Box"])
    x = np.arange(ox+r, ox+lx, 2*r)
    y = np.arange(oy+r, oy+ly, 2*r)
    for i in range(x.shape[0]) :
       for j in range(y.shape[0]) :
          rr=r*(0.8+random.random()*0.2)
          p.addParticle((x[i], y[j]), rr, rr**2 * np.pi * rhop);
    p.write(filename)


outputdir = "output-initpressure"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0

r=1e-4*3
p = scontact2.ParticleProblem()

#physical parameters
g =  [0, -9.81]
rho = 1200
rhop = 1500 #2400
nu = 1e-6
tEnd = 20

#numerical parameters
#h = 2e-4 # approx r*100 but should match the mesh size

c =  20 #h/dt* 7./50.
dt = 1e-4 #0.5*r/V
epsilon = 1e-5 #2e-2*c*h # ?? not sure ??

print("c = %g, eps = %g" % (c, epsilon))

y0part = 0.051
lypart = 0.025
lxpart = 0.015
genInitialPosition(outputdir + "/part-00000", r, 0., y0part, lxpart, lypart, rhop)
p = scontact2Interface.ScontactInterface(outputdir+"/part-00000")

deltamass = p.volume().sum()*(rhop-rho)
deltapinit = deltamass*-g[1]/lxpart

print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 100 #000

fluid = pyfefp.FluidSolver("mesh.msh", dt, rhop=rhop, rho=rho, g = g, mu=nu*rho, c=c*c, epsilon = epsilon,timeScheme="Explicit", dragFactor=0.3, TypeEMu=TriangleP1, TypeEMp=TriangleP1)

# impermeability
fluid.add_boundary_condition("Top", 1, 0.)
fluid.add_boundary_condition("Top", 0, 0.)
fluid.add_boundary_condition("Box", 1, 0.)
fluid.add_boundary_condition("Box", 0, 0.)

fluid.complete()
y = fluid.mesh_position()[...,1].transpose()
pinit = (y0part+lypart-y)*deltapinit/lypart
pinit[pinit<0] = 0
pinit[pinit>deltapinit] = deltapinit
fluid._dof.setField(2, pinit, fluid.solution())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

tic = time.clock()
while t < tEnd :
    forces = fluid.solve() + g * p.mass()
    NS = 10
    if ii%NS == 0 :
        p.iterate(dt*NS, forces)
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        #p.write_vtk(outputdir, ioutput, t)
        p.write(outputdir, ioutput, t)
        #p.write_boundary_vtk(outputdir, ioutput, t)
        fluid.write_solution(outputdir, ioutput, t, "msh")
        #fluid.write_solution(outputdir, ioutput, t, "vtk")
    ii += 1
    print("%i : %.2g/%.2g (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))