Solving Stokes equation with NewtonSolver

dolfinx
4

There are several issues with your implementation

  1. The sign of the second term in the variational formulation. You have integrated both terms by part, so it should have the opposite sign
  2. Setting of boundary conditions. You are setting them on a collapsed sub space, not the mixed space.
    Here is a functional version of the code:
import numpy as np
import ufl
from ufl import inner, grad, dx, div, VectorElement, FiniteElement, TestFunctions
import dolfinx
from dolfinx import default_scalar_type
from dolfinx import mesh, fem, io
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc

# %% mesh and function spaces
domain = mesh.create_rectangle(MPI.COMM_WORLD, [(0.0, 0.0), (1.0, 1.0)],
                               (15, 15), mesh.CellType.quadrilateral)

# combined function space, taylor hood elements with 2nd order for momentum
velocity_element = VectorElement("Lagrange", domain.ufl_cell(), degree=2)
pressure_element = FiniteElement("Lagrange", domain.ufl_cell(), degree=1)
space = fem.functionspace(domain, velocity_element*pressure_element)

v, q = TestFunctions(space)

# w should hold the solution
w = fem.Function(space)
u, p = ufl.split(w)

# %% boundary conditions
U, _ = space.sub(0).collapse()
Q = space.sub(1)


def left(x):
    return np.isclose(x[0], 0.0)


def right(x):
    return np.isclose(x[0], 1.0)


def top(x):
    return np.isclose(x[1], 1.0)


def bottom(x):
    return np.isclose(x[1], 0.0)


def corner(x):
    return np.logical_and(np.isclose(x[1], 0.0), np.isclose(x[0], 0.0))


bc = []  # list of Dirichlet boundary conditions

# velocities for cavity flow
noslip = dolfinx.fem.Function(U)
noslip.x.array[:] = 0
fdim = domain.topology.dim-1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
dofs = fem.locate_dofs_topological((space.sub(0), U), fdim, left_facets)
bc.append(fem.dirichletbc(noslip, dofs, space.sub(0)))
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
dofs = fem.locate_dofs_topological((space.sub(0), U), fdim, right_facets)
bc.append(fem.dirichletbc(noslip, dofs, space.sub(0)))
bottom_facets = mesh.locate_entities_boundary(domain, fdim, bottom)
dofs = fem.locate_dofs_topological((space.sub(0), U), fdim, bottom_facets)
bc.append(fem.dirichletbc(noslip, dofs, space.sub(0)))
top_facets = mesh.locate_entities_boundary(domain, fdim, top)
dofs = fem.locate_dofs_topological((space.sub(0), U), fdim, top_facets)
slip = dolfinx.fem.Function(U)
slip.interpolate(lambda x: (np.ones(x.shape[1]), np.zeros(x.shape[1])))
bc.append(fem.dirichletbc(slip, dofs, space.sub(0)))

# fix pressure in one corner
corner_vertex = mesh.locate_entities_boundary(domain, 0, corner)
dofs = fem.locate_dofs_topological(Q, 0,  corner_vertex)
bc.append(fem.dirichletbc(fem.Constant(domain, 0.0), dofs, Q))

# %% solve the problem
# residual of weak form of Stokes

F = inner(grad(u), grad(v))*dx - p*div(v)*dx - div(u)*q*dx
problem = NonlinearProblem(F, w, bc)
solver = NewtonSolver(domain.comm, problem)

solver.rtol = 1e-6
solver.report = True
solver.error_on_nonconvergence = False
solver.max_it = 1
info = solver.solve(w)
print(info)

with dolfinx.io.VTXWriter(domain.comm, "u.bp", [w.sub(0).collapse()], engine="BP4") as vtx:
    vtx.write(0.0)

image

3 Likes