Dear community,
I tried to adapt this Navier-Stokes example to a similar scenario with an L-shaped domain. The velocities are very turbulent and many velocities have an extreme value compared to their neighbours (see Fig. 1). Later on in the simulation all velocities are NaN (see Fig. 2). The boundary conditions are trivial, which means it can’t be the Gibbs phenomenon mentioned in this thread.
I don’t know what causes that bug and thus clueless how to fix it.
Many thanks in advance for your response!
Code to reproduce:
import pyvista
import os
import numpy as np
from petsc4py import PETSc
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, create_vector, create_matrix, set_bc)
from dolfinx.mesh import locate_entities_boundary
from dolfinx.plot import create_vtk_mesh
from ufl import (FiniteElement, TestFunction, TrialFunction, VectorElement, div, dot, dx, inner, lhs, grad, nabla_grad, rhs)
from dolfinx import io
from mpi4py import MPI
import pygmsh
import meshio
with pygmsh.geo.Geometry() as geom:
geom.add_polygon(
[
[0.0, 0.0],
[0.0, 1.0],
[0.5, 1.0],
[0.5, 0.5],
[1.0, 0.5],
[1.0, 0.0],
],
mesh_size=0.1,
)
mesh = geom.generate_mesh()
cells = [("triangle", mesh.cells_dict['triangle'])]
mesh = meshio.Mesh(mesh.points[:, :2], cells)
meshio.write("lshape.xdmf", mesh)
with io.XDMFFile(MPI.COMM_WORLD, "lshape.xdmf", 'r') as infile:
mesh = infile.read_mesh(name="Grid")
m_grid = pyvista.UnstructuredGrid(*create_vtk_mesh(mesh))
plotter = pyvista.Plotter()
plotter.add_mesh(m_grid, show_edges=True, opacity=0.5)
plotter.show_bounds(
grid='front',
location='outer',
all_edges=True)
plotter.view_xy()
plotter.show()
# Problem parameters
t = 0
T = 3 # Final time
dt = 0.1 # Time step size
num_steps = int(T / dt)
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(0.001)) # Dynamic viscosity
rho = Constant(mesh, PETSc.ScalarType(1)) # Density
v_cg2 = VectorElement("CG", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)
fdim = mesh.topology.dim - 1
gdim = 2
# Define boundary conditions
class InletPressure():
def __init__(self, t):
self.t = t
def __call__(self, x):
values = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = np.sin(3.0 * self.t)
return values
DOLFIN_EPS = 0.000001
def inflow(x):
return x[1] > 1 - DOLFIN_EPS
def outflow(x):
return x[0] > 1 - DOLFIN_EPS
def boundary(x):
return np.logical_or(x[0] < DOLFIN_EPS,
x[1] < DOLFIN_EPS,
np.logical_and(x[0] > 0.5 - DOLFIN_EPS, x[1] > 0.5 - DOLFIN_EPS))
# Inlet
p_inlet = Function(Q)
inlet_pressure = InletPressure(t)
p_inlet.interpolate(inlet_pressure)
inflow_facets = locate_entities_boundary(mesh, fdim, inflow)
bcu_inflow = dirichletbc(p_inlet, locate_dofs_topological(Q, fdim, inflow_facets))
# Walls
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bndry_facets = locate_entities_boundary(mesh, fdim, boundary)
bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, bndry_facets), V)
bcu = [bcu_inflow, bcu_walls]
# Outlet
outflow_facets = locate_entities_boundary(mesh, fdim, outflow)
bcp_outlet = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(Q, fdim, outflow_facets), Q)
bcp = [bcp_outlet]
# Define variables
u = TrialFunction(V)
v = TestFunction(V)
u_ = Function(V)
u_.name = "u"
u_s = Function(V)
u_n = Function(V)
u_n1 = Function(V)
p = TrialFunction(Q)
q = TestFunction(Q)
p_ = Function(Q)
p_.name = "p"
phi = Function(Q)
f = Constant(mesh, PETSc.ScalarType((0, 0)))
F1 = rho / k * dot(u - u_n, v) * dx
F1 += inner(dot(1.5 * u_n - 0.5 * u_n1, 0.5 * nabla_grad(u + u_n)), v) * dx
F1 += 0.5 * mu * inner(grad(u + u_n), grad(v)) * dx - dot(p_, div(v)) * dx
F1 += dot(f, v) * dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
A1 = create_matrix(a1)
b1 = create_vector(L1)
a2 = form(dot(grad(p), grad(q)) * dx)
L2 = form(-1 / k * dot(div(u_s), q) * dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)
a3 = form(dot(u, v) * dx)
L3 = form(dot(u_s, v) * dx - k * dot(nabla_grad(phi), v) * dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)
# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.JACOBI)
# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.MINRES)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")
# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)
a3 = form(dot(u, v) * dx)
L3 = form(dot(u_s, v) * dx - k * dot(nabla_grad(phi), v) * dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)
u_out = './u-shapeL'
p_out = './p-shapeL'
if not os.path.exists(u_out):
os.makedirs(u_out)
if not os.path.exists(p_out):
os.makedirs(p_out)
topology, cell_types, geometry = create_vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
p_grid = pyvista.UnstructuredGrid(*create_vtk_mesh(Q))
for i in range(num_steps):
# Update current time step
t += dt
# Update inlet velocity
inlet_pressure.t = t
p_inlet.interpolate(inlet_pressure)
# Step 1: Tentative velocity step
A1.zeroEntries()
assemble_matrix(A1, a1, bcs=bcu)
A1.assemble()
with b1.localForm() as loc:
loc.set(0)
assemble_vector(b1, L1)
apply_lifting(b1, [a1], [bcu])
b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b1, bcu)
solver1.solve(b1, u_s.vector)
u_s.x.scatter_forward()
# Step 2: Pressure corrrection step
with b2.localForm() as loc:
loc.set(0)
assemble_vector(b2, L2)
apply_lifting(b2, [a2], [bcp])
b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b2, bcp)
solver2.solve(b2, phi.vector)
phi.x.scatter_forward()
p_.vector.axpy(1, phi.vector)
p_.x.scatter_forward()
# Step 3: Velocity correction step
with b3.localForm() as loc:
loc.set(0)
assemble_vector(b3, L3)
b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
solver3.solve(b3, u_.vector)
u_.x.scatter_forward()
# Create plotter
u_plotter = pyvista.Plotter(off_screen=True)
p_plotter = pyvista.Plotter(off_screen=True)
p_plotter.add_title("Frame {}".format(i))
p_plotter.add_mesh(p_grid, style="wireframe", color="k")
p_grid["p"] = p_.x.array.real
p_plotter.add_mesh(p_grid, cmap='viridis')
p_plotter.show_bounds(
grid='front',
location='outer',
all_edges=True)
values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
values[:, :len(u_)] = u_.x.array.real.reshape((geometry.shape[0], len(u_)))
u_grid["u"] = values
glyphs = u_grid.glyph(orient="u", scale=False, factor=0.1)
u_plotter.add_title("Frame {}".format(i))
u_plotter.add_mesh(glyphs, cmap='jet')
u_plotter.show_bounds(
grid='front',
location='outer',
all_edges=True)
u_plotter.view_xy()
u_plotter.screenshot("{}/{}.png".format(u_out, i))
p_plotter.view_xy()
p_plotter.screenshot("{}/{}.png".format(p_out, i))
# Update variable with solution form this time step
with u_.vector.localForm() as loc_, u_n.vector.localForm() as loc_n, u_n1.vector.localForm() as loc_n1:
loc_n.copy(loc_n1)
loc_.copy(loc_n)
Figure 1:
Figure 2:
