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

all test cases article with paraview state with relative path

parent 5acff405
from math import *
import matplotlib
from matplotlib.tri import Triangulation
from matplotlib import pyplot
from matplotlib.mlab import griddata
import matplotlib.cm as cm
#!/usr/bin/env python
from marblesbag import scontact3
import numpy as np
import os
import time
import shutil
import random
p = scontact3.ParticleProblem()
p.read_vtk("outputIm",20)
i = 70
R = 2.7e-3
r = []
x = []
y = []
b1 = False
b2 = False
rhop=2400
n = len(p.mass())
print(n)
print(p.position()[0,1])
fig, ax = pyplot.subplots()
for i in range(n) :
if i<n/2.:
if p.position()[i,2]<.00002 and p.position()[i,2]>-.00002:
circle = pyplot.Circle((p.position()[i,0]/R, p.position()[i,1]/R), 2*(3*p.mass()[i]/(4*np.pi*rhop))**(1./3.)/R, edgecolor='r', facecolor=(0,0,1),linewidth=0.1)
ax.add_artist(circle)
else:
if p.position()[i,2]<.00002 and p.position()[i,2]>-.00002:
circle = pyplot.Circle((p.position()[i,0]/R, p.position()[i,1]/R), 2*(3*p.mass()[i]/(4*np.pi*rhop))**(1./3.)/R, edgecolor='b', facecolor=(1,0,0),linewidth=0.1)
ax.add_artist(circle)
ax.set_xlim((-1.5, 3.5))
ax.set_ylim((.17/R, .18/R))
ax.set_aspect('equal')
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
fig.savefig('state_3d_0.png',bbox_inches='tight', dpi=900)
from math import *
import matplotlib
from matplotlib.tri import Triangulation
from matplotlib import pyplot
from matplotlib.mlab import griddata
import matplotlib.cm as cm
#!/usr/bin/env python
from marblesbag import scontact3
import numpy as np
import os
import time
import shutil
import random
p = scontact3.ParticleProblem()
p.read_vtk("outputIm",50)
i = 70
R = 2e-3
r = []
x = []
y = []
b1 = False
b2 = False
rhop=2400
n = len(p.mass())
print(n)
print(p.position()[0,1])
fig, ax = pyplot.subplots()
for i in range(n) :
if i<n/2.:
if p.position()[i,2]<.0001 and p.position()[i,2]>-.0001:
circle = pyplot.Circle((p.position()[i,0]/R, p.position()[i,1]/R), 2*(3*p.mass()[i]/(4*np.pi*rhop))**(1./3.)/R, edgecolor='r', facecolor=(0,0,1),linewidth=0.1)
ax.add_artist(circle)
else:
if p.position()[i,2]<.0001 and p.position()[i,2]>-.0001:
circle = pyplot.Circle((p.position()[i,0]/R, p.position()[i,1]/R), 2*(3*p.mass()[i]/(4*np.pi*rhop))**(1./3.)/R, edgecolor='b', facecolor=(1,0,0),linewidth=0.1)
ax.add_artist(circle)
ax.set_xlim((-2.5, 2.5))
ax.set_ylim((.181/R, .191/R))
ax.set_aspect('equal')
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
fig.savefig('stateVert_3d_0.png',bbox_inches='tight', dpi=900)
# Marblesbag - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of Marblesbag: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (Marblesbag) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
#!/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
#Physical parameters for the drops are the ones presented by Machu et al. (2001) "Coalescence, torus formation and breakup of sedimenting drops: experiments and computer simulations"
#grains creation and placement + physical boundaries definition
'''filename is the name of the output file
r is the radius of the grains
rout is the radius of the drop
rhop is the particles density
compacity is the solid volume fraction inside the drop'''
def genInitialPosition(filename, r, rout, rhop, compacity) :
p = scontact2.ParticleProblem()
#Loading of the mesh.msh file specifying physical boundaries name
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral"])
#Space between the grains to obtain the expected compacity
N = compacity*4*rout**2/(np.pi*r**2)
e = 2*rout/(N)**(1./2.)
#Definition of the points where the grains are located
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]
#First drop
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);
#Second drop
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);
#Vertical shift of the grains to the top of the box
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
#Enable computation of the initial velocity
initialize = 1
#physical parameters
g = -9.81 # gravity
rhop = 2400 # grains density
nu = 1e-4 # kinematic viscosity
drho = 35.4 # density differences rhop-rho
compacity = 0.03 # solid volume fraction in the drop
rho = rhop-drho/compacity # fluid density
r=25e-6 # grains radii
R = 2.7e-3 # drop radius
mu = nu*rho # dynamic viscosity
tEnd = 1000 # final time
print("RHOP = %g" % rhop)
#numerical parameters
lcmin = 0.0008 # mesh size
dt = 5e-2 # time step
alpha = 1e-3 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#Object particles creation
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]))
outf = 1 # number of iterations between output files
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
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
#Initial fluid velocities are computed based on the forces applied by the grains without moving the grains
if initialize:
for jj in range(4):
dt1 = dt/100
t = 0
if jj!=0:
#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
fluid.adapt_mesh(5e-3,8e-4,50000)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd1 :
#Fluid solver
niter = fluid.implicit_euler(dt1)
forces = fluid.compute_node_force(dt1)
#Changes the velocities of the grains without moving the grains
p.velocity()[:] += dt1*forces/p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#Increase time step if convergence
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()
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
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 solver
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
#Computation of the new velocities
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
#number of sub time step
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()))
#Contact solver
for i in range(nsub) :
p.iterate(dt/nsub, forces)
t += dt
#Output files writting
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))
This diff is collapsed.
# Marblesbag - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of Marblesbag: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (Marblesbag) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
#!/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
#Physical parameters for the drops are the ones presented by Machu et al. (2001) "Coalescence, torus formation and breakup of sedimenting drops: experiments and computer simulations"
#grains creation and placement + physical boundaries definition
'''filename is the name of the output file
r is the radius of the grains
rout is the radius of the drop
rhop is the particles density
compacity is the solid volume fraction inside the drop'''
def genInitialPosition(filename, r, rout, rhop, compacity) :
p = scontact2.ParticleProblem()
#Loading of the mesh.msh file specifying physical boundaries name
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral"])
#Space between the grains to obtain the expected compacity
N = compacity*4*rout**2/(np.pi*r**2)
e = 2*rout/(N)**(1./2.)
#Definition of the points where the grains are located
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]
#First drop
for i in range(x.shape[0]) :
p.add_particle((x[i]+random.uniform(-e/2+r,e/2-r), y[i]+random.uniform(-e/2+r,e/2-r)-2*rout), r, r**2 * np.pi * rhop);
#Second drop
for i in range(x.shape[0]) :
p.add_particle((x[i]+random.uniform(-e/2+r,e/2-r), y[i]+random.uniform(-e/2+r,e/2-r)+2*rout), r, r**2 * np.pi * rhop);
#Vertical shift of the grains to the top of the box
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
#Enable computation of the initial velocity
initialize = 1
#physical parameters
g = -9.81 # gravity
rhop = 2400 # grains density
nu = 1e-4 # kinematic viscosity
compacity = 0.013 # solid volume fraction in the drop
drho = 15.6 # density differences rhop-rho
rho = rhop-drho/compacity # fluid density
r=25e-6 # grains radii
mu = nu*rho # dynamic viscosity
R = 2e-3 # drop radius
tEnd = 1000 # final time
#numerical parameters
lcmin = 0.0004 # mesh size
dt = 5e-2 # time step
alpha = 1e-3 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#Object particles creation
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]))
outf = 1 # number of iterations between output files
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
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
#Initial fluid velocities are computed based on the forces applied by the grains without moving the grains
if initialize:
for jj in range(4):
dt1 = dt/100
t = 0
if jj!=0:
#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
fluid.adapt_mesh(5e-3,1e-4,50000)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd1 :
#Fluid solver
niter = fluid.implicit_euler(dt1)
forces = fluid.compute_node_force(dt1)
#Changes the velocities of the grains without moving the grains
p.velocity()[:] += dt1*forces/p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#Increase time step if convergence
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()
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
if (ii%5==0 and ii != 0):
fluid.adapt_mesh(5e-3,1e-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 solver
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
#Computation of the new velocities
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
#number of sub time step
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()))
#Contact solver
for i in range(nsub) :
p.iterate(dt/nsub, forces)
t += dt
#Output files writting
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))
......@@ -90,7 +90,7 @@ tEnd = 1000 # final time
#numerical parameters
lcmin = 0.0004 # mesh size
dt = 5e-2 # time step
alpha = 1e-3 # stabilization coefficient
alpha = 5e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
......
#!/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
#grains creation and placement + physical boundaries definition
'''filename is the name of the output file
r is the radius of the grains
rout is the radius of the drop
rhop is the particles density
compacity is the solid volume fraction inside the drop'''
def genInitialPosition(filename, r, rout, rhop1, rhop2, compacity) :
p = scontact3.ParticleProblem()
#Loading of the mesh.msh file specifying physical boundaries name
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","X","Z"])
#Space between the grains to obtain the expected compacity
N = compacity*rout**3/(r**3)
e = 2*rout/(6*N/np.pi)**(1./3.)
#Definition of the points where the grains are located
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
#First drop
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)-6*rout, z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop2);
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
#Second drop
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)+6*rout, z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop1);
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
#Second drop
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)+2*rout, z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop1);
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
#Second drop
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)-2*rout, z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop1);
print("Number of grains=%d" % len(p.position()[:, 1] ))
p.write_vtk(filename,0,0)