cas test 3d

testcases/periodicTest/3D/depot.py

+# 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
+
+
+# TESTCASE DESCRIPTION
+# Deposit of 3D particles in a box full of fluid.
+
+from migflow import fluid
+from migflow import scontact
+from migflow import time_integration
+
+import numpy as np
+import os
+import time
+import shutil
+
+def genInitialPosition(p, filename, rhop) :
+    """Set all the particles centre positions and create the particles objects to add in the computing structure.
+    For this testcase, radius are contained in a file radius.csv created with a special distribution
+
+    Keyword arguments:
+    filename -- name of the output file
+    rhop -- particles density
+    """
+    # Reading the file containing the radii
+    myRadius = np.genfromtxt('radius.csv', delimiter=',')*2
+    r2 = np.amax(myRadius)*8e-6
+    i = 0
+    # Positioning the particles on a regular grid
+    for y in np.arange(0, -.05, -2*r2):
+        for z in np.arange(.01675-r2, -.01675, -2*r2):
+            for x in np.arange(0.014-r2, -.014+r2, -2*r2):
+                p.add_particle((x, y, z), myRadius[i%500000]*8e-6, 4./3.* (myRadius[i%500000]*8e-6)**3 * np.pi * rhop);
+                i += 1
+    # Translation of the particles to the top of the box
+    state = p.state()
+    state.x[:, 1] += .07/2
+    p.set_state(state)
+    print("Number of particles=%d" % len(p.state().x[:, 1]))
+    # p.write_vtk(filename, 0, 0)
+
+
+# Define output directory
+outputdir = "output"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+# Physical parameters
+g = np.array([0,-9.81,0])                        # gravity
+rho = 1000                                       # fluid density
+rhop = 2640                                      # particles density
+nu = 1e-6                                        # kinematic viscosity
+
+# Numerical parameters
+outf = 50                                        # number of iterations between output files
+dt = 1e-4                                        # time step
+tEnd = 100                                       # final time
+
+#
+# PARTICLE PROBLEM
+#
+# Particles structure creation and setting particles properties
+
+p = scontact.ParticleProblem(3)
+p.load_msh_boundaries("mesh.msh", ["Bottom", "Top", "X", "Z"])
+genInitialPosition(p,outputdir, rhop)
+print(p.n_particles())
+p.write_vtk(outputdir,0,0)
+
+# Initial time and iteration
+t         =  0
+ii        =  0
+
+#
+# FLUID PROBLEM
+#
+
+fluid = fluid.FluidProblem(3,g,[nu*rho],[rho])
+#Set the mesh geometry for the fluid computation
+fluid.load_msh("mesh.msh")
+fluid.set_wall_boundary("Top",pressure=0)
+fluid.set_wall_boundary("Bottom")
+fluid.set_wall_boundary("X")
+fluid.set_wall_boundary("Z")
+#Set location of the particles in the mesh and compute the porosity in each computation cell
+fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces(),init=True)
+
+tic = time.time()
+fluid.export_vtk(outputdir,0,0)
+G = np.zeros((p.n_particles(),p.dim()))
+G[:,1] = p.mass()[:,0]*g[1]
+#
+# COMPUTATION LOOP
+#
+while t < tEnd :
+    time_integration.iterate(fluid, p, dt,min_nsub=20, contact_tol = 1e-8,external_particles_forces=G)
+    t += dt
+    # Output files writing
+    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.time() - tic))

testcases/periodicTest/3D/mesh.geo

+SetFactory("OpenCASCADE");
+L = 0.014;
+H = 0.036;
+P = 0.01675;
+y = 0;
+lc = 1e-2;//2e-3;
+
+SetFactory("OpenCASCADE");
+Box(1) = {-L,-H,-P,2*L,2*H,2*P};
+Periodic Surface {1} = {2};
+Volume("Domain") = {1};
+Physical Surface("Top") = {4};
+Physical Surface("Bottom") = {3};
+Physical Surface("X") = {1,2};
+Physical Surface("Z") = {5,6};
+Physical Point("PtFix") = {5};
+Characteristic Length{ PointsOf{ Volume{:}; } } = lc;
+
+Mesh.MshFileVersion = 2;

testcases/periodicTest/depot.py

+# 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 mbfluid
+from migflow import time_integration
+from migflow import scontact
+
+import numpy as np
+import os
+import subprocess
+import time
+import shutil
+import random
+import unittest
+
+dir_path = os.path.dirname(os.path.realpath(__file__))
+os.chdir(dir_path)
+# Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
+
+outputdir = "outputDepotNP"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","1"])
+
+t = 0
+ii = 0
+
+
+#physical parameters
+# g = np.array([0.001,-0.001])
+# g =  np.array([0.1,-0.5])                                      # gravity
+g =  np.array([0.,-9.81])                                      # gravity
+rho = 1000                                      # fluid density
+rhop = 1500                                     # particle density
+nu = 1e-6                                   # kinematic viscosity
+mu = nu*rho                                     # dynamic viscosity
+tEnd = 5                                     # final time
+r = 1e-3                                       #particle radius
+L = 0.4
+H = 0.6
+
+#numerical parameters
+lcmin = .1                                  # mesh size
+dt = 2.5e-3                                       # time step
+alpha = 1e-4                                    # stabilization coefficient
+epsilon = alpha*lcmin**2 /nu                    # stabilization parametre
+print('epsilon',epsilon)
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+
+p = scontact.ParticleProblem(2)
+p.load_msh_boundaries("mesh.msh", ["Bottom", "Top", "Right", "Left"])
+for x in np.arange(2*r,L-2*r,2.1*r):
+   for y in np.arange(H-0.2,H-r,2.1*r):
+       R = r*(1-np.random.random()/4)
+       p.add_particle((x+R,y),R,R**2*np.pi*rhop)
+
+p.write_vtk(outputdir,0,0)
+print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
+print("RHOP = %g" % rhop)
+
+outf = 5                                       # 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)
+fluid = mbfluid.FluidProblem(2,g,nu*rho,rho,petsc_solver_type="-pc_type lu")
+fluid.load_msh("mesh.msh")
+#fluid.set_open_boundary("Right",pressure = 0)
+#fluid.set_open_boundary("Left", pressure = 0)
+fluid.set_wall_boundary("Bottom",velocity=[0,0])
+fluid.set_wall_boundary("Top",velocity=[0,0])
+
+# fluid.set_strong_boundary("Bottom",0,0)
+# if strong boundary on periodic line, it should be forced on both sides
+fluid.set_strong_boundary("Right",2,0)
+fluid.set_strong_boundary("Left",2,0)
+
+
+ii = 0
+t = 0
+
+#set initial_condition
+fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces())
+fluid.export_vtk(outputdir,0,0)
+
+G = np.zeros((p.n_particles(),2))
+G[:,1] = p.mass()[:,0]*g[1]
+
+ii = 0
+tic = time.time()
+while t < tEnd :
+    #Fluid solver
+    oldstate = np.copy(p.state())
+    time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=G)
+
+    state = p.state()
+    state.x[:,0] = np.fmod(state.x[:,0],L)
+    ind = np.where(state.x[:,0] <= 0)
+    state.x[ind,0] += L
+    p.set_state(state)
+
+    fluid.set_particles(p.mass(), p.volume(), oldstate,p.contact_forces(),init=False)
+    fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces(),init=False)
+    #Output files writting
+    t += dt
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        fluid.export_vtk(outputdir, t, ioutput)
+        p.write_vtk(outputdir, ioutput, t)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
+
+s = fluid.solution()
+x = fluid.coordinates()