Commit 124e9065 authored by Michel Henry's avatar Michel Henry
Browse files

Add DieSwell

parent 01d619af
Pipeline #9838 failed with stage
in 1 minute and 52 seconds
from migflow import fluid as mbfluid
from migflow import time_integration
import freeSurface as fs
import numpy as np
import os
import time
import shutil
import subprocess
import random
import math
# def apply_tension_weak(x):
# lapl_h = fs.get_lapl_fs(fluid, 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"tension de surface : {max(abs(tension))}")
# return tension
# Define output directory
outputdir = "output_MichelALE"
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
# Numerical parameters
outf = 2 # number of iterations between output files
dt = 0.025 # time step
tEnd = 50 # final time
# Geometrical parameters
H = 1 # domain height
L = 13.5 # domain width
Ox = 0
Oy = 0
# 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 = 0)
# fluid.set_strong_boundary("TopRight",2,0)
fluid.set_wall_boundary("BottomLeft")
fluid.set_wall_boundary("BottomRight")
fluid.set_open_boundary("Right", pressure = 0)
fluid.set_strong_boundary("Right", 1,0)
newx = fluid.coordinates().copy()
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)
fluid.export_vtk(outputdir,0,0)
tic = time.time()
#
# COMPUTATION LOOP
while t < tEnd :
time_integration.iterate(fluid,None,dt)
# if ii >= 100 :
# 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()[:] = (newx[:,:2] - fluid.coordinates()[:,:2])*fluid.porosity()[:]/dt
# # fluid.coordinates()[:] = newx[:]
# fluid.set_coordinates(newx)
# gonflement = (newx[fs_nodes[-1],1])/H
# print(f"\tGonflement = {gonflement}")
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))
from migflow import fluid
import numpy as np
import os
import time
import shutil
import random
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------- PART 1 : Data structure-----------------------------------------
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------
# Return all the nodes which belong to the boundary identified by tag
def get_nodes(fluid,tag):
edges = get_boundary(fluid,tag)
nodes = edges.reshape(np.size(edges))
nodes = np.unique(nodes)
nodes = sorted(nodes,key=lambda node: fluid.coordinates()[node,0])
return nodes
# Return the edges which belong to the boundary identified by tag
def get_boundary(fluid,tag):
bnds = fluid._mesh_boundaries()
edges = bnds.get(tag.encode())
return(edges)
def get_nodes_struct(fluid,top_nodes, bottom_nodes, Oy = 0, Ox = 0):
eps = 1e-4
x = fluid.coordinates()[:,0]
y = fluid.coordinates()[:,1]
size = np.size(x)
nodes = np.arange(size)
nodes = np.setdiff1d(nodes, top_nodes)
nodes = np.setdiff1d(nodes, bottom_nodes)
min_top_bottom = top_nodes[0] if x[top_nodes[0]] < x[bottom_nodes[0]] else bottom_nodes[0]
nodes_to_delete = [node for node in nodes if x[node] < x[min_top_bottom]]
nodes = np.setdiff1d(nodes,nodes_to_delete)
nodes = sorted(nodes, key = lambda x: fluid.coordinates()[x,0])
# print(top_nodes)
# print(bottom_nodes)
L = []
for node in nodes:
if y[node] > Oy-eps and x[node] > Ox-eps:
# We need to find the 2 closest top_nodes and the 2 closest bottom_nodes
node_struct = dict(node_index = node)
n = np.size(top_nodes)
for i in range(0,n):
current_top_node = top_nodes[i]
next_top_node = top_nodes[i+1]
if x[node] >= x[current_top_node]-eps and x[node] < x[next_top_node]+eps:
node_struct["upper_nodes"] = [top_nodes[i], top_nodes[i+1]]
break
m = np.size(bottom_nodes)
for i in range(0,m):
current_bottom_node = bottom_nodes[i]
next_bottom_node = bottom_nodes[i+1]
if x[node] >= x[current_bottom_node] - eps and x[node] < x[next_bottom_node] + eps:
node_struct["lower_nodes"] = [bottom_nodes[i], bottom_nodes[i+1]]
break
# print(node_struct)
L.append(node_struct)
# print(L)
return L
def get_alpha(fluid,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
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------- PART 2 : Free surface Computation ------------------------------
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------
#
# Compute the free surface in the case of an extrusion problem where the velocity has a positive x direction.
# To do so, we solve :
# Dh/Dt = dh/dt + ux dh/dx - uy = 0,
# Where h = y - Y, y is the initial ordinate and Y is the new ordinate.
def get_h(fluid,nodes,dt,newx,vbox = 0):
x = fluid.coordinates()[nodes,0]
y = fluid.coordinates()[nodes,1]
u = fluid.solution()[nodes,:2]/fluid.porosity()[nodes]
dx = x[1:] - x[:-1]
dhdx = (y[1:]- y[:-1])/dx
# ux = (u[1:,0] + u[:-1,0])/2
ux = u[:-1,0]
# uy = (u[1:,1] + u[:-1,1])/2
uy = u[:-1,1]
rhs = uy - ux*dhdx
# print(f"dhdx = {dhdx[0]}\t v/u = {uy[0]/ux[0]}")
# print(f"angle = {np.rad2deg(np.arctan(uy[0]/ux[0]))}")
newx[nodes[1:],1] += rhs*dt
def get_h_Kruskal(fluid,nodes,dt,newx):
porosity = fluid.porosity()[nodes]
x = fluid.coordinates()[nodes,0]
y = fluid.coordinates()[nodes,1]
u = fluid.solution()[nodes,:2]/fluid.porosity()[nodes]
dx = x[2:] - x[:-2]
dhdx = (y[2:] - y[:-2])/dx
ux = (u[1:-1,0] + u[2:,0] + u[:-2,0])/3
uy = (u[1:-1,1] + u[2:,1] + u[:-2,1])/3
C = 1/(12*dx);
diffusion = C[:]*(abs(u[2:,0]+u[1:-1,0])/2 * (y[2:] - y[1:-1]) - abs(u[1:-1,0]+u[:-2,0])/2 * (y[1:-1] - y[:-2]))
rhs = np.zeros_like(x)
rhs[1:-1] = uy - ux*dhdx- diffusion
rhs[0] = 0
rhs[-1] = u[-1,1]
newx[nodes,1] += rhs*dt
def get_h_conservative(fluid,nodes,dt,newx):
x = fluid.coordinates()[:,0]
y = fluid.coordinates()[:,1]
u = fluid.solution()[:,:2]/fluid.porosity()[:]
dx = x[nodes[1:]]-x[nodes[:-1]]
dhdx = (y[nodes[1:]]-y[nodes[:-1]])/dx
uel = (u[nodes[1:],:]+u[nodes[:-1],:])/2
rhsel = (uel[:,1]-uel[:,0]*dhdx)*dx
lnode = np.zeros([len(nodes)])
rhsnode = np.zeros([len(nodes)])
lnode[:-1] += dx/2
lnode[1:] += dx/2
rhsnode[:-1] += rhsel/2
rhsnode[1:] += rhsel/2
rhsnode /= lnode
newx[nodes,1] += rhsnode*dt
def get_h_decentered(fluid,nodes,dt,newx):
porosity = fluid.porosity()[nodes]
x = fluid.coordinates()[nodes,0]
y = fluid.coordinates()[nodes,1]
u = fluid.solution()[nodes,:2]/fluid.porosity()[nodes]
u_up = np.maximum(u[1:-1,0],0)
u_down = -np.maximum(-u[1:-1,0],0)
dx_up = x[1:-1] - x[:-2]
dx_down = x[2:] - x[1:-1]
dhdx_up = (y[1:-1] - y[:-2])/dx_up
dhdx_down = (y[2:] - y[1:-1])/dx_down
u_dhdx = u_up*dhdx_up + u_down*dhdx_down
rhs = np.zeros_like(x)
rhs[1:-1] = u[1:-1,1] - u_dhdx
rhs[0] = u[0,1] #dhdx = 0 sur les bords
rhs[-1] = u[-1,1] #dhdx = 0 sur les bords
newx[nodes,1] += rhs*dt
def get_h_centered(fluid,nodes,dt,newx):
porosity = fluid.porosity()[nodes]
x = fluid.coordinates()[nodes,0]
y = fluid.coordinates()[nodes,1]
u = fluid.solution()[nodes,:2]/fluid.porosity()[nodes]
dx = x[2:] - x[:-2]
dhdx = (y[2:] - y[:-2])/dx
ux = u[1:-1,0]
uy = u[1:-1,1]
rhs = np.zeros_like(x)
rhs[1:-1] = uy - ux*dhdx
rhs[0] = u[0,1]
rhs[-1] = u[-1,1]
newx[nodes,1] += rhs*dt
def get_h_fe(fluid,nodes,dt,newx):
n = len(nodes)-1
enodes = np.column_stack([np.arange(0,n),np.arange(1,n+1)])
x = fluid.coordinates()[nodes,0]
y = fluid.coordinates()[nodes,1]
u = fluid.solution()[nodes,:2]
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
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------- PART 3 : Mesh regularisation -----------------------------------
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------
def update_inside_points(fluid, nodes_struct, newx, alpha):
x = newx[:,0]
h = newx[:,1]
h_old = fluid.coordinates()[:,1]
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]]
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] = alpha[node]*(H_new - O_new) + O_new
def arbitrary_inside_update(fluid,newx, inside_nodes, A,w,t,Ox = 0, L = 1):
eps = 1e-4
for col in inside_nodes:
for node in col:
x = fluid.coordinates()[node,0]
if x > Ox + eps and x < Ox + L - eps :
newx[node,0] += A*np.cos(w*t)
newx[node,1] += A*np.sin(w*t)
class UnIntegrator() :
def __init__(self,fluid,tag) :
self.bnd = fluid._mesh_boundaries()[tag.encode()]
self.fluid = fluid
def compute(self) :
x = self.fluid.coordinates()[self.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()[self.bnd,:2],axis=1)
un = u[:,0]*n[:,0]+u[:,1]*n[:,1]
return np.sum(un*l)
def corr_un(fluid,un_integrator,nodes,vbox=0):
integral_un = un_integrator.compute()
print(f"Integral un before correction: {integral_un-vbox}")
fluid.solution()[nodes,1] += (vbox-integral_un)
def integral_h(fluid,fs_nodes,Oy,H):
h = fluid.coordinates()[fs_nodes,1] - Oy - H
x = fluid.coordinates()[fs_nodes,0]
integral = np.trapz(h,x)
print(f"Integral de h sur le contour : {integral}")
return(integral)
def corr_h(fs_nodes, integral_h,newx,L):
corr = integral_h/L
newx[fs_nodes,1] -= corr
# Return the height of a gaussian function
def gaussian(x,a=1,b=0,c=1):
y = a*np.exp(-(x-b)*(x-b)/(2*c*c))
return y
# Initialise the height as a gaussian curve
# A set_coordinates is still required after this initialisation
def init_h(fluid,fs_nodes,a,b,c,newx):
x = fluid.coordinates()[:,:2]
for node in fs_nodes:
newx[node,1] += gaussian(x[node,0],a,b,c)
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------- Old Code Version -----------------------------------------------
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------
# Return the indices of all the inside points in a nodes_struct which will be used in the update function
def get_inside_points(fluid, top_nodes, Oy, newx):
eps = 1e-4
x = newx[:,0]
y = newx[:,1]
n = np.size(x)
nodes = np.arange(n)
nodes = np.setdiff1d(nodes,top_nodes)
nodes = sorted(nodes, key= lambda x: fluid.coordinates()[x,0])
nodes_struct = []
i = 0
L = []
for node in nodes:
if y[node] > (Oy + eps): # to exclude the Bottom points
current_top_node = top_nodes[i]
next_top_node = top_nodes[i+1]
if x[node] >= x[current_top_node]-eps and x[node] < x[next_top_node]+eps:
L.append(node)
if x[node] >= x[next_top_node]+eps:
i += 1
nodes_struct.append(L.copy())
L.clear()
L.append(node)
nodes_struct.append(L.copy())
return nodes_struct
def update_inside_points_old(fluid, fs_nodes, inside_nodes, Origin_old, dt, newx, vbox = 0,eps = 1e-8):
i = 0
x = newx[:,0]
h = newx[:,1]
h_old = fluid.coordinates()[:,1]
Origin_new = Origin_old + vbox*dt
for col in inside_nodes:
nodeLeft = fs_nodes[i]
nodeRight = fs_nodes[i+1]
dx = x[nodeRight] - x[nodeLeft]
for node in col:
wRight = (x[node] - x[nodeLeft])/dx
wLeft = 1 - wRight
H_old = wLeft*h_old[nodeLeft] + wRight*h_old[nodeRight]
h_node_old = h_old[node]
alpha = (h_node_old - Origin_old) / (H_old - Origin_old)
# if alpha <= eps : alpha = 0
H_new = wLeft*h[nodeLeft] + wRight*h[nodeRight]
rhs = alpha*(H_new - Origin_new) + Origin_old
newx[node,1] = rhs
i +=1
def get_h_MacCormack(fluid, nodes, dt, newx):
porosity = fluid.porosity()[:]
x = fluid.coordinates()[:,0]
h = fluid.coordinates()[:,1]
u = fluid.solution()[:,:2]
u[:,0] -= fluid.mesh_velocity()[:,0]
rhs = np.zeros_like(h)
n = np.size(nodes)
for i in range(1,n-1):
nodeLeft = nodes[i-1]
node = nodes[i]
nodeRight = nodes[i+1]
phi = porosity[node]
ux = u[node,0]
uy = u[node,1]
uxLeft = u[nodeLeft,0]
uyLeft = u[nodeLeft,1]
dx = x[nodeRight] - x[node]
dxLeft = x[nodeLeft] - x[node]
pred = h[node] - ux*dt/dx*(h[nodeRight] - h[node]) + uy*dt
pred_left = h[nodeLeft] - uxLeft*dt/dx*(h[node] - h[nodeLeft]) + uyLeft*dt
h_node_half_step = (h[node]+pred)/2
corr = h_node_half_step - ux*dt/(2*dx)*(pred - pred_left) + uy*dt/2
rhs[node] = corr
# Leftmost point
i = 0
node = nodes[i]
nodeRight = nodes[i+1]
phi = (porosity[node] + porosity[nodeRight])/2
ux = u[node,0]
uy = u[node,1]
dx = x[nodeRight] - x[node]
# dhdx = (y[nodeRight] - y[node])/dx
dhdx = 0
rhs[node] = h[node] - (ux*dhdx - uy)*dt/phi
#Rightmost point
i = n-1
node = nodes[i]
nodeLeft = nodes[i-1]
phi = (porosity[node] + porosity[nodeLeft])/2
ux = u[node,0]
uy = u[node,1]
dx = x[node] - x[nodeLeft]
# dhdx = (y[node] - y[nodeLeft])/dx
dhdx = 0
rhs[node] = h[node] - (ux*dhdx - uy)*dt/phi
for i in range(0,n):
newx[nodes[i],1] = rhs[nodes[i]]
lc = 0.2;
L1 = 3.5;
L2 = 10;
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};