Skip to content

GitLab

Commit 40249b20 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

start injection

parent 984e7949
Hide whitespace changes
Inline Side-by-side
testcases/injectionWGrains/dep.py → testcases/injectionWGrains/depot.py
View file @ 40249b20
#!/usr/bin/env python
from migflow import fluid as fluid
from migflow import scontact2
from migflow import fluid
from migflow import scontact
import numpy as np
import os
......@@ -9,7 +9,7 @@ import shutil
import random
def genInitialPosition(filename, r, N, rhop) :
p = scontact2.ParticleProblem()
p = scontact.ParticleProblem(2)
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral","Injection"])
dr = r/10
x = np.arange(r+dr, .222-r-dr, 2*(r+dr))
......@@ -36,30 +36,21 @@ ii = 0
r = 5e-4
ly = 5e-2
p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
p = scontact.ParticleProblem(2)
#physical parameters
# physical parameters
g = -9.81
rho = 1000
rhop = 2500
nu = 1e-6
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 2.5
#numerical parameters
lcmin = 0.05 # approx r*100 but should match the mesh size
# numerical parameters
dt = 1e-4
alpha = 2.5e-5
epsilon = alpha*lcmin**2 /nu
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
N = 40000
genInitialPosition(outputdir, r, N, rhop)
p = scontact2.ParticleProblem()
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
......@@ -69,10 +60,6 @@ outf1 = 100000
ii = 0
tic = time.clock()
forces = np.zeros_like(p.velocity())
print(forces[:,1].shape)
......
testcases/injectionWGrains/inject.py
View file @ 40249b20
......@@ -20,16 +20,16 @@
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
from migflow import fluid as fluid
from migflow import scontact2
# TESTCASE DESCRIPTION
# Fluid injection jet in same fluid
from migflow import fluid
import numpy as np
import os
import time
import shutil
import random
#Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
outputdir = "output"
if not os.path.isdir(outputdir) :
......@@ -44,22 +44,20 @@ g = 0 # gravity
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
mu = nu*rho # dynamic viscosity
tEnd = 10 # final time
#numerical parameters
lcmin = .1 # mesh size
tEnd = 10 # final time
dt = .01 # time step
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
outf = 1 # number of iterations between output files
def outerBndV(x) :
print(.4*max(np.sin(t*np.pi*2./1),0))
return 0.4*max(np.sin(t*np.pi*2./1),0)
strong_boundaries = [("Top",2,2,0.),("Injection",1,1,outerBndV)]#,("Lateral",0,0,0),("Bottom",1,1,0)]
fluid = fluid.fluid_problem("mesh.msh",g,[nu*rho],[rho],strong_boundaries,0)
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Injection",1,outerBndV)
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Outflow")
......@@ -67,12 +65,8 @@ fluid.set_weak_boundary("Injection","Inflow")
ii = 0
t = 0
#set initial_condition
fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.clock()
while t < tEnd :
#Fluid solver
......
testcases/injectionWGrains/inject1fGrains.py
View file @ 40249b20
......@@ -20,18 +20,20 @@
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
from migflow import fluid as fluid
from migflow import scontact2
# TESTCASE DESCRIPTION
# Injection of fluid in granular matrix immersed in the same fluid
from migflow import fluid
from migflow import scontact
import numpy as np
import os
import time
import shutil
import random
#Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
def genInitialPosition(filename, r, N, rhop) :
p = scontact2.ParticleProblem()
p = scontact.ParticleProblem(2)
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral","Injection"])
x = np.arange(-.111+r, .111-r, 2*r)
y = np.arange(-.08+r, .08-r, 2*r)
......@@ -53,23 +55,19 @@ t = 0
ii = 0
#physical parameters
# physical parameters
g = -9.81 # gravity
rho0 = 900 # fluid density
rho1 = 1000
nu1 = 1e-6 # kinematic viscosity
nu0 = nu1
tEnd = 10 # final time
rho = 900 # fluid density
nu = 1e-6
r = 5e-4
N = 10000
rhop = 1500
#numerical parameters
lcmin = .1 # mesh size
# numerical parameters
tEnd = 10 # final time
dt = .01 # time step
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
genInitialPosition(outputdir, r, N, rhop)
p = scontact2.ParticleProblem()
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
outf = 1 # number of iterations between output files
......@@ -78,7 +76,10 @@ def outerBndV(x) :
return 0.4*max(np.sin(t*np.pi*2./1),0)
strong_boundaries = [("Top",2,2,0.),("Injection",1,1,outerBndV)]#,("Lateral",0,0,0),("Bottom",1,1,0)]
fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho0,nu1*rho1],[rho0,rho1],strong_boundaries,0)
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Injection",1,outerBndV)
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Outflow")
......@@ -87,19 +88,10 @@ ii = 0
t = 0
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#set initial_condition
fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.clock()
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(1e-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)
......@@ -119,4 +111,4 @@ while t < tEnd :
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))
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
\ No newline at end of file
testcases/injectionWGrains/testReload/inject.py deleted 100644 → 0
View file @ 984e7949
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (MigFlow) 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 migflow import fluid as fluid
from migflow import scontact2
import numpy as np
import os
import time
import shutil
import random
#Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
#physical parameters
g = -9.81 # gravity
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
mu = nu*rho # dynamic viscosity
tEnd = 10 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = .05 # time step
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
outf = 1 # number of iterations between output files
def outerBndV(x) :
print(-0.01/(0.004**2)* (x[:,0]-0.004)*(x[:,0]+0.004))
return -0.01/(0.004**2)* (x[:,0]-0.004)*(x[:,0]+0.004)
fluid = fluid.fluid_problem(g,[nu*rho],[rho])
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Injection",1,outerBndV )
fluid.set_strong_boundary("Injection",0,0)
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Outflow")
fluid.set_weak_boundary("Injection","Inflow")
ii = 0
t = 0
#set initial_condition
fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.clock()
file = open("log.txt","w")
while t < tEnd :
if t>4.95:
file.write("TIME=%g, ITER=%d, IOUTPUT=%d\n" % (t,ii,int(ii/outf)+1))
for i in range(fluid.solution().shape[0]):
file.write("Node %d: u=%g, v=%g, p=%g\n"%(i,fluid.solution()[i,0],fluid.solution()[i,1],fluid.solution()[i,2]))
#Fluid solver
fluid.implicit_euler(dt)
t += dt
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
#print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
file.close()
\ No newline at end of file
testcases/injectionWGrains/testReload/inject2f.py deleted 100644 → 0
View file @ 984e7949
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (MigFlow) 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 migflow import fluid as fluid
from migflow import scontact2
import numpy as np
import os
import time
import shutil
import random
#Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
#physical parameters
g = -9.81 # gravity
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
mu = nu*rho # dynamic viscosity
tEnd = 10 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = .005 # time step
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
outf = 1 # number of iterations between output files
def outerBndV(x) :
print(-0.01/(0.004**2)* (x[:,0]-0.004)*(x[:,0]+0.004))
return -0.01/(0.004**2)* (x[:,0]-0.004)*(x[:,0]+0.004)
fluid = fluid.fluid_problem(g,[nu*rho/2, nu*rho],[rho/2, rho])
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Injection",1,outerBndV )
fluid.set_strong_boundary("Injection",0,0)
fluid.set_strong_boundary("Injection",3,1)
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Outflow")
fluid.set_weak_boundary("Injection","Inflow")
ii = 0
t = 0
#set initial_condition
fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.clock()
file = open("log2f.txt","w")
while t < tEnd :
if t>5-dt:
file.write("TIME=%g, ITER=%d, IOUTPUT=%d\n" % (t,ii,int(ii/outf)+1))
for i in range(fluid.solution().shape[0]):
file.write("Node %d: u=%g, v=%g, p=%g\n"%(i,fluid.solution()[i,0],fluid.solution()[i,1],fluid.solution()[i,2]))
#Fluid solver
fluid.implicit_euler(dt)
t += dt
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
file.close()
\ No newline at end of file
testcases/injectionWGrains/testReload/inject2fGrains.py deleted 100644 → 0
View file @ 984e7949
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (MigFlow) 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 migflow import fluid as fluid
from migflow import scontact2
import numpy as np
import os
import time
import shutil
import random
#Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
def genInitialPosition(filename, r, N, rhop) :
p = scontact2.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral","Injection"])
x = np.arange(-.01+r, .01-r, 2*r)
i = 0
for xl in x:
p.add_particle((xl, -.01+r), r, r**2 * np.pi * rhop);
for xl in x:
p.add_particle((xl, -.01+3*r), r, r**2 * np.pi * rhop);
for xl in x:
p.add_particle((xl, -.01+5*r), r, r**2 * np.pi * rhop);
for xl in x:
p.add_particle((xl, -.01+6*r), r, r**2 * np.pi * rhop);
print(len(p.mass()))
p.write_vtk(filename,0,0)
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
#physical parameters
g = -9.81 # gravity
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
mu = nu*rho # dynamic viscosity
tEnd = 10 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = .005 # time step
genInitialPosition(outputdir, 0.0001, 1000, 2600)
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
outf = 1 # number of iterations between output files
def outerBndV(x) :
print(-0.01/(0.004**2)* (x[:,0]-0.004)*(x[:,0]+0.004))
return -0.01/(0.004**2)* (x[:,0]-0.004)*(x[:,0]+0.004)
fluid = fluid.fluid_problem(g,[nu*rho/2, nu*rho],[rho/2, rho])
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Injection",1,outerBndV )
fluid.set_strong_boundary("Injection",0,0)
fluid.set_strong_boundary("Injection",3,1)
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Outflow")
fluid.set_weak_boundary("Injection","Inflow")
ii = 0
t = 0
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#set initial_condition
fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.clock()
file = open("log2fG.txt","w")
while t < tEnd :
if t>5-dt:
file.write("TIME=%g, ITER=%d, IOUTPUT=%d\n" % (t,ii,int(ii/outf)+1))
for i in range(fluid.solution().shape[0]):
file.write("Node %d: u=%g, v=%g, p=%g\n"%(i,fluid.solution()[i,0],fluid.solution()[i,1],fluid.solution()[i,2]))
for i in range(p.mass().shape[0]):
file.write("Grain %d: up=%g, vp=%g, xp=%g, yp=%g\n"%(i,p.velocity()[i,0],p.velocity()[i,1],p.position()[i,0],p.position()[i,1]))
#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
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#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))
file.close()
\ No newline at end of file
testcases/injectionWGrains/testReload/inject2fGrainsReload.py deleted 100644 → 0
View file @ 984e7949
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (MigFlow) 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 migflow import fluid as fluid
from migflow import scontact2
import numpy as np
import os
import time
import shutil
import random
#Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
def genInitialPosition(filename, r, N, rhop) :
p = scontact2.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral","Injection"])
x = np.arange(-.01+r, .01-r, 2*r)
i = 0
for xl in x:
p.add_particle((xl, -.01+r), r, r**2 * np.pi * rhop);
for xl in x:
p.add_particle((xl, -.01+3*r), r, r**2 * np.pi * rhop);