Commit 86e8b17c authored by Nathan Coppin's avatar Nathan Coppin
Browse files

updating testcases part 1

parent 9ccc26a3
Pipeline #8908 passed with stages
in 5 minutes and 22 seconds
......@@ -63,7 +63,7 @@ fluid.set_open_boundary("Bottom",velocity=[0,outerBndV])
ii = 0
t = 0
fluid.export_vtk(outputdir,0,0)
fluid.write_vtk(outputdir,0,0)
tic = time.clock()
while t < tEnd :
......@@ -73,6 +73,6 @@ while t < tEnd :
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
......@@ -25,6 +25,7 @@
# Injection of fluid in granular matrix immersed in the same fluid
from migflow import fluid
from migflow import scontact
from migflow import time_integration
import numpy as np
import os
......@@ -56,10 +57,10 @@ ii = 0
# physical parameters
g = -9.81 # gravity
g = np.array([0,-9.81]) # gravity
rho = 1000
nu = 1e-6 # kinematic viscosity
r = 5e-4
r = 1e-3
N = 20000
lx = .222
ly = .16
......@@ -88,28 +89,16 @@ fluid.set_open_boundary("Bottom",velocity=[0,outerBndV])
ii = 0
t = 0
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export_vtk(outputdir,0,0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
tic = time.clock()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
fluid.write_vtk(outputdir,0,0)
tic = time.time()
while t < tEnd :
#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(2, 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)
time_integration.iterate(fluid,p,dt,min_nsub=5,external_particles_forces=g*p.mass())
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)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......@@ -24,6 +24,7 @@
# TESTCASE DESCRIPTION
# Injection of a light fluid in a heavier fluid
from migflow import fluid
from migflow import time_integration
import numpy as np
import os
......@@ -74,8 +75,6 @@ fluid.set_wall_boundary("Lateral")
fluid.set_open_boundary("Top",pressure=0,concentration=0)
fluid.set_open_boundary("Bottom",velocity=[0,outerBndV],concentration=1)
#fluid.set_concentration_cg(np.where(fluid.coordinates()[:,1]>0.32,1,0))
ii = 0
t = 0
......@@ -83,11 +82,11 @@ fluid.export_vtk(outputdir,0,0)
tic = time.time()
while t < tEnd :
#Fluid solver
fluid.implicit_euler(dt)
time_integration.iterate(fluid,None,dt)
t += dt
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......@@ -27,6 +27,7 @@
from migflow import fluid
from migflow import scontact
from pylmgc90 import pre
from migflow import time_integration
import numpy as np
import os
import time
......@@ -54,13 +55,13 @@ if not os.path.isdir(outputdir) :
# physical parameters
alpha = 0*np.pi/4.
g = -9.81*np.cos(alpha) # gravity
g = np.array([0,-9.81]) # gravity
rho0 = 1.117 # fluid density
rho1 = 1000
nu1 = 1.57e-5 # kinematic viscosity
nu0 = 1e-6
tEnd = 10 # final time
r = 5e-4
r = 1e-3
N = 10000*3
lx = .222
ly = .14*2
......@@ -69,8 +70,6 @@ rhop = 2500
dt = .0005 # time step
genInitialPosition(outputdir, N, r, lx, ly, rhop)
friction=0.3
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
# number of iterations between output files
......@@ -87,42 +86,21 @@ fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
fluid.set_open_boundary("Top",pressure=0, concentration=1)
fluid.set_open_boundary("Injection",velocity=[0, outerBndV], concentration=1)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_open_boundary("Injection",velocity=[0, outerBndV(1)], concentration=1)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
#set initial_condition
print(fluid.solution().shape, fluid.coordinates().shape)
fluid.solution()[fluid.coordinates()[:,1]>0.3,3] = 1
fluid.export_vtk(outputdir,0,0)
fluid.write_vtk(outputdir,0,0)
ii = 0
tic = time.time()
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
#Fluid solver
fluid.implicit_euler(dt, newton_max_it=20)
if (ii%50==0 and ii != 0):
fluid.adapt_mesh(8e-2,3e-3,10000)
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)
time_integration.iterate(fluid,p,dt,external_particles_forces=g*p.mass())
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)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......@@ -38,7 +38,7 @@ if not os.path.isdir(outputdir) :
#physical parameters
alpha = 0*np.pi/4.
g = -9.81*np.cos(alpha) # gravity
g = np.array([0,-9.81]) # gravity
rho0 = 1.117 # fluid density
rho1 = 1000
nu1 = 1.57e-5 # kinematic viscosity
......@@ -75,37 +75,20 @@ fluid.import_vtk(outputrel+"/fluid_%05d.vtu"%iReload)
rel = 1
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),rel)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
ii = (iReload-1)*outf+1
t=ii*dt
tic = time.clock()
tic = time.time()
outf1 = 25
ii1 = 0
while t < tEnd :
#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
#Fluid solver
fluid.implicit_euler(dt)
if (ii%65==0 and ii != 0):
fluid.adapt_mesh(8e-2,3e-3,10000)
rel = 0
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)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
time_integration.iterate(fluid,p,dt,external_particles_forces=g*p.mass())
#Output files writting
if ii %outf1 == 0 :
ioutput = int(ii1/outf1) + 1
print("IOUTPUT",ioutput,"TIME",t,"II", ii)
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
t += dt
ii += 1
ii1 += 1
......
......@@ -79,7 +79,7 @@ fluid.set_open_boundary("Right", pressure = 0)
#
fs = time_integration.FreeSurface(fluid,"TopRight","BottomRight",isCentered=False)
fluid.export_vtk(outputdir,0,0)
fluid.write_vtk(outputdir,0,0)
tic = time.time()
#
# COMPUTATION LOOP
......@@ -94,7 +94,7 @@ while t < tEnd :
# Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
print(f"\tSwelling = {swelling}")
\ No newline at end of file
print(f"\tSwelling = {swelling}")
......@@ -37,11 +37,11 @@ def genInitialPosition(filename, r, H, ly, lx, rin, rhop) :
rhop -- particles density
"""
# Particles structure builder
p = scontact.ParticleProblem(2)
p = scontact.ParticleProblem(2,True)
# Load mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh.msh", ["Inner","Left","Right","Bottom"],material="Steel")
x = np.arange(-lx+r, lx-r, 2*r)
y = np.arange(-ly+r, H-r, 2*r)
x = np.arange(-lx+r, lx-r, 2.1*r)
y = np.arange(-ly+r, H-r, 2.1*r)
#Definition of the points where the particles are located
x, y = np.meshgrid(x, y)
......@@ -49,7 +49,7 @@ def genInitialPosition(filename, r, H, ly, lx, rin, rhop) :
y = y.flat
# Add a grain at each centre position
for i in range(len(x)) :
r1 = np.random.uniform(0.4*r,r)
r1 = np.random.uniform(0.8*r,r)
if x[i]**2 + y[i]**2 > (rin+r1)**2:
p.add_particle((x[i], y[i]), r1, r**2 * np.pi * rhop,"Sand")
p.write_vtk(filename,0,0)
......@@ -70,16 +70,17 @@ tEnd = 10 # final time
# Geometrical parameters
ly = 0.7 # particles area height
lx = 0.15 # particles area widht
H = 2 # domain height
H = 1.3 # domain height
#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r, H, ly, lx, rin, rhop)
p = scontact.ParticleProblem(2,True,True)
p.set_friction_coefficient(0.3,"Sand","Sand")
p.set_friction_coefficient(0.3,"Sand","Steel")
p = scontact.ParticleProblem(2,True)
p.read_vtk(outputdir,0)
p.set_friction_coefficient(0.2,"Sand","Sand")
p.set_friction_coefficient(0.2,"Sand","Steel")
# Initial time and iteration
t = 0
ii = 0
......@@ -87,7 +88,7 @@ ii = 0
# COMPUTATION LOOP
#
while t < tEnd :
time_integration.iterate(None, p, dt, min_nsub=1,contact_tol = 5e-7,external_particles_forces=g*p.mass())
time_integration.iterate(None, p, dt, min_nsub=1,external_particles_forces=g*p.mass())
t += dt
# Output files writting
if ii %outf == 0 :
......
......@@ -40,4 +40,4 @@ Line Loop(20) = {20, 21, 22, 23};
Plane Surface(1) = {20, 10};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {13};
Mesh.MshFileVersion = 2;
Mesh.MshFileVersion = 4;
......@@ -43,7 +43,7 @@ dt = 1e-3 #time step
t = 0 #initial time
tEnd = 10 #final time
ii = 0 #iteration number
outf = 1000 #iterations between data frames
outf = 10 #iterations between data frames
# Define output directory
outputdir = "output/" + vitesse[index]
if not os.path.isdir(outputdir) :
......@@ -68,15 +68,15 @@ p2.write_vtk(outputdir, 0, 0)
#
# FLUID PROBLEM
#
#fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=0)
#fluid.load_msh("mesh.msh")
#fluid.set_wall_boundary("Inner")
#fluid.set_wall_boundary("Left", velocity = [0,0])
#fluid.set_wall_boundary("Right", velocity = [0,0])
#fluid.set_open_boundary("Bottom", velocity = [0,vit])
#fluid.set_open_boundary("Top",pressure = 0)
#fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces(),init=True)
#fluid.export_vtk(outputdir, 0.,0)
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=0)
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Inner")
fluid.set_wall_boundary("Left", velocity = [0,0])
fluid.set_wall_boundary("Right", velocity = [0,0])
fluid.set_open_boundary("Bottom", velocity = [0,vit])
fluid.set_open_boundary("Top",pressure = 0)
fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces())
fluid.write_vtk(outputdir, 0,0)
#
# COMPUTATION LOOP
#
......@@ -90,16 +90,16 @@ while t < tEnd :
p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.remove_particles_flag( (p2.position()[:,1] + p2.r()[:,0] >-0.55))
F = np.zeros(2)
#time_integration.particle_changed(fluid,p2)
time_integration.iterate(None,p2,dt,1,contact_tol=5e-7,external_particles_forces=g*p2.mass(),after_sub_iter=lambda n_divide : accumulate(F,n_divide))
fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces())
time_integration.iterate(fluid,p2,dt,1,contact_tol=5e-7,external_particles_forces=g*p2.mass(),after_sub_iter=lambda n_divide : accumulate(F,n_divide))
with open(filename,"a") as file1:
#F += fluid.boundary_force()[0,:]
F += fluid.boundary_force()[0,:]
file1.write( str(F[0])+";"+str(F[1])+";"+str(t)+"\n")
t += dt
# Output files writing
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p2.write_vtk(outputdir, ioutput, t)
#fluid.export_vtk(outputdir+"/Flow", t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g" % (ii, t, tEnd))
......@@ -40,4 +40,4 @@ Line Loop(20) = {20, 21, 22, 23};
Plane Surface(1) = {20, 10};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {13};
Mesh.MshFileVersion = 2;
Mesh.MshFileVersion = 4;
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment