Hi everybody, I’m trying to solve the lid driven cavity problem with dolfinx version 0.7.3.
In some sense it works but I wonder if someone can solve some doubts I have.
- How can I plot at different time steps (for example every 20 time steps)? I tried to delete the plot after visualization with plotter.delete() and make him wait but it can’t read following lines until the plot is closed.
- When I run it, greater values of magnitude of velocity are enormous (from 10^5 to 10^6) and I can’t understand why (the driven velocity has magnitude 1).
- I can’t see the typical vortex at the center of the plot. I tried to reduce meshing size (maybe not enough). Are there any other problems?
I attach the code for completeness.
"""
Equations:
Momentum: du/dt + (u nabla)u = -(1/rho)(nabla p) + nu (laplacian u)
Incompressibility: nabla u = 0
Notations:
u: the velocity field from the previous iteration
u_: the velocity field after a tentative iteration step
u__: the velocity field after incompressibility correction (= after projection)
p_: the pressure field at the next step
dt: the time step size
v: the test function for u
q: the test function for p
Theoretical tips:
I recall that inner product in FEM strategy means integrate the product of the functions over the domain.
Because of the strategy used (Taylor-Hood elements) we chose the degree of the velocity function space
one order higher than the degree of the pressure function space.
The chioce of the time step size has to be careful because of the hybrid explicit-implicit method used.
Solution strategy:
1. Solve the weak form of the momentum equation without the pressure gradient.
Treat the convection explicity and the diffusion implicity. Use an Euler time step.
Solve for u_
0 = 1/dt <(u_ - u), v> + <(nabla u)*u, v> + nu <(nabla u_), (nabla v)>
2. Solve the weak form of the pressure poisson equation with the divergence of the tentative
velocity field as the right-hand side. Solve for p_
<(nabla p_), (nabla q)> = - rho/dt <(nabla u_), q>
3. Solve the weak form of the incompressibility correction. Solve for u__.
rho*<u__, v> = rho*<u_, v> - (1/dt) <(nabla p_), v>
"""
import time
from mpi4py import MPI
from petsc4py import PETSc
import matplotlib as plt
import numpy as np
from tqdm import tqdm
import pyvista
import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem
from dolfinx.fem import (
Constant,
Function,
dirichletbc,
extract_function_spaces,
form,
functionspace,
locate_dofs_topological,
)
from dolfinx.fem.petsc import (assemble_matrix_block, assemble_vector_block, LinearProblem, NonlinearProblem,
assemble_matrix, assemble_vector, create_vector, apply_lifting, set_bc)
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner, TestFunction, TrialFunction, lhs, rhs, nabla_grad
from dolfinx.plot import vtk_mesh
from dolfinx.io import VTXWriter
"""Create the mesh"""
msh = create_rectangle(
MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [20, 20], CellType.triangle
)
"""Taylor-Hood elements and solution fields"""
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)
V, Q = functionspace(msh, P2), functionspace(msh, P1)
u_trial = TrialFunction(V)
p_trial = TrialFunction(Q)
u_prev = Function(V)
u_tent = Function(V)
u_next = Function(V)
p_next = Function(Q)
v_test = TestFunction(V)
p_test = TestFunction(Q)
t = 0
T = 10
num_step = 1000
dt = T/num_step
nu = 1
rho = 1
"""Define the boundary conditions"""
def noslip_boundary(x):
return np.logical_or(
np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)), np.isclose(x[1], 0.0)
)
def lid(x):
return np.isclose(x[1], 1.0)
def lid_velocity_expression(x):
return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType) # type: ignore
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)
lid_velocity = Function(V)
lid_velocity.interpolate(lid_velocity_expression)
facets = locate_entities_boundary(msh, 1, lid)
bc1 = dirichletbc(lid_velocity, locate_dofs_topological(V, 1, facets))
velocity_boundary_conditions = [bc0, bc1]
"""Define the variational problem """
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0))) # type: ignore
"""First step"""
Step1 = (1/dt) * inner(u_trial - u_prev, v_test) * dx
Step1 += inner(nabla_grad(u_prev) * u_prev, v_test) * dx
Step1 += nu * inner(nabla_grad(u_trial), nabla_grad(v_test)) * dx
momentum_weak_lhs = form(lhs(Step1))
momentum_weak_rhs = form(rhs(Step1))
A1 = assemble_matrix(momentum_weak_lhs, velocity_boundary_conditions)
A1.assemble()
B1 = create_vector(momentum_weak_rhs)
"""Second step"""
Step2 = inner(nabla_grad(p_trial), nabla_grad(p_test)) * dx
Step2 += - (rho/dt) * inner(div(u_tent), p_test) * dx
pressure_poisson_lhs = form(lhs(Step2))
pressure_poisson_rhs = form(rhs(Step2))
A2 = assemble_matrix(pressure_poisson_lhs)
A2.assemble()
B2 = create_vector(pressure_poisson_rhs)
"""Third step"""
Step3 = inner(u_trial, v_test) * dx
Step3 += inner(u_tent, v_test) * dx
Step3 += - (1/(rho*dt)) * inner(nabla_grad(p_next), v_test) * dx
velocity_update_weak_lhs = form(lhs(Step3))
velocity_update_weak_rhs = form(rhs(Step3))
A3 = assemble_matrix(velocity_update_weak_lhs)
A3.assemble()
B3 = create_vector(velocity_update_weak_rhs)
"""Create solvers"""
solver1 = PETSc.KSP().create(msh.comm) #create the solver
solver1.setOperators(A1) #tell the solver what matrix it's going to solve
solver1.setType(PETSc.KSP.Type.BCGS) #the next 4 lines set the preconditioner
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")
solver2 = PETSc.KSP().create(msh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.CG)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.SOR)
solver3 = PETSc.KSP().create(msh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)
from pathlib import Path #make a folder where we are going to write
# other lines are just some consequences
folder = Path("results")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(msh.comm, folder / "poiseuille_u.bp", u_prev, engine="BP4")
vtx_p = VTXWriter(msh.comm, folder / "poiseuille_p.bp", u_prev, engine="BP4")
vtx_u.write(t)
vtx_p.write(t)
topology, cell_types, geometry = vtk_mesh(V) #computing the mesh on which vectors will live
values = np.zeros((geometry.shape[0], 3), dtype=np.float64) #making values tab and
#inizializing it to 0
function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(*vtk_mesh(msh, msh.topology.dim))
plotter.add_mesh(grid, style="wireframe", color="k")
plotter.view_xy()
for t in tqdm(range(num_step)):
"""(1) solve for the tentative velocity"""
with B1.localForm() as loc_1:
loc_1.set(0)
assemble_vector(B1, momentum_weak_rhs)
apply_lifting(B1, [momentum_weak_lhs], [velocity_boundary_conditions])
B1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(B1, velocity_boundary_conditions)
solver1.solve(B1, u_tent.vector) #in u_tent the solution is memorized
#u_tent.x.scatter_forward() #I can't understand what it does
"""(2) solve for the pressure"""
with B2.localForm() as loc_2:
loc_2.set(0)
assemble_vector(B2, pressure_poisson_rhs)
B2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
solver2.solve(B2, p_next.vector)
"""(3) solve for """
with B3.localForm() as loc_3:
loc_3.set(0)
assemble_vector(B3, velocity_update_weak_rhs)
B3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
solver3.solve(B3, u_next.vector)
"""Update variable with solution form this time step"""
u_prev = u_next
"""Interactive visualization"""
values[:, :len(u_next)] = u_next.x.array.real.reshape((geometry.shape[0], len(u_next)))
#selecting first len(u_n) columns and assigning u_next values
function_grid["u"] = values * pow(10, -3.5)
"""Advancing in time"""
t += dt
glyphs = function_grid.glyph(orient="u", factor=0.1)
plotter.add_mesh(glyphs)
if not pyvista.OFF_SCREEN:
plotter.show()
else:
fig_as_array = plotter.screenshot("glyphs.png")
Thanks in advance.