Commit d0bc9bac authored by Matthieu Constant's avatar Matthieu Constant
Browse files

modif testcases

parent c761e07d
......@@ -24,7 +24,7 @@ def genInitialPosition(filename, rhop) :
p.position()[:, 1] += .15/2
print(len(p.position()[:, 1]))
p.write(filename)
p.write_vtk(filename, 0, 0)
outputdir = "outputFunnelG"
......@@ -49,69 +49,48 @@ tEnd = 100
#numerical parameters
lcmin = 0.001 # approx r*100 but should match the mesh size
dt = 1e-3
alpha = 2.5e-5
alpha = 2.5e-4
epsilon = alpha*lcmin**2 /nu#2e-2*c*h*2#2e-2*c*h # ?? not sure ??1e-5
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
p.write(outputdir+"/part-00000")
genInitialPosition(outputdir + "/part-00000", rhop)
p = scontact3.ParticleProblem(outputdir+"/part-00000")
genInitialPosition(outputdir, rhop)
p = scontact3.ParticleProblem()
p.read_vtk(outputdir,0)
#print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 100
strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)]#("Top",3,0.),("Top",1,0.),("Top",0,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries, 1)
print("Initialisation fluid done!\n")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
print("Initialisation particles done!\n")
fluid.export(outputdir,0,0)
ii = 0
tic = time.clock()
#forces = g*p.mass()
while t < tEnd and ii < 50:
ticf = time.clock()
fluid.export_vtk(outputdir,0,0)
forces = g*p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
tocf = time.clock()
print("Time to compute the force : %.6e\n" % (tocf-ticf))
#print(forces)
ticp = time.clock()
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max((vn[:, 0]**2+ vn[:, 1]**2+vn[:,2]**2)**0.5)
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r()))))
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
for i in range(nsub) :
p.iterate(dt/nsub, forces)
tocp = time.clock()
print("Time to move the grains : %.6e\n" % (tocp-ticp))
ticfp = time.clock()
p.iterate(dt/nsub, forces)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
tocfp = time.clock()
print("Time to create the particle fluid : %.6e\n" % (tocfp-ticfp))
tice = time.clock()
fluid.implicit_euler(dt)
toce = time.clock()
print("Time to solve the fluid : %.6e\n" % (toce-tice))
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
p.write(outputdir+"/part-%05d" % ioutput)
p.write_boundary_vtk(outputdir, ioutput, t)
fluid.export(outputdir, t, ioutput)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
if ii==400:
break
L = 0.014;
H = 0.076;
y = 0;
lc = 0.001;
Point(1) = {-L, H, 0, lc};
Point(2) = {-L, -H, 0, lc};
Point(3) = {L, -H, 0, lc};
Point(4) = {L, H, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(1) = {1:4};
Plane Surface(1) = {1};
Physical Line("Box") = {1,2,3};
Physical Line("Top") = {4};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
import pyfefp
import scontact2
from pyfefp import scontact2Interface
import numpy as np
import time
import os
from pyfefp.legendre import *
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)
x = x.flat
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.005
tEnd = 2
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
rhop=2500
d = 0.005
mu = 0.001
rho = 1000
ii = 0
filename = outputdir + "/part-00000"
genInitialPosition(filename, d, 2, 1, rhop)
## scontact
p = scontact2Interface.ScontactInterface(outputdir+"/part-%05d"%ii)
## lmgc
#from pyfefp import lmgc2Interface
#lmgc2Interface.scontactToLmgc90(filename)
#p = lmgc2Interface.LmgcInterface()
g = [0, -9.81]
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
fluid = pyfefp.FluidSolver("mesh.msh", dt, rhop=rhop, rho=rho, mu=mu, g = g, epsilon=1e-02, imex=False, TypeEMu=TriangleP1, TypeEMp=TriangleP1)
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
p.write_vtk(outputdir, ii, t)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.write_solution(outputdir, ii, t, "vtk")
tic = time.clock()
while t < tEnd :
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)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
for i in range(nsub) :
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
p.write_vtk(outputdir, ioutput, t)
p.write(outputdir, ioutput, t)
fluid.write_solution(outputdir, ioutput, t, "vtk")
ii += 1
print("%i : %.2g/%.2g (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))
import pyfefp
import scontact2
from pyfefp import scontact2Interface
import numpy as np
import time
import os
from pyfefp.legendre import *
def genInitialPosition(filename, r, lx, ly, rhop) :
p = scontact2.ParticleProblem()
scontact2Interface.MeshLoader(p, "mesh1.msh", ["Top", "Box"])
x = np.arange(lx-r, -lx + r*0.999, -2.1*r)
y = np.arange(ly-r, -ly + r*0.999, -2.1*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);
p.position()[:, 1] += 0.07542
p.write(filename)
dt = 0.005
tEnd = 110
outputdir = "outputarttestpl1"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
rhop=2500
d = 0.000045
mu = 0.001
rho = 1000
ii = 0
filename = outputdir + "/part-00000"
genInitialPosition(filename, d/2, 0.014, 0.00055, rhop)
## scontact
p = scontact2Interface.ScontactInterface(outputdir+"/part-%05d"%ii)
## lmgc
#from pyfefp import lmgc2Interface
#lmgc2Interface.scontactToLmgc90(filename)
#p = lmgc2Interface.LmgcInterface()
g = [0, -9.81]
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 = 10
fluid = pyfefp.FluidSolver("mesh1.msh", dt, rhop=rhop, rho=rho, g = g, mu=mu, epsilon=1e-04, imex=False, TypeEMu=TriangleP1, TypeEMp=TriangleP1)
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
p.write_vtk(outputdir, ii, t)
p.write(outputdir, ii, t)
p.write_boundary_vtk(outputdir, ii, t)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.write_solution(outputdir, ii, t, "vtk")
tic = time.clock()
while t < tEnd :
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)/min(p.r()))))
print("NSUB", nsub,"v", vn[:, 1],"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
for i in range(nsub) :
p.iterate(dt/nsub, forces)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
t += dt
#dt=min(0.05,d/vmax)
print("dt",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, "vtk")
ii += 1
print("%i : %.2g/%.2g (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))
alpha 0
beta 0
scale 2.120488091862423
shift 0 1.400000000000001
visible 3 3 3 3 3 3 3 3 3 3
This diff is collapsed.
L = 0.05;
H = 0.6;
P = 0.02;
y = 0;
lc = 0.01;
Point(1) = {-L, H, -P, lc};
Point(2) = {-L, 0, -P, lc};
Point(3) = {L, 0, -P, lc};
Point(4) = {L, H, -P, lc};
Point(5) = {-L, H, P, lc};
Point(6) = {-L, 0, P, lc};
Point(7) = {L, 0, P, lc};
Point(8) = {L, H, P, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Line(9) = {5, 1};
Line(10) = {4, 8};
Line(11) = {2, 6};
Line(12) = {3, 7};
Line Loop(1) = {1:4};
Line Loop(2) = {5:8};
Line Loop(3) = {9,-4,10,8};
Line Loop(4) = {2,12,-6,-11};
Line Loop(5) = {12,7,-10,-3};
Line Loop(6) = {1,11,-5,9};
Plane Surface(1) = {1};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
Plane Surface(4) = {4};
Plane Surface(5) = {5};
Plane Surface(6) = {6};
Surface Loop(1) = {1,-3,-2,4,1,-5,-2,6};
Volume(1) = {1};
Physical Surface("Z") = {1,2};
Physical Surface("X") = {5,6};
Physical Surface("Top") = {3};
Physical Surface("Bottom") = {4};
Physical Volume("Domain") = {1};
Physical Point("PtFix") = {1};
Mesh.Algorithm = 5;
Merge "lc.pos";
Field[1] = PostView;
Field[1].IView = 0;
Background Field = 1;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.CharacteristicLengthFromPoints = 0;
#!/usr/bin/env python
from marblesbag import fluid3 as fluid
from marblesbag import scontact3
import numpy as np
import os
import time
import shutil
import random
def genInitialPosition(r, rout, rhop, compacity) :
p = scontact3.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","X","Z"])
N = compacity*rout**3/(r**3)
e = 2*rout/(6*N/np.pi)**(1./3.)
for x in np.arange(rout, -rout, -e):
for y in np.arange(rout, -rout, -e):
for z in np.arange(rout, -rout, -e):
R2 = x**2 + y**2 + z**2
if R2<(rout-r)**2:
p.add_particle((x+random.uniform(-e/2+r,e/2-r), y+random.uniform(-e/2+r,e/2-r), z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop);
p.position()[:, 1] += 0.52
print(len(p.position()[:, 1] ))
p.write_vtk(outputdir,0,0)
outputdir = "MetzgerB"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
r=154e-6
R = 3.3e-3
compacity = .20
drho = 35.4
p = scontact3.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
#physical parameters
g = -9.81
rhop = 2450
nu = 1.17/1030.
rho = 1030#rhop-drho/compacity
V = 0.5 # todo : estimate V base on limit velocity
tEnd = 1000
#numerical parameters
lcmin = 0.0004 # approx r*100 but should match the mesh size
dt = 5e-5*1000
alpha = 1e-3
epsilon = alpha*lcmin**2/nu
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
p.write_vtk(outputdir,0,0)
mu = nu*rho
genInitialPosition(r, R, rhop, compacity)
p = scontact3.ParticleProblem()
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 10
outf1 = 100000
#strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)]
#strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Bottom",0,.0),("Bottom",1,0.),("Bottom",2,0.),("X",0,.0),("X",1,0.),("X",2,0.),("Z",0,.0),("Z",1,0.),("Z",2,0.),("Top",3,0.)]
strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries, 1)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
ii = 0
jj = 0
tEnd1 = 100
#p.velocity()[:,1] = 2./9. *(compacity)* R**2*g*(rhop-rho)/mu
print(p.velocity()[1,1])
print(1-compacity)
for jj in range(4):
dt1 = dt
t = 0
if jj!=0:
fluid.adapt_mesh(1e-3,1e-4,300000, epsilon, 1)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd1 :
niter = fluid.implicit_euler(dt1)
forces = fluid.compute_node_force(dt1)
vn = p.velocity() + forces * dt1 / p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
if niter < 5 and dt1 < 20:
dt1 *= 3
elif niter > 5 and dt1 > dt:
dt1 /= 3
print(dt1)
t += dt1
# if ii %outf == 0 :
# ioutput = int(ii/outf) + 1
# p.write_vtk(outputdir, ioutput, t)
# fluid.export_vtk(outputdir, t, ioutput)
ii += 1
ii = 0
t = 0
tic = time.clock()
fluid.export_vtk(outputdir,0,0)
forces = g*p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
if (ii%5==0 and ii != 0):
fluid.adapt_mesh(1e-3,1e-4,300000, epsilon, 1)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
vn = p.velocity() + forces*dt/p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
for i in range(nsub) :
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
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
L = 0.05;
H = 0.6;
P = 0.02;
y = 0;
lc = 0.01;
Point(1) = {-L, H, -P, lc};
Point(2) = {-L, -H, -P, lc};
Point(3) = {L, -H, -P, lc};
Point(4) = {L, H, -P, lc};
Point(5) = {-L, H, P, lc};
Point(6) = {-L, -H, P, lc};
Point(7) = {L, -H, P, lc};
Point(8) = {L, H, P, lc};
Point(9) = {0,H,0,lc};
Point(10) = {0,-H,0,lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Line(9) = {5, 1};
Line(10) = {4, 8};
Line(11) = {2, 6};
Line(12) = {3, 7};
Line Loop(1) = {1:4};
Line Loop(2) = {5:8};
Line Loop(3) = {9,-4,10,8};
Line Loop(4) = {2,12,-6,-11};
Line Loop(5) = {12,7,-10,-3};
Line Loop(6) = {1,11,-5,9};
Plane Surface(1) = {1};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
Plane Surface(4) = {4};
Plane Surface(5) = {5};
Plane Surface(6) = {6};
Surface Loop(1) = {1,-3,-2,4,1,-5,-2,6};
Volume(1) = {1};
Physical Surface("Z") = {1,2};
Physical Surface("X") = {5,6};
Physical Surface("Top") = {3};
Physical Surface("Bottom") = {4};
Physical Volume("Domain") = {1};
Physical Point("PtFix") = {1};
Point(15) = {0,.52,0};
Field[1] = Attractor;
Field[1].NodesList = {15};
Field[2] = Threshold;
Field[2].DistMax = 0.06;
Field[2].DistMin = 0.01;
Field[2].LcMax = 0.008;
Field[2].LcMin = 0.0004;
Field[2].IField = 1;
Background Field = 2;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
L = 0.02;
H = 0.6;
P = 0.02;
y = 0;
lc = 0.01;
Point(1) = {-L, H, -P, lc};
Point(2) = {-L, -H, -P, lc};
Point(3) = {L, -H, -P, lc};
Point(4) = {L, H, -P, lc};
Point(5) = {-L, H, P, lc};
Point(6) = {-L, -H, P, lc};
Point(7) = {L, -H, P, lc};
Point(8) = {L, H, P, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Line(9) = {5, 1};
Line(10) = {4, 8};
Line(11) = {2, 6};
Line(12) = {3, 7};
Line Loop(1) = {1:4};
Line Loop(2) = {5:8};
Line Loop(3) = {9,-4,10,8};
Line Loop(4) = {2,12,-6,-11};