Commit acb5bd0d authored by Michel Henry's avatar Michel Henry
Browse files

Ajout dieSwell validation + free surface into time_integration

parent 3fd73132
......@@ -19,6 +19,8 @@
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve
def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_sub_iter=None,max_split=1,nsubt=[]) :
"""Contact solver for the grains
......@@ -159,3 +161,202 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
fluid.move_particles(particles.position(), particles.velocity(),-fluid.compute_node_force(dt))
else:
predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,after_sub_iter=after_sub_iter,max_split=max_split)
class FreeSurface():
def __init__(self, fluid, tagFS, tagOFS, isCentered=True, correctVelocity=False, origin = (0,0)):
if fluid._dim != 2:
print("Free Surface not available for this dimension !")
self.fluid = fluid
self.tagFS = tagFS
self.tafOFS = tagOFS
self.fs_nodes = self.get_nodes(tagFS)
self.ofs_nodes = self.get_nodes(tagOFS)
self.scheme = self.compute_h_centered if isCentered else self.compute_h_decentered
self.correctVelocity = correctVelocity
self.nodes_struct = self.get_nodes_struct(self.fs_nodes, self.ofs_nodes, origin[0], origin[1])
self.alpha = self.get_alpha()
def iterate(self,dt,vbox = [0,0]):
dim = self.fluid._dim
vbox = np.array(vbox, dtype=np.float)
newx = np.copy(self.fluid.coordinates())
if self.correctVelocity :
self.corr_un(vbox)
self.scheme(dt,newx)
newx[self.ofs_nodes,:dim] += vbox[None,:dim]*dt
self.update_inside_points(newx)
self.fluid.mesh_velocity()[:,:dim] = (newx[:,:dim] - self.fluid.coordinates()[:,:dim])*self.fluid.porosity()[:]
self.fluid.coordinates()[:] = newx[:]
def get_nodes_struct(self, fs_nodes, ofs_nodes, Ox, Oy):
eps = 1e-4
x = self.fluid.coordinates()[:,0]
y = self.fluid.coordinates()[:,1]
nodes = np.arange(np.size(x))
nodes = np.setdiff1d(nodes, fs_nodes)
nodes = np.setdiff1d(nodes, ofs_nodes)
min_node = fs_nodes[0] if x[fs_nodes[0]] < x[ofs_nodes[0]] else ofs_nodes[0]
nodes_to_delete = [node for node in nodes if x[node] < x[min_node]]
nodes = np.setdiff1d(nodes,nodes_to_delete)
nodes = sorted(nodes, key = lambda x: self.fluid.coordinates()[x,0])
L = []
for node in nodes:
if y[node] > Oy-eps and x[node] > Ox-eps:
node_struct = dict(node_index = node)
n = np.size(fs_nodes)
for i in range(0,n):
current_top_node = fs_nodes[i]
next_top_node = fs_nodes[i+1]
if x[node] >= x[current_top_node]-eps and x[node] < x[next_top_node]+eps:
node_struct["upper_nodes"] = [fs_nodes[i], fs_nodes[i+1]]
break
m = np.size(ofs_nodes)
for i in range(0,m):
current_bottom_node = ofs_nodes[i]
next_bottom_node = ofs_nodes[i+1]
if x[node] >= x[current_bottom_node] - eps and x[node] < x[next_bottom_node] + eps:
node_struct["lower_nodes"] = [ofs_nodes[i], ofs_nodes[i+1]]
break
L.append(node_struct)
return L
def get_nodes(self, tag):
nodes = (self.fluid._mesh_boundaries())[tag.encode()]
nodes = np.unique(nodes.reshape(np.size(nodes)))
nodes = sorted(nodes, key=lambda node: self.fluid.coordinates()[node,0])
return nodes
def get_alpha(self):
fluid = self.fluid
nodes_struct = self.nodes_struct
x = fluid.coordinates()[:,0]
h_old = fluid.coordinates()[:,1]
alpha = np.zeros_like(x)
for node_struct in nodes_struct:
node = node_struct.get("node_index")
upper_nodes = node_struct.get("upper_nodes")
lower_nodes = node_struct.get("lower_nodes")
upper_dx = x[upper_nodes[1]] - x[upper_nodes[0]]
lower_dx = x[lower_nodes[1]] - x[lower_nodes[0]]
upper_wRight = (x[node] - x[upper_nodes[0]])/upper_dx
upper_wLeft = 1 - upper_wRight
lower_wRight = (x[node] - x[lower_nodes[0]])/lower_dx
lower_wLeft = 1 - lower_wRight
H_old = upper_wLeft*h_old[upper_nodes[0]] + upper_wRight*h_old[upper_nodes[1]]
O_old = lower_wLeft*h_old[lower_nodes[0]] + lower_wRight*h_old[lower_nodes[1]]
alpha[node] = (h_old[node] - O_old)/(H_old - O_old)
return alpha
def update_inside_points(self, newx):
x, h = newx[:,0], newx[:,1]
h_old = self.fluid.coordinates()[:,1]
for node_struct in self.nodes_struct:
node = node_struct.get("node_index")
upper_nodes = node_struct.get("upper_nodes")
lower_nodes = node_struct.get("lower_nodes")
upper_dx = x[upper_nodes[1]] - x[upper_nodes[0]]
lower_dx = x[lower_nodes[1]] - x[lower_nodes[0]]
upper_wRight = (x[node] - x[upper_nodes[0]])/upper_dx
upper_wLeft = 1 - upper_wRight
lower_wRight = (x[node] - x[lower_nodes[0]])/lower_dx
lower_wLeft = 1 - lower_wRight
# H_old = upper_wLeft*h_old[upper_nodes[0]] + upper_wRight*h_old[upper_nodes[1]]
# O_old = lower_wLeft*h_old[lower_nodes[0]] + lower_wRight*h_old[lower_nodes[1]]
H_new = upper_wLeft*h[upper_nodes[0]] + upper_wRight*h[upper_nodes[1]]
O_new = lower_wLeft*h[lower_nodes[0]] + lower_wRight*h[lower_nodes[1]]
newx[node,1] = (self.alpha)[node]*(H_new - O_new) + O_new
def compute_h_centered(self,dt,newx):
nodes = self.fs_nodes
n = len(nodes)-1
enodes = np.column_stack([np.arange(0,n),np.arange(1,n+1)])
x = self.fluid.coordinates()[nodes,0]
y = self.fluid.coordinates()[nodes,1]
u = self.fluid.solution()[nodes,:2]/self.fluid.porosity()[:]
l = (x[1:]-x[:-1])
J = l/2
#nu = 1e-5
dphidx = np.column_stack([-1/l,1/l])
xi = np.array([-1/3,1/3])
phi = np.column_stack([(1-xi)/2,(1+xi)/2])
qw = np.array([1,1])
yqp = y[enodes]@phi
uqp = u[enodes,0]@phi
vqp = u[enodes,1]@phi
dudx = np.sum(u[enodes,0]*dphidx,axis=1)
dydx = np.sum(y[enodes]*dphidx,axis=1)
# rhsqp : element, test function, qp
rhsqp = J[:,None,None]*qw[None,None,:]*(
vqp[:,None,:]*phi[None,:,:]
- uqp[:,None,:]*phi[None,:,:]*dydx[:,None,None]
#+ nu*dphidx[:,:,None]*dudx[:,None,None]]
)
rhslocal = np.sum(rhsqp,axis=2)
rhs = np.zeros((n+1))
np.add.at(rhs,enodes.flat,rhslocal.flat)
massparent = np.array([[2/3,1/3],[1/3,2/3]])
masslocal = J[:,None,None]*massparent[None,:,:]
rows = np.repeat(enodes[:,:,None],2,axis=2)
cols = np.repeat(enodes[:,None,:],2,axis=1)
massglobal = coo_matrix((masslocal.flat,(rows.flat,cols.flat))).tocsr()
dhdt = spsolve(massglobal,rhs)
newx[nodes,1] += dhdt*dt
def compute_h_decentered(self,dt,newx):
nodes = self.fs_nodes
x = self.fluid.coordinates()[nodes,0]
y = self.fluid.coordinates()[nodes,1]
u = self.fluid.solution()[nodes,:2]/self.fluid.porosity()[nodes]
dx = x[1:] - x[:-1]
dhdx = (y[1:]- y[:-1])/dx
ux = (u[1:,0] + u[:-1,0])/2
uy = (u[1:,1] + u[:-1,1])/2
rhs = uy - ux*dhdx
newx[nodes[1:],1] += rhs*dt
def compute_lapl_fs(self,nodes):
n = len(nodes)-1
enodes = np.column_stack([np.arange(0,n), np.arange(1,n+1)])
x = self.fluid.coordinates()[nodes,0]
h = self.fluid.coordinates()[nodes,1]
l = (x[1:] - x[:-1])
J = l/2
dphidx = np.column_stack([-1/l,1/l])
xi = np.array([-1/3,1/3])
phi = np.column_stack([(1-xi)/2, (1+xi)/2])
qw = np.array([1,1])
dhqp = np.sum(h[enodes]*dphidx,axis=1)
rhsqp = -J[:,None,None]*qw[None,None,:]*(dphidx[:,:,None]*dhqp[:,None,None])
rhslocal= np.sum(rhsqp,axis=2)
rhs = np.zeros((n+1))
np.add.at(rhs,enodes.flat,rhslocal.flat)
massparent = np.array([[2/3,1/3],[1/3,2/3]])
masslocal = J[:,None,None]*massparent[None,:,:]
rows = np.repeat(enodes[:,:,None],2,axis=2)
cols = np.repeat(enodes[:,None,:],2,axis=1)
massglobal= coo_matrix((masslocal.flat, (rows.flat, cols.flat))).tocsr()
ddhdx = spsolve(massglobal,rhs)
return ddhdx
def compute_integral_un(self):
bnd = self.fluid._mesh_boundaries()[(self.tagFS).encode()]
x = self.fluid.coordinates()[bnd]
dx = x[:,1,:]-x[:,0,:]
l = np.linalg.norm(dx,axis=1)
n = np.column_stack([dx[:,1],-dx[:,0]])/l[:,None]
u = np.mean(self.fluid.solution()[bnd,:2],axis=1)
un = u[:,0]*n[:,0]+u[:,1]*n[:,1]
return np.sum(un*l)
def corr_un(self, vbox):
integral_un = self.compute_integral_un()
print(f"Integral un before correction: {integral_un - vbox[-1]}")
self.fluid.solution()[self.fs_nodes,1] += (vbox[-1]-integral_un)
\ No newline at end of file
from migflow import fluid as mbfluid
from migflow import time_integration
import freeSurface as fs
import numpy as np
import os
import time
......@@ -12,7 +11,8 @@ import math
# Function used to apply surface tension at the free surface
def apply_tension_weak(x):
lapl_h = fs.get_lapl_fs(fluid, fs_nodes)
fs_nodes = fs.fs_nodes
lapl_h = fs.compute_lapl_fs(fs_nodes)
n_e = len(lapl_h)-1
kappa = np.zeros_like(x[:,0])
X = fluid.coordinates()[fs_nodes,0]
......@@ -27,31 +27,31 @@ def apply_tension_weak(x):
phi = np.column_stack([(1-xi)/2, (1+xi)/2])
kappa[ind_i] = (phi[0,:]*lapl_h[i] + phi[1,:]*lapl_h[i+1])*dl[i]
tension = -gamma*kappa[:]
print(f"\tSurface tension : {max(abs(tension))}")
# print(f"\tSurface tension : {max(abs(tension))}")
return tension
# Define output directory
outputdir = "output"
outputdir = "outputRE11_8"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","1"])
# Physical parameters
g = 0 # gravity
rho = 1 # fluid density
nu = 1 # kinematic viscosity
mu = 1 # dynamic viscosity
gamma = 0 # superficial tension coefficient
rho = 5.9 # fluid density
mu = 1 # dynamic viscosity
nu = mu/rho # kinematic viscosity
gamma = 1 # superficial tension coefficient
# Numerical parameters
outf = 5 # number of iterations between output files
dt = 0.05 # time step
tEnd = 50 # final time
dt = 0.1 # time step
tEnd = 75 # final time
# Geometrical parameters
H = 1 # domain height
L = 13.5 # domain width
L = 20.5 # domain width
Ox, Oy = 0, 0 # leftmost and bottommost point
......@@ -66,38 +66,30 @@ fluid = mbfluid.FluidProblem(2,[0,g],[mu],[rho])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_open_boundary("Left", velocity = [lambda x : 2*(1-x[:,1]**2),0])
# fluid.set_strong_boundary("Left",0, lambda x : 2*(1-x[:,1]**2))
# fluid.set_strong_boundary("Left",1,0)
fluid.set_strong_boundary("Left",0, lambda x : 2*(1-x[:,1]**2))
fluid.set_strong_boundary("Left",1,0)
fluid.set_wall_boundary("TopLeft", velocity = [0,0])
fluid.set_open_boundary("TopRight", pressure = apply_tension_weak, compute_viscous_term = 0)
fluid.set_symmetry_boundary("BottomLeft")
fluid.set_symmetry_boundary("BottomRight")
fluid.set_open_boundary("Right", pressure = 0)
fluid.set_strong_boundary("Right", 1, 0)
#
# FREE SURFACE STRUCTURES
#
fs_nodes = fs.get_nodes(fluid, tag = "TopRight")
bottom_right = fs.get_nodes(fluid,tag = "BottomRight")
nodes_struct = fs.get_nodes_struct(fluid, fs_nodes,bottom_right)
alpha = fs.get_alpha(fluid, nodes_struct)
fs = time_integration.FreeSurface(fluid,"TopRight","BottomRight",isCentered=False)
fluid.export_vtk(outputdir,0,0)
tic = time.time()
#
# COMPUTATION LOOP
#
swelling = 0
while t < tEnd :
time_integration.iterate(fluid,None,dt)
if ii >= 20:
newx = np.copy(fluid.coordinates())
fs.get_h(fluid,fs_nodes,dt,newx)
fs.update_inside_points(fluid, nodes_struct, newx, alpha)
fluid.mesh_velocity()[:,:2] = (newx[:,:2] - fluid.coordinates()[:,:2])*fluid.porosity()[:]/dt
fluid.coordinates()[:] = newx[:]
gonflement = newx[fs_nodes[-1],1]/H
print(f"\tSwelling = {gonflement}")
if t >= 50.0:
fs.iterate(dt)
swelling = fluid.coordinates()[fs.fs_nodes[-1],1]/H
t += dt
# Output files writting
if ii %outf == 0 :
......@@ -105,3 +97,4 @@ while t < tEnd :
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
print(f"\tSwelling = {swelling}")
\ No newline at end of file
lc = 0.2;
L1 = 3.5;
L2 = 10;
L2 = 17;
H = 1;
// Génération du rectangle
......@@ -37,4 +37,4 @@ Physical Line("Right") = {6};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
Mesh.MshFileVersion=2;
// Mesh.MshFileVersion=2;
# MigFlow - Copyright (C) <2010-2020>
# <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
import numpy as np
import os
import subprocess
import time
import shutil
import random
import unittest
class DieSwell(unittest.TestCase):
def runTest(self):
# Function used to apply surface tension at the free surface
def apply_tension_weak(x):
fs_nodes = fs.fs_nodes
lapl_h = fs.compute_lapl_fs(fs_nodes)
n_e = len(lapl_h)-1
kappa = np.zeros_like(x[:,0])
X = fluid.coordinates()[fs_nodes,0]
h = fluid.coordinates()[fs_nodes,1]
dX = X[1:] - X[:-1]
dhdx = (h[1:] - h[:-1])/dX[:]
dl = np.sqrt(dhdx[:]*dhdx[:] + 1)**(-3)
ind = np.argsort(x[:,0])
for i in range(0,n_e):
ind_i = ind[2*i:2*i+2]
xi = (2*x[ind_i,0] - (X[i+1] + X[i]))/(X[i+1]- X[i])
phi = np.column_stack([(1-xi)/2, (1+xi)/2])
kappa[ind_i] = (phi[0,:]*lapl_h[i] + phi[1,:]*lapl_h[i+1])*dl[i]
tension = -gamma*kappa[:]
# print(f"\tSurface tension : {max(abs(tension))}")
return tension
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","1"])
# Physical parameters
g = 0 # gravity
rho = 5.9 # fluid density
mu = 1 # dynamic viscosity
nu = mu/rho # kinematic viscosity
gamma = 1 # superficial tension coefficient
# Numerical parameters
outf = 5 # number of iterations between output files
dt = 0.05 # time step
tEnd = 75 # final time
# Geometrical parameters
H = 1 # domain height
L = 20.5 # domain width
Ox, Oy = 0, 0 # leftmost and bottommost point
# Initial time and iteration
t = 0
ii = 0
#
# FLUID PROBLEM
#
fluid = mbfluid.FluidProblem(2,[0,g],[mu],[rho])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_open_boundary("Left", velocity = [lambda x : 2*(1-x[:,1]**2),0])
fluid.set_strong_boundary("Left",0, lambda x : 2*(1-x[:,1]**2))
fluid.set_strong_boundary("Left",1,0)
fluid.set_wall_boundary("TopLeft", velocity = [0,0])
fluid.set_open_boundary("TopRight", pressure = apply_tension_weak, compute_viscous_term = 0)
fluid.set_symmetry_boundary("BottomLeft")
fluid.set_symmetry_boundary("BottomRight")
fluid.set_open_boundary("Right", pressure = 0)
#
# FREE SURFACE STRUCTURES
#
fs = time_integration.FreeSurface(fluid,"TopRight","BottomRight",isCentered=False)
fluid.export_vtk(outputdir,0,0)
tic = time.time()
#
# COMPUTATION LOOP
#
swelling = 0
while t < tEnd :
time_integration.iterate(fluid,None,dt)
if t >= 50.0:
fs.iterate(dt)
swelling = fluid.coordinates()[fs.fs_nodes[-1],1]/H
t += dt
# Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
swelling_experiment = 1.02 # Value from "Goren, S. L., & Wronski, S. (1966).
# The shape of low-speed capillary jets of Newtonian liquids.
# Journal of Fluid Mechanics, 25(1), 185-198."
error = np.abs(swelling - swelling_experiment)
self.assertLess(error, 2*10-2, "error is too large in dieSwell")
\ No newline at end of file
lc = 0.2;
L1 = 3.5;
L2 = 17;
H = 1;
// Génération du rectangle
Point(1) = {L1+L2,H,0,lc};
Point(2) = {L1,H,0,lc/3};
Point(3) = {0,H,0,lc};
Point(4) = {0,0,0,lc};
Point(5) = {L1,0,0,lc};
Point(6) = {L2+L1,0,0,lc};
// Génération des lignes
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,5};
Line(5) = {5,6};
Line(6) = {6,1};
// Génération des contours
Line Loop(1) = {1:6};
// Génération de la surface
Plane Surface(1) = {1};
// On pose les conditions physiques
Physical Line("TopRight") = {1};
Physical Line("TopLeft") = {2};
Physical Line("Left") = {3};
Physical Line("BottomLeft") = {4};
Physical Line("BottomRight") = {5};
Physical Line("Right") = {6};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
Mesh.MshFileVersion=2;
......@@ -8,7 +8,8 @@ all_tests = ["weight.weight",
"oneGrain.oneGrain",
"darcy.darcy",
"periodic2d.periodic2d",
"periodic3d.periodic3d",]
"periodic3d.periodic3d",
"dieSwell.dieSwell"]
suite = unittest.TestSuite()
for t in all_tests :
suite.addTest(unittest.defaultTestLoader.loadTestsFromName(t))
......
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