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

updatetestcases

parent a1e0eaec
......@@ -371,6 +371,7 @@ HXTStatus hxtLinearSystemPETScMapFromVec(HXTLinearSystemPETSc *system, Vec vec,
int nf = system->nFields;
for (int i=0; i< system->nNodes; ++i) {
uint32_t n= system->nodeMap[i];
if(n==-1) continue;
for (int j = 0; j < nf; ++j)
v[i*nf+j] = vv[n*nf+j];
}
......
......@@ -10,7 +10,7 @@
#define n_fields (DIMENSION+1)
#define newton_max_it 10
#define newton_atol 2.5e-6
#define newton_atol 5e-7
#define newton_rtol 1e-5
#if DIMENSION==2
......@@ -449,7 +449,12 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
problem->porosity[el[j]] += sf[j]*volume[i];
}
for (int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = (1-problem->porosity[i]/problem->node_volume[i]);
if(problem->node_volume[i]==0.){
problem->porosity[i] = 1.;
}
else{
problem->porosity[i] = (1-problem->porosity[i]/problem->node_volume[i] /4.);
}
}
}
......
#!/usr/bin/env python
use_lmgc = 0
use_lmgc = 1
from marblesbag import fluid as fluid
......@@ -38,17 +38,17 @@ def genInitialPosition(filename, r, rout, rhop) :
x = x[keep]
y = y[keep]
for i in range(x.shape[0]) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
# x = np.arange(rout, -rout, 2 * -r/3)
# y = np.arange(-rout/2.+r, -5*rout/12., 2 * r/3)
# x, y = np.meshgrid(x, y)
# R2 = x**2 + y**2
# keep = R2 < (rout - r/3)**2
# x = x[keep]
# y = y[keep]
# for i in range(x.shape[0]) :
# p.add_particle((x[i], y[i]), r/3, r**2 /9 * np.pi * rhop);
p.write(filename)
p.add_particle((x[i], y[i]), r/2., r**2/4. * np.pi * rhop);
x = np.arange(rout, -rout, 2 * -r/3)
y = np.arange(-rout/2.+r, -5*rout/12., 2 * r/3)
x, y = np.meshgrid(x, y)
R2 = x**2 + y**2
keep = R2 < (rout - r/3)**2
x = x[keep]
y = y[keep]
for i in range(x.shape[0]) :
p.add_particle((x[i], y[i]), r/3, r**2 /9 * np.pi * rhop);
p.write_vtk(filename,0,0)
t = 0
......@@ -63,7 +63,7 @@ tEnd = 100
#numerical parameters
h = 0.002 # approx r*100 but should match the mesh size
dt = 5e-4#h/V*0.1
dt = 5e-3#h/V*0.1
epsilon = 1e-4 # ?? not sure ??import pyfefp
......@@ -74,7 +74,7 @@ compacity3d = 0.18
N = int((rout**2 - rin**2) * compacity3d * 3./2 / r**2)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
genInitialPosition(outputdir + "/part-00000", r, rout, rhop)
genInitialPosition(outputdir, r, rout, rhop)
if use_lmgc :
......@@ -82,9 +82,10 @@ if use_lmgc :
lmgc2Interface.scontactToLmgc90(outputdir+"/part-00000", friction)
p = lmgc2Interface.LmgcInterface()
else :
p = scontact2.ParticleProblem(outputdir+"/part-00000")
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
outf = 10
outf = 5
strong_boundaries = [("Outer",0,lambda x : -4*x[:, 1]),("Outer",1,lambda x : 4*x[:, 0]),("Inner",0,0.),("Inner",1,0.),("PtFix",2,0.)]
#g is define in fluid_problem
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
......
......@@ -26,7 +26,7 @@ def genInitialPosition(filename, r, rout, rin, rhop) :
y = y[keep]
for i in range(x.shape[0]) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
p.write(filename)
p.write_vtk(filename,0,0)
t = 0
......@@ -41,7 +41,7 @@ tEnd = 50
#numerical parameters
h = 0.002 # approx r*100 but should match the mesh size
dt = 1e-3#h/V*0.1
dt = 1e-2#h/V*0.1
epsilon = 1e-5 # ?? not sure ??import pyfefp
rout = 0.0254
......@@ -51,12 +51,12 @@ compacity3d = 0.18
N = int((rout**2 - rin**2) * compacity3d * 3./2 / r**2)
p = scontact2.ParticleProblem()
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
p.write(outputdir+"/part-00000")
genInitialPosition(outputdir + "/part-00000", r, rout, rin, rhop)
p = scontact2.ParticleProblem(outputdir+"/part-00000")
genInitialPosition(outputdir, r, rout, rin, rhop)
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
outf = 50
outf = 5
strong_boundaries = [("Outer",0,lambda x : -x[:, 1]),("Outer",1,lambda x : x[:, 0]),("Inner",0,0.),("Inner",1,0.),("PtFix",2,0.)]
#g is define in fluid_problem
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
......@@ -65,21 +65,21 @@ fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)
tic = time.clock()
while t < tEnd :
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
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)
p.iterate(dt/nsub, forces,1e-6)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.implicit_euler(dt)
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
p.write(outputdir+"/part-%05d" % ioutput)
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))
......@@ -17,12 +17,12 @@ def genInitialPosition(filename, r, ly, rhop) :
x = x.flat
y = y.flat
for i in range(len(x)) :
rhop = random.choice([1100,1100,1100])
rhop = random.choice([1100,1200,1300])
p.add_particle((x[i], y[i]), r/2, r**2/4 * np.pi * rhop);
p.write(filename)
p.write_vtk(filename,0,0)
outputdir = "outputdep"
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
......@@ -45,35 +45,25 @@ tEnd = 100
#numerical parameters
lcmin = 0.001 # approx r*100 but should match the mesh size
dt = 1e-4
dt = 1e-2
alpha = 2.5e-5
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(outputdir+"/part-00000")
genInitialPosition(outputdir + "/part-00000", r, ly, rhop)
genInitialPosition(outputdir, r, ly, rhop)
p = scontact2.ParticleProblem(outputdir+"/part-00000")
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 200
outf = 20
outf1 = 100000
ii = 0
def outerBndV(x) :
bndLoop = 12000
if ii%bndLoop<11000 and ii%bndLoop>9500:
print(min(((ii-9500)%bndLoop)*5e-5,.05))
return min(((ii-9500)%bndLoop)*5e-5,.05)
else:
print(.05/((ii-11000)%bndLoop + 1))
return .05/((ii-11000)%bndLoop + 1)
strong_boundaries = [("Top",2,0.),("TopOut",1,0),("Top",1,0),("BottomOut",1,0),("Bottom",1,0),("Lateral",0,0.),("Lateral",1,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
......@@ -83,17 +73,17 @@ fluid.export(outputdir,0,0)
tic = time.clock()
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)
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)
p.iterate(dt/nsub, forces,1e-6)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.implicit_euler(dt)
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
......
......@@ -8,75 +8,75 @@ import time
import shutil
import random
outputdir = "MetzgerBReload"
def genInitialPosition(filename, rhop) :
p = scontact3.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
myRadius = np.genfromtxt('radius.csv', delimiter=',')
r1=np.amin(myRadius)*1e-6
r2=np.amax(myRadius)*1e-6
i = 0
for y in np.arange(0, -.004, -2*r2):
for z in np.arange(.01675-r2, -.01675, -2*r2):
for x in np.arange(0.014-r2, -.014+r2, -2*r2):
p.add_particle((x, y, z), myRadius[i%500000]*1e-6, 4./3.* (myRadius[i%500000]*1e-6)**3 * np.pi * rhop);
i += 1
p.position()[:, 1] += .15/2
print(len(p.position()[:, 1]))
p.write_vtk(filename, 0, 0)
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
iReload = 326
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
rho = 1000
rhop = 2640
nu = 1e-6
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 1000
tEnd = 100
#numerical parameters
lcmin = 0.0004 # approx r*100 but should match the mesh size
dt = 5e-2
alpha = 1e-3
epsilon = alpha*lcmin**2 /nu
lcmin = 0.001 # approx r*100 but should match the mesh size
dt = 1e-3
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"))
mu = nu*rho
genInitialPosition(outputdir, rhop)
p = scontact3.ParticleProblem()
p.read_vtk("MetzgerB1",iReload)
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
#print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 1
outf1 = 100000
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.)]
#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.)]
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)
print("Import start!!!")
fluid.import_vtk("MetzgerB"+"/fluid_%5g.vtu"%iReload)
print("Import done!!!")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
print("Initialisation fluid done!\n")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
print("Initialisation particles done!\n")
tEnd1 = 100
ii = iReaload*outf+1
t = ii*dt
tic = time.clock()
forces = g*p.mass()
fluid.export_vtk(outputdir,0,0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
if (ii%5==0 and ii != 0):
fluid.adapt_mesh(5e-3,8e-4,600000)
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()
......
L = 0.014;
H = 0.076;
P = 0.01675;
y = 0;
lc = 0.004;
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(13) = {9,10};
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("Box") = {1,2,4,5,6};
Physical Surface("Top") = {3};
Physical Volume("Domain") = {1};
Physical Point("PtFix") = {1};
This diff is collapsed.
#!/usr/bin/env python
from marblesbag import fluid as fluid
from marblesbag import scontact2
import numpy as np
import os
import time
import shutil
import random
def genInitialPosition(filename, r, rout, rhop, compacity) :
p = scontact2.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral"])
N = compacity*4*rout**2/(np.pi*r**2)
e = 2*rout/(N)**(1./2.)
x = np.arange(rout, -rout, -e)
x, y = np.meshgrid(x, x)
R2 = x**2 + y**2
keep = R2 < (rout - r)**2
x = x[keep]
y = y[keep]
for i in range(x.shape[0]) :
p.add_particle((x[i]+random.uniform(-e/2+r,e/2-r)+rout, y[i]+random.uniform(-e/2+r,e/2-r)-2*rout), r, r**2 * np.pi * rhop);
for i in range(x.shape[0]) :
p.add_particle((x[i]+random.uniform(-e/2+r,e/2-r)-rout, y[i]+random.uniform(-e/2+r,e/2-r)+2*rout), r, r**2 * np.pi * rhop);
p.position()[:, 1] += 0.22
p.write_vtk(filename,0,0)
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
initialize = 1
r=25e-6
R = 2.7e-3
compacity = 0.03
drho = 35.4
p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
#physical parameters
g = -9.81
rhop = 2400
nu = 1e-4
rho = rhop-drho/compacity
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 1000
#numerical parameters
lcmin = 0.0008 # approx r*100 but should match the mesh size
dt = 5e-2
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"))
mu = nu*rho
genInitialPosition(outputdir, r, R, rhop,compacity)
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 1
outf1 = 100000
strong_boundaries = [("Top",2,0.),("Top",1,0.),("Bottom",1,0.),("Lateral",0,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
ii = 0
jj = 0
tEnd1 = .2
if initialize:
for jj in range(4):
dt1 = dt/100
t = 0
if jj!=0:
fluid.adapt_mesh(5e-3,8e-4,50000)
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)
p.velocity()[:] += dt1*forces/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/100:
dt1 /= 3
print(dt1)
t += dt1
ii += 1
ii = 0
t = 0
fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.clock()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
if (ii%5==0 and ii != 0):
fluid.adapt_mesh(5e-3,8e-4,50000)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
else :
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(10*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)
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.025;
H = 0.25;
y = 0;
lc = 0.001;
Point(1) = {-L, H, 0};
Point(2) = {-L, -H, 0};
Point(3) = {L, -H, 0};
Point(4) = {L, H, 0};
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("Bottom") = {2};
Physical Line("Lateral") = {1,3};
Physical Line("Top") = {4};
Physical Surface("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;
L = 0.025;
H = 0.25;
y = 0;
lc = 0.001;
Point(1) = {-L, H, 0};
Point(2) = {-L, -H, 0};
Point(3) = {L, -H, 0};