example of reload

testcases/drop-2d/SimpleDrop/dropReload.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
+# This tescase presents the fall of a cloud made of particles in a viscous fluid (Stokes cloud).
+# Physical parameters for the drops are the ones presented by Metzger et al. (2007) 
+# "Falling clouds of particles in viscous fluids"
+
+from migflow import fluid
+from migflow import scontact
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+def genInitialPosition(filename, r, rout, rhop, compacity) :
+    """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
+    rout -- outer radius of the cloud
+    rhop -- particles density        
+    compacity -- initial compacity in the cloud
+    """
+    # Particles structure builder
+    p = scontact.ParticleProblem(2)
+    # Load mesh.msh file specifying physical boundaries names
+    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral"])
+
+    # Space between the particles to obtain the expected compacity
+    N = compacity*4*rout**2/(np.pi*r**2)
+    e = 2*rout/(N)**(1./2.)
+
+    # Definition of the points where the particles are located
+    x = np.arange(rout, -rout, -e)
+    x, y = np.meshgrid(x, x)
+    R2 = x**2 + y**2
+    keep = R2 < (rout - r)**2
+    x = x[keep]
+    y = y[keep]
+
+    for i in range(x.shape[0]) :
+        # Add particles object with properties defined above
+        p.add_particle((x[i]+2*random.uniform(-e/2+r,e/2-r), y[i]+2*random.uniform(-e/2+r,e/2-r)), r, r**2 * np.pi * rhop)
+    
+    # Vertical shift of the particles to the top of the box
+    p.position()[:, 1] += 0.52
+    p.write_vtk(filename,0,0)
+
+# Define output directory
+outputdir = "output1"
+outputrel = "output"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+# Physical parameters
+g = -9.81                                       # gravity
+rhop = 2450                                     # particles density
+r = 154e-6                                      # particles radii
+compacity = 0.20                                # solid volume fraction in the drop
+rho = 1030                                      # fluid density
+nu = 1.17/rho                                   # kinematic viscosity
+rout = 3.3e-3                                   # cloud radius
+mu = nu*rho                                     # dynamic viscosity
+print("RHOP = %g, r = %g, RHO = %g" % (rhop,r,rho))
+
+# Numerical parameters
+outf = 1                                        # number of iterations between output files
+dt = 5e-2                                       # time step
+tEnd = 100                                      # final time
+
+#
+# PARTICLE PROBLEM
+#
+# Reload particles
+iReload = 10                                    # iteration file reloaded
+p = scontact.ParticleProblem(2)
+p.read_vtk(outputrel,iReload)
+
+# Initial time and iteration
+t = 0
+ii = 0
+jj = 0
+tEnd1 = 0.2
+
+#
+# FLUID PROBLEM
+#
+fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
+# Set the mesh geometry for the fluid computation
+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("Lateral",0,0)
+fluid.set_weak_boundary("Top","Wall")
+fluid.set_weak_boundary("Bottom","Wall")
+fluid.set_weak_boundary("Lateral","Wall")
+fluid.import_vtk(outputrel+"/fluid_%05d.vtu"%iReload)
+
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),reload=1)
+ii = (iReload-1)*outf+1
+t=ii*dt
+
+#
+# COMPUTATION LOOP
+#
+tic = time.time()
+
+while t < tEnd : 
+    # Fluid solver
+    fluid.implicit_euler(dt)
+
+    # Adaptation of the mesh.
+    if (ii%15==0 and ii != 0):
+       fluid.adapt_mesh(5e-3,8e-4,50000)
+    
+    # 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()))
+    # NLGS iterations
+    for i in range(nsub) :
+        p.iterate(dt/nsub, forces)  
+
+    t += dt
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+
+    # 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))
\ No newline at end of file