Commit 6bf0ea3e authored by Matthieu Constant's avatar Matthieu Constant
Browse files

update sablier testcase

parent 476aceaa
......@@ -20,8 +20,11 @@
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
from migflow import fluid as fluid
from migflow import scontact2
# TESTCASE DESCRIPTION
# Falling particles in a sandglass
from migflow import fluid
from migflow import scontact
from migflow import lmgc2Interface
import numpy as np
......@@ -29,9 +32,10 @@ import os
import time
import shutil
import random
#grains creation and placement + physical boundaries definition
# grains creation and placement + physical boundaries definition
def genInitialPosition(filename, r, ox, oy, lx, ly, rhop) :
p = scontact2.ParticleProblem()
p = scontact.ParticleProblem(2)
#Loading of the mesh.msh file specifying physical boundaries name
p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
#Definition of the points where the grains are located
......@@ -65,13 +69,10 @@ lxpart = 0.015 # semi-width of the grains
#numerical parameters
dt = 1e-3 # time step
epsilon = 1e-4 # stabilization parameter
print("eps = %g" % ( epsilon))
#Enable the use of lmgc90
use_lmgc = True
#Object particles creation
genInitialPosition(outputdir, r, 0., y0part, lxpart, lypart, rhop)
if use_lmgc :
......@@ -79,19 +80,22 @@ if use_lmgc :
lmgc2Interface.scontactToLmgc90(outputdir, 0, friction)
p = lmgc2Interface.LmgcInterface()
else :
p = scontact2.ParticleProblem()
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
outf = 20 # 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 = [("Box",0,0),("Box",1,0),("Top",0,0.),("Top",1,0.),("Top",2,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid = fluid.FluidProblem(2,g,nu*rho,rho)
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Box",0,0)
fluid.set_strong_boundary("Box",1,0)
fluid.set_strong_boundary("Box",0,0)
fluid.set_strong_boundary("Box",1,0)
fluid.set_strong_boundary("Box",2,0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)
fluid.export_vtk(outputdir,0,0)
tic = time.time()
#Computation loop
......@@ -108,7 +112,7 @@ while t < tEnd :
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
#Contact solver
for i in range(nsub) :
tol = 1e-6
tol = 1e-8
p.iterate(dt/nsub, forces, tol)
#Localisation of the grains in the fluid
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
......@@ -119,5 +123,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.time() - tic))
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
\ No newline at end of file
Supports Markdown
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