Cavity (2D)
This notebook implements a steady incompressible Navier-Stokes solver for the cavity problem.
The problem is strong form reads: \begin{equation} \left\{ \begin{array}{ll} \nabla \cdot \mathbf{u} =0& in\;\Omega\\ \displaystyle \left(\mathbf{u}\cdot \nabla\right)\mathbf{u}= \nu\Delta \mathbf{u}-\nabla p & in\;\Omega\\ & \\ \mathbf{u} = \mathbf{i} & on\;\Gamma_{lid}\\ \mathbf{u} = \mathbf{0} & on\;\partial\Omega\setminus\Gamma_{lid} \end{array} \right. \end{equation}
[1]:
import numpy as np
import pandas as pd
# Mesh generation
import gmsh
import dolfinx
from dolfinx.io import gmshio
from mpi4py import MPI
from dolfinx.mesh import (CellType, GhostMode, create_rectangle, locate_entities_boundary)
from dolfinx.fem import Function, locate_dofs_geometrical, dirichletbc
import matplotlib.pyplot as plt
from matplotlib import cm
import sys
mesh_path = '../../../mesh/'
benchmark_path = '../../../BenchmarkData/'
sys.path.append('../../../models/fenicsx')
Mesh generation from gmsh
[2]:
gdim = 2
model_rank = 0
mesh_comm = MPI.COMM_WORLD
mesh_factor = .01
# Initialize the gmsh module
gmsh.initialize()
# Load the .geo file
gmsh.merge(mesh_path+'cavity.geo')
gmsh.model.geo.synchronize()
gmsh.option.setNumber("Mesh.MeshSizeFactor", mesh_factor)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")
# Domain
domain, ct, ft = gmshio.model_to_mesh(gmsh.model, comm = mesh_comm, rank = model_rank, gdim = gdim )
gmsh.finalize()
bound_markers = {'lid': 10,
'walls': 20}
domain_marker = 100
Info : Reading '../../../mesh/cavity.geo'...
Info : Done reading '../../../mesh/cavity.geo'
Info : Meshing 1D...
Info : [ 0%] Meshing curve 1 (Line)
Info : [ 30%] Meshing curve 2 (Line)
Info : [ 50%] Meshing curve 3 (Line)
Info : [ 80%] Meshing curve 4 (Line)
Info : Done meshing 1D (Wall 0.00200965s, CPU 0.000675s)
Info : Meshing 2D...
Info : Meshing surface 1 (Transfinite)
Info : Done meshing 2D (Wall 0.00728691s, CPU 0.003989s)
Info : 10201 nodes 20404 elements
Info : Optimizing mesh (Netgen)...
Info : Done optimizing mesh (Wall 1.98001e-06s, CPU 2e-06s)
The problem we want to face is non-linear, whose weak formulation reads: \begin{equation} \int_\Omega \left(\mathbf{u}\cdot \nabla\right)\mathbf{u}\cdot \mathbf{v}\,d\Omega + \nu \int_\Omega\nabla \mathbf{u}\cdot \nabla \mathbf{v}\,d\Omega -\int_\Omega p(\nabla\cdot\mathbf{v})\,d\Omega -\int_\Omega q(\nabla\cdot\mathbf{u})\,d\Omega=0 \end{equation}
[5]:
from fluid_dynamics.steady_ns import ns_steady_nl
cavity = ns_steady_nl(domain, ct = ct, ft = ft, bound_markers=bound_markers)
# ReVec = np.array([100, 400, 1000, 2500, 5000, 7500])
ReVec = np.array([100, 400, 1000])
u_list = []
boundary_type = {'lid': 0,
'walls': 1}
boundary_value = {'lid': np.array([1, 0, 0])}
for ii in range(len(ReVec)):
cavity.parameters(1/ReVec[ii])
cavity.set_bc(boundary_type=boundary_type, boundary_value=boundary_value)
dofs_p = locate_dofs_geometrical((cavity.W.sub(1), cavity.W.sub(1).collapse()[0]),
lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
cavity.bcs.append( dirichletbc(Function(cavity.W.sub(1).collapse()[0]),
dofs_p, cavity.W.sub(1)) )
# cavity.set_cavity_bc()
cavity.assemble(maxIter=20, verbose = False)
up_sol = cavity.solve()
(u_sol, p_sol) = (up_sol.sub(0).collapse(), up_sol.sub(1).collapse())
u_list.append(u_sol)
Below a function to extract data from a function is implemented
[6]:
def extract2D(problem, N, u_sol):
grid = np.linspace(0, 1, N)
ux = np.zeros((N,N))
uy = np.zeros((N,N))
for ii in range(N):
points = np.zeros((3, N))
points[0, :] = grid[ii]
points[1, :] = grid
bb_tree = dolfinx.geometry.BoundingBoxTree(problem.domain, problem.domain.topology.dim)
cells = []
points_on_proc = []
cell_candidates = dolfinx.geometry.compute_collisions(bb_tree, points.T)
colliding_cells = dolfinx.geometry.compute_colliding_cells(problem.domain, cell_candidates, points.T)
for i, point in enumerate(points.T):
if len(colliding_cells.links(i))>0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
xPlot = np.array(points_on_proc, dtype=np.float64)
ux[ii, :] = u_sol.sub(0).eval(xPlot, cells).flatten()
uy[ii, :] = u_sol.sub(1).eval(xPlot, cells).flatten()
return ux.T, uy.T, grid
Let’s make some 2D streamline plots
[7]:
nrows = 1
ncols = 3
fontsize = 30
labelsize = 20
fig, axs = plt.subplots(nrows = nrows, ncols = ncols, sharex=True, sharey=True,
figsize = (6 * ncols, 6 * nrows))
axs = axs.flatten()
N = 100
for ii in range(nrows * ncols):
ux, uy, grid = extract2D(cavity, N, u_list[ii])
X, Y = np.meshgrid(grid, grid)
color = np.sqrt(((uy+2)/2)*2 + ((ux+2)/2)*2)
axs[ii].streamplot(X, Y, ux, uy, color=np.sqrt(ux**2+uy**2), linewidth=1, cmap='jet')
axs[ii].set_title(r'$Re = {:.0f}'.format(ReVec[ii])+'$', fontsize = fontsize)
axs[ii].set_xlim(0,1)
axs[ii].set_ylim(0,1)
axs[ii].set_xticks(np.arange(0.25, 1.01, 0.25))
axs[ii].set_yticks(np.arange(0., 1.01, 0.25))
axs[ii].tick_params(axis='both', labelsize=labelsize)
for ii in range(3, len(ReVec)):
axs[ii].set_xlabel(r'$x$', fontsize = fontsize)
axs[0].set_ylabel(r'$y$', fontsize = fontsize)
# axs[3].set_ylabel(r'$y$', fontsize = fontsize)
fig.subplots_adjust(hspace=0.1, wspace = 0)
Comparison with benchmark data
Benchmark data are taken from https://www.acenumerics.com/the-benchmarks.html
[8]:
def extract1D_y(problem, N, u_sol, y_lines: list, component = 0):
grid = np.linspace(0, 1, N)
u_extracted = np.zeros((len(y_lines), N))
for idx_y in range(len(y_lines)):
points = np.zeros((3, N))
points[0, :] = grid
points[1, :] = y_lines[idx_y]
bb_tree = dolfinx.geometry.BoundingBoxTree(problem.domain, problem.domain.topology.dim)
cells = []
points_on_proc = []
cell_candidates = dolfinx.geometry.compute_collisions(bb_tree, points.T)
colliding_cells = dolfinx.geometry.compute_colliding_cells(problem.domain, cell_candidates, points.T)
for i, point in enumerate(points.T):
if len(colliding_cells.links(i))>0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
xPlot = np.array(points_on_proc, dtype=np.float64)
if np.isclose(component, 0):
u_extracted[idx_y, :] = u_sol.sub(0).eval(xPlot, cells).flatten()
else:
u_extracted[idx_y, :] = u_sol.sub(1).eval(xPlot, cells).flatten()
return u_extracted, xPlot[:,0]
def extract1D_x(problem, N, u_sol, x_lines: list, component = 0):
grid = np.linspace(0, 1, N)
u_extracted = np.zeros((len(x_lines), N))
for idx_x in range(len(x_lines)):
points = np.zeros((3, N))
points[0, :] = x_lines[idx_x]
points[1, :] = grid
bb_tree = dolfinx.geometry.BoundingBoxTree(problem.domain, problem.domain.topology.dim)
cells = []
points_on_proc = []
cell_candidates = dolfinx.geometry.compute_collisions(bb_tree, points.T)
colliding_cells = dolfinx.geometry.compute_colliding_cells(problem.domain, cell_candidates, points.T)
for i, point in enumerate(points.T):
if len(colliding_cells.links(i))>0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
xPlot = np.array(points_on_proc, dtype=np.float64)
if np.isclose(component, 0):
u_extracted[idx_x, :] = u_sol.sub(0).eval(xPlot, cells).flatten()
else:
u_extracted[idx_x, :] = u_sol.sub(1).eval(xPlot, cells).flatten()
return u_extracted, xPlot[:,1]
Vertical velocity across various horizontal planes
[11]:
idx_Re = 2
df = pd.read_excel(benchmark_path+'fluid_dynamics/cavity/re_'+str(ReVec[idx_Re])+'_driven_cavity_benchmark_results.xlsx',
sheet_name='Table 6')
y_lines = [0.05, 0.1, 0.5, 0.9, 0.99]
df = np.asarray(df.to_numpy())
x_bench = df[:,0].flatten()
uy_bench = df[:,1:]
del df
uy, x_plot = extract1D_y(cavity, N, u_list[idx_Re], y_lines, component = 1)
###############################################################################################
fontsize = 30
labelsize = 20
nrows, ncols = 2, 3
fig, axs = plt.subplots(nrows = nrows, ncols = ncols, sharex=True,
figsize = (6 * ncols, 5 * nrows))
axs = axs.flatten()
for ii in range(len(y_lines)):
axs[ii+1].plot(x_plot, uy[ii], '-', color=cm.jet(0.1), label='FEniCSx')
axs[ii+1].plot(x_bench, uy_bench[:, ii], 'o', color=cm.jet(0.9), label='Benchmark')
axs[ii+1].grid()
axs[ii+1].set_xlim(0,1)
axs[ii+1].set_ylabel(r'$u_y$ at $y={:.2f}'.format(y_lines[ii])+'$', fontsize=fontsize)
axs[ii+1].tick_params(axis='both', labelsize = labelsize)
axs[ii+1].set_xticks(np.arange(0, 1.01, 0.25))
axs[0].axis('off')
axs[3].set_xlabel(r'$x$', fontsize=fontsize)
axs[4].set_xlabel(r'$x$', fontsize=fontsize)
axs[5].set_xlabel(r'$x$', fontsize=fontsize)
fig.text(0.16, 0.8, s=r'$Re='+str(ReVec[idx_Re])+'$', fontsize=fontsize)
Line, Label = axs[1].get_legend_handles_labels()
fig.legend(Line, Label, framealpha = 1, fontsize=labelsize, loc=(0.1, 0.75), ncols=1)
fig.subplots_adjust(hspace=0.075, wspace=0.4)
Horizontal velocity across various vertical planes
[12]:
idx_Re = 2
df = pd.read_excel(benchmark_path+'fluid_dynamics/cavity/re_'+str(ReVec[idx_Re])+'_driven_cavity_benchmark_results.xlsx',
sheet_name='Table 9')
x_lines = [0.05, 0.1, 0.5, 0.9, 0.95]
df = np.asarray(df.to_numpy())
y_bench = df[:,0].flatten()
ux_bench = df[:,1:]
del df
ux, y_plot = extract1D_x(cavity, N, u_list[idx_Re], x_lines, component = 0)
###############################################################################################
fontsize = 30
labelsize = 20
nrows, ncols = 2, 3
fig, axs = plt.subplots(nrows = nrows, ncols = ncols, sharex=True,
figsize = (6 * ncols, 5 * nrows))
axs = axs.flatten()
for ii in range(len(x_lines)):
axs[ii+1].plot(y_plot, ux[ii], '-', color=cm.jet(0.1), label='FEniCSx')
axs[ii+1].plot(y_bench, ux_bench[:, ii], 'o', color=cm.jet(0.9), label='Benchmark')
axs[ii+1].grid()
axs[ii+1].set_xlim(0,1)
axs[ii+1].set_ylabel(r'$u_x$ at $x={:.2f}'.format(x_lines[ii])+'$', fontsize=fontsize)
axs[ii+1].tick_params(axis='both', labelsize = labelsize)
axs[0].axis('off')
axs[3].set_xlabel(r'$y$', fontsize=fontsize)
axs[4].set_xlabel(r'$y$', fontsize=fontsize)
fig.text(0.175, 0.8, s=r'$Re='+str(ReVec[idx_Re])+'$', fontsize=fontsize)
Line, Label = axs[1].get_legend_handles_labels()
fig.legend(Line, Label, framealpha = 1, fontsize=labelsize, loc=(0.1, 0.75), ncols=1)
fig.subplots_adjust(hspace=0.0, wspace=0.35)
