Seg error with Netwon Solverq

I am trying to solve the tutorial problem of
“Flow past a cylinder (DFG 2D-3 benchmark)”
using Newton solver to get solution of steady state Navier Stokes equation.
But I am getting segmentation error upon invoking the solve function.
The geometry generation is same as the tutorial i.e.

``````# +
# # Modified
# the Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark)  by Jørgen S. Dokken
# to solve steady state by Newtons iteration
import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import tqdm.autonotebook

from mpi4py import MPI
from petsc4py import PETSc

from basix.ufl import (element,mixed_element)
from dolfinx.io import (gmshio)
from dolfinx.fem import (Constant, Function, functionspace,dirichletbc,locate_dofs_topological,form)

from dolfinx.cpp.mesh import to_type, cell_entity_type

from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
create_vector, create_matrix, set_bc)
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from dolfinx.mesh import create_mesh, meshtags_from_entities

#for nonlinear solution
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx import log

gmsh.initialize()

L = 2.2
H = 0.41
c_x = c_y = 0.2
r = 0.05
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
if mesh_comm.rank == model_rank:
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)
# -

# The next step is to subtract the obstacle from the channel, such that we do not mesh the interior of the circle.
#

if mesh_comm.rank == model_rank:
fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
gmsh.model.occ.synchronize()

# To get GMSH to mesh the fluid, we add a physical volume marker
#

fluid_marker = 1
if mesh_comm.rank == model_rank:
volumes = gmsh.model.getEntities(dim=gdim)
assert (len(volumes) == 1)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")

# To tag the different surfaces of the mesh, we tag the inflow (left hand side) with marker 2, the outflow (right hand side) with marker 3 and the fluid walls with 4 and obstacle with 5. We will do this by computing the center of mass for each geometrical entity.
#

inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
if mesh_comm.rank == model_rank:
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
if np.allclose(center_of_mass, [0, H / 2, 0]):
inflow.append(boundary[1])
elif np.allclose(center_of_mass, [L, H / 2, 0]):
outflow.append(boundary[1])
elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]):
walls.append(boundary[1])
else:
obstacle.append(boundary[1])
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")

# In our previous meshes, we have used uniform mesh sizes. In this example, we will have variable mesh sizes to resolve the flow solution in the area of interest; close to the circular obstacle. To do this, we use GMSH Fields.
#

# Create distance field from obstacle.
# Add threshold of mesh sizes based on the distance field
# LcMax -                  /--------
#                      /
# LcMin -o---------/
#        |         |       |
#       Point    DistMin DistMax
res_min = r / 3
if mesh_comm.rank == model_rank:
gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle)
gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.25 * H)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2 * H)
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)

# ## Generating the mesh
#
# We are now ready to generate the mesh. However, we have to decide if our mesh should consist of triangles or quadrilaterals. In this demo, to match the DFG 2D-3 benchmark, we use second order quadrilateral elements.
#

if mesh_comm.rank == model_rank:
gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")

#
# As we have generated the mesh, we now need to load the mesh and corresponding facet markers into DOLFINx.
# To load the mesh, we follow the same structure as in [Deflection of a membrane](./../chapter1/membrane_code.ipynb), with the difference being that we will load in facet markers as well. To learn more about the specifics of the function below, see [A GMSH tutorial for DOLFINx](https://jsdokken.com/src/tutorial_gmsh.html).
#

mesh, _, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
ft.name = "Facet markers"

# ## Physical and discretization parameters
#
# Following the DGF-2 benchmark, we define our problem specific parameters
#

mu = Constant(mesh, PETSc.ScalarType(0.001))  # Dynamic viscosity
rho = Constant(mesh, PETSc.ScalarType(1))     # Density

##### Added boundary conditions and Newton solver to solve steady state equations

# Domain
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank=0
#mesh, _, ft =gmshio.read_from_msh("cylinder.msh", mesh_comm, model_rank, gdim)
mesh=mesh
ft=ft
ft.name = "Facet markers"
fdim = mesh.topology.dim - 1

# Functional Spaces

P2 = element("Lagrange", mesh.topology.cell_name(), 2,shape=(mesh.geometry.dim, ))
P1=s_cg1 = element("Lagrange", mesh.topology.cell_name(), 1)
mixEl = mixed_element([P2, P1])
#combined space
W = functionspace(mesh, mixEl)
# Test and trial functions: monolithic
(v, q) = TestFunctions(W)
dup = TrialFunction(W)
up  = Function(W)
(u, p) = split(up)

# Computing the viscosity
# ## Physical and discretization parameters
Re=20
InvReNo=1.0/Re
rho = Constant(mesh, PETSc.ScalarType(1))     # Density
InvRe=Constant(mesh, PETSc.ScalarType(InvReNo))
print('Re = ',format(Re) )
``````

until here the code is almost same as tutorial. I just added different boundary conditions and solution method using Newton solver as follows:

``````
# Setting Boundary Conditions
##nonslip on obstacles. no penetration at walls normal gradient of velocity at wall and outlet is zero
##pressure is zero at outlet
# inlet velocity is uniform and =(1,0)
#(**markers**)
inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
#inlet
u_inlet=np.array([1.0,0.0], dtype=PETSc.ScalarType)
dofs_inlet=locate_dofs_topological(W.sub(0), fdim, ft.find(inlet_marker))
bcu_inflow = dirichletbc(u_inlet,dofs_inlet ,W.sub(0).collapse()[0])
# Walls
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
dofs_walls=locate_dofs_topological(W.sub(0), fdim, ft.find(wall_marker))
bcu_walls = dirichletbc(u_nonslip, dofs_walls, W.sub(0).collapse()[0])
# Obstacle
dofs_obstacle=locate_dofs_topological(W.sub(0), fdim, ft.find(obstacle_marker))
bcu_obstacle = dirichletbc(u_nonslip, dofs_obstacle, W.sub(0).collapse()[0])
# Outlet pressure bc
dofs_out=locate_dofs_topological(W.sub(1), fdim, ft.find(outlet_marker))
bcp_outlet = dirichletbc(PETSc.ScalarType(0), dofs_out, W.sub(1).collapse()[0])

#combined bc
bcs1=[bcu_inflow, bcu_obstacle, bcu_walls, bcp_outlet ]
##########################
# Variational forms
F = dot(dot(u, nabla_grad(u)), v) * dx
- div(v)*p * dx  #pressure grad term
+q*div(u) * dx ##continuity equation

_F = form(F)

# Solver settings
problem = NonlinearProblem(F, up, bcs=bcs1)

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-3
solver.atol= 1e-3
solver.report = True

# We can modify the linear solver in each Newton iteration by accessing the underlying `PETSc` object.

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"

ksp.setFromOptions()

# We are now ready to solve the non-linear problem. We assert that the solver has converged and print the number of iterations.
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(up)
assert (converged)
print(f"Number of interations: {n:d}")
``````

Any help regarding this is appreciated