màj cas tests

testcases/couette-2d/couette.py

@@ -103,14 +103,14 @@ r = 397e-6/2                                   # grains radius
 # Object particles creation
 genInitialPosition(outputdir, r, rout, rin, rhop)
 if use_lmgc90 :
-    friction = 0.1                                    # friction coefficient
+    friction = 0.1                             # friction coefficient
     lmgc90Interface.scontactTolmgc90(outputdir,2,0,friction)
     p = lmgc90Interface.ParticleProblem(2)
 else :
   p = scontact.ParticleProblem(2,True)
   p.read_vtk(outputdir,0)
-  p.set_friction_coefficient(0.1,"Sand","Sand")#Particle-Particle
-  p.set_friction_coefficient(0.1,"Sand","Steel")#Particle-Wall
+  p.set_friction_coefficient(0.1,"Sand","Sand") # Particle-Particle
+  p.set_friction_coefficient(0.1,"Sand","Steel")# Particle-Wall
 
 # Initial time and iteration
 t = 0

testcases/depot-2d/depot.py

@@ -76,7 +76,7 @@ nu = 1e-6                                       # kinematic viscosity
 
 # Numerical parameters
 outf = 5                                        # number of iterations between output files
-dt = 0.25e-2                                       # time step
+dt = 2.5e-3                                     # time step
 tEnd = 100                                      # final time
 
 # Geometrical parameters

testcases/depot-2d/mesh.geo

-L = 0.1;
+L = 0.2;
 H = 0.3;
 y = 0;
 lc = 0.01;
 
 Point(1) = {-L, H, 0,lc};
-Point(2) = {-L, -H, 0,lc};
-Point(3) = {L, -H, 0,lc};
+Point(2) = {-L, .2, 0,lc};
+Point(3) = {L, .2, 0,lc};
 Point(4) = {L, H, 0,lc};
 Line(1) = {1, 2};
 Line(2) = {2, 3};

testcases/drop-2d/InteractingDrops/Diag/2DiagDrops.py

@@ -95,7 +95,7 @@ H = 0.5
 print("RHOP = %g, r = %g, RHO = %g" % (rhop,r,rho))
 
 # Numerical parameters
-outf = 10                                        # number of iterations between output files
+outf = 10                                       # number of iterations between output files
 dt = 5e-3                                       # time step
 tEnd = 100                                      # final time
 

testcases/drop-2d/InteractingDrops/Vert/2VertDrops.py

@@ -90,12 +90,12 @@ nu = 1e-4                                        # kinematic viscosity
 drho = 15.6                                      # density differences rhop-rho
 rho = rhop-drho/compacity                        # fluid density
 rout = 2e-3                                      # cloud radius
-mu = nu*rho    
-H = 0.5                                  # dynamic viscosity
+mu = nu*rho                                      # dynamic viscosity
+H = 0.5                                  
 print("RHOP = %g, r = %g, RHO = %g" % (rhop,r,rho))
 
 # Numerical parameters
-outf = 10                                         # number of iterations between output files
+outf = 10                                        # number of iterations between output files
 dt = 5e-3                                        # time step
 tEnd = 100                                       # final time
 

testcases/drop-2d/SimpleDrop/drop.py

@@ -79,7 +79,7 @@ if not os.path.isdir(outputdir) :
 # Physical parameters
 g = -9.81                                       # gravity
 rhop = 2450                                     # particles density
-r = 154e-6/2.5                                      # particles radii
+r = 154e-6/2.5                                  # particles radii
 compacity = 0.20                                # solid volume fraction in the drop
 rho = 1030                                      # fluid density
 nu = 1.17/rho                                   # kinematic viscosity
@@ -113,7 +113,6 @@ fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=1)
 # Set the mesh geometry for the fluid computation
 fluid.load_msh("mesh.msh")
 fluid.set_wall_boundary("Top", pressure=0)
-#fluid.set_wall_boundary("Top")
 fluid.set_wall_boundary("Bottom")
 fluid.set_wall_boundary("Lateral")
 # Set location of the particles in the mesh and compute the porosity in each computation cell

testcases/mill-2d/mill.py

@@ -33,12 +33,17 @@ outputdir = "output"
 if not os.path.isdir(outputdir) :
     os.makedirs(outputdir)
 
-#grains creation and placement + physical boundaries definition
-'''r is the radius of the grains
-rout is the outer radius of the geometry
-rin is the inner radius of the geometry
-rhop is the particles density'''
+
 def genInitialPosition(filename, r, rout, rhop) :
+    """Set all the particles centre positions and create the particles objects to add in the computing structure
+
+    Keyword arguments:
+    filename -- name of the output file
+    r -- radius of the particles
+    rout -- the outer radius of the geometry
+    rin -- the inner radius of the geometry
+    rhop -- the particles density
+    """
     p = scontact.ParticleProblem(2)
     #Loading of the mesh.msh file specifying physical boundaries name and material
     p.load_msh_boundaries("mesh.msh", ["Outer","Top",],material="PVC")
@@ -73,22 +78,26 @@ dt = 1e-3                               # time step
 #geometry parameters
 rout = 0.05                             # outer radius
 r  = 1e-3                               # grains radius
-Fr = 8e-3				        #Froude number
-v = 0.0715          #Corresponding tangent velocity
+Fr = 8e-3				                # Froude number
+v = 0.0715                              # Corresponding tangent velocity
 
 shutil.copy("mesh.msh", outputdir +"/mesh.msh")
 
-#Object particles creation
-#genInitialPosition(outputdir, r, rout, rhop)
+#
+# PARTICLE PROBLEM
+#
+# Initialise particles
+genInitialPosition(outputdir, r, rout, rhop)
 p = scontact.ParticleProblem(2,True)
-p.read_vtk("depot/output",4000)
-outf = 250                                # number of iterations between output files
+p.read_vtk(outputdir,0)
+outf = 250                                       # number of iterations between output files
 #Taking friction into account
-p.set_friction_coefficient(.4,"Glass","Glass")#Between particles
-p.set_friction_coefficient(.5,"Glass","PVC")#Between particles and boundaries
+p.set_friction_coefficient(.4,"Glass","Glass")   # between particles
+p.set_friction_coefficient(.5,"Glass","PVC")     # between particles and boundaries
+
+p.set_tangent_boundary_velocity("Outer",-v)      # setting the tangent velocity of the rotating mill
+p.set_tangent_boundary_velocity("Top",-v)        # setting the tangent velocity of the rotating mill
 
-p.set_tangent_boundary_velocity("Outer",-v)#Setting the tangent velocity of the rotating mill
-p.set_tangent_boundary_velocity("Top",-v)#Setting the tangent velocity of the rotating mill
 #
 # FLUID PROBLEM
 #
@@ -106,7 +115,10 @@ fluid.export_vtk(outputdir, 0.,0)
 
 G = np.zeros_like(p.velocity())
 G[:,1] = p.mass()[:,0]*g
-#Computation loop
+
+#
+# COMPUTATION LOOP
+#
 while t < tEnd :
     time_integration.iterate(fluid, p, dt, min_nsub=15, external_particles_forces=G)
     t += dt

testcases/shaker/shaker.py

@@ -25,6 +25,7 @@
 # Sort particles with respect to their mass by jigging the granular matrix
 from migflow import fluid
 from migflow import scontact
+from migflow import time_integration
 
 import numpy as np
 import os
@@ -33,18 +34,26 @@ import shutil
 import random
 
 def genInitialPosition(filename, r, ly, rhop) :
+    """Set all the particles centre positions and create the particles objects to add in the computing structure
+    
+    Keyword arguments:    
+    filename -- name of the output file
+    r -- radius of the particles
+    ly - particles area height
+    rhop -- particles density        
+    """
     p = scontact.ParticleProblem(2)
-    #Loading of the mesh.msh file specifying physical boundaries name
+    # Loading of the mesh.msh file specifying physical boundaries name
     p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom","BottomOut","TopOut"])
-    #Definition of the points where the grains are located
-    x = np.arange(-0.5+r, 0.5-r, 2*r)
-    y = np.arange(-ly/2+r, ly/2-r, 2*r)
+    # Definition of the points where the grains are located
+    x = np.arange(-0.5+r, 0.5-r+1e-10, 2*r)
+    y = np.arange(-ly/2+r, 0.5*ly/2-r, 2*r)
     x, y = np.meshgrid(x, y)
     x = x.flat
     y = y.flat
     for i in range(len(x)) :
         rhop = random.choice([1100,1350,1600])
-            #Addition of an particle object at each point
+        # Addition of a particle object at each point
         p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
     p.write_vtk(filename,0,0)
 
@@ -52,65 +61,63 @@ outputdir = "output"
 if not os.path.isdir(outputdir) :
     os.makedirs(outputdir)
 
-#physical parameters
+# physical parameters
 g =  -9.81                                  # gravity
 rho = 1000                                  # fluid density
 rhop = 1500                                 # grains density
 nu = 1e-6                                   # kinematic viscosity
-tEnd = 100                                  # final time
 r = 5e-3                                    # grains radius
 ly = 1                                      # box height
 
-#numerical parameters
-dt = 1e-3                                   # time step
+# numerical parameters
+dt = 5e-3                                   # time step
+tEnd = 100                                  # final time
+outf = 10                                   # number of iterations between output files
 
-#Object particles creation
+#
+# PARTICLE PROBLEM
+#
+# Initialise particles
 genInitialPosition(outputdir, r, ly, rhop)
 p = scontact.ParticleProblem(2)
 p.read_vtk(outputdir,0)
 
-print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
-outf = 10                                   # number of iterations between output files
+print("r = %g, m = %g" %  (p.r()[0], p.mass()[0]))
+print("RHOP = %g" % rhop)
+print("number of grains = %d"%len(p.r()))
 
-ii = 0
-
-# periodic boundary condition function
-def outerBndV(x) :
-    return 0.3 * max(np.sin(-t*np.pi*2./5-4*np.pi/5),0)
+# Initial time and iteration
+t         =  0
+ii        =  0
 
+#
+# FLUID PROBLEM
+#
 fluid = fluid.FluidProblem(2,g,nu*rho,rho)
+# Set the mesh geometry for the fluid computation
 fluid.load_msh("mesh.msh")
+# periodic boundary condition function
+def outerBndV(x) :
+    return 0.15*max(np.sin(-t*np.pi*2./5-4*np.pi/5),0)
 fluid.set_wall_boundary("Lateral")
 fluid.set_open_boundary("Bottom",velocity=[0, outerBndV])
 fluid.set_open_boundary("BottomOut",velocity=[0, outerBndV])
 fluid.set_open_boundary("TopOut",pressure=0)
 fluid.set_open_boundary("Top",pressure=0)
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+# Set location of the particles in the mesh and compute the porosity in each computation cell
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces(), init=True)
 fluid.export_vtk(outputdir,0,0)
 
-t = 0
-ii = 0
 tic = time.time()
-#Computation loop
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+G = np.zeros_like(p.velocity())
+G[:,1] = p.mass()[:,0]*g
+#
+# COMPUTATION LOOP
+#
 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(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) :
-        tol = 1e-7
-        p.iterate(dt/nsub, forces, tol)
-    #Localisation of the grains in the fluid
-    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    time_integration.iterate(fluid, p, dt, min_nsub=25, external_particles_forces=G)
     t += dt
-    #Output files writting
+    # Output files writting
     if ii %outf == 0 :
         ioutput = int(ii/outf) + 1
         p.write_vtk(outputdir, ioutput, t)