Hi,
we are experiencing an unexpected issue in implementing a solution for a benchmark problem in the context of Stokes Equations (flow in a driven cavity, slipping condition on the upper boundary with u=1).
We built the weak formulation and the solver following indications in the Taylor-Hood tutorial and imposed the value of pressure on (0,0) in order to make the problem well-posed. Unfortunately, the results we are having is still nan for both pressure and velocities.
Could you help us detect the issue in our code?
#importing the packages
import numpy as np
# import ufl
from dolfinx import fem, io, mesh, plot
from dolfinx.fem import (Function, FunctionSpace, dirichletbc, form,
locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.mesh import (CellType, GhostMode, create_rectangle,
locate_entities_boundary)
from ufl import ds, dx, SpatialCoordinate, \
sin, cos, dot, \
grad, inner, div, FacetNormal, \
TrialFunctions, TestFunctions, as_vector, as_tensor,\
VectorElement, FiniteElement, CellDiameter
import matplotlib.pyplot as plt
from mpi4py import MPI
from petsc4py import PETSc
import gmsh
from dolfinx.io import gmshio
######################################
# Lid-Driven Cavity Scenario:
# ------>>>>> u_top
# border 3
# 1 +-------------------------------------------------+
# | |
# | * * |
# | * * * * |
# 0.8 | |
# | * |
# | * * |
# | * * |
# 0.6 | * |
# u = 0 | * * | u = 0
# v = 0 | * | v = 0
# border 4 | * | border 2
# | * * * |
# 0.4 | |
# | |
# | * * * |
# | * * |
# 0.2 | * * |
# | * |
# | * * * * * |
# | * |
# 0 +-------------------------------------------------+
# 0 0.2 0.4 0.6 0.8 1
# u = 0
# v = 0
# border 1
###########################
########### Creation of the mesh ##################
# DEFINITION OF THE MESH - type 1
meshsize= 0.1
gmsh.initialize()
gmsh.option.setNumber("Mesh.Algorithm", 1)
gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize)
gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize)
gmsh.model.occ.addRectangle(0, 0, 0, 1, 1, tag=1)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], tag=2)
gmsh.model.mesh.generate(2)
msh, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, MPI.COMM_SELF, 0, gdim=2)
x = SpatialCoordinate(msh)
#####################################
########### Definition of the functional spaces ###############
P2 = VectorElement("CG", msh.ufl_cell(), 2)
P1 = FiniteElement("CG",msh.ufl_cell(), 1)
V0, Q0 = FunctionSpace(msh, P2), FunctionSpace(msh, P1)
MIXED_SPACE = fem.FunctionSpace(msh, P2*P1)
V, _ = MIXED_SPACE.sub(0).collapse()
Q, _ = MIXED_SPACE.sub(1).collapse()
# Definition of fixed boundaries
def fixed_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))
# Definition of the lid
def lid(x):
return np.isclose(x[1], 1.0)
####### Definition of the boundary conditions #########
# Lid boundary condition
class LidVelocity():
def __call__(self, x):
values = np.zeros((2, x.shape[1]),dtype=PETSc.ScalarType)
values[0] = np.ones((1,x.shape[1]), dtype=PETSc.ScalarType)
values[1]= np.zeros((1,x.shape[1]), dtype=PETSc.ScalarType)
return values
lid_velocity= LidVelocity()
# interpolation of the given velocity in the FEM space
lid_u= fem.Function(V)
lid_u.interpolate(lid_velocity)
facets_lid = mesh.locate_entities_boundary(msh, 1, lid)
# set the non-homogeneous dirichlet boundary condition
dofs_lid = fem.locate_dofs_topological((MIXED_SPACE.sub(0), V0), 1, facets_lid)
bc1 = fem.dirichletbc(lid_u, dofs_lid, MIXED_SPACE.sub(0))
# No slip boundary condition
fixed = fem.Function(V) # QUESTO IN AUTOMATICO TI METTE NOSLIP = 0?
facets = mesh.locate_entities_boundary(msh, 1, fixed_boundary)
dofs = fem.locate_dofs_topological((MIXED_SPACE.sub(0), V0), 1, facets)
bc0 = fem.dirichletbc(fixed, dofs, MIXED_SPACE.sub(0))
# Since for this problem the pressure is only determined up to a
# constant, we pin the pressure at the point (0, 0)
zero = Function(Q)
with zero.vector.localForm() as zero_local:
zero_local.set(0.0)
dofs = locate_dofs_geometrical((MIXED_SPACE.sub(1), Q),
lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = dirichletbc(zero, dofs, MIXED_SPACE.sub(1))
#define the final boundary conditions
bcs = [bc0, bc1, bc2]
################ Definition of the variational problem #####################
(u, p) = TrialFunctions(MIXED_SPACE)
(v, q) = TestFunctions(MIXED_SPACE)
zero = fem.Constant(msh, 0.0)
f = fem.Function(V)
# f.x.array[:] = 0.0
f= fem.Constant(msh, (0.0, 0.0))
a = fem.form((inner(grad(u), grad(v)) - div(v)*p - q*div(u)) * dx)
L = fem.form((inner(f, v))*dx)
# Assemble LHS matrix and RHS vector
A = fem.petsc.assemble_matrix(a, bcs=bcs)
A.assemble()
b = fem.petsc.assemble_vector(L)
# Implement boundary conditions in the right hand side
fem.petsc.apply_lifting(b, [a], bcs=[bcs]) # inhomgeneous dirichlet boundary
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) #update right hand side
fem.petsc.set_bc(b, bcs) # Set Dirichlet boundary condition values in the RHS
# Create and configure solver
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("superlu_dist")
# compute the solution and split it
U = fem.Function(MIXED_SPACE)
ksp.solve(b, U.vector)
# Split the mixed solution and collapse
u_h = U.sub(0).collapse()
p_h = U.sub(1).collapse()
uh_x, uh_y = u_h.split()
# Compute norms
norm_u_3 = u_h.vector.norm()
norm_p_3 = p_h.vector.norm()
if MPI.COMM_WORLD.rank == 0:
print("(D) Norm of velocity coefficient vector (monolithic, direct): {}".format(norm_u_3))
print("(D) Norm of pressure coefficient vector (monolithic, direct): {}".format(norm_p_3))
assert np.isclose(norm_u_3, norm_p_3)
################ PLOTTING THE PRESSURE ###################################################
import pyvista
cells, types, geometry = plot.create_vtk_mesh(Q0)
grid = pyvista.UnstructuredGrid(cells, types, geometry)
grid.point_data["p"] = p_h.x.array.real
grid.set_active_scalars("p")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=False)
plotter.view_xy()
plotter.show(screenshot='pressure.png')
###########################################################################################
##### SAVING THE DATA
# Plotting the mesh and solution, copy+paste for the exercises!
with io.XDMFFile(msh.comm, "pressure_out.xdmf", "w") as file:
file.write_mesh(msh)
file.write_function(p_h)
with io.XDMFFile(msh.comm, "velocity_out.xdmf", "w") as file:
file.write_mesh(msh)
file.write_function(u_h)
print('Computation is over! Success!')