Skip to content
GitLab
Commit 3fd73132 authored by Michel Henry's avatar Michel Henry
Browse files

Symmetry Boundary & Die Swell

parent 2cc2e4b4
Pipeline #8729 passed with stages
in 3 minutes and 58 seconds
Hide whitespace changes
Inline Side-by-side
fluid/fluid_problem.c
View file @ 3fd73132
......@@ -217,15 +217,17 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
for (int jD = 0; jD < D; ++jD)
c_du_o_c[iD][jD] = ds[(U+iD)*D+jD] -u[iD]/c*dc[jD];
}
for (int id = 0; id < D; ++id) {
f0[U+id] += vid<0?0:sigma*(u[id]-data[vid+id]+s_c*n[id]);
f00[(U+id)*n_fields+U+id] += (vid<0?0:sigma);
for (int jd = 0; jd <D; ++jd) {
f0[U+id] -= mu*(c_du_o_c[id][jd]+c_du_o_c[jd][id])*n[jd];
f00[(U+id)*n_fields+U+id] += (vid<0?0:n[id]*n[jd]*sigma) + mu/c*dc[jd]*n[jd];
f00[(U+id)*n_fields+U+jd] += mu/c*dc[id]*n[jd];
f01[(U+id)*n_fields*D+(U+id)*D+jd] -= mu*n[jd];
f01[(U+id)*n_fields*D+(U+jd)*D+id] -= mu*n[jd];
if(wbnd->type != BND_SYMMETRY && wbnd->compute_viscous_term == 1){
for (int id = 0; id < D; ++id) {
f0[U+id] += vid<0?0:sigma*(u[id]-data[vid+id]+s_c*n[id]);
f00[(U+id)*n_fields+U+id] += (vid<0?0:sigma);
for (int jd = 0; jd <D; ++jd) {
f0[U+id] -= mu*(c_du_o_c[id][jd]+c_du_o_c[jd][id])*n[jd];
f00[(U+id)*n_fields+U+id] += (vid<0?0:n[id]*n[jd]*sigma) + mu/c*dc[jd]*n[jd];
f00[(U+id)*n_fields+U+jd] += mu/c*dc[id]*n[jd];
f01[(U+id)*n_fields*D+(U+id)*D+jd] -= mu*n[jd];
f01[(U+id)*n_fields*D+(U+jd)*D+id] -= mu*n[jd];
}
}
}
if (problem->advection) {
......@@ -1916,7 +1918,7 @@ void fluid_problem_set_strong_boundary(FluidProblem *problem, const char *tag, i
problem->strong_boundaries[i].row = row;
}
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, BoundaryType type, BoundaryCallback *cb, int vid, int pid, int cid, int aid)
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, BoundaryType type, BoundaryCallback *cb, int vid, int pid, int cid, int aid, int compute_viscous_term)
{
for (int i = 0; i < p->n_weak_boundaries; ++i) {
if (strcmp(p->weak_boundaries[i].tag, tag) == 0){
......@@ -1932,6 +1934,7 @@ void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, BoundaryT
bnd->pid = pid;
bnd->cid = cid;
bnd->aid = aid;
bnd->compute_viscous_term = compute_viscous_term;
bnd->type = type;
bnd->field_cb = cb;
}
......
fluid/fluid_problem.h
View file @ 3fd73132
......@@ -44,12 +44,12 @@ typedef struct {
} StrongBoundary;
typedef struct FluidProblem FluidProblem;
typedef enum {BND_WALL=0,BND_OPEN=1} BoundaryType;
typedef enum {BND_WALL=0,BND_OPEN=1,BND_SYMMETRY=2} BoundaryType;
typedef struct {
char *tag;
BoundaryType type;
BoundaryCallback *field_cb;
int vid,pid,cid,aid;
int vid,pid,cid,aid,compute_viscous_term;
} WeakBoundary;
......@@ -119,7 +119,7 @@ void fluid_problem_after_import(FluidProblem *problem);
int fluid_problem_n_fields(const FluidProblem *p);
double *fluid_problem_solution(const FluidProblem *p);
double *fluid_problem_mesh_velocity(const FluidProblem *p);
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, BoundaryType type, BoundaryCallback *cb, int vid, int pid, int cid, int aid);
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, BoundaryType type, BoundaryCallback *cb, int vid, int pid, int cid, int aid, int compute_viscous_term);
void fluid_problem_set_strong_boundary(FluidProblem *p, const char *tag, int field, int row, BoundaryCallback *apply);
double *fluid_problem_old_porosity(const FluidProblem *p);
double *fluid_problem_porosity(const FluidProblem *p);
......
python/fluid.py
View file @ 3fd73132
......@@ -186,7 +186,7 @@ class FluidProblem :
self.sys = None
gmsh.model.remove()
def set_wall_boundary(self, tag, pressure=None, velocity=None) :
def set_wall_boundary(self, tag, pressure=None, velocity=None, compute_viscous_term=1) :
"""Sets the weak boundaries (=normal fluxes) for the fluid problem.
Keyword arguments:
......@@ -213,9 +213,31 @@ class FluidProblem :
pid = self._dim
bndcb = BNDCB(_Bnd(cb_or_value,self._dim).apply)
self.weak_cb_ref.append(bndcb)
self._lib.fluid_problem_set_weak_boundary(self._fp,tag.encode(), c_int(0), bndcb, c_int(vid), c_int(pid), c_int(-1), c_int(-1))
self._lib.fluid_problem_set_weak_boundary(self._fp,tag.encode(), c_int(0), bndcb, c_int(vid), c_int(pid), c_int(-1), c_int(-1), c_int(int(compute_viscous_term)))
def set_open_boundary(self, tag, velocity=None, pressure=None, porosity=1, concentration=1):
def set_symmetry_boundary(self, tag, pressure=None):
"""Sets the symmetry boundaries (=normal fluxes) for the fluid problem.
Keyword arguments:
tag -- Physical tag (or list of tags), set in the mesh.geo file, of the symmetry boundaries
pressure -- The pressure value if imposed (callback or number)
velocity -- The velocity value if imposed (callback or number)
"""
if not _is_string(tag):
for t in tag:
self.set_symmetry_boundary(t,pressure)
return
pid = -1
vid = -1
cb_or_value = None
if pressure is not None:
cb_or_value = [pressure]
pid = 0
bndcb = BNDCB(_Bnd(cb_or_value, self._dim).apply)
self.weak_cb_ref.append(bndcb)
self._lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(), c_int(0), bndcb, c_int(vid), c_int(pid), c_int(-1), c_int(-1), c_int(-1))
def set_open_boundary(self, tag, velocity=None, pressure=None, porosity=1, concentration=1, compute_viscous_term=1):
"""Sets the weak boundaries (=normal fluxes) for the fluid problem.
Keyword arguments:
......@@ -256,7 +278,7 @@ class FluidProblem :
bndcb = BNDCB(_Bnd(cb_or_value, self._dim).apply)
self.weak_cb_ref.append(bndcb)
self._lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(),c_int(1),bndcb,c_int(vid),c_int(pid),c_int(cid),c_int(aid))
self._lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(),c_int(1),bndcb,c_int(vid),c_int(pid),c_int(cid),c_int(aid),c_int(int(compute_viscous_term)))
def set_strong_boundary(self, tag, field,callback_or_value, row=None) :
"""Sets the strong boundaries (=constrains) for the fluid problem.
......
testcases/dieSwell/dieSwell.py
View file @ 3fd73132
......@@ -10,6 +10,7 @@ import subprocess
import random
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)
n_e = len(lapl_h)-1
......@@ -26,11 +27,11 @@ 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"tension de surface : {max(abs(tension))}")
print(f"\tSurface tension : {max(abs(tension))}")
return tension
# Define output directory
outputdir = "output_MichelALE"
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","1"])
......@@ -40,19 +41,18 @@ g = 0 # gravity
rho = 1 # fluid density
nu = 1 # kinematic viscosity
mu = 1 # dynamic viscosity
gamma = 0
gamma = 0 # superficial tension coefficient
# Numerical parameters
outf = 2 # number of iterations between output files
dt = 0.025 # time step
tEnd = 50 # final time
outf = 5 # number of iterations between output files
dt = 0.05 # time step
tEnd = 50 # final time
# Geometrical parameters
H = 1 # domain height
L = 13.5 # domain width
Ox = 0
Oy = 0
L = 13.5 # domain width
Ox, Oy = 0, 0 # leftmost and bottommost point
# Initial time and iteration
......@@ -66,38 +66,38 @@ 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)
# fluid.set_strong_boundary("TopRight",2,0)
fluid.set_wall_boundary("BottomLeft")
fluid.set_wall_boundary("BottomRight")
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)
fluid.set_strong_boundary("Right", 1, 0)
newx = fluid.coordinates().copy()
#
# 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)
fluid.export_vtk(outputdir,0,0)
tic = time.time()
#
# COMPUTATION LOOP
while t < tEnd :
time_integration.iterate(fluid,None,dt)
if ii >= 100 :
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()[:] = (newx[:,:2] - fluid.coordinates()[:,:2])*fluid.porosity()[:]/dt
fluid.mesh_velocity()[:,:2] = (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}")
gonflement = newx[fs_nodes[-1],1]/H
print(f"\tSwelling = {gonflement}")
t += dt
# Output files writting
if ii %outf == 0 :
......
testcases/dieSwell/freeSurface.py
View file @ 3fd73132
......@@ -42,6 +42,10 @@ def get_nodes_struct(fluid,top_nodes, bottom_nodes, Oy = 0, Ox = 0):
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:
......@@ -61,7 +65,9 @@ def get_nodes_struct(fluid,top_nodes, bottom_nodes, Oy = 0, Ox = 0):
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
......@@ -88,30 +94,7 @@ def get_alpha(fluid,nodes_struct):
alpha[node] = (h_old[node] - O_old)/(H_old - O_old)
return alpha
def get_lapl_fs(fluid,nodes):
n = len(nodes)-1
enodes = np.column_stack([np.arange(0,n), np.arange(1,n+1)])
x = fluid.coordinates()[nodes,0]
h = 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
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------- PART 2 : Free surface Computation ------------------------------
......@@ -128,36 +111,10 @@ def get_h(fluid,nodes,dt,newx,vbox = 0):
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]
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
# w = 0.5
# h = newx[nodes[1:],1] + rhs*dt
# newx[nodes[1:],1] += w*(h - newx[nodes[1:],1])
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]
......@@ -251,6 +208,29 @@ def get_h_fe(fluid,nodes,dt,newx):
dhdt = spsolve(massglobal,rhs)
newx[nodes,1] += dhdt*dt
def get_lapl_fs(fluid,nodes):
n = len(nodes)-1
enodes = np.column_stack([np.arange(0,n), np.arange(1,n+1)])
x = fluid.coordinates()[nodes,0]
h = 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
# ----------------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------
......@@ -283,14 +263,14 @@ def update_inside_points(fluid, nodes_struct, newx, alpha):
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):
def arbitrary_inside_update(fluid,newx, nodes_struct, 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)
for node_struct in nodes_struct:
node = node_struct.get("node_index")
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) :
......@@ -334,6 +314,10 @@ def init_h(fluid,fs_nodes,a,b,c,newx):
for node in fs_nodes:
newx[node,1] += gaussian(x[node,0],a,b,c)
def sine(fluid,fs_nodes,A,w, newx):
x = fluid.coordinates()[:,:]
newx[fs_nodes,1] -= A*np.sin(w*x[fs_nodes,0] - w/2)
# newx[fs_nodes,1] -= A*np.sin(w*x[fs_nodes,0])
......
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