add depot-small

testcases/depot-2d/depot-small.py

+# MigFlow - Copyright (C) <2010-2020>
+# <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
+
+# TESTCASE DESCRIPTION
+# Bidimensional particles sedimentation in fluid
+
+
+from migflow import fluid
+from migflow import scontact
+from migflow import time_integration
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+def genInitialPosition(filename, r, H, ly, lx, 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 -- max radius of the particles
+    H -- domain height
+    ly - particles area height
+    lx -- particles area width
+    rhop -- particles density
+    """
+    # Particles structure builder
+    p = scontact.ParticleProblem(2)
+    # Load mesh.msh file specifying physical boundaries names
+    p.load_msh_boundaries("mesh-small.msh", ["Top", "Lateral","Bottom"])
+
+    #Definition of the points where the particles are located
+    x = np.arange(-lx/2+r-1e-8, lx/2-r+1e-8, 2*r)
+    y = np.arange(H/2-r, H/2-ly+r, -2*r)
+    x, y = np.meshgrid(x, y)
+    x = x.flat
+    y = y.flat
+    # Add a grain at each centre position
+    for i in range(len(x)) :
+        p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop)
+    p.write_vtk(filename,0,0)
+
+# Define output directory
+outputdir = "output-small"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+# Physical parameters
+g = np.array([0,-9.81])                         # gravity
+r = 1.5e-3                                        # particles radius
+rhop = 1500                                     # particles density
+rho = 1000                                      # fluid density
+nu = 1e-6                                       # kinematic viscosity
+
+# Numerical parameters
+outf = 100                                        # number of iterations between output files
+dt = 1e-3                                     # time step
+tEnd = 100                                      # final time
+
+# Geometrical parameters
+ly = 1e-1                                       # particles area height
+lx = 2e-1                                       # particles area widht
+H = 0.2                                         # domain height
+
+#
+# PARTICLE PROBLEM
+#
+# Initialise particles
+genInitialPosition(outputdir, r, H, ly, lx, rhop)
+p = scontact.ParticleProblem(2)
+p.read_vtk(outputdir,0)
+
+print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
+print("RHOP = %g" % rhop)
+
+# Initial time and iteration
+t         =  0
+ii        =  0
+
+#
+# FLUID PROBLEM
+#
+fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=1,solver="petsc4py",solver_options="-pc_type lu",usolid=True)
+# Set the mesh geometry for the fluid computation
+fluid.load_msh("mesh-small.msh")
+fluid.set_wall_boundary("Bottom")
+fluid.set_wall_boundary("Lateral")
+fluid.set_wall_boundary("Top",pressure=0)
+# 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())
+fluid.write_vtk(outputdir,0,0)
+p.write_vtk(outputdir, 0, 0)
+
+tic = time.time()
+#
+# COMPUTATION LOOP
+#
+while t < tEnd :
+    time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=g*p.mass(), use_predictor_corrector=True)
+    t += dt
+
+    # Output files writting
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        p.write_vtk(outputdir, ioutput, t)
+        fluid.write_vtk(outputdir, ioutput, t)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))

testcases/depot-2d/mesh-small.geo

+L = 0.1;
+H = 0.1;
+y = 0;
+lc = 0.01;
+
+Point(1) = {-L, H, 0,lc};
+Point(2) = {-L, -H, 0,lc};
+Point(3) = {L, -H, 0,lc};
+Point(4) = {L, H, 0,lc};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Point(5) = {0,0.52,0};
+
+Line Loop(1) = {1:4};
+
+Plane Surface(1) = {1};
+Physical Line("Bottom") = {2};
+Physical Line("Lateral") = {1,3};
+Physical Line("Top") = {4};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+
+Mesh.CharacteristicLengthFromPoints = 1;
+Mesh.MshFileVersion = 2;