Commit 6f127124 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

update testcases

parent b79d70fa
......@@ -15,6 +15,7 @@ class FluidSolver :
self._pf = ParticleFluid(self._jac)
self._imex = imex
self._dt = dt
self._stabFactor = 1
law = Brinkman(
mu = mu,
pf = self._pf,
......@@ -54,19 +55,22 @@ class FluidSolver :
tic = time.clock()
if self._imex :
v0 = self._pf.velocity.copy()
self._law.stab = self._dt *0.5
self._law.stab = self._dt *0.5 * self._stabFactor
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._pf.velocity[:] = v0
else :
self._law.stab = self._dt
self._law.stab = self._dt * self._stabFactor
self._sol[:] = self._solver.solve()
print("cpu : ", time.clock() - tic)
return self.forces
def set_stab_factor(self, s) :
self._stabFactor = s
def write_position(self, odir, oiter, time, filetype="msh"):
if filetype == "msh" :
basename = "%s/f_%05i" % (odir, oiter)
......
......@@ -83,9 +83,11 @@ class LmgcInterface(object):
self._externalF= np.zeros([self._nbDisk,3], 'd')
self._volume = np.zeros([self._nbDisk,1], 'd')
self._mass = np.zeros([self._nbDisk,1], 'd')
self._r = np.zeros([self._nbDisk,1], 'd')
for i in range(self._nbDisk):
self._volume[i] = np.pi * chipy.DISKx_GetContactorRadius(i+1)**2 * 2./3
self._r[i] = chipy.DISKx_GetContactorRadius(i+1)
self._mass[i] = chipy.RBDY2_GetBodyMass(int(self._d2bMap[i,0]))
self._volume = np.pi * self._r**2 * 2./3
self._tag2bnd = {}
self._tag2id = {}
for i in range(chipy.RBDY2_GetNbRBDY2()):
......@@ -154,6 +156,10 @@ class LmgcInterface(object):
""" Return the mass of the particles """
return self._mass
def r(self):
""" Return the radius of the particles """
return self._r
def position(self):
""" Get current position of contactors and return it """
for i in range(self._nbDisk):
......
import numpy as np
import shutil
import scontact2
from . import mesh
class ScontactInterface :
def __init__(self, inputfile) :
import scontact2
p = scontact2.ParticleProblem()
self._p = p
p.load(inputfile)
p.velocity()[:, :] = 0;
self._r = np.mean(p.particles()[:, 0])
self._r = p.particles()[:, 0]
self._volume = np.pi * p.particles()[:, [0]]**2 * 2./3 ## 2./3 = 2d-3d correction
self._mass = p.particles()[:, [1]]
......@@ -17,6 +18,9 @@ class ScontactInterface :
def volume(self):
return self._volume
def r(self) :
return self._r
def mass(self):
return self._mass
......@@ -34,7 +38,7 @@ class ScontactInterface :
def iterate(self, dt, forces) :
self._p.externalForces()[:] = forces
self._p.iterate(self._r, dt, 1e-5)
self._p.iterate(np.min(self._r), dt, 1e-5)
def write(self, odir, i, t) :
outputname = "%s/part-%05i" % (odir, i)
......@@ -43,3 +47,28 @@ class ScontactInterface :
class MeshLoader(object) :
def __init__(self, p, filename, tags = None, shift = [0, 0]) :
self._p = p
self._m = mesh.Mesh(filename)
self._addv = set()
if tags is not None :
self.load(tags, shift)
def load(self, tags, shift= [0, 0]) :
pnum = [self._m.getPhysicalNumber(1, i) for i in tags]
for ent in self._m.entities :
if ent.dimension != 1 or not [i for i in pnum if i in ent.physicals]:
continue
tag = 0 if not ent.physicals else ent.physicals[0]
stag = self._m.getPhysicalName(ent.dimension, tag) or str(tag)
for i, el in enumerate(ent.elements) :
for v in el.vertices :
if v[3] in self._addv :
continue
self._addv.add(v[3])
self._p.addBoundaryDisk((v[0] + shift[0], v[1] + shift[1]), 0.000, stag)
v0 = el.vertices[1]
v1 = el.vertices[0]
self._p.addBoundarySegment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), stag)
import numpy as np
import time
import pyfefp
class Hydro :
def __init__(self, mesh, dt, r, rho = 1, mu = 1, g = [0, -1], rhop = 1) :
self._jac = pyfefp.MeshJacobian(mesh, 2)
self._pf = pyfefp.ParticleFluid(self._jac)
self._dt = dt
#todo add vmesh in law
law = pyfefp.Brinkman(
mu = mu,
pf = self._pf,
g = g,
ns = True,
rho = rho,
stab = dt/rho * 100
)
d = pyfefp.DofManager(mesh)
d.add_field(pyfefp.TriangleP1Mini)
d.add_field(pyfefp.TriangleP1Mini)
d.add_field(pyfefp.TriangleP1)
d.add_boundary_condition("Inner", 0, cb = lambda x : -x[:, 1])
d.add_boundary_condition("Inner", 1, cb = lambda x : x[:, 0])
d.add_boundary_condition("Outer", 0, 0.)
d.add_boundary_condition("Outer", 1, 0.)
d.add_boundary_condition("PtFix", 2, 0.)
d.complete()
self._dof = d
self._sol = d.new_vector()
self._solver = pyfefp.Steady(law, d, self._jac)
self.law = law
def set_particles(self, pos, vol, vel, forces) :
self._pf.set_particles(pos, vol, vel, forces)
def solve(self) :
tic = time.clock()
self._sol[:] = self._solver.solve(atol = 1e-10)
print("cpu : ", time.clock() - tic)
def write_solution(self, odir, oiter, time, filetype="msh") :
if filetype == "msh" :
basename = "%s/f_%05i" % (odir, oiter)
self._jac.writeScalar(basename + "-porosity.msh", "porosity", oiter, time, self._jac.dofManager.getField(0, self._pf.porosity))
self._jac.writeScalar(basename + "-pressure.msh", "pressure", oiter, time, self._dof.getField(2, self._sol))
X = np.zeros((9, self._jac.x.shape[0]))
X[0::3, :] = self._dof.getField(0, self._sol)[:-1,:]
X[1::3, :] = self._dof.getField(1, self._sol)[:-1,:]
self._jac.writeScalar(basename + "-velocity.msh", "velocity", oiter, time, X)
else :
print("unknown export format : \"%s\"." % filetype)
return
from pyfefp.mesh import Mesh
from hydro import Hydro
import scontact2
import pyfefp
import numpy as np
import shutil
import time
import os
def writep(p, odir, i) :
outputname = "%s/part-%05i" % (odir, i)
tmpoutputname = "%s/part-%05i_tmp" % (odir, i)
p.write(tmpoutputname)
shutil.move(tmpoutputname, outputname)
mu = 500e-3 #100e-3#588e-3
mu = 100e-3 #100e-3#588e-3
outputdir = "output/mu%.10g" % mu
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
p = scontact2.ParticleProblem()
p.load("deposited")
## scontact
from pyfefp import scontact2Interface
p = scontact2Interface.ScontactInterface("deposited")
tEnd = 200
dt = 0.02
predictor_corrector = False
g = [0, 0]
r = np.mean(p.particles()[:, 0])
r = np.min(p.r())
print("th dt : %g dt %g" % (r / 0.0064, dt))
volume = np.pi * p.particles()[:, [0]]**2* 2./3 ## 2./3 = 2d-3d correction
rhop = (p.particles()[:, [1]] / (volume * 3./2))[0, 0]
rhop = np.mean(p.mass()/p.volume())
rho = 1.253e3
print(rho, rhop)
outf = 1
m = Mesh("disk.msh")
hydro = Hydro(m, dt, r, rho = rho, g=g, mu = mu, rhop= rhop)
fluid = pyfefp.FluidSolver("disk.msh", dt, rhop=rhop, rho=rho, g = g, mu=mu, imex=True)
fluid.add_boundary_condition("Inner", 0, cb = lambda x : -x[:, 1])
fluid.add_boundary_condition("Inner", 1, cb = lambda x : x[:, 0])
fluid.add_boundary_condition("Outer", 0, 0.)
fluid.add_boundary_condition("Outer", 1, 0.)
fluid.add_boundary_condition("PtFix", 2, 0.)
#this has NO justification, this testcase does not really work !!
fluid.set_stab_factor(100)
fluid.complete()
t = 0
ii = 0
outputname = "%s/part-%05i" % (outputdir, 0)
p.velocity()[:, :] = 0;
p.write(outputname)
writep(p, outputdir, 0)
p.write(outputdir, 0, 0.)
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
hydro.write_solution(outputdir, 0, 0, filetype="msh")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.write_solution(outputdir, 0, 0.)
tic = time.clock()
while t < tEnd :
if predictor_corrector :
print("predictor")
pos0 = p.position().copy()
vel0 = p.velocity().copy()
hydro.solve()
p.iterate(r, dt/2, 1e-5)
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
print("corrector")
hydro.solve()
p.position()[:] = pos0
p.velocity()[:] = vel0
p.iterate(r, dt, 1e-5)
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
else :
hydro.solve()
p.iterate(r, dt, 1e-6)
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
p._p.externalForces()[:] = fluid.solve()
p.iterate(dt, p._p.externalForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
writep(p, outputdir, ioutput)
hydro.write_solution(outputdir, ioutput, t, filetype="msh")
p.write(outputdir, ioutput, t)
fluid.write_solution(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))
import numpy as np
import time
import pyfefp
class Hydro :
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
#todo add vmesh in law
law = pyfefp.Brinkman(
mu = mu,
pf = self._pf,
ns = False,
g = self.g,
rho = rho,
rhop = rhop,
stab = dt)
d = pyfefp.DofManager(mesh)
d.add_field(pyfefp.TriangleP1Mini)
d.add_field(pyfefp.TriangleP1Mini)
d.add_field(pyfefp.TriangleP1)
d.add_boundary_condition("Box", 0, 0.)
d.add_boundary_condition("Box", 1, 0.)
d.add_boundary_condition("Top", 0, 0.)
d.add_boundary_condition("Top", 1, 0.)
d.add_boundary_condition("Top", 2, 0.)
d.complete()
self.rhop = rhop
self._dof = d
self.dt = dt
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 *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._pf.velocity[:] = v0
print("cpu : ", time.clock() - tic)
def write_solution(self, odir, oiter, time, filetype="msh") :
print(oiter, time)
if filetype == "msh" :
basename = "%s/f_%05i" % (odir, oiter)
self._jac.writeScalar(basename + "-porosity.msh", "porosity", oiter, time, self._jac.dofManager.getField(0, self._pf.porosity))
self._jac.writeScalar(basename + "-pressure.msh", "pressure", oiter, time, self._dof.getField(2, self._sol))
X = np.zeros((9, self._jac.x.shape[0]))
X[0::3, :] = self._dof.getField(0, self._sol)[:-1,:]
X[1::3, :] = self._dof.getField(1, self._sol)[:-1,:]
self._jac.writeScalar(basename + "-velocity.msh", "velocity", oiter, time, X)
else :
print("unknown export format : \"%s\"." % filetype)
return
from pyfefp.mesh import Mesh
from hydro import Hydro
import scontact2
import numpy as np
import shutil
import time
import os
import sys
def writep(p, odir, i) :
outputname = "%s/part-%05i" % (odir, i)
tmpoutputname = "%s/part-%05i_tmp" % (odir, i)
p.write(tmpoutputname)
shutil.move(tmpoutputname, outputname)
def loadMeshBoundary(p, filename, tags, shift= [0, 0]) :
addv = set()
m = Mesh(filename)
pnum = [m.getPhysicalNumber(1, i) for i in tags]
for ent in m.entities :
if ent.dimension != 1 or not [i for i in pnum if i in ent.physicals]:
continue
tag = 0 if not ent.physicals else ent.physicals[0]
for i, el in enumerate(ent.elements) :
for v in el.vertices :
if v[3] in addv :
continue
addv.add(v[3])
p.addBoundaryDisk((v[0] + shift[0], v[1] + shift[1]), 0.000, tag)
v0 = el.vertices[1]
v1 = el.vertices[0]
p.addBoundarySegment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), tag)
def genInitialPosition(p, r, lx, ly, rhop) :
x = np.arange(lx-r, -lx, 2 * -r)
y = np.arange(ly-r, -ly, 2 * -r)
x, y = np.meshgrid(x, y)
x = x.flat
y = y.flat
for i in range(len(x)) :
p.addParticle((x[i], y[i]), r, r**2 * np.pi * rhop);
dt = 0.01
tEnd = 1000
outputdir = "output/fine_water"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
p = scontact2.ParticleProblem()
rhop=1100
genInitialPosition(p, 0.01, 2, 1, rhop)
loadMeshBoundary(p, "mesh.msh", ["Top", "Box"])
p.position()[:, 1] += 4
g = [0, -9.81]
r = np.mean(p.particles()[:, 0])
volume = np.pi * p.particles()[:, [0]]**2 * 2./3 ## 2./3 = 2d-3d correction
#rhop = np.mean(p.particles()[:, 1] / volume[:, 0])
print("RHOP = %g" % rhop)
outf = 5
m = Mesh("mesh.msh")
mu = 1e-3
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
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())
hydro.write_solution(outputdir, 0, 0, filetype="msh")
tic = time.clock()
while t < tEnd :
hydro.solve()
p.iterate(r, dt, 1e-5)
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
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 (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))
from pyfefp.mesh import Mesh
from hydro import Hydro
import pyfefp
import scontact2
from pyfefp import scontact2Interface
import numpy as np
import shutil
import time
import os
import sys
def writep(p, odir, i) :
outputname = "%s/part-%05i" % (odir, i)
tmpoutputname = "%s/part-%05i_tmp" % (odir, i)
p.write(tmpoutputname)
shutil.move(tmpoutputname, outputname)
def loadMeshBoundary(p, filename, tags, shift= [0, 0]) :
addv = set()
m = Mesh(filename)
pnum = [m.getPhysicalNumber(1, i) for i in tags]
for ent in m.entities :
if ent.dimension != 1 or not [i for i in pnum if i in ent.physicals]:
continue
tag = 0 if not ent.physicals else ent.physicals[0]
for i, el in enumerate(ent.elements) :
for v in el.vertices :
if v[3] in addv :
continue
addv.add(v[3])
p.addBoundaryDisk((v[0] + shift[0], v[1] + shift[1]), 0.000, tag)
v0 = el.vertices[1]
v1 = el.vertices[0]
p.addBoundarySegment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), tag)
def genInitialPosition(p, r, lx, ly, rhop) :
def genInitialPosition(filename, r, lx, ly, rhop) :
p = scontact2.ParticleProblem()
scontact2Interface.MeshLoader(p, "mesh.msh", ["Top", "Box"])
x = np.arange(lx-r, -lx + r*0.999, 2 * -r)
y = np.arange(ly-r, -ly + r*0.999, 2 * -r)
x, y = np.meshgrid(x, y)
......@@ -39,67 +15,66 @@ def genInitialPosition(p, r, lx, ly, rhop) :
y = y.flat
for i in range(len(x)) :
p.addParticle((x[i], y[i]), r, r**2 * np.pi * rhop);
p.position()[:, 1] += 4
p.write(filename)
dt = 0.0025
tEnd = 15000
outputdir = "output/fine5-imex-dt-00025-outf5"
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
p = scontact2.ParticleProblem()
rhop=1100
d = 0.01
mu = 0.001
rho = 1000
ii = 0
if ii == 0 :
genInitialPosition(p, d, 2, 1, rhop)
loadMeshBoundary(p, "mesh.msh", ["Top", "Box"])
p.position()[:, 1] += 4
p.velocity()[:, :] = 0;
else :
p.load(outputdir+"/part-%05d"%ii)
filename = outputdir + "/part-00000"
genInitialPosition(filename, d, 2, 1, rhop)
mass = p.particles()[:, [1]]
## scontact
#p = scontact2Interface.ScontactInterface(outputdir+"/part-%05d"%ii)
## lmgc
from pyfefp import lmgc2Interface
lmgc2Interface.scontactToLmgc90(filename)
p = lmgc2Interface.LmgcInterface()
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("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("dtmax = %g Re_p/u = %g" % (p.mass()[0] / (150 * mu * d), d * rho * 0.3/mu))
print("RHOP = %g" % rhop)
outf = 5
m = Mesh("mesh.msh")
hydro = Hydro(m, dt, r, rhop=rhop, rho = rho, g=g, mu = mu)
fluid = pyfefp.FluidSolver("mesh.msh", dt, rhop=rhop, rho=rho, g = g, mu=mu, imex=True)
fluid.add_boundary_condition("Box", 0, 0.)
fluid.add_boundary_condition("Box", 1, 0.)
fluid.add_boundary_condition("Top", 0, 0.)
fluid.add_boundary_condition("Top", 1, 0.)
fluid.add_boundary_condition("Top", 2, 0.)
fluid.complete()
t = 0
outputname = "%s/part-%05i" % (outputdir, 0)
p.write(outputname)
p.r = r
p.dt = dt
writep(p, outputdir, 0)
pileseg = np.where(p.segmentTag() == 3)[0]
piledisk = np.where(p.diskTag() == 3)[0]
p.write(outputdir, ii, t)
hydro._pf.set_particles(p.position(), volume, p.velocity(), p.externalForces(), mass)
hydro.write_solution(outputdir, 0, 0, filetype="msh")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.write_solution(outputdir, ii, t)
tic = time.clock()
while t < tEnd :
hydro.solve()
nsub = 20
vn = p.velocity() + p.externalForces() * dt / mass
forces = fluid.solve()
vn = p.velocity() + forces * dt / p.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)
nsub = max(1, int(np.ceil((vmax * dt * 2)/min(p.r()))))
print("NSUB", nsub, "VMAX * dt", vmax * dt, "r", min(p.r()))
for i in range(nsub) :
p.iterate(r, dt/nsub, 1e-5)
hydro._pf.set_particles(p.position(), volume, p.velocity(), p.externalForces(), mass)
p.iterate(dt/nsub, forces)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
writep(p, outputdir, ioutput)
hydro.write_solution(outputdir, ioutput, t, filetype="msh")
p.write(outputdir, ioutput, t)
fluid.write_solution(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))
This diff is collapsed.
import numpy as np
import time
import pyfefp
class Hydro :
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
#todo add vmesh in law
law = pyfefp.Brinkman(
mu = mu,
pf = self._pf,
ns = False,
g = self.g, #np.array([0, 0]),
rho = rho,
rhop = rhop,
stab = dt)
d = pyfefp.DofManager(mesh)
d.add_field(pyfefp.TriangleP1Mini)
d.add_field(pyfefp.TriangleP1Mini)
d.add_field(pyfefp.TriangleP1)
d.add_boundary_condition("Box", 0, 0.)
d.add_boundary_condition("Box", 1, 0.)
d.add_boundary_condition("Top", 0, 0.)
d.add_boundary_condition("Top", 1, 0.)
d.add_boundary_condition("Top", 2, 0.)
d.complete()
self.rhop = rhop
self._dof = d
self.dt = dt
self.law = law
self._sol = d.new_vector()
self._solver = pyfefp.Steady(law, d, self._jac, ilu=False)
self._solverEx = pyfefp.Steady(law