Commit da30176e authored by Margaux Boxho's avatar Margaux Boxho
Browse files

Commit for merge with master branch (2)

parent b84f02eb
Pipeline #5370 passed with stage
in 27 seconds
# 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
# Example of reload
# This tescase presents the fall of a cloud made of particles in a viscous fluid (Stokes cloud) reloaded from drop.py test case.
# 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
# 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],petsc_solver_type="-pc_type ilu -ksp_max_it 30")
# 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("Lateral")
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
# 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, shift) :
"""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
shift -- vertical shift of the particles to the top of the box
"""
# 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], y[i]), r, r**2 * np.pi * rhop)
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] += shift
p.write_vtk(filename,0,0)
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Enable computation of the initial velocity
init = False
# Physical parameters
g = -9.81 # gravity
rhop = 2450 # particles density
r = 154e-6 # particles radii
compacity = 0.25 # solid volume fraction in the drop
rho = 1030 # fluid density
nu = 1.17/rho # kinematic viscosity
rout = 3.3e-3 #1.1e-2 #1.1e-2; # cloud radius
mu = nu*rho # dynamic viscosity
shift = -0.05 # vertical shift of the particles to the top (depend on the mesh.geo file)
print("RHOP = %g, r = %g, RHO = %g" % (rhop,r,rho))
e_target = 4e-8
H = 0.3
L = 0.05
# Numerical parameters
outf = 5 # number of iterations between output files
dt = 5e-3 # time step
tEnd = 20 # final time
#
# PARTICLE PROBLEM
#
# Initialise particles
#genInitialPosition(outputdir, r, rout, rhop, compacity, shift)
p = scontact.ParticleProblem(2)
p.read_vtk("/Users/margauxboxho/migflow/testcases/drop-2d/Reference_solution/output_refined_mesh",0)
list_of_position = p.position()
N_0 = len(list_of_position[:,1])
print(" ======= The initial number of grains is %i [-]======= " %N_0)
Vol_p = 4*np.pi*r**3/3
F = Vol_p*(rhop-rho)*np.abs(g)
U_0 = F/(5*np.pi*mu*rout)*N_0
print(" ======= The initial velocity is %.5f [m/s]=======" %U_0)
P_0 = mu*U_0/rout
print(" ======= The initial pressure is %.5f [Pa]=======" %P_0)
Re = rhop*U_0*rout/mu
print(" ======= The Reynolds number is %.3f [Pa]=======" %Re)
# Initial time and iteration
t = 0
ii = 0
jj = 0
tEnd1 = 0.2
delta = 5
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,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("Lateral")
# 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())
# Initial fluid velocities are computed based on the forces applied by the particles without moving the particles
# The initial mesh is obtained by adapting for times the mesh on the converged fields
if init:
for jj in range(4):
dt1 = dt/100
t = 0
if jj!=0:
# Adaptation of the mesh.
#fluid.adapt_mesh(5e-3,8e-4,50000)
fluid.adapt_gen_mesh(30*r, 10*r, 6602*2, "adapt/lc.pos")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
while t < tEnd1 :
# Fluid solver
niter = fluid.implicit_euler(dt1)
forces = fluid.compute_node_force(dt1)
# Changes the velocities of the particles without moving the particles
p.velocity()[:] += dt1*forces/p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
# Increase time step if convergence
if niter < 5 and dt1 < 20:
dt1 *= 3
elif niter > 5 and dt1 > dt/100:
dt1 /= 3
print(dt1)
t += dt1
t = 0
fluid.export_vtk(outputdir,0,0)
tic = time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
#
# COMPUTATION LOOP
#
while t < tEnd :
# Save the state of the fluid at certain iteration
if (ii==jj*delta):
print("Save state at %i" %(ii))
fluid.export_vtk(outputdir,0,0)
p.write_vtk(outputdir,0,0)
# Fluid solver
fluid.implicit_euler(dt)
if (ii==0):
# initialise the h_target array
fluid.adapt_mesh(60*r, 10*r, 6602*2, e_target, 1, U_0, P_0)
if (ii%(delta*(jj+1))!=0):
# update the h_target array using the recovery gradient estimator
fluid.adapt_mesh(60*r, 10*r, 6602*2, e_target, 0, U_0, P_0)
# Adaptation of the mesh and loading the past save state
if (ii%(delta*(jj+1))==0 and ii != 0):
print("Adapt mesh at %i" %(ii))
# Adapt the mesh by creating the lc.pos file
fluid.adapt_gen_mesh(60*r, 10*r, 6602*2, e_target, U_0, P_0)
ii = ii-delta
t = t-dt*delta
jj += 1
print("ii=%i, t=%i, jj=%i" %(ii,t,jj))
# load the save past state
fluid.import_vtk("output/fluid_00000.vtu")
p.read_vtk(outputdir,0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
# create the .msh on the adapt/lc.pos and transfer the solution
fluid.adapt_load_mesh("adapt/mesh.msh")
#fluid.implicit_euler(dt)
print("Restart computation %i iterations back at %i" %(delta, ii))
# Initialize the min_h_target attribut of the Mesh structure
print('P_0 %.5f ' %P_0)
fluid.adapt_mesh(60*r, 10*r, 6602*2, e_target, 1, U_0, P_0)
# 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(),p.contact_forces())
# 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))
This source diff could not be displayed because it is too large. You can view the blob instead.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment