You have not set any solver options:
replacing this with
problem = dolfinx_mpc.LinearProblem(
a,
L,
mpc,
bcs=bcw,
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"ksp_monitor": None,
"ksp_error_if_not_converged": True,
},
)
yields a solution, but not what you would expect.
This is due to a bug in:
which you can replace with:
pc = MultiPointConstraint(W)
mpc.create_periodic_constraint_topological(
W.sub(0).sub(0), ft, 2, periodic_relation, bcu
)
mpc.create_periodic_constraint_topological(
W.sub(0).sub(1), ft, 2, periodic_relation, bcu
)
mpc.finalize()
i.e. constrain each velocity component individually, which as a complete code yields:
import ufl
import dolfinx
from mpi4py import MPI
from dolfinx import default_real_type, default_scalar_type
from dolfinx.fem import (
functionspace,
Function,
Constant,
dirichletbc,
locate_dofs_topological,
assemble_scalar,
form,
)
import numpy as np
import basix.ufl
from ufl import dot, inner, grad, div, dx, ds, nabla_grad, lhs, rhs
import dolfinx_mpc
from dolfinx_mpc import MultiPointConstraint
mesh_comm = MPI.COMM_WORLD
model_rank = 0
# create mesh
height = 1.0
width = 1.0
nx = 50
ny = 100
def meshing(lx, ly, Nx, Ny):
m = dolfinx.mesh.create_unit_square(mesh_comm, Nx, Ny)
x = m.geometry.x
# Refine on top and bottom sides
x[:, 1] = 0.5 * (np.arctan(2 * np.pi * (x[:, 1] - 0.5)) / np.arctan(np.pi) + 1)
# Scale
x[:, 0] = x[:, 0] * lx
x[:, 1] = x[:, 1] * ly
return m
mesh = meshing(width, height, nx, ny)
tdim = mesh.geometry.dim
fdim = tdim - 1
# Defining the domain markers (that will be used later for the BC):
facet_indices, facet_markers = [], []
mesh.topology.create_connectivity(fdim, tdim)
def walls_marker(x):
return np.isclose(x[1], 0.0) | np.isclose(x[1], height)
walls_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, walls_marker)
facet_indices.append(walls_facets)
facet_markers.append(np.full_like(walls_facets, 3)) # 3-->walls marker
def inl_marker(x):
return np.isclose(x[0], 0.0)
inlet_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, inl_marker)
facet_indices.append(inlet_facets)
facet_markers.append(np.full_like(inlet_facets, 1)) # 1-->inlet marker
def out_marker(x):
return np.isclose(x[0], width)
out_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, out_marker)
facet_indices.append(out_facets)
facet_markers.append(np.full_like(out_facets, 2)) # 2-->outlet marker
# Creating the meshargs object for easier BC imposition
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
ft = dolfinx.mesh.meshtags(
mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets]
) # Create a MeshTags object that associates data with a subset of mesh entities.
# ft.find(1) <-- inlet facets
# ft.find(2) <-- outlet facets
# ft.find(3) <-- walls facets
normal = ufl.FacetNormal(mesh)
# Construct function spaces
P2 = basix.ufl.element(
"Lagrange",
mesh.topology.cell_name(),
2,
shape=(mesh.geometry.dim,),
dtype=default_real_type,
)
P1 = basix.ufl.element(
"Lagrange", mesh.topology.cell_name(), 1, dtype=default_real_type
)
# 2. Crea l'elemento misto
Mel = basix.ufl.mixed_element([P2, P1])
# 3. Crea il FunctionSpace misto DOLFINx
W = functionspace(mesh, Mel) # Ora W è un dolfinx.fem.FunctionSpace
# Nota: Non ti servono V e Q separati qui sopra per creare W.
# Se ti servono per definire le funzioni (es. u_no_slip), li ottieni collassando i sottospazi di W:
V, V_to_W_map = W.sub(0).collapse()
Q, Q_to_W_map = W.sub(1).collapse()
# V = dolfinx.fem.functionspace(mesh, P2)
# Q = dolfinx.fem.functionspace(mesh, P1)
# W = MixedFunctionSpace(V, Q)
# Costruct the boundary conditions:
bcu = []
bcp = []
bcw = [] # bcu+bcp
def periodic_relation(x):
out_x = np.zeros(x.shape)
out_x[0] = x[0] - width
out_x[1] = x[1]
return out_x
# ---------------------- DIREHLET BC -----------------------#
# Velocity (ux, uy)
# Top or bottom walls
u_no_slip = Function(V)
u_no_slip.x.array[:] = 0
u_no_slip.x.scatter_forward() # <-- updates ghost values
# Pressure:
p_in = Constant(mesh, dolfinx.default_scalar_type(2.0))
p_out = Constant(mesh, dolfinx.default_scalar_type(0.0))
# ---------------------- DIRECT BC -----------------------#
# --- BC Velocità ---
u_no_slip = Function(V) # Definito su V (collassato)
u_no_slip.x.array[:] = 0.0
dofs_walls_u = locate_dofs_topological(
(W.sub(0), V), fdim, ft.find(3)
) # Nota la tupla (W.sub, V)
bc_walls_u = dirichletbc(u_no_slip, dofs_walls_u, W.sub(0)) # Applicato su W.sub(0)
bcu.append(bc_walls_u)
# Pressure(p)
dofs_in_p = locate_dofs_topological(W.sub(1), fdim, ft.find(1))
dofs_out_p = locate_dofs_topological(W.sub(1), fdim, ft.find(2))
bc_in_p = dirichletbc(p_in, dofs_in_p, W.sub(1))
bc_out_p = dirichletbc(p_out, dofs_out_p, W.sub(1))
bcp.append(bc_in_p)
bcp.append(bc_out_p)
bcw = np.concatenate((bcu, bcp))
with dolfinx.io.XDMFFile(mesh_comm, "mesh_mwe_10.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(ft, mesh.geometry)
mpc = MultiPointConstraint(W)
mpc.create_periodic_constraint_topological(
W.sub(0).sub(0), ft, 2, periodic_relation, bcu
)
mpc.create_periodic_constraint_topological(
W.sub(0).sub(1), ft, 2, periodic_relation, bcu
)
mpc.finalize()
# Initialize constants and expressions
Re = 10
nu = Constant(mesh, default_scalar_type(1 / Re))
# --- Set forcing ---
alpha = 0 # iniziale ( U_bulk=..... con alpha= -1)
force = Constant(mesh, default_scalar_type((alpha, 0)))
# Initialise the functions
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
u0 = Function(V)
u1 = Function(V)
p0 = Function(Q)
p1 = Function(Q)
FW = (
dot(dot(u0, nabla_grad(u)), v) * dx
+ (nu) * inner(nabla_grad(u), nabla_grad(v)) * dx
- div(v) * p * dx
- div(u) * q * dx
- dot(force, v) * dx
)
# Surface terms non considered:
# FW+= dot(p1*normal, v)*ds - dot((nu)*nabla_grad(u1)*normal, v)*ds
a = form(lhs(FW))
L = form(rhs(FW))
problem = dolfinx_mpc.LinearProblem(
a,
L,
mpc,
bcs=bcw,
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"ksp_monitor": None,
"ksp_error_if_not_converged": True,
},
)
wh = problem.solve()
u1 = wh.sub(0).collapse()
p1 = wh.sub(1).collapse()
with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u1]) as vtx:
vtx.write(0.0)
with dolfinx.io.VTXWriter(mesh.comm, "p.bp", [p1]) as vtx:
vtx.write(0.0)
which yields
