Commit f8aad78e authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

Merge branch 'oneFluid' into debugOneFluid

parents 598361bb 3aa1bb9a
......@@ -25,6 +25,7 @@ mbtests :
- mkdir build
- cd build/
- cmake .. -DENABLE_PETSC=1 -DPETSC_DIR=/usr/lib/petscdir/petsc3.9/ -DPETSC_ARCH=x86_64-linux-gnu-real
- PPATH=$(pwd)
- make -j4
- cd ../validation
- ../build/migflow.sh mbtests.py
- PYTHONPATH=$PPATH:$PYTHONPATH python3 mbtests.py
......@@ -40,7 +40,7 @@ if(ENABLE_PETSC)
find_package(PETSc)
if(PETSC_FOUND)
add_definitions("-DHXT_HAVE_PETSC ${PETSC_CFLAGS}")
link_directories("$PETSC_LIBDIR")
link_directories("${PETSC_LIBDIR}")
list(APPEND EXTERNAL_LIBRARIES ${PETSC_LIBS})
list(APPEND EXTERNAL_INCLUDES ${PETSC_INC})
endif(PETSC_FOUND)
......
......@@ -70,17 +70,15 @@ if(PETSC_VARIABLES)
${PETSC_PACKAGES_INCLUDES})
endif(PETSC_PACKAGES_INCLUDES)
if(PETSC_PACKAGES_INCLUDES)
if(PETSC_PACKAGES_INCLUDES)
string(REPLACE "-I" "" PETSC_PACKAGES_INCLUDES ${PETSC_PACKAGES_INCLUDES})
string(REPLACE " " ";" PETSC_PACKAGES_INCLUDES ${PETSC_PACKAGES_INCLUDES})
foreach(VAR ${PETSC_PACKAGES_INCLUDES})
# seem to include unexisting directories (/usr/include/lib64)
# check to avoid warnings
if(EXISTS ${VAR})
list(APPEND PETSC_INC ${VAR})
# check to avoid warnings
if(EXISTS ${VAR})
list(APPEND PETSC_INC ${VAR})
endif(EXISTS ${VAR})
endforeach(VAR)
endif(PETSC_PACKAGES_INCLUDES)
endif(PETSC_PACKAGES_INCLUDES)
string(REGEX MATCHALL "\nCC = [^\n\r]*" PETSC_CC ${PETSC_VARIABLES})
if(PETSC_CC)
......
......@@ -53,7 +53,7 @@ class FluidProblem :
g -- gravity (scalar)
viscosity -- list of fluid phases viscosities (this should be a vector whom dimension specify if there is one or two fluid)
density -- list of fluid phases densities (this should be a vector whom dimension specify if there is one or two fluid)
coeff_stab -- optional argument used as coefficient for axtra diffusivity added to stabilise the advection of concentration (only for two fluid problem)
coeff_stab -- optional argument used as coefficient for extra diffusivity added to stabilise the advection of concentration (only for two fluid problem)
petsc_solver_type -- optional argument to specify solver option (only when petsc is used)
Raises:
......
......@@ -23,8 +23,9 @@ from pylmgc90 import chipy
import numpy as np
import os
import math
import sys
class Lmgc90Interface(object):
class ParticleProblem(object):
"""
A class definition holding data of lmgc90
in a compatible form to use coupling with
......@@ -51,7 +52,7 @@ class Lmgc90Interface(object):
def __init__(self,dim,period=0,xp=0,yp=0) :
"""
Initialize LMGC90
Initialize lmgc90
"""
self._dim=dim
......@@ -143,9 +144,9 @@ class Lmgc90Interface(object):
tag = chipy.RBDY3_GetContactorColor(i+1,1).rstrip('_')
self._tag2id.setdefault(tag, []).append(i)
else :
chipy.utilities_logMes('dim must be 2 or 3')
sys.exit()
else :
chipy.utilities_logMes('dim must be 2 or 3')
sys.exit()
self._freq_detect = 1
self._solver_params = { 'type' :'Stored_Delassus_Loops ',
......@@ -167,7 +168,7 @@ class Lmgc90Interface(object):
def iterate(self, dt, forces, tol=None, gsit1=None, gsit2=None):
"""
Do one step of a LMGC90 computation.
Do one step of a lmgc90 computation.
"""
if tol is not None : self._solver_params['conv'] = tol
......@@ -269,7 +270,7 @@ class Lmgc90Interface(object):
chipy.CloseDisplayFiles()
chipy.Finalize()
def scontactToLmgc90(dirname, dim, it=0, fric = 0., assume_box=False):
def scontactTolmgc90(dirname, dim, it=0, fric = 0., assume_box=False):
from migflow import scontact
from pylmgc90 import pre
datbox_path = 'DATBOX'
......
......@@ -29,7 +29,7 @@
from migflow import fluid
from migflow import scontact
from migflow import lmgc2Interface
from migflow import lmgc90Interface
from pylmgc90 import pre
import numpy as np
import os
......@@ -103,8 +103,8 @@ tEnd = 10 #final time
# Initialise particles
genInitialPosition(outputdir, N, r, lx, ly, rhop)
friction = 0.3 #friction coefficient
lmgc2Interface.scontactToLmgc90(outputdir, 0, friction)
p = lmgc2Interface.LmgcInterface()
lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
p = lmgc90Interface.ParticleProblem(2)
# Initial time and iteration
t = 0
......
# 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
# Mixing of grains in a rotating drum.
# This example shows how to set boundary conditions as a function of some parametres.
# A boolean parametre give the possibility to use lmgc90 to solve contacts or not.
#variable to use lmgc or not
use_lmgc = 1
from migflow import fluid
from migflow import scontact
if use_lmgc :
from migflow import lmgc90Interface
import numpy as np
import os
import time
import shutil
def genInitialPosition(filename, r, rout, rhop) :
"""Create a container to initialize particles and write particles in an initial output file.
Keyword arguments:
filename -- directory to write the output file
r -- max particle radius
rout -- outer radius of the drum
rhop -- particle density
"""
# Grains structure builder
p = scontact.ParticleProblem(2)
# Load mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
# First layer of grains
# Definition of the points where the grains are located
x = np.arange(rout, -rout, 2 * -r)
y = np.arange(-rout, -2*rout/3., 2 * r)
x, y = np.meshgrid(x, y)
R2 = x**2 + y**2
# Condition to be inside the outer boundary
keep = R2 < (rout-r)**2
x = x[keep]
y = y[keep]
# Add a grain at each centre position
for i in range(x.shape[0]) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
# Second layer of grains
# Definition of the points where the grains are located
x = np.arange(rout, -rout, -r)
y = np.arange(-2*rout/3.+r, -rout/2., r)
x, y = np.meshgrid(x, y)
# Condition to be inside the outer boundary
R2 = x**2 + y**2
keep = R2 < (rout-r)**2
x = x[keep]
y = y[keep]
# Add a grain at each centre position
for i in range(x.shape[0]) :
p.add_particle((x[i], y[i]), r/2., r**2/4. * np.pi * rhop);
# Third layer of grains
# Definition of the points where the grains are located
x = np.arange(rout, -rout, 2 * -r/3)
y = np.arange(-rout/2.+r, -5*rout/12., 2 * r/3)
x, y = np.meshgrid(x, y)
R2 = x**2 + y**2
# Condition to be inside the outer boundary
keep = R2 < (rout-r/3)**2
x = x[keep]
y = y[keep]
# Add a grain at each centre position
for i in range(x.shape[0]) :
p.add_particle((x[i], y[i]), r/3, r**2 /9*np.pi*rhop);
p.write_vtk(filename,0,0)
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = -9.81 # gravity
rhop = 1600 # grains density
r = 397e-6/2 # grains radius
rho = 1e3 # fluid density
nu = 1e-4 # kinematic viscosity
# Numerical parameters
outf = 1 # number of iterations between output files
dt = 5e-3 # time step
tEnd = 100 # final time
# Geometrical parametres
rout = 0.0254 # outer boundary
#
# PARTICLE PROBLEM
#
# Initialise grains
genInitialPosition(outputdir, r, rout, rhop)
if use_lmgc :
friction = 0.1 # friction coefficient
lmgc90Interface.scontactToLmgc90(outputdir,2,0,friction)
p = lmgc90Interface.Lmgc90Interface(2)
else :
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
# Initial time and iteration
t = 0
ii = 0
#
# 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_strong_boundary("Outer",0,lambda x : -4*x[:, 1])
fluid.set_strong_boundary("Outer",1,lambda x : 4*x[:, 0])
fluid.set_strong_boundary("Inner",0,0)
fluid.set_strong_boundary("Inner",1,0)
fluid.set_strong_boundary("PtFix",2,0)
fluid.set_weak_boundary("Inner","Wall")
fluid.set_weak_boundary("Outer","Wall")
#Set location of the grains in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export_vtk(outputdir,0,0)
tic = time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#
# COMPUTATION LOOP
#
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 contact 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) :
tol = 1e-6 #numerical tolerance for grains intersection
p.iterate(dt/nsub, forces, tol)
# Localisation of the grains in the fluid
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))
......@@ -22,11 +22,16 @@
#!/usr/bin/env python
# TESTCASE DESCRIPTION
# Grains in circular shear flow.
# This example shows how to set boundary conditions as a function of some parametres.
# Grains in circular shear flow generated by a rotating drum.
# This example shows how to set boundary conditions as a function of some parameters.
# A boolean parameter gives the possibility to use lmgc90 to solve contacts instead of scontact.
use_lmgc90 = 0
from migflow import fluid
from migflow import scontact
if use_lmgc90 :
from migflow import lmgc90Interface
import numpy as np
import os
......@@ -34,11 +39,6 @@ import time
import shutil
import random
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
def genInitialPosition(filename, r, rout, rin, rhop) :
"""Create a container to initialize particles and write particles in an initial output file.
......@@ -73,6 +73,13 @@ def genInitialPosition(filename, r, rout, rin, rhop) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop1)
p.write_vtk(filename,0,0)
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = 0 # gravity
rho = 1.253e3 # fluid density
......@@ -94,8 +101,13 @@ r = 397e-6/2 # grains radius
#
# Object particles creation
genInitialPosition(outputdir, r, rout, rin, rhop)
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
if use_lmgc90 :
friction = 0.1 # friction coefficient
lmgc90Interface.scontactTolmgc90(outputdir,2,0,friction)
p = lmgc90Interface.ParticleProblem(2)
else :
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
# Initial time and iteration
t = 0
......
......@@ -3,6 +3,7 @@ L = 0.014;
H = 0.076;
P = 0.01675;
y = 0;
lc= 2e-3;
SetFactory("OpenCASCADE");
Box(1) = {-L,-H,-P,2*L,2*H,2*P};
......
......@@ -20,9 +20,9 @@
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
from marblesbag import fluid3 as fluid
from marblesbag import scontact3
from marblesbag import lmgc3Interface
from migflow import fluid
from migflow import scontact
from migflow import lmgc90Interface
import numpy as np
import os
......@@ -31,7 +31,7 @@ import shutil
import random
def genInitialPosition(filename, rM, rm, ly, lxz, rhop) :
p = scontact3.ParticleProblem()
p = scontact.ParticleProblem(3)
p.load_msh_boundaries("mesh.msh", ["Top","Bottom","Z","X"])
for x in np.arange(-lxz+rM, lxz-rM, 2*rM):
for y in np.arange(0.015+rM, 0.015+ly-rM, 2*rM):
......@@ -50,7 +50,7 @@ g = -9.81
rho = 1000
rhop = 2500
nu = 1e-8
V = 0.5 # todo : estimate V base on limit velocity
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 5.
......@@ -66,22 +66,28 @@ lxz = 0.025
genInitialPosition(outputdir, rM, rm, ly, lxz, rhop)
friction = 0.1
lmgc3Interface.scontactToLmgc90(outputdir,0,friction,assume_box=True)
p = lmgc3Interface.LmgcInterface()
lmgc90Interface.scontactTolmgc90(outputdir,3,0,friction,assume_box=True)
p = lmgc90Interface.ParticleProblem(3)
#numerical parameters
lcmin = 0.001 # mesh size
dt = 1e-4 # time step
alpha = 2.5e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
fluid = fluid.FluidProblem(3,g,[nu*rho],[rho])
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",0,0.)
fluid.set_strong_boundary("Top",1,0.)
fluid.set_strong_boundary("Top",2,0.)
fluid.set_strong_boundary("X",1,0.)
fluid.set_strong_boundary("X",2,0.)
fluid.set_strong_boundary("X",0,0.)
fluid.set_strong_boundary("Z",1,0.)
fluid.set_strong_boundary("Z",2,0.)
fluid.set_strong_boundary("Z",0,0.)
fluid.set_strong_boundary("Bottom",0,0.)
fluid.set_strong_boundary("Bottom",1,0.)
fluid.set_strong_boundary("Bottom",2,0.)
fluid.set_strong_boundary("PtFix",3,0.)
strong_boundaries = [("Top",0,0.),("Top",1,0.),("Top",2,0.),("PtFix",3,0.),
("X",0,0.),("X",1,0.),("X",2,0.),
("Z",0,0.),("Z",1,0.),("Z",2,0.),
("Bottom",0,0.),("Bottom",1,0.),("Bottom",2,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export_vtk(outputdir,0,0)
......
......@@ -26,7 +26,6 @@
from migflow import fluid
from migflow import scontact
from migflow import lmgc2Interface
from pylmgc90 import pre
import numpy as np
import os
......@@ -34,9 +33,6 @@ import time
import shutil
import random
outputdir = "outputVid"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
def genInitialPosition(filename, N, r, lx, ly, rhop) :
p = scontact.ParticleProblem(2)
......@@ -51,23 +47,29 @@ def genInitialPosition(filename, N, r, lx, ly, rhop) :
p.add_particle(coors[2*i:2*i+2], radii[i], radii[i]**2 * np.pi * rhop);
p.write_vtk(filename,0,0)
outputdir = "outputVid"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# physical parameters
alpha = 0*np.pi/4.
g = -9.81*np.cos(alpha) # gravity
rho0 = 1.117 # fluid density
g = -9.81*np.cos(alpha) # gravity
rho0 = 1.117 # fluid density
rho1 = 1000
nu1 = 1.57e-5 # kinematic viscosity
nu1 = 1.57e-5 # kinematic viscosity
nu0 = 1e-6
tEnd = 10 # final time
tEnd = 10 # final time
r = 5e-4
N = 10000*3
lx = .222
ly = .14*2
rhop = 2500
# numerical parameters
dt = .0001 # time step
dt = .0001 # time step
genInitialPosition(outputdir, N, r, lx, ly, rhop)
friction=0.3
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
......
#!/usr/bin/env python
use_lmgc = False #True
use_lmgc90 = False #True
from migflow import fluid
from migflow import scontact
if use_lmgc :
from migflow import lmgc2Interface
if use_lmgc90 :
from migflow import lmgc90Interface
import numpy as np
import os
......@@ -51,10 +51,11 @@ lypart = 150*r
lxpart = 50*r
genInitialPosition(outputdir, r, x0part, y0part, lxpart, lypart, rhop)
if use_lmgc :
friction = 0.1 # friction coefficient
lmgc2Interface.scontactToLmgc90(outputdir,0,friction)
p = lmgc2Interface.LmgcInterface()
if use_lmgc90 :
# friction coefficient
friction = 0.1
lmgc90Interface.scontactTolmgc90(outputdir,2,0,friction)
p = lmgc90Interface.ParticleProblem(2)
else :
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
......@@ -66,8 +67,6 @@ print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 20 #000
fluid = fluid.FluidProblem(2,g,nu*rho,rho)
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Box",0,0)
......@@ -75,7 +74,7 @@ fluid.set_strong_boundary("Box",1,0)
fluid.set_strong_boundary("Left",0,0.001)
fluid.set_strong_boundary("Right",2,0)
fluid.set_weak_boundary("Left","Inflow")
fluid.set_weak_boundary("Right","Outlow")
fluid.set_weak_boundary("Right","Outflow")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export_vtk(outputdir,0,0)
......@@ -95,8 +94,7 @@ while t < tEnd :
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
if not use_lmgc : p.write_vtk(outputdir, ioutput, t)
#if not use_lmgc : p.write_boundary_vtk(outputdir, ioutput, t)
if not use_lmgc90 : 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))
......
......@@ -25,7 +25,7 @@
# Falling particles in a sandglass
from migflow import fluid
from migflow import scontact
from migflow import lmgc2Interface
from migflow import lmgc90Interface
import numpy as np
import os
......@@ -71,14 +71,14 @@ lxpart = 0.015 # semi-width of the grains
dt = 1e-3 # time step
#Enable the use of lmgc90
use_lmgc = True
use_lmgc90 = True
#Object particles creation
genInitialPosition(outputdir, r, 0., y0part, lxpart, lypart, rhop)
if use_lmgc :
if use_lmgc90 :
friction=0.1 # friction coefficient
lmgc2Interface.scontactToLmgc90(outputdir, 0, friction)
p = lmgc2Interface.LmgcInterface()
lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
p = lmgc90Interface.ParticleProblem(2)
else :
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
......@@ -123,4 +123,4 @@ while t < tEnd :
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
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......@@ -22,7 +22,7 @@
#!/usr/bin/env python
from migflow import fluid
from migflow import scontact
from migflow import lmgc2Interface
from migflow import lmgc90Interface
import numpy as np
import os
......@@ -58,12 +58,12 @@ for i in range(len(p1.position()[:,0])):
p.add_particle((p1.position()[i,0], p1.position()[i,1]), p1.r()[i], p1.r()[i]**2 * np.pi * rhop)
p.write_vtk(outputdir, 0, 0)