Problem in mixed-function space boundary conditions for stokes equation

General
,
1

Hello Everyone,

I am trying to understand the mixed function spaces implementation in the context of stokes problem. I am using Taylor-Hood elements. I have followed the Mixed formulation for the Poisson equation demo from the fenicsx documentation for setting up spaces, bcs etc. However, I am getting inf solution for velocity and pressure everywhere. It should atleast be equal to the applied boundary condition values at the Dirichlet dofs. I checked my identified degrees of freedom for velocity and pressure, which seems correct. Can someone point out what I might be missing?

from dolfinx import fem, mesh, io, geometry
import numpy as np
import ufl
from mpi4py import MPI
from basix.ufl import element, mixed_element
from petsc4py import PETSc

from dolfinx.fem.petsc import LinearProblem

# Create mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 1, 1)

# Taylor-Hood elements (P2-P1)
P2 = element("Lagrange", domain.basix_cell(), 2, shape=(2,))  # Vector element for velocity
P1 = element("Lagrange", domain.basix_cell(), 1)              # Scalar element for pressure
TH = mixed_element([P2, P1])

W = fem.functionspace(domain, TH)

# Get subspace of V
V0 = W.sub(0)
Q0 = W.sub(1)

V, V_to_W = W.sub(0).collapse()
Q, Q_to_W = W.sub(1).collapse()

# Define trial and test functions
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)  # Test functions

# Define manufactured solution
x = ufl.SpatialCoordinate(domain)
u_exact = ufl.as_vector([
    ufl.sin(2 * np.pi * x[0]) * ufl.cos(2 * np.pi * x[1]),
    -ufl.cos(2 * np.pi * x[0]) * ufl.sin(2 * np.pi * x[1])])
p_exact = ufl.sin(2 * np.pi * x[0]) * ufl.sin(2 * np.pi * x[1])

# Create a Function to hold the exact solution
u_exact_func = fem.Function(V)
u_exact_func.interpolate(lambda x: np.array([
    np.sin(2 * np.pi * x[0]) * np.cos(2 * np.pi * x[1]),
    -np.cos(2 * np.pi * x[0]) * np.sin(2 * np.pi * x[1])
], dtype=np.float64))

# Apply boundary conditions
facets = mesh.locate_entities_boundary(domain, 1, lambda x: np.full(x.shape[1], True))
dofs_V = fem.locate_dofs_topological((V0, V), 1, facets)
print("velocity essential dofs", dofs_V)
bc = fem.dirichletbc(u_exact_func, dofs_V, V0)

# Fix pressure at one point to remove null space
# Find a point (e.g., first vertex)
def is_corner(x, tol=1e-5):
    return np.logical_and(np.isclose(x[0], 0.0, atol=tol),
                          np.isclose(x[1], 0.0, atol=tol))
corner = mesh.locate_entities_boundary(domain, 0, is_corner)
pressure_dof = fem.locate_dofs_topological((Q0, Q), 0, corner)

print("pressure essential dofs", pressure_dof)

# Create pressure constraint
zero = fem.Function(Q)
zero.interpolate(lambda x: np.array([0.0], dtype=np.float64))
bc_p = fem.dirichletbc(zero, pressure_dof, Q0)

# Combine BCs
bcs = [bc, bc_p]

# Compute f = -νΔu + ∇p
ν = 1.0
f = -ν * ufl.div(ufl.grad(u_exact)) + ufl.grad(p_exact)

a = ν * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx - p * ufl.div(v) * ufl.dx - q * ufl.div(u) * ufl.dx
L = ufl.inner(f, v) * ufl.dx


# Solve
wh = fem.Function(W)
problem = fem.petsc.LinearProblem(a, L, bcs= bcs, 
                                petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
wh = problem.solve()
uh, ph = wh.sub(0), wh.sub(1)  # Split solution

print(uh.x.array)
print(ph.x.array)

## Verify solution
print("Velocity min/max:", np.min(uh.x.array), np.max(uh.x.array))
print("Pressure min/max:", np.min(ph.x.array), np.max(ph.x.array))

# Compute errors
error_u = fem.assemble_scalar(fem.form(ufl.inner(uh - u_exact, uh - u_exact) * ufl.dx)) ** 0.5
error_p = fem.assemble_scalar(fem.form((ph - p_exact)**2 * ufl.dx)) ** 0.5
print(f"Taylor-Hood errors: u_L2 = {error_u:.2e}, p_L2 = {error_p:.2e}")

The output looks like this.

velocity dofs [array([ 0,  1,  2,  3,  4,  5,  6,  7, 10, 11, 15, 16, 17, 18, 19, 20],
      dtype=int32), array([ 6,  7,  0,  1,  8,  9,  2,  3,  4,  5, 12, 13, 14, 15, 16, 17],
      dtype=int32)]
pressure dofs [array([12], dtype=int32), array([0], dtype=int32)]
[inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf]
[inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf]
Velocity min/max: inf inf
Pressure min/max: inf inf
Taylor-Hood errors: u_L2 = nan, p_L2 = inf
2

First, I would say: Stokes equations using Taylor-Hood elements — DOLFINx 0.9.0 documentation
would give you some hints.

The first is tuning the LU-solver.
I would use:

problem = fem.petsc.LinearProblem(
    a,
    L,
    bcs=bcs,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "ksp_error_if_not_converged": True,
        "pc_factor_mat_solver_type": "mumps",
        "mat_mumps_icntl_24": 1,
        "mat_mumps_icntl_25": 0,
    },
)

I would also increase the mesh-size from 1x1 to say 10x10.

Finally, for looking at the values of the solution, I would do:

wh = problem.solve()
uh, ph = wh.sub(0).collapse(), wh.sub(1).collapse()  # Split solution

## Verify solution
print("Velocity min/max:", np.min(uh.x.array), np.max(uh.x.array))
print("Pressure min/max:", np.min(ph.x.array), np.max(ph.x.array))

Complete code:

from dolfinx import fem, mesh, io, geometry
import numpy as np
import ufl
from mpi4py import MPI
from basix.ufl import element, mixed_element
from petsc4py import PETSc

from dolfinx.fem.petsc import LinearProblem

# Create mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

# Taylor-Hood elements (P2-P1)
P2 = element(
    "Lagrange", domain.basix_cell(), 2, shape=(2,)
)  # Vector element for velocity
P1 = element("Lagrange", domain.basix_cell(), 1)  # Scalar element for pressure
TH = mixed_element([P2, P1])

W = fem.functionspace(domain, TH)

# Get subspace of V
V0 = W.sub(0)
Q0 = W.sub(1)

V, V_to_W = W.sub(0).collapse()
Q, Q_to_W = W.sub(1).collapse()

# Define trial and test functions
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)  # Test functions

# Define manufactured solution
x = ufl.SpatialCoordinate(domain)
u_exact = ufl.as_vector(
    [
        ufl.sin(2 * np.pi * x[0]) * ufl.cos(2 * np.pi * x[1]),
        -ufl.cos(2 * np.pi * x[0]) * ufl.sin(2 * np.pi * x[1]),
    ]
)
p_exact = ufl.sin(2 * np.pi * x[0]) * ufl.sin(2 * np.pi * x[1])

# Create a Function to hold the exact solution
u_exact_func = fem.Function(V)
u_exact_func.interpolate(
    lambda x: np.array(
        [
            np.sin(2 * np.pi * x[0]) * np.cos(2 * np.pi * x[1]),
            -np.cos(2 * np.pi * x[0]) * np.sin(2 * np.pi * x[1]),
        ],
        dtype=np.float64,
    )
)

# Apply boundary conditions
facets = mesh.locate_entities_boundary(domain, 1, lambda x: np.full(x.shape[1], True))
dofs_V = fem.locate_dofs_topological((V0, V), 1, facets)
print("velocity essential dofs", dofs_V)
bc = fem.dirichletbc(u_exact_func, dofs_V, V0)


# Fix pressure at one point to remove null space
# Find a point (e.g., first vertex)
def is_corner(x, tol=1e-5):
    return np.logical_and(
        np.isclose(x[0], 0.0, atol=tol), np.isclose(x[1], 0.0, atol=tol)
    )


corner = mesh.locate_entities_boundary(domain, 0, is_corner)
pressure_dof = fem.locate_dofs_topological((Q0, Q), 0, corner)

print("pressure essential dofs", pressure_dof)

# Create pressure constraint
zero = fem.Function(Q)
zero.interpolate(lambda x: np.array([0.0], dtype=np.float64))
bc_p = fem.dirichletbc(zero, pressure_dof, Q0)

# Combine BCs
bcs = [bc, bc_p]

# Compute f = -νΔu + ∇p
ν = 1.0
f = -ν * ufl.div(ufl.grad(u_exact)) + ufl.grad(p_exact)

a = (
    ν * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
    - p * ufl.div(v) * ufl.dx
    - q * ufl.div(u) * ufl.dx
)
L = ufl.inner(f, v) * ufl.dx


# Solve
wh = fem.Function(W)
problem = fem.petsc.LinearProblem(
    a,
    L,
    bcs=bcs,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "ksp_error_if_not_converged": True,
        "pc_factor_mat_solver_type": "mumps",
        "mat_mumps_icntl_24": 1,
        "mat_mumps_icntl_25": 0,
    },
)
wh = problem.solve()
uh, ph = wh.sub(0).collapse(), wh.sub(1).collapse()  # Split solution

print(uh.x.array)
print(ph.x.array)

## Verify solution
print("Velocity min/max:", np.min(uh.x.array), np.max(uh.x.array))
print("Pressure min/max:", np.min(ph.x.array), np.max(ph.x.array))

# Compute errors
error_u = (
    fem.assemble_scalar(fem.form(ufl.inner(uh - u_exact, uh - u_exact) * ufl.dx)) ** 0.5
)
error_p = fem.assemble_scalar(fem.form((ph - p_exact) ** 2 * ufl.dx)) ** 0.5
print(f"Taylor-Hood errors: u_L2 = {error_u:.2e}, p_L2 = {error_p:.2e}")

import dolfinx

with dolfinx.io.VTXWriter(
    MPI.COMM_WORLD,
    "output.bp",
    [uh],
) as vtk:
    vtk.write(0.0)

with output

Velocity min/max: -1.0 1.0
Pressure min/max: -0.9241271481311483 1.020048274066652
Taylor-Hood errors: u_L2 = 3.15e-03, p_L2 = 5.05e-02
1 Like
3

Hi @dokken, thanks for your help. It works.

However, something really weird is happening. I ran the code once, and I get the same results as yours. When, I ran the same code again, the results are different and near inf. Does it imply faulty installation?

4

Are you now talking about the pressure, or both pressure and velocity?

5

Both pressure and velocity. The output is as follows:

Velocity min/max: -3.086480136725926e+154 8.394637962972945e+153
Pressure min/max: 5.262859580023891e+170 5.262859580023954e+170
Taylor-Hood errors: u_L2 = inf, p_L2 = inf

If I restart the computer and run the code, the solution is correct for the first time. If I rerun it, I am getting the above garbage solution.

6

Are you running with all the options:

        "ksp_type": "preonly",
        "pc_type": "lu",
        "ksp_error_if_not_converged": True,
        "pc_factor_mat_solver_type": "mumps",
        "mat_mumps_icntl_24": 1,
        "mat_mumps_icntl_25": 0,

from

as I would expect PETSc to complain if you hit these limits.

Edit:
I also occupationally see “mumps” struggle with this system.
What you could do is to remove the “pressure” bc, and set the constant pressure nullspace:

from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, mesh
import dolfinx.fem.petsc
import numpy as np
import ufl

from basix.ufl import element, mixed_element


# Create mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

# Taylor-Hood elements (P2-P1)
P2 = element(
    "Lagrange", domain.basix_cell(), 2, shape=(2,)
)  # Vector element for velocity
P1 = element("Lagrange", domain.basix_cell(), 1)  # Scalar element for pressure
TH = mixed_element([P2, P1])

W = fem.functionspace(domain, TH)

# Get subspace of V
V0 = W.sub(0)
Q0 = W.sub(1)

V, V_to_W = W.sub(0).collapse()
Q, Q_to_W = W.sub(1).collapse()

# Define trial and test functions
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)  # Test functions

# Define manufactured solution
x = ufl.SpatialCoordinate(domain)
u_exact = ufl.as_vector(
    [
        ufl.sin(2 * np.pi * x[0]) * ufl.cos(2 * np.pi * x[1]),
        -ufl.cos(2 * np.pi * x[0]) * ufl.sin(2 * np.pi * x[1]),
    ]
)
p_exact = ufl.sin(2 * np.pi * x[0]) * ufl.sin(2 * np.pi * x[1])

# Create a Function to hold the exact solution
u_exact_func = fem.Function(V)
u_exact_func.interpolate(
    lambda x: np.array(
        [
            np.sin(2 * np.pi * x[0]) * np.cos(2 * np.pi * x[1]),
            -np.cos(2 * np.pi * x[0]) * np.sin(2 * np.pi * x[1]),
        ],
        dtype=np.float64,
    )
)

# Apply boundary conditions
facets = mesh.locate_entities_boundary(domain, 1, lambda x: np.full(x.shape[1], True))
dofs_V = fem.locate_dofs_topological((V0, V), 1, facets)
print("velocity essential dofs", dofs_V)
bc = fem.dirichletbc(u_exact_func, dofs_V, V0)

# Combine BCs
bcs = [bc]

# Compute f = -νΔu + ∇p
ν = 1.0
f = -ν * ufl.div(ufl.grad(u_exact)) + ufl.grad(p_exact)

a = (
    ν * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
    - p * ufl.div(v) * ufl.dx
    - q * ufl.div(u) * ufl.dx
)
L = ufl.inner(f, v) * ufl.dx


# Solve
wh = fem.Function(W)
problem = fem.petsc.LinearProblem(
    a,
    L,
    bcs=bcs,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "ksp_error_if_not_converged": True,
        "pc_factor_mat_solver_type": "mumps",
        "mat_mumps_icntl_24": 1,
        "mat_mumps_icntl_25": 0,
    },
)
nsp = PETSc.NullSpace().create(constant=True)
problem.A.setNullSpace(nsp)
wh = problem.solve()
uh, ph = wh.sub(0).collapse(), wh.sub(1).collapse()  # Split solution

print(uh.x.array)
print(ph.x.array)

## Verify solution
print("Velocity min/max:", np.min(uh.x.array), np.max(uh.x.array))
print("Pressure min/max:", np.min(ph.x.array), np.max(ph.x.array))

# Compute errors
error_u = (
    fem.assemble_scalar(fem.form(ufl.inner(uh - u_exact, uh - u_exact) * ufl.dx)) ** 0.5
)
error_p = fem.assemble_scalar(fem.form((ph - p_exact) ** 2 * ufl.dx)) ** 0.5
print(f"Taylor-Hood errors: u_L2 = {error_u:.2e}, p_L2 = {error_p:.2e}")

import dolfinx

with dolfinx.io.VTXWriter(
    MPI.COMM_WORLD,
    "output.bp",
    [uh],
) as vtk:
    vtk.write(0.0)

Alternatively, change the options to;

  petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "ksp_error_if_not_converged": True,
        "pc_factor_mat_solver_type": "mumps",
    },

which gives me a constant solution (here you would need to normalize the pressure after the fact).

7

Setting the constant pressure nullspace works. However, I was avoiding it as I plan to extend this to NS equations, where I will use nonlinear solver. I am not sure how nullspace works with nonlinear problems, if it works at all.

Still, the unusual behavior of mumps is puzzling.

8

Inspecting the output with:

        "mat_mumps_icntl_4": 3,

I observe that the major difference between a working and failing simulation is:


 ****** SOLVE & CHECK STEP ********

 GLOBAL STATISTICS PRIOR SOLVE PHASE ...........
 Number of right-hand-sides                    =           1
 Blocking factor for multiple rhs              =           1
 ICNTL (9)                                     =           1
  --- (10)                                     =           0
  --- (11)                                     =           0
  --- (20)                                     =           0
  --- (21)                                     =           0
  --- (30)                                     =           0
  --- (35)                                     =           0
  --- (48) (effective)                         =           0


 Vector solution for column            1
 RHS
  -9.510565D-01  8.365346D-18 -2.449294D-16  2.707083D-17 -7.568733D-17
  -9.510565D-01 -1.981520D-16 -5.877853D-01  1.130382D+38 -1.651398D+38

and

Entering DMUMPS 5.7.3 from C interface with JOB, N, NNZ =   3         278           7234
      executing #MPI =      1, without OMP



 ****** SOLVE & CHECK STEP ********

 GLOBAL STATISTICS PRIOR SOLVE PHASE ...........
 Number of right-hand-sides                    =           1
 Blocking factor for multiple rhs              =           1
 ICNTL (9)                                     =           1
  --- (10)                                     =           0
  --- (11)                                     =           0
  --- (20)                                     =           0
  --- (21)                                     =           0
  --- (30)                                     =           0
  --- (35)                                     =           0
  --- (48) (effective)                         =           0


 Vector solution for column            1
 RHS
  -9.510565D-01  8.365346D-18 -2.449294D-16  2.707083D-17 -7.568733D-17
  -9.510565D-01 -1.981520D-16 -5.877853D-01 -4.669933D-01 -4.669933D-01
 ** Space in MBYTES used for solve                        :         0

which can be traced back to a difference in the RHS.
This is due to:

# Create pressure constraint
zero = fem.Function(Q)
zero.interpolate(lambda x: np.array([0.0], dtype=np.float64))
bc_p = fem.dirichletbc(zero, pressure_dof, Q0)

which is not correct.
It should be:

zero.interpolate(lambda x: np.zeros_like(x[1], dtype=np.float64))

I am suprised the code above actually doesn’t throw an error.

2 Likes
9

Thanks @dokken for your help.