Commit ffe8de84 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

update drop2d for temporal convergence

parent f1a866a4
......@@ -4,23 +4,21 @@ import pyfefp
class Hydro :
def __init__(self, mesh, dt, r, rhop, rho, mu, g) :
def __init__(self, mesh, dt, r, rhop = 1, rho = 1, mu = 1, g = [0, -1]) :
self._jac = pyfefp.MeshJacobian(mesh, 2)
self._pf = pyfefp.ParticleFluid(self._jac)
self._dt = dt
self.g = np.array(g)
self.rho = rho
self.rhop = rhop
#todo add vmesh in law
self.law = pyfefp.Brinkman(
law = pyfefp.Brinkman(
mu = mu,
pf = self._pf,
ns = False,
g = self.g,
g = self.g, #np.array([0, 0]),
rho = rho,
rhop = rhop,
stab = 0
)
stab = dt)
d = pyfefp.DofManager(mesh)
d.add_field(pyfefp.TriangleP1Mini)
d.add_field(pyfefp.TriangleP1Mini)
......@@ -31,25 +29,23 @@ class Hydro :
d.add_boundary_condition("Top", 1, 0.)
d.add_boundary_condition("Top", 2, 0.)
d.complete()
self.rhop = rhop
self._dof = d
self._sol = d.new_vector()
self._solver = pyfefp.Steady(self.law, d, self._jac)
self._solverex = pyfefp.Steady(self.law, d, self._jac)
self.dt = dt
def set_particles(self, pos, vol, vel, forces, mass) :
self._pf.set_particles(pos, vol, vel, forces, mass)
self.law = law
self._sol = d.new_vector()
self._solver = pyfefp.Steady(law, d, self._jac, ilu=False)
self._solverEx = pyfefp.Steady(law, d, self._jac, ilu=False)
def solve(self) :
tic = time.clock()
v0 = self._pf.velocity.copy()
#self.law.stab = self.dt
#self._solver.solve()
#self._pf.velocity[:] += self._pf.forces * self.dt/self._pf.mass
self.law.stab = self.dt *0.5
self._sol[:] = self._solver.solve()
self._pf.velocity += self._pf.forces * self.dt*0.5 / self._pf.mass
self.law.stab = 0
self._sol[:] = self._solverex.solve()
self._sol[:] = self._solverEx.solve()
self._pf.velocity[:] = v0
#self._pf.forces[:] += self._pf.volume * self.rhop * self.g
print("cpu : ", time.clock() - tic)
def write_solution(self, odir, oiter, time, filetype="msh") :
......
......@@ -41,32 +41,41 @@ def genInitialPosition(p, r, rout, rhop) :
for i in range(x.shape[0]) :
p.addParticle((x[i], y[i]), r, r**2 * np.pi * rhop);
dt = 0.005
tEnd = 10
dt = 0.025
tEnd = 15
outputdir = "output_run_%g" % dt
outputdir = "output/"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
p = scontact2.ParticleProblem()
rhop=1100
genInitialPosition(p, 0.005, 0.4, rhop)
loadMeshBoundary(p, "mesh.msh", ["Top", "Box"])
p.position()[:, 1] += 2.5
d = 0.005
mu = 0.001
rho = 1000
ii = 0
if ii == 0 :
genInitialPosition(p, d, 0.4, rhop)
loadMeshBoundary(p, "mesh.msh", ["Top", "Box"])
p.position()[:, 1] += 2.5
p.velocity()[:, :] = 0;
else :
p.load(outputdir+"/part-%05d"%ii)
mass = p.particles()[:, [1]]
g = [0, -9.81]
r = np.mean(p.particles()[:, 0])
volume = np.pi * p.particles()[:, [0]]**2 * 2./3 ## 2./3 = 2d-3d correction
print("r = %g, m = %g\n" % (r, p.particles()[0, 1]))
print("dtmax = %g Re_p/u = %g" % (p.particles()[0, 1] / (150 * mu * d), d * rho * 0.3/mu))
#rhop = np.mean(p.particles()[:, 1] / volume[:, 0])
print("RHOP = %g" % rhop)
outf = 1
m = Mesh("mesh.msh")
mu = 1
rho = 1000
hydro = Hydro(m, dt, r, rhop=rhop, rho = rho, g=g, mu = mu)
t = 0
ii = 0
outputname = "%s/part-%05i" % (outputdir, 0)
p.velocity()[:, :] = 0;
p.write(outputname)
p.r = r
p.dt = dt
......@@ -74,23 +83,24 @@ writep(p, outputdir, 0)
pileseg = np.where(p.segmentTag() == 3)[0]
piledisk = np.where(p.diskTag() == 3)[0]
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces(), p.particles()[:, [1]])
hydro._pf.set_particles(p.position(), volume, p.velocity(), p.externalForces(), mass)
hydro.write_solution(outputdir, 0, 0, filetype="msh")
tic = time.clock()
v0 = 0
while t < tEnd :
hydro.solve()
p.velocity()[:] += p.externalForces() * dt/ p.particles()[:, [1]]
nsub = 1
nsub = 20
vn = p.velocity() + p.externalForces() * dt / mass
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 2)/r)))
print("NSUB", nsub, "VMAX * dt", vmax * dt, "r", r)
for i in range(nsub) :
p.solve(dt/nsub, 1e-5)
p.position()[:] += p.velocity() * dt/nsub
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces(), p.particles()[:, [1]])
p.iterate(r, dt/nsub, 1e-5)
hydro._pf.set_particles(p.position(), volume, p.velocity(), p.externalForces(), mass)
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
writep(p, outputdir, ioutput)
hydro.write_solution(outputdir, ioutput, t, filetype="msh")
ii += 1
print("%i : %.2g/%.2g" % (ii, t, tEnd))
print("%i : %.2g/%.2g (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment