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

testcases ready to use + gravity in args

parent bb4517f6
......@@ -15,7 +15,7 @@ class StrongBoundary(Structure):
class fluid_problem :
def __init__(self, mesh_file_name, mu, rho, alpha, autoEpsilon, strong_boundaries):
def __init__(self, mesh_file_name, g, mu, rho, alpha, autoEpsilon, strong_boundaries):
nsb = len(strong_boundaries)
class Bnd :
def __init__(self, b) :
......@@ -31,7 +31,7 @@ class fluid_problem :
asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],BNDCB(Bnd(i[2]).apply)) for i in strong_boundaries])
self.asb = asb
lib.fluid_problem_new.restype = c_void_p
self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(mu), c_double(rho), c_double(alpha), c_bool(autoEpsilon), nsb, asb))
self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), c_double(mu), c_double(rho), c_double(alpha), c_bool(autoEpsilon), nsb, asb))
if self._fp == None :
raise NameError("cannot create fluid problem")
......
......@@ -67,7 +67,7 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
double vol = problem->particle_volume[i];
double mass = problem->particle_mass[i];
particle_force[i*2+0] = 0;
particle_force[i*2+1] = -9.81*(mass-problem->rho*vol);
particle_force[i*2+1] = problem->g*(mass-problem->rho*vol);
if(iel < 0){
continue;
}
......@@ -385,12 +385,13 @@ static void computeEpsilon(FluidProblem *problem, bool autoEpsilon)
printf("Epsilon Min = %g\n",epsMin);
}
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double alpha, bool autoEpsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double alpha, bool autoEpsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
{
FluidProblem *problem = malloc(sizeof(FluidProblem));
problem->rho = rho;
problem->alpha = alpha;
problem->mu = mu;
problem->g = g;
Mesh *mesh = mesh_load_msh(mesh_file_name);
problem->epsilon = malloc(sizeof(double)*mesh->n_triangles);
if (!mesh)
......
......@@ -16,6 +16,7 @@ typedef struct {
double rho;
double mu;
double alpha;
double g;
Mesh *mesh;
double *porosity;
double *old_porosity;
......@@ -39,7 +40,7 @@ int fluid_problem_export(const FluidProblem *problem, const char *output_dir, do
// complete force on the particle (including gravity)
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double alpha, bool autoEpsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double alpha, bool autoEpsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity);
void fluid_problem_free(FluidProblem *problem);
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin, bool autoEpsilon);
......
#!/usr/bin/env python
from marblesbag import fluid as fluid
from marblesbag import scontact2
import numpy as np
import os
import time
import shutil
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
def genInitialPosition(filename, r, rout, rhop) :
p = scontact2.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
x = np.arange(rout, -rout, 2 * -r)
y = np.arange(-rout, -2*rout/3., 2 * r)
x, y = np.meshgrid(x, y)
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], y[i]), r, r**2 * np.pi * rhop);
x = np.arange(rout, -rout, 2 * -r/2)
y = np.arange(-2*rout/3.+r, -rout/2., 2 * r/2)
x, y = np.meshgrid(x, y)
R2 = x**2 + y**2
keep = R2 < (rout - r/2)**2
x = x[keep]
y = y[keep]
for i in range(x.shape[0]) :
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(filename)
t = 0
ii = 0
#physical parameters
g = -9.81
rho = 1e3
rhop = 1600#1.18e3
nu = 1e-6
tEnd = 100
#numerical parameters
h = 0.002 # approx r*100 but should match the mesh size
dt = 1e-4#h/V*0.1
epsilon = 1e-4 # ?? not sure ??import pyfefp
rout = 0.0254
rin = 0.0064
r = 397e-6/2
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, rhop)
p = scontact2.ParticleProblem(outputdir+"/part-00000")
outf = 50
strong_boundaries = [("Outer",0,lambda x : -2*x[:, 1]),("Outer",1,lambda x : 2*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,False,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)
tic = time.clock()
while t < tEnd :
#if ii%100==0 and ii != 0:
# fluid.adapt_mesh(0.1,10,8e-4,False)
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)
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)
p.write_boundary_vtk(outputdir, ioutput, t)
fluid.export(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
#!/usr/bin/env python
import scontact2
import numpy as np
import mesh
import os
meshBoundary = True
outputdir = "output_depot"
rout = 0.0254
rin = 0.0064
r = 397e-6/2
compacity3d = 0.18
def genInitialPosition(p, N, r, rout, rin, rhop) :
x = np.arange(rout, -rout, 2 * -r)
x, y = np.meshgrid(x, x)
R2 = x**2 + y**2
keep = np.logical_and(R2 < (rout - r)**2, R2 > (rin + r) **2)
x = x[keep][:N]
y = y[keep][:N]
for i in range(N) :
p.addParticle((x[i], y[i]), r, r**2 * np.pi * rhop);
addv = set()
def loadMeshBoundary(p, filename) :
m = mesh.mesh(filename)
for ent in m.entities :
if ent.dimension != 1 :
continue
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], v[1]), 0.000, 0)
v0 = el.vertices[1]
v1 = el.vertices[0]
p.addBoundarySegment((v0[0], v0[1]), (v1[0], v1[1]), 0)
N = int((rout**2 - rin**2) * compacity3d * 3./2 / r**2)
surfacedense = N * 2 * np.sqrt(3) * r**2
surfacebulk = 0.34 * np.pi * (rout**2 - rin**2)
p = scontact2.ParticleProblem()
if meshBoundary :
loadMeshBoundary(p, "disk.msh")
else :
p.addBoundaryDisk((0, 0), rin, 0)
p.addBoundaryDisk((0, 0), -rout, 0)
genInitialPosition(p, N, r, rout, rin, rhop=1.18e3)
scale = (surfacebulk/surfacedense)**0.5
print("total compacity 3d (0.2) = %g\n" % (N * r**2 * 2./3 / (rout**2 - rin**2)))
print("scale = %g (%g %g)" % (scale, surfacedense, surfacebulk))
r *= scale
g = 1
tEnd = 3 * (4 * 0.0254 / g)**0.5;
t = 0
i = 0
ioutput = 1
os.makedirs(outputdir, exist_ok=True)
while t < tEnd :
p.externalForces()[:, :] = - 3 * p.velocity()
p.externalForces()[:, 1] += -g
p.externalForces()[:, :] *= p.particles()[:, [1]]
dt = min(0.001, p.maxdt(r)/2)
p.iterate(r, dt, tol=1e-5)
t += dt
i += 1
print("%.2g %.2g" % (t, tEnd))
if i %100 == 0 :
p.write("%s/part-%05i" % (outputdir, ioutput))
ioutput += 1
p.write("deposited")
#!/usr/bin/env python
from marblesbag import fluid as fluid
from marblesbag import scontact2
import numpy as np
import os
import time
import shutil
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
def genInitialPosition(filename, r, rout, rhop) :
p = scontact2.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
x = np.arange(rout, -rout, 2 * -r)
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]) :
if x[i]>-y[i]+0.005:
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
p.write(filename)
t = 0
ii = 0
#physical parameters
g = 0
rho = 1.253e3
rhop = 1000#1.18e3
nu = 1e-4
tEnd = 50
#numerical parameters
h = 0.002 # approx r*100 but should match the mesh size
dt = 1e-3#h/V*0.1
epsilon = 1e-5 # ?? not sure ??import pyfefp
rout = 0.0254
rin = 0.0064
r = 397e-6/2
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, rhop)
p = scontact2.ParticleProblem(outputdir+"/part-00000")
outf = 50
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,False,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)
tic = time.clock()
while t < tEnd :
#if ii%100==0 and ii != 0:
# fluid.adapt_mesh(0.1,10,8e-4,False)
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)
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)
p.write_boundary_vtk(outputdir, ioutput, t)
fluid.export(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
import pyfefp
import numpy as np
import time
import os
mu = 100e-3 #100e-3#588e-3
outputdir = "output/mu%.10g" % mu
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
## scontact
from pyfefp import scontact2Interface
p = scontact2Interface.ScontactInterface("deposited")
tEnd = 200
dt = 0.02
g = [0, 0]
r = np.min(p.r())
print("th dt : %g dt %g" % (r / 0.0064, dt))
rhop = np.mean(p.mass()/p.volume())
rho = 1.253e3
print(rho, rhop)
outf = 1
fluid = pyfefp.FluidSolver("disk.msh", dt, rhop=rhop, rho=rho, g = g, mu=mu, imex=False)
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.write(outputdir, 0, 0.)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.write_solution(outputdir, 0, 0.)
tic = time.clock()
while t < tEnd :
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
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))
......@@ -11,60 +11,60 @@ import shutil
def genInitialPosition(filename, r, rout, rhop) :
p = scontact2.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
x = np.arange(rout, -rout, 2 * -r)
x = np.arange(rout, -rout, -2*r)#-1.1441e-4/ 2**0.5)
x, y = np.meshgrid(x, x)
R2 = x**2 + y**2
keep = R2 < (rout - r)**2
x = x[keep]
y = y[keep]
x1 = np.arange(rout, -rout, 2 * -r)
x1 = np.arange(rout, -rout, -2*r)#-1.1441e-4/ 2**0.5)
x1, y1 = np.meshgrid(x1, x1)
R2 = x1**2 + y1**2
keep = R2 < (rout - r)**2
x1 = x1[keep]
y1 = y1[keep]
for i in range(x.shape[0]) :
p.add_particle((x[i]-rout, y[i]), r, r**2 * np.pi * rhop);
p.add_particle((x[i]-rout, y[i]+4*rout), r, r**2 * np.pi * rhop);
for i in range(x1.shape[0]) :
p.add_particle((x1[i]+rout, y1[i]-4*rout), r, r**2 * np.pi * rhop);
p.add_particle((x1[i]+rout, y1[i]), r, r**2 * np.pi * rhop);
print("n_particles", x.shape[0]+x1.shape[0])
p.position()[:, 1] += 0.13
p.position()[:, 1] += 0.17
p.write(filename)
outputdir = "output-2drops"
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
r=25e-6/(2)
r=25e-6
p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
#physical parameters
g = [0, -9.81]
g = -9.81
rho = 1200
rhop = 2400
nu = 1e-4
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 10
tEnd = 20
R = 1e-3
#numerical parameters
h = 0.001 # approx r*100 but should match the mesh size
dt = 1e-4
c = h/dt * 7/50/2#50
#dt= 1e-6#0.5*h/c
epsilon = 1e-6#2e-2*c*h*2#2e-2*c*h # ?? not sure ??1e-5
print("c = %g, eps = %g" % (c, epsilon))
lcmin = 5e-4 # approx r*100 but should match the mesh size
dt = 5e-4
alpha = 2.5e-3
epsilon = alpha*lcmin**2 /nu#2e-2*c*h*2#2e-2*c*h # ?? not sure ??1e-5
autoEpsilon = False
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, 1e-3, rhop)
genInitialPosition(outputdir + "/part-00000", r, R, rhop)
p = scontact2.ParticleProblem(outputdir+"/part-00000")
......@@ -74,7 +74,7 @@ outf = 100
outf1 = 100000
strong_boundaries = [("Top",2,0.),("Top",1,0.),("Top",0,0.),("Box",1,0.),("Box",0,0.)]
fluid = fluid.fluid_problem("mesh.msh",nu*rho,rho,epsilon,False,strong_boundaries)
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,autoEpsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)
......@@ -82,8 +82,8 @@ fluid.export(outputdir,0,0)
tic = time.clock()
forces = g*p.mass()
while t < tEnd :
if ii%50==0 and ii != 0:
fluid.adapt_mesh(0.05,5,5e-4,False)
if (ii%200==0 and ii != 0) or ii==100:
fluid.adapt_mesh(0.01,500,1e-4,autoEpsilon)
forces = fluid.compute_node_force(dt)
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
......
L = 0.025;
H = 0.15;
H = 0.5;
y = 0;
lc = 0.001;
lc = 0.1;
Point(1) = {-L, H, 0, lc};
Point(2) = {-L, 0.0, 0, lc};
Point(3) = {L, 0.0, 0, lc};
Point(2) = {-L, 0, 0, lc};
Point(3) = {L, 0, 0, lc};
Point(4) = {L, H, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
......
#!/usr/bin/env python
from marblesbag import fluid as Fluid
from marblesbag import scontact2Interface
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) :
D = 1.2e-4/2**0.5
p = scontact2.ParticleProblem()
scontact2Interface.MeshLoader(p, "mesh.msh", ["Top", "Box"])
x = np.arange(rout, -rout, -2*r)
p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
x = np.arange(rout, -rout, -D)
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]) :
if y[i]<rout/2 :
p.addParticle((x[i], y[i]), r, r**2 * np.pi * rhop);
if y[i]<0 :
p.add_particle((x[i]-D/2.1+2*random.uniform(0,D/4), y[i]-D/2.1+2*random.uniform(0,D/4)), r, r**2 * np.pi * rhop);
x = np.arange(rout*(3**0.5)/2, -rout*(3**0.5)/2, -2*r)
y = np.arange(0.5*rout, 2*rout, 2*r)
x = np.arange(rout, -rout, -D)
y = np.arange(2*r, 4*rout, D)
x, y = np.meshgrid(x, y)
x = x.flat
y = y.flat
T = rout
x1 = -4/5*rout
x2 = 4/5*rout
for i in range(len(x)) :
if y[i] < -4/(3**0.5)*x[i]+5*rout/2 and y[i] < 4/(3**0.5)*x[i]+5*rout/2 :
p.addParticle((x[i], y[i]), r, r**2 * np.pi * rhop);
p.position()[:, 1] += 0.14
if y[i] < -T/x1*x[i]+T and y[i] < -T/x2*x[i]+T :
p.add_particle((x[i]-D/2.1+2*random.uniform(0,D/4.1), y[i]-D/2.1+2*random.uniform(0,D/4.1)), r, r**2 * np.pi * rhop);
p.position()[:, 1] += 0.45
p.write(filename)
......@@ -41,50 +44,53 @@ if not os.path.isdir(outputdir) :
t = 0
ii = 0
r=25e-6/(2)
r=25e-6/2
p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
#physical parameters
g = [0, -9.81]
rho = 1200
g = -9.81
rho = 561
rhop = 2400
nu = 1e-4
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 1