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

wip

parent 320a6429
......@@ -32,6 +32,7 @@ set(SRC
scipylsys.py
VTK.py
dof_numbering.py
free_surface.py
)
if(BUILD_MUMPS)
......
......@@ -373,8 +373,10 @@ class FluidProblem :
"""
da = self.concentration_dg_grad()
cell_data = []
v = self.velocity()
if self._dim == 2 :
da = np.insert(da,self._dim,0,axis=1)
v = np.insert(v,self._dim,0,axis=1)
data = [
("pressure",self.pressure().reshape(-1,1)),
("porosity",self.porosity()),
......@@ -392,11 +394,11 @@ class FluidProblem :
offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0])).astype(np.int32)
isp2p1 = np.max(self.order_fields) == 2
if isp2p1 :
field_data.append((b"velocity", self.velocity()))
field_data.append((b"velocity", v))
idx = self.get_numbering(self.elements(),[2])
cell_data.append(("velocity_map", idx.reshape(self.n_elements(),-1)))
else:
data.append(("velocity", self.velocity()))
data.append(("velocity", v))
VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data,cell_data)
def read_vtk(self, dirname, i):
......@@ -430,7 +432,7 @@ class FluidProblem :
sol = self.solution()
isp2p1 = np.max(self.order_fields) == 2
velocity = fdata['velocity'] if isp2p1 else data["velocity"]
sol[id_v] = velocity.reshape([-1])
sol[id_v] = velocity[:,:self._dim].reshape([-1])
sol[id_p] = data["pressure"].reshape([-1])
if self._n_fluids == 2 :
self.concentration_dg()[:] = cdata["concentration"]
......@@ -574,6 +576,7 @@ class FluidProblem :
id_v = np.r_[np.repeat(np.arange(self._dim)[None,:],self.n_nodes(),axis=0).reshape(-1)
+ self.n_fields()*np.repeat(np.arange(self.n_nodes()), self._dim),
np.arange(self.n_fields()*self.n_nodes(), ndof)]
return self.solution()[id_v].reshape(-1, self._dim)
def mesh_velocity(self) :
......
......@@ -19,8 +19,6 @@
# 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,max_nsub=None,after_sub_iter=None,n_divide=None) :
"""Contact solver for the grains
......@@ -165,202 +163,3 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
else:
predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,
after_sub_iter=after_sub_iter,max_nsub=max_nsub,check_residual_norm=check_residual_norm)
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.float64)
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()[:]/dt
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.velocity()[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.velocity()[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.velocity()[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.velocity()[self.fs_nodes,1] += (vbox[-1]-integral_un)
......@@ -4,34 +4,56 @@ import time
from migflow import fluid
from migflow import time_integration
from migflow import free_surface
outputdir = "output"
outputdir = "output_2750"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
outf = 5
outf = 2
h = 0.03
r = 0.001
xc = np.array([0.0,0])
# xc = np.array([0,-h/3+r])
vc = np.array([0.0,0.0])
g = np.array([0, -9.81])
rho = 1000
nu = 1e-6
dt = 1e-2
tEnd = 100
rho = 1200
rho_b = 2750#1200#2500#2750
m_b = rho_b*np.pi*r**2
print(f"m : {m_b} ---- r = {r}")
nu = 1.85e-5
dt = 1e-3
tEnd = 5
fluid = fluid.FluidProblem(2,g,nu*rho,rho, order=[2,2,1], solver="petsc4py", solver_options="-pc_type lu")
fluid.load_msh("mesh.msh")
c = np.array([0.5,0.5])
r = 0.1
v_max = 0.1
def v_bubble(x, axis=0):
r = x[:,:]-xc[None,:]
r = r[:,:]/np.linalg.norm(r,axis=1)[:,None]
return v_max*r[:,axis]
# r = x[:,axis]-xc[None,axis]
# r /= np.linalg.norm(r,axis=0)
# return 0.1*r
fluid.set_open_boundary("Lateral", velocity = [0,0])
fluid.set_open_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom", velocity=[0,0])
fluid.set_wall_boundary("Bubble", velocity=[lambda x : (x[:,0]-c[0])/r, lambda x : (x[:,1]-c[1])/r])
fluid.set_strong_boundary("Bubble", 0, lambda x : x[:,0]-c[0])
fluid.set_strong_boundary("Bubble", 1, lambda x : x[:,1]-c[1])
fluid.set_open_boundary("Bubble", velocity=[lambda x: v_bubble(x,0), lambda x : v_bubble(x,1)])
fluid.set_strong_boundary("Bubble", 0, lambda x : v_bubble(x,0))
fluid.set_strong_boundary("Bubble", 1, lambda x : v_bubble(x,1))
nodes_bubble= np.unique(fluid.mesh_boundaries()[b"Bubble"])
nodes_boundary = np.unique(np.concatenate(list(fluid.mesh_boundaries().values()), axis=0))
x = fluid.coordinates()
t = 0
ii = 0
......@@ -39,12 +61,24 @@ ii = 0
fluid.write_vtk(outputdir, t, ii)
tic = time.time()
while(t < tEnd):
fluid.implicit_euler(dt, 1e-2)
fboundary = fluid.boundary_force()[3,:]
print(fboundary)
fluid.implicit_euler(dt, 1e-1)
if t>0.01:
fbubble = fluid.boundary_force()[3,:] + m_b*g
print("\t fb = ", fluid.boundary_force()[3,1], "-----", m_b*g[1])
print("\t fbubble = ", fbubble)
vc = vc + fbubble/m_b*dt
dxc = vc*dt
xold = x.copy()
R = min(0.75*(h/2-xc[1]), 0.75*(h/2+xc[1]))
free_surface.regularize_point(x, R, xc, r, dxc, nodes_boundary)
xc = xc+dxc
x[nodes_bubble,:2] = x[nodes_bubble,:2] + dxc[None,:2]
fluid.mesh_velocity()[:,:2] = (x[:,:2]-xold[:,:2])/dt
t += dt
if ii%outf == 0:
ioutput = int(ii/outf)+1
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
print(f"xc : {xc}")
\ No newline at end of file
H = 1;
L = 1;
Cx = 0.5;
Cy = 0.5;
r = 0.1;
lc = 0.1;
lc_circle = lc/10;
H = 0.03;
L = 0.03;
r = 0.001;
Cx = 0.0;
// Cy = -H/3+r;
Cy = 0;
lc = 0.004;
lc_circle = 0.0005;
Point(1) = {0,0,0,lc};
Point(2) = {L,0,0,lc};
Point(3) = {L,H,0,lc};
Point(4) = {0,H,0,lc};
Point(1) = {-L/2,-H/2,0,lc};
Point(2) = {L/2,-H/2,0,lc};
Point(3) = {L/2,H/2,0,lc};
Point(4) = {-L/2,H/2,0,lc};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
......@@ -34,6 +35,7 @@ Physical Line("Bottom") = {1};
Physical Line("Lateral") = {2,4};
Physical Line("Top") = {3};
Physical Line("Bubble") = {5:8};
Physical Surface("Domain") = {1};
Mesh.MshFileVersion = 4.1;
......@@ -76,6 +76,7 @@ class Cavity(unittest.TestCase) :
fluid.set_wall_boundary("Bottom",velocity=[0,0])
fluid.set_wall_boundary("Top",velocity=[vUp,0])
fluid.set_wall_boundary("Right",velocity=[0,0])
ii = 0
t = 0
......@@ -84,7 +85,6 @@ class Cavity(unittest.TestCase) :
s = fluid.solution()
#s[:,3] = -rho*g*(1-x[:,1])
fluid.write_vtk(outputdir,0,0)
ii = 0
......@@ -99,7 +99,6 @@ class Cavity(unittest.TestCase) :
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
s = fluid.velocity()
x = fluid.coordinates(order=2)
......
......@@ -104,7 +104,7 @@ class Darcy(unittest.TestCase) :
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
ii = 0
fluid = mbfluid.FluidProblem(2,g,mu,rho,order=[2,2,1], solver="petsc4py", solver_options="-pc_type lu", stability=False)#, stability=True, drag_in_stab=True)
fluid = mbfluid.FluidProblem(2,g,mu,rho,order=[2,2,1], solver="petsc4py", solver_options="-pc_type lu")#, stability=True, drag_in_stab=True)
fluid.load_msh("mesh.msh")
U = 0.00001 # averaged fluid velocity
V = U/(1-compacity) # fluid velocity
......@@ -123,8 +123,7 @@ class Darcy(unittest.TestCase) :
#Computation loop
while t < tEnd :
#Fluid solver
time_integration.iterate(fluid,p,dt,fixed_grains=True,check_residual_norm=1e-4)
time_integration.iterate(fluid,p,dt,fixed_grains=True,check_residual_norm=1)
t += dt
#Output files writting
if ii %outf == 0 :
......
......@@ -22,6 +22,7 @@
#!/usr/bin/env python
from migflow import fluid as mbfluid
from migflow import time_integration
from migflow import free_surface
import numpy as np
import os
......@@ -33,8 +34,7 @@ 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)
lapl_h = free_surface.compute_lapl_fs(fluid.coordinates(), fs_nodes)
n_e = len(lapl_h)-1
kappa = np.zeros_like(x[:,0])
X = fluid.coordinates()[fs_nodes,0]
......@@ -86,7 +86,7 @@ class DieSwell(unittest.TestCase):
#
# FLUID PROBLEM
#
fluid = mbfluid.FluidProblem(2,[0,g],mu,rho, order=[2,2,1], solver_options="-pc_type lu")
fluid = mbfluid.FluidProblem(2,[0,g],mu,rho, order=[1,1,1], solver_options="-pc_type lu")
# 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])
......@@ -99,7 +99,11 @@ class DieSwell(unittest.TestCase):
#
# FREE SURFACE STRUCTURES
#
fs = time_integration.FreeSurface(fluid,"TopRight","BottomRight",isCentered=False)
fs_nodes = np.unique(fluid.mesh_boundaries()[b"TopRight"])
x = fluid.coordinates()
# fs = time_integration.FreeSurface(fluid,"TopRight","BottomRight",isCentered=False)
fluid.write_vtk(outputdir,0,0)
tic = time.time()
......@@ -108,10 +112,13 @@ class DieSwell(unittest.TestCase):
#
swelling = 0
while t < tEnd :
time_integration.iterate(fluid,None,dt,check_residual_norm=1e-6)
if t >= 50.0:
fs.iterate(dt)
swelling = fluid.coordinates()[fs.fs_nodes[-1],1]/H
time_integration.iterate(fluid,None,dt,check_residual_norm=1e-2)
if t >= 10:
xold = np.copy(fluid.coordinates())
x[fs_nodes,1] = free_surface.compute_h_upwind(x[fs_nodes], fluid.velocity()[fs_nodes], dt)
fluid.mesh_velocity()[:,:2] = (x[:,:2]-xold[:,:2])/dt
swelling = x[fs_nodes[-1],1]/H
print("gonflement = ", swelling)
t += dt
# Output files writting
if ii %outf == 0 :
......
......@@ -64,7 +64,7 @@ class Poiseuille(unittest.TestCase) :
outf = 1 # number of iterations between output files
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
fluid = mbfluid.FluidProblem(2,g,nu*rho,rho, order=[2,2,1], solver="petsc4py", solver_options="-pc_type lu", stability=False)
fluid = mbfluid.FluidProblem(2,g,nu*rho,rho, order=[2,2,1], solver="petsc4py", solver_options="-pc_type lu")
fluid.load_msh("mesh.msh")
fluid.set_open_boundary("LeftUp",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
fluid.set_open_boundary("LeftDown",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),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