Incompressible Navier Stokes with dolfinx_mpc unconverged

General
1

Hello I am solving Navier Stokes equation using dolfinx_mpc. The goal is to do it under periodic boundary condition for the top and bottom wall. A sticky obstacle is placed in the domain to reproduce Karman vortex street.

I first solve this problem under normal no_slip condition at the top, bottom and obstacle using averaged Euler integration. The inlet velocity and the outlet pressure are specified. However, when I do it using dolfinx, I can get convergent results over a period of time. When I switch to dolfinx_mpc, it doesnt converge anymore. Can you take a look to see if there is any mistake? Thanks!

Below is the code for dolfinx_mpc. In comparison, the code for dolfinx is at the end.

BTW the text editor keeps complaining I didnt use preformatted text for the code which I did used and to quote every code…

import basix.ufl
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.nls.petsc
import dolfinx.log
import dolfinx.plot
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex
import matplotlib.pyplot as plt
from petsc4py.PETSc import ScalarType
import pyvista
import sys
import scipy.sparse
import os
from pathlib import Path
import shutil
import dolfinx_mpc

def create_rectangle_with_a_hole(l_left, l_right, u_right, u_left, hole_center, hole_radius, lc):
    gmsh.initialize()
    gmsh.model.add("mesh")

    p0 = gmsh.model.geo.addPoint(l_left[0], l_left[1], 0, lc)
    p1 = gmsh.model.geo.addPoint(l_right[0], l_right[1], 0, lc)
    p2 = gmsh.model.geo.addPoint(u_right[0], u_right[1], 0, lc)
    p3 = gmsh.model.geo.addPoint(u_left[0], u_left[1], 0, lc)

    p0p1 = gmsh.model.geo.addLine(p0, p1)
    p1p2 = gmsh.model.geo.addLine(p1, p2)
    p2p3 = gmsh.model.geo.addLine(p2, p3)
    p3p0 = gmsh.model.geo.addLine(p3, p0)

    hole_c = gmsh.model.geo.addPoint(hole_center[0], hole_center[1], 0, hole_radius/5)
    hole_0 = gmsh.model.geo.addPoint(hole_center[0] + hole_radius, hole_center[1], 0, hole_radius/5)
    hole_1 = gmsh.model.geo.addPoint(hole_center[0], hole_center[1] + hole_radius, 0, hole_radius/5)
    hole_2 = gmsh.model.geo.addPoint(hole_center[0] - hole_radius, hole_center[1], 0, hole_radius/5)
    hole_3 = gmsh.model.geo.addPoint(hole_center[0], hole_center[1] - hole_radius, 0, hole_radius/5)

    hole_01 = gmsh.model.geo.addCircleArc(hole_0, hole_c, hole_1)
    hole_12 = gmsh.model.geo.addCircleArc(hole_1, hole_c, hole_2)
    hole_23 = gmsh.model.geo.addCircleArc(hole_2, hole_c, hole_3)
    hole_30 = gmsh.model.geo.addCircleArc(hole_3, hole_c, hole_0)

    hole = gmsh.model.geo.addCurveLoop([hole_01, hole_12, hole_23, hole_30])
    boundary = gmsh.model.geo.addCurveLoop([p0p1, p1p2, p2p3, p3p0])
    domain = gmsh.model.geo.addPlaneSurface([boundary, hole])

    gmsh.model.geo.synchronize()
    
    gmsh.model.addPhysicalGroup(2, [domain], 1)
    gmsh.model.addPhysicalGroup(1, [p3p0], 2)
    gmsh.model.addPhysicalGroup(1, [p1p2], 3)
    gmsh.model.addPhysicalGroup(1, [hole_01, hole_12, hole_23, hole_30], 4)
    gmsh.model.addPhysicalGroup(1, [p0p1], 5)
    gmsh.model.addPhysicalGroup(1, [p2p3], 6)

    gmsh.model.mesh.generate(2)

    mesh_data = dolfinx.io.gmsh.model_to_mesh(gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
    
    gmsh.finalize()

    # viskex.dolfinx.plot_mesh(mesh, line_width=4)
    # viskex.dolfinx.plot_mesh_tags(mesh, subdomains, 'subdomains', line_width=4)
    # viskex.dolfinx.plot_mesh_tags(mesh, boundaries, 'boundaries', line_width=6)

    return mesh_data

# l_left, l_right, u_right, u_left, hole_center, hole_radius, lc = [-15,-8], [25,-8], [25,8], [-15,8], [0.,0.], 0.5, 1.4
l_left, l_right, u_right, u_left, hole_center, hole_radius, lc = [0.,0.], [2.2,0], [2.2,0.41], [0.,0.41], [0.2,0.2], 0.05, 0.04
mesh_data = create_rectangle_with_a_hole(l_left, l_right, u_right, u_left, hole_center, hole_radius, lc)
mesh = mesh_data.mesh
boundaries = mesh_data.facet_tags
mdim = mesh.topology.dim
fdim = mdim - 1

t = 0.0
T = 8.0
num_steps = 12800
dt = T / num_steps
delta_t = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(dt))
nu = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0.001))
rho = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(1))
theta = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(.5))
from_top_to_bot = np.array(l_left) - np.array(u_left)

def periodic_relation_top_bottom(x):
    out_x = np.zeros(x.shape)
    out_x[0] = x[0] + from_top_to_bot[0]
    out_x[1] = x[1] + from_top_to_bot[1]
    out_x[2] = x[2]
    return out_x

Gamma_in, Gamma_out, Gamma_hole, Gamma_down, Gamma_up = 2, 3, 4, 5, 6

element_family = basix.ElementFamily.P
cellname = mesh.basix_cell()
dtype = dolfinx.default_real_type
Ve = basix.ufl.element(element_family, cellname, 2, shape=(mdim,), dtype=dtype)
Qe = basix.ufl.element(element_family, cellname, 1, dtype=dtype)

V = dolfinx.fem.functionspace(mesh, Ve)
Q = dolfinx.fem.functionspace(mesh, Qe)
W = ufl.MixedFunctionSpace(V, Q)

u_no_slip = dolfinx.fem.Function(V)
u_no_slip.x.array[:] = 0
dofs_no_slip = dolfinx.fem.locate_dofs_topological(V, fdim, np.hstack((np.hstack((boundaries.find(Gamma_up), boundaries.find(Gamma_down))), boundaries.find(Gamma_hole))))
# dofs_no_slip = dolfinx.fem.locate_dofs_topological(V, fdim, boundaries.find(Gamma_hole))
bc_hole = dolfinx.fem.dirichletbc(u_no_slip, dofs_no_slip)

u_in = dolfinx.fem.Function(V)
# u_in.interpolate(lambda x: [1. + x[1] * 0, 0 * x[0]])
u_in.interpolate(lambda x: [4. * 1 * x[1] * (u_left[1] - x[1]) / (u_left[1])**2, 0 * x[0]])
dofs_in = dolfinx.fem.locate_dofs_topological(V, fdim, boundaries.find(Gamma_in))
bc_in = dolfinx.fem.dirichletbc(u_in, dofs_in)

p_out = dolfinx.fem.Function(Q)
p_out.x.array[:] = 0
dofs_p_out = dolfinx.fem.locate_dofs_topological(Q, fdim, boundaries.find(Gamma_out))
bc_out = dolfinx.fem.dirichletbc(p_out, dofs_p_out)

bcs = [bc_hole, bc_in, bc_out]
bcs_u = [bc_hole, bc_in]
bcs_p = [bc_out]

mpc_V = dolfinx_mpc.MultiPointConstraint(V)
# mpc_V.create_periodic_constraint_topological(V, boundaries, Gamma_up, periodic_relation_top_bottom, bcs=bcs_u)
mpc_V.finalize()

mpc_Q = dolfinx_mpc.MultiPointConstraint(Q)
# mpc_Q.create_periodic_constraint_topological(Q, boundaries, Gamma_up, periodic_relation_top_bottom, bcs=bcs_p)
mpc_Q.finalize()

u, u_old = dolfinx.fem.Function(mpc_V.function_space), dolfinx.fem.Function(mpc_V.function_space)
p, p_old = dolfinx.fem.Function(mpc_Q.function_space), dolfinx.fem.Function(mpc_Q.function_space)
v, q = ufl.TestFunctions(W)

f = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0, 0)))
f_old = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0, 0)))

n = ufl.FacetNormal(mesh)

F = (ufl.dot(u - u_old, v) / delta_t 
     + theta * nu * ufl.inner(ufl.grad(u), ufl.grad(v)) + (1 - theta) * nu * ufl.inner(ufl.grad(u_old), ufl.grad(v))
     + theta * ufl.dot(ufl.dot(ufl.grad(u), u), v) + (1 - theta) * ufl.dot(ufl.dot(ufl.grad(u_old), u_old), v)
     - theta * ufl.inner(p, ufl.div(v)) - (1 - theta) * ufl.inner(p_old, ufl.div(v))
     - theta * ufl.inner(f, v) - (1 - theta) * ufl.inner(f_old, v)
     + ufl.inner(ufl.div(u), q)) * ufl.dx \
     + (theta * ufl.dot(p * n, v) + (1 - theta) * ufl.dot(p_old * n, v)) * (ufl.ds(Gamma_hole) + ufl.ds(Gamma_up) + ufl.ds(Gamma_down)) \
     + (theta * ufl.dot(p_out * n - nu * ufl.dot(ufl.grad(u), n), v) + (1 - theta) * ufl.dot(p_out * n - nu * ufl.dot(ufl.grad(u_old), n), v)) * ufl.ds(Gamma_out) \
     + (theta * ufl.dot(p * n - nu * ufl.dot(ufl.grad(u_in), n), v) + (1 - theta) * ufl.dot(p_old * n - nu * ufl.dot(ufl.grad(u_in), n), v)) * ufl.ds(Gamma_in)

tol = 1e-4
petsc_options = {
    "snes_type": "newtonls",
    "snes_atol": tol,
    "snes_rtol": tol,
    "snes_error_if_not_converged": True,
    "snes_linesearch_type": "none",
    "ksp_error_if_not_converged": True,
    "ksp_type": "minres",
    "ksp_rtol": 1e-4,
    "pc_type": "fieldsplit",
    "pc_fieldsplit_type": "additive",
    "monitorCancel": True
}
problem = dolfinx_mpc.NonlinearProblem(ufl.extract_blocks(F), [u, p], mpc=[mpc_V, mpc_Q], bcs=bcs, kind='nest', petsc_options=petsc_options)
nested_IS = problem.A.getNestISs()
ksp = problem.solver.getKSP()
pc = ksp.getPC()
pc.setFieldSplitIS(("u", nested_IS[0][0]), ("p", nested_IS[0][1]))
ksp_u, ksp_p = pc.getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")
# ksp.setMonitor(
#     lambda ctx, it, r: petsc4py.PETSc.Sys.Print(  # type: ignore
#         f"Iteration: {it:>4d}, |r| = {r:.3e}"
#     )
# )

export_path = Path('sol')
if export_path.exists():
    shutil.rmtree(export_path)
export_path.mkdir(exist_ok=True, parents=True)
vtx_u = dolfinx.io.VTXWriter(mesh.comm, export_path / "velocity.bp", [u], engine="BP4")
vtx_p = dolfinx.io.VTXWriter(mesh.comm, export_path / "pressure.bp", [p], engine="BP4")
vtx_u.write(t)
vtx_p.write(t)

t = 0.
step = 0
while t < T:
    print(t)
    _ = problem.solve()
    converged_reason = problem.solver.getConvergedReason()
    assert converged_reason > 0
    vtx_u.write(t)
    vtx_p.write(t)
    t += dt
    step += 1
    if step % 100 ==0:
        print(f"Converged = {converged_reason} in {problem.solver.getIterationNumber()} iterations at time {t}")
    u_old.x.array[:], p_old.x.array[:] = u.x.array, p.x.array
vtx_u.close()
vtx_p.close()

In dolfinx:

import basix.ufl
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.nls.petsc
import dolfinx.log
import dolfinx.plot
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex
import matplotlib.pyplot as plt
from petsc4py.PETSc import ScalarType
import pyvista
import sys
import scipy.sparse
import os
from pathlib import Path
import shutil
import copy

def create_rectangle_with_a_hole(l_left, l_right, u_right, u_left, hole_center, hole_radius, lc):
    gmsh.initialize()
    gmsh.model.add("mesh")

    p0 = gmsh.model.geo.addPoint(l_left[0], l_left[1], 0, lc)
    p1 = gmsh.model.geo.addPoint(l_right[0], l_right[1], 0, lc)
    p2 = gmsh.model.geo.addPoint(u_right[0], u_right[1], 0, lc)
    p3 = gmsh.model.geo.addPoint(u_left[0], u_left[1], 0, lc)

    p0p1 = gmsh.model.geo.addLine(p0, p1)
    p1p2 = gmsh.model.geo.addLine(p1, p2)
    p2p3 = gmsh.model.geo.addLine(p2, p3)
    p3p0 = gmsh.model.geo.addLine(p3, p0)

    hole_c = gmsh.model.geo.addPoint(hole_center[0], hole_center[1], 0, hole_radius/5)
    hole_0 = gmsh.model.geo.addPoint(hole_center[0] + hole_radius, hole_center[1], 0, hole_radius/5)
    hole_1 = gmsh.model.geo.addPoint(hole_center[0], hole_center[1] + hole_radius, 0, hole_radius/5)
    hole_2 = gmsh.model.geo.addPoint(hole_center[0] - hole_radius, hole_center[1], 0, hole_radius/5)
    hole_3 = gmsh.model.geo.addPoint(hole_center[0], hole_center[1] - hole_radius, 0, hole_radius/5)

    hole_01 = gmsh.model.geo.addCircleArc(hole_0, hole_c, hole_1)
    hole_12 = gmsh.model.geo.addCircleArc(hole_1, hole_c, hole_2)
    hole_23 = gmsh.model.geo.addCircleArc(hole_2, hole_c, hole_3)
    hole_30 = gmsh.model.geo.addCircleArc(hole_3, hole_c, hole_0)

    hole = gmsh.model.geo.addCurveLoop([hole_01, hole_12, hole_23, hole_30])
    boundary = gmsh.model.geo.addCurveLoop([p0p1, p1p2, p2p3, p3p0])
    domain = gmsh.model.geo.addPlaneSurface([boundary, hole])

    gmsh.model.geo.synchronize()
    
    gmsh.model.addPhysicalGroup(2, [domain], 1)
    gmsh.model.addPhysicalGroup(1, [p3p0], 2)
    gmsh.model.addPhysicalGroup(1, [p1p2], 3)
    gmsh.model.addPhysicalGroup(1, [p0p1, p2p3, hole_01, hole_12, hole_23, hole_30], 4)

    gmsh.model.mesh.generate(2)

    mesh_data = dolfinx.io.gmsh.model_to_mesh(gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
    
    gmsh.finalize()

    # viskex.dolfinx.plot_mesh(mesh, line_width=4)
    # viskex.dolfinx.plot_mesh_tags(mesh, subdomains, 'subdomains', line_width=4)
    # viskex.dolfinx.plot_mesh_tags(mesh, boundaries, 'boundaries', line_width=6)

    return mesh_data

l_left, l_right, u_right, u_left, hole_center, hole_radius, lc = [0.,0.], [2.2,0], [2.2,0.41], [0.,0.41], [0.2,0.2], 0.05, 0.04

mesh_data = create_rectangle_with_a_hole(l_left, l_right, u_right, u_left, hole_center, hole_radius, lc)
mesh = mesh_data.mesh
boundaries = mesh_data.facet_tags
mdim = mesh.topology.dim
fdim = mdim - 1

t = 0.0
T = 8.0
num_steps = 12800
dt = T / num_steps
delta_t = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(dt))
nu = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0.001))
U = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(1.))
rho = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(1))
theta = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(.5))
from_top_to_bot = np.array(l_left) - np.array(u_left)

Gamma_in, Gamma_out, Gamma_no_slip = 2, 3, 4

el_u = basix.ufl.element('Lagrange', mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))
el_p = basix.ufl.element('Lagrange', mesh.basix_cell(), 1)
el_mixed = basix.ufl.mixed_element([el_u, el_p])

W = dolfinx.fem.functionspace(mesh, el_mixed)
w = dolfinx.fem.Function(W)
u, p = ufl.split(w)
w_old = dolfinx.fem.Function(W)
u_old, p_old = ufl.split(w_old)
v, q = ufl.TestFunctions(W)

g = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0, 0)))
f = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0, 0)))
f_old = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0, 0)))

W0 = W.sub(0)
V, V_to_W0 = W0.collapse()

u_no_slip = dolfinx.fem.Function(V)
u_no_slip.x.array[:] = 0
dofs_no_slip = dolfinx.fem.locate_dofs_topological((W0, V), mesh.topology.dim-1, boundaries.find(Gamma_no_slip))
bc_no_slip = dolfinx.fem.dirichletbc(u_no_slip, dofs_no_slip, W0)

u_in = dolfinx.fem.Function(V)
u_in.interpolate(lambda x: [4 * U * x[1] * (u_left[1] - x[1]) / (u_left[1])**2, 0 * x[0]])
dofs_in = dolfinx.fem.locate_dofs_topological((W0, V), mesh.topology.dim-1, boundaries.find(Gamma_in))
bc_in = dolfinx.fem.dirichletbc(u_in, dofs_in, W0)

W1 = W.sub(1)
Q, Q_to_W1 = W1.collapse()

p_out = dolfinx.fem.Function(Q)
p_out.x.array[:] = 0
dofs_out = dolfinx.fem.locate_dofs_topological((W1, Q), mesh.topology.dim-1, boundaries.find(Gamma_out))
bc_out = dolfinx.fem.dirichletbc(p_out, dofs_out, W1)

bcs = [bc_no_slip, bc_in, bc_out]

n = ufl.FacetNormal(mesh)

F = (ufl.dot(u - u_old, v) / delta_t 
     + theta * nu * ufl.inner(ufl.grad(u), ufl.grad(v)) + (1 - theta) * nu * ufl.inner(ufl.grad(u_old), ufl.grad(v))
     + theta * ufl.dot(ufl.dot(ufl.grad(u), u), v) + (1 - theta) * ufl.dot(ufl.dot(ufl.grad(u_old), u_old), v)
     - theta * ufl.inner(p, ufl.div(v)) - (1 - theta) * ufl.inner(p_old, ufl.div(v))
     - theta * ufl.inner(f, v) - (1 - theta) * ufl.inner(f_old, v)
     + ufl.inner(ufl.div(u), q)) * ufl.dx \
     + (theta * ufl.dot(p * n - nu * ufl.dot(ufl.grad(u), n), v) + (1 - theta) * ufl.dot(p_old * n - nu * ufl.dot(ufl.grad(u_old), n), v)) * ufl.ds(Gamma_no_slip) \
     + (theta * ufl.dot(p_out * n - nu * ufl.dot(ufl.grad(u), n), v) + (1 - theta) * ufl.dot(p_out * n - nu * ufl.dot(ufl.grad(u_old), n), v)) * ufl.ds(Gamma_out) \
     + (theta * ufl.dot(p * n - nu * ufl.dot(ufl.grad(u), n), v) + (1 - theta) * ufl.dot(p_old * n - nu * ufl.dot(ufl.grad(u_old), n), v)) * ufl.ds(Gamma_in)

use_superlu = petsc4py.PETSc.IntType == np.int64  # or PETSc.ScalarType == np.complex64
sys = petsc4py.PETSc.Sys()  # type: ignore
if sys.hasExternalPackage("mumps") and not use_superlu:
    linear_solver = "mumps"
elif sys.hasExternalPackage("superlu_dist"):
    linear_solver = "superlu_dist"
else:
    linear_solver = "petsc"
petsc_options = {
    "snes_type": "newtonls",
    "snes_linesearch_type": "none",
    "snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps) * 1e-2,
    "snes_atol": 0,
    "snes_rtol": 0,
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": linear_solver,
    "snes_monitor_cancel": True,
    "ksp_monitor_cancel": True
}
problem = dolfinx.fem.petsc.NonlinearProblem(
    F, w, bcs=bcs, petsc_options_prefix="my_prefix", petsc_options=petsc_options
)

uh, ph = w.sub(0).collapse(), w.sub(1).collapse()
export_path = Path('sol_test')
if export_path.exists():
    shutil.rmtree(export_path)
export_path.mkdir(exist_ok=True, parents=True)
vtx_u = dolfinx.io.VTXWriter(mesh.comm, export_path / "velocity.bp", [uh], engine="BP4")
vtx_p = dolfinx.io.VTXWriter(mesh.comm, export_path / "pressure.bp", [ph], engine="BP4")
vtx_u.write(t)
vtx_p.write(t)

t = 0.
step = 0
while t < T:
    _ = problem.solve()
    converged_reason = problem.solver.getConvergedReason()
    num_iterations = problem.solver.getIterationNumber()
    if converged_reason <= 0:
        print([converged_reason, num_iterations, t])
    vtx_u.write(t)
    vtx_p.write(t)
    t += dt
    step += 1
    uh.x.array[:], ph.x.array[:] = w.sub(0).collapse().x.array, w.sub(1).collapse().x.array
    if step % 100 ==0:
        print(f"Converged = {converged_reason} in {num_iterations} iterations at time {t}")
    w_old.x.array[:] = w.x.array
vtx_u.close()
vtx_p.close()
2

I’ll start by providing a code that does run with periodic conditions and nest matrices, without using fieldsplit, instead using mumps:

import basix.ufl
import dolfinx.fem
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
from pathlib import Path
import shutil
import dolfinx_mpc


def create_rectangle_with_a_hole(
    l_left, l_right, u_right, u_left, hole_center, hole_radius, lc
):
    gmsh.initialize()
    gmsh.model.add("mesh")

    p0 = gmsh.model.geo.addPoint(l_left[0], l_left[1], 0, lc)
    p1 = gmsh.model.geo.addPoint(l_right[0], l_right[1], 0, lc)
    p2 = gmsh.model.geo.addPoint(u_right[0], u_right[1], 0, lc)
    p3 = gmsh.model.geo.addPoint(u_left[0], u_left[1], 0, lc)

    p0p1 = gmsh.model.geo.addLine(p0, p1)
    p1p2 = gmsh.model.geo.addLine(p1, p2)
    p2p3 = gmsh.model.geo.addLine(p2, p3)
    p3p0 = gmsh.model.geo.addLine(p3, p0)

    hole_c = gmsh.model.geo.addPoint(hole_center[0], hole_center[1], 0, hole_radius / 5)
    hole_0 = gmsh.model.geo.addPoint(
        hole_center[0] + hole_radius, hole_center[1], 0, hole_radius / 5
    )
    hole_1 = gmsh.model.geo.addPoint(
        hole_center[0], hole_center[1] + hole_radius, 0, hole_radius / 5
    )
    hole_2 = gmsh.model.geo.addPoint(
        hole_center[0] - hole_radius, hole_center[1], 0, hole_radius / 5
    )
    hole_3 = gmsh.model.geo.addPoint(
        hole_center[0], hole_center[1] - hole_radius, 0, hole_radius / 5
    )

    hole_01 = gmsh.model.geo.addCircleArc(hole_0, hole_c, hole_1)
    hole_12 = gmsh.model.geo.addCircleArc(hole_1, hole_c, hole_2)
    hole_23 = gmsh.model.geo.addCircleArc(hole_2, hole_c, hole_3)
    hole_30 = gmsh.model.geo.addCircleArc(hole_3, hole_c, hole_0)

    hole = gmsh.model.geo.addCurveLoop([hole_01, hole_12, hole_23, hole_30])
    boundary = gmsh.model.geo.addCurveLoop([p0p1, p1p2, p2p3, p3p0])
    domain = gmsh.model.geo.addPlaneSurface([boundary, hole])

    gmsh.model.geo.synchronize()

    gmsh.model.addPhysicalGroup(2, [domain], 1)
    gmsh.model.addPhysicalGroup(1, [p3p0], 2)
    gmsh.model.addPhysicalGroup(1, [p1p2], 3)
    gmsh.model.addPhysicalGroup(1, [hole_01, hole_12, hole_23, hole_30], 4)
    gmsh.model.addPhysicalGroup(1, [p0p1], 5)
    gmsh.model.addPhysicalGroup(1, [p2p3], 6)

    gmsh.model.mesh.generate(2)

    mesh_data = dolfinx.io.gmsh.model_to_mesh(
        gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2
    )

    gmsh.finalize()

    # viskex.dolfinx.plot_mesh(mesh, line_width=4)
    # viskex.dolfinx.plot_mesh_tags(mesh, subdomains, 'subdomains', line_width=4)
    # viskex.dolfinx.plot_mesh_tags(mesh, boundaries, 'boundaries', line_width=6)

    return mesh_data


# l_left, l_right, u_right, u_left, hole_center, hole_radius, lc = [-15,-8], [25,-8], [25,8], [-15,8], [0.,0.], 0.5, 1.4
l_left, l_right, u_right, u_left, hole_center, hole_radius, lc = (
    [0.0, 0.0],
    [2.2, 0],
    [2.2, 0.41],
    [0.0, 0.41],
    [0.2, 0.2],
    0.05,
    0.04,
)
mesh_data = create_rectangle_with_a_hole(
    l_left, l_right, u_right, u_left, hole_center, hole_radius, lc
)
mesh = mesh_data.mesh
boundaries = mesh_data.facet_tags
mdim = mesh.topology.dim
fdim = mdim - 1

t = 0.0
T = 8.0
num_steps = 12800
dt = T / num_steps
delta_t = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(dt))
nu = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0.001))
rho = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(1))
theta = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.5))
from_top_to_bot = np.array(l_left) - np.array(u_left)


def periodic_relation_top_bottom(x):
    out_x = np.zeros(x.shape)
    out_x[0] = x[0] + from_top_to_bot[0]
    out_x[1] = x[1] + from_top_to_bot[1]
    out_x[2] = x[2]
    return out_x


Gamma_in, Gamma_out, Gamma_hole, Gamma_down, Gamma_up = 2, 3, 4, 5, 6

element_family = basix.ElementFamily.P
cellname = mesh.basix_cell()
dtype = dolfinx.default_real_type
Ve = basix.ufl.element(element_family, cellname, 2, shape=(mdim,), dtype=dtype)
Qe = basix.ufl.element(element_family, cellname, 1, dtype=dtype)

V = dolfinx.fem.functionspace(mesh, Ve)
Q = dolfinx.fem.functionspace(mesh, Qe)
W = ufl.MixedFunctionSpace(V, Q)

u_no_slip = dolfinx.fem.Function(V)
u_no_slip.x.array[:] = 0
dofs_no_slip = dolfinx.fem.locate_dofs_topological(V, fdim, boundaries.find(Gamma_hole))
bc_hole = dolfinx.fem.dirichletbc(u_no_slip, dofs_no_slip)

u_in = dolfinx.fem.Function(V)
u_in.interpolate(lambda x: [1.0 + x[1] * 0.1, 10 * x[1]])
dofs_in = dolfinx.fem.locate_dofs_topological(V, fdim, boundaries.find(Gamma_in))
bc_in = dolfinx.fem.dirichletbc(u_in, dofs_in)

p_out = dolfinx.fem.Function(Q)
p_out.x.array[:] = 0
dofs_p_out = dolfinx.fem.locate_dofs_topological(Q, fdim, boundaries.find(Gamma_out))
bc_out = dolfinx.fem.dirichletbc(p_out, dofs_p_out)

bcs = [bc_hole, bc_in, bc_out]
bcs_u = [bc_hole, bc_in]
bcs_p = [bc_out]

mpc_V = dolfinx_mpc.MultiPointConstraint(V)
mpc_V.create_periodic_constraint_topological(
    V, boundaries, Gamma_up, periodic_relation_top_bottom, bcs=bcs_u
)
mpc_V.finalize()
mpc_Q = dolfinx_mpc.MultiPointConstraint(Q)
# mpc_Q.create_periodic_constraint_topological(
#     Q, boundaries, Gamma_up, periodic_relation_top_bottom, bcs=bcs_p
# )
mpc_Q.finalize()

u, u_old = (
    dolfinx.fem.Function(mpc_V.function_space),
    dolfinx.fem.Function(mpc_V.function_space),
)
p, p_old = (
    dolfinx.fem.Function(mpc_Q.function_space),
    dolfinx.fem.Function(mpc_Q.function_space),
)
v, q = ufl.TestFunctions(W)

f = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0, 0)))
f_old = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0, 0)))

n = ufl.FacetNormal(mesh)

F = (
    (
        ufl.dot(u - u_old, v) / delta_t
        + theta * nu * ufl.inner(ufl.grad(u), ufl.grad(v))
        + (1 - theta) * nu * ufl.inner(ufl.grad(u_old), ufl.grad(v))
        + theta * ufl.dot(ufl.dot(ufl.grad(u), u), v)
        + (1 - theta) * ufl.dot(ufl.dot(ufl.grad(u_old), u_old), v)
        - theta * ufl.inner(p, ufl.div(v))
        - (1 - theta) * ufl.inner(p_old, ufl.div(v))
        - theta * ufl.inner(f, v)
        - (1 - theta) * ufl.inner(f_old, v)
        + ufl.inner(ufl.div(u), q)
    )
    * ufl.dx
    + (theta * ufl.dot(p * n, v) + (1 - theta) * ufl.dot(p_old * n, v))
    * (ufl.ds(Gamma_hole) + ufl.ds(Gamma_up) + ufl.ds(Gamma_down))
    + (
        theta * ufl.dot(p_out * n - nu * ufl.dot(ufl.grad(u), n), v)
        + (1 - theta) * ufl.dot(p_out * n - nu * ufl.dot(ufl.grad(u_old), n), v)
    )
    * ufl.ds(Gamma_out)
    + (
        theta * ufl.dot(p * n - nu * ufl.dot(ufl.grad(u_in), n), v)
        + (1 - theta) * ufl.dot(p_old * n - nu * ufl.dot(ufl.grad(u_in), n), v)
    )
    * ufl.ds(Gamma_in)
)

tol = 1e-4
petsc_options = {
    "snes_monitor": None,
    "snes_type": "newtonls",
    "snes_atol": tol,
    "snes_rtol": tol,
    "snes_error_if_not_converged": True,
    "snes_linesearch_type": "none",
    "ksp_error_if_not_converged": True,
    "ksp_rtol": 1e-4,
    "snes_max_iterations": 250,
    "ksp_type": "preonly",
    "pc_type": "lu",
    "mat_mumps_icntl_24": 1,
    "mat_mumps_icntl_25": 0,
    "pc_mat_factor_type": "mumps",
    # "ksp_type": "minres",
    # "pc_type": "fieldsplit",
    # "pc_fieldsplit_type": "additive",
    # "monitorCancel": True,
}
problem = dolfinx_mpc.NonlinearProblem(
    ufl.extract_blocks(F),
    [u, p],
    mpc=[mpc_V, mpc_Q],
    bcs=bcs,
    kind="nest",
    petsc_options=petsc_options,
)
nested_IS = problem.A.getNestISs()
ksp = problem.solver.getKSP()
pc = ksp.getPC()
if pc.getType() == "fieldsplit":
    pc.setFieldSplitIS(("u", nested_IS[0][0]), ("p", nested_IS[0][1]))
    ksp_u, ksp_p = pc.getFieldSplitSubKSP()
    ksp_u.setType("preonly")
    ksp_u.getPC().setType("gamg")
    ksp_p.setType("preonly")
    ksp_p.getPC().setType("jacobi")
ksp.setMonitor(
    lambda ctx, it, r: petsc4py.PETSc.Sys.Print(  # type: ignore
        f"Iteration: {it:>4d}, |r| = {r:.3e}"
    )
)

export_path = Path("sol")
if export_path.exists():
    shutil.rmtree(export_path)
export_path.mkdir(exist_ok=True, parents=True)
vtx_u = dolfinx.io.VTXWriter(mesh.comm, export_path / "velocity.bp", [u], engine="BP4")
vtx_p = dolfinx.io.VTXWriter(mesh.comm, export_path / "pressure.bp", [p], engine="BP4")
vtx_u.write(t)
vtx_p.write(t)

t = 0.0
step = 0
while t < T:
    print(t)
    _ = problem.solve()
    converged_reason = problem.solver.getConvergedReason()
    assert converged_reason > 0
    vtx_u.write(t)
    vtx_p.write(t)
    t += dt
    step += 1
    if step % 100 == 0:
        print(
            f"Converged = {converged_reason} in {problem.solver.getIterationNumber()} iterations at time {t}"
        )
    u_old.x.array[:], p_old.x.array[:] = u.x.array, p.x.array

vtx_u.close()
vtx_p.close()

If you run this for a while (450 time steps), you get something like:

image
image719×428 39.2 KB

So the question becomes, why doesn’t it work with fieldsplit:

3

Though I dont understand the solver interface or the necessity of nest matrices, I have one question.

# mpc_Q.create_periodic_constraint_topological(
#     Q, boundaries, Gamma_up, periodic_relation_top_bottom, bcs=bcs_p
# )

This line was quoted. How to tell this line is not needed?

4

You don’t need to use nest to solve these problems. You can apply the MPC to a subspace and pass the form to nonoinearproblem ( dolfinx_mpc — DOLFINx-MPC: An extension to DOLFINx for multi point constraints)

Fieldsplit is used to efficiently solve blocked system, but relies on mathematical knowledge of what solver is efficient for each block.

Here you need to think about the physics. A pressure drop drives flow. If the pressure is periodic, there will be no drive.
This is mentioned in page 80 of: https://web.stanford.edu/class/me469b/handouts/incompressible.pdf

5

I talked to my mathematical wizard (Espen Sande) and he told me that since the pressure is in L^2 having a periodic condition in pressure doesn’t directly make mathematical sense without changing the variational formulation and corresponding function spaces.

Pressure boundary conditions are in general tricky, see for instance: Redirecting for a discussion regarding periodic conditions in a simple case.

6

This link is not working!!!

Now it works on campus network.

This one seems to be my unresolved question. I dont really know how to directly use the nonlinear solver from dolfinx_mpc. This code doesnt seem right.

use_superlu = petsc4py.PETSc.IntType == np.int64  # or PETSc.ScalarType == np.complex64
sys = petsc4py.PETSc.Sys()  # type: ignore
if sys.hasExternalPackage("mumps") and not use_superlu:
    linear_solver = "mumps"
elif sys.hasExternalPackage("superlu_dist"):
    linear_solver = "superlu_dist"
else:
    linear_solver = "petsc"
petsc_options = {
    "snes_type": "newtonls",
    "snes_linesearch_type": "none",
    "snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps) * 1e-2,
    "snes_atol": 0,
    "snes_rtol": 0,
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": linear_solver,
    "snes_monitor_cancel": True,
    "ksp_monitor_cancel": True
}
problem = dolfinx_mpc.NonlinearProblem(
    ufl.extract_blocks(F),
    [u, p],
    mpc=[mpc_V, mpc_Q],
    bcs=bcs,
    kind='None',
    petsc_options=petsc_options,
)
problem.solver.cancelMonitor()
7

I think a pressure gradient is from the left to the right overall. The pbc is imposed from top to bottom. btw I am in physics background, so most of the math is abstract for me.

First of all, I tested by imposing periodic condition in pressure space. The results appear to be the same, then nothing wrong with it and no effect with this imposition.

My understanding is that, in the weak form, the pressure can appear as itself acting on \nabla \cdot v through integration by part. No derivative of pressure needs to exist. Thus, the pressure doesnt need to be differentiable: L^2 space is sufficient. Mathematically, elements in this space doesnt need to have a well-defined boundary value for reasons I dont understand. Due to this fact, it is not meaningful to impose a pbc for pressure unless changing the variational form to let pressure live in other space where the boundary value is well-defined. In my problem set up, the pressure is fixed at one boundary. It doesnt make any difference whether to have the pbc or not for the pressure.

Another fact is that the pressure can differ by a constant and has no impact on the velocity field. This behaves like a gauge field. The easiest way to see is the in the strong form the pressure only appears in gradient \nabla p. Physically the driving form is the gradient as you mentioned.

8

Hello dokken could you take a look at this code why it’s not working by just passing the form to nonlinearproblem in dolfinx_mpc? Thank you!

use_superlu = petsc4py.PETSc.IntType == np.int64  # or PETSc.ScalarType == np.complex64
sys = petsc4py.PETSc.Sys()  # type: ignore
if sys.hasExternalPackage("mumps") and not use_superlu:
    linear_solver = "mumps"
elif sys.hasExternalPackage("superlu_dist"):
    linear_solver = "superlu_dist"
else:
    linear_solver = "petsc"
petsc_options = {
    "snes_type": "newtonls",
    "snes_linesearch_type": "none",
    "snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps) * 1e-2,
    "snes_atol": 0,
    "snes_rtol": 0,
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": linear_solver,
    "snes_monitor_cancel": True,
    "ksp_monitor_cancel": True
}
problem = dolfinx_mpc.NonlinearProblem(
    ufl.extract_blocks(F),
    [u, p],
    mpc=[mpc_V, mpc_Q],
    bcs=bcs,
    kind='None',
    petsc_options=petsc_options,
)
problem.solver.cancelMonitor()
9

I don’t know if you are still interested in the matter but you cannot L^2 “functions” do not have a determined value at the boundary (nor at any single point) because they are not functions but classes of functions.

In L^2 you are interested on functions as they act on other functions by integration.

That is, if you have two functions f and g that only differ at a single point (or at a set with zero measure), then when you compute

\int f(x)v(x)dx

you get exactly the same as if you compute

\int g(x)v(x)dx

for any v. So as they “act” the same, they belong to the same class of functions (they are the same element of L^2) and as such, you could not say what is the value of the functions of that class at a given point, it differs from function to function.

10

Thank you for your explanation. I think it explains

Can you also tell how to see that pressure in this problems belong to L^2 space?

Meanwhile I really want my question of how to implement this get answered by someone…

11

DOLFINx_mpc does not support kind=None.
You have to use nest matrices, and mumps if you want to use a direct solver with blocked problems.

It would have been useful to see the error traceback here, as I could have told you the answer earlier then (as I would have to inspect the code, but the error message).

12

Ok thank you so much! Since the nest matrices now work well, I will ask again when I solve some other problems and nest no longer works anymore.