testcase to develop periodicity

testcases/shear/mesh.geo

+L = .2;
+H = .5;
+y = 0;
+lc = 0.03;
+f = 4;
+
+Point(1) = {0, 0, 0,lc};
+Point(2) = {L, 0, 0,lc};
+Point(3) = {L, H, 0,lc};
+Point(4) = {0, H, 0,lc};
+
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+
+Line Loop(1) = {1:4};
+
+Plane Surface(1) = {1};
+
+Physical Line("Bottom") = {1};
+Physical Line("Left") = {4};
+Physical Line("Top") = {3};
+Physical Line("Right") = {2};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+
+Mesh.MshFileVersion = 2;
\ No newline at end of file

testcases/shear/shear.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
+from migflow import scontact
+from migflow import lmgc2Interface
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+# TESTCASE TO DEVELOP PERIODIC: COMING SOON
+def genInitialPosition(filename, r, lx, ly, rhop) :
+    p = scontact.ParticleProblem(2)
+    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom"])
+    dr = r/10.
+    dr = r/10.
+    x = np.arange(r+dr, lx-(r+dr), 10*(r+dr))
+    y = np.arange(r+dr, 0.66*ly, 10*(r+dr))
+    # Problem with lots of particles. It doesn't work with small time step
+    x, y = np.meshgrid(x, y)
+    x = x.flat
+    y = y.flat
+    for i in range(len(x)) :
+       d=np.random.random()
+       p.add_particle((x[i]+(d*dr), y[i]+(d*dr)), r+(d*dr), (r+(d*dr))**2 * np.pi * rhop);
+    p.write_vtk(filename,0,0)
+
+outputdir = "output"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+r = 0.002
+lx= 0.2
+ly= 0.5
+
+#physical parameters
+g =  -9.81
+rho = 1000
+nu = 1e-6
+rhop = 2500
+
+#numerical parameters
+outf = 10
+dt = 1e-3
+tEnd = 10
+
+genInitialPosition(outputdir, r, lx, ly, rhop)
+friction=0.3
+lmgc2Interface.scontactToLmgc90(outputdir, 0, friction)
+p = lmgc2Interface.LmgcInterface(period=1,xp=.2)
+
+
+t = 0
+ii = 0
+
+def shear(x) :
+    if t < 0.5 :
+        return 0.
+    else:
+        print(np.max(.05*min(t-.5,4.)*x[:,1]/ly))
+        return .05*min(t-.5,4.)*x[:,1]/ly
+
+
+fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
+fluid.load_msh("mesh.msh")
+
+fluid.set_strong_boundary("Top",2,0)
+fluid.set_strong_boundary("Top",1,0)
+fluid.set_strong_boundary("Bottom",1,0)
+fluid.set_strong_boundary("Left",0,shear)
+fluid.set_strong_boundary("Right",0,shear)
+fluid.set_strong_boundary("Right",1,0)
+fluid.set_strong_boundary("Left",1,0)
+
+fluid.set_weak_boundary("Bottom","Wall")
+fluid.set_weak_boundary("Left","Inflow")
+fluid.set_weak_boundary("Top","Wall")
+fluid.set_weak_boundary("Right","Inflow")
+
+ii = 0
+t = 0
+
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+#set initial_condition
+
+fluid.export_vtk(outputdir,0,0)
+
+ii = 0
+tic = time.time()
+while t < tEnd : 
+    #Fluid solver
+    fluid.implicit_euler(dt)
+    
+    #Computation of the new velocities
+    forces = fluid.compute_node_force(dt)
+    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())    
+    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)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))