Hi,
I’m trying to simulate the transport of a quantity n flowing through a channel. I have successfully solved the problem in serial. When
I execute the same code in parallel, I obtained the same results as in a serial run for \mathrm{np} \leq 21, for \mathrm{np} \gt 21 the results for n present some strange behavior close to the inlet. Is that an issue with the DoFs ordering?
python3 stokes.py
mpirun -np 24 python3 stokes.py
stokes.py
import os
os.environ["OMP_NUM_THREADS"] = "1"
import datetime
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.io import XDMFFile
import dolfinx
from dolfinx.cpp.mesh import CellType
from dolfinx import (RectangleMesh, DirichletBC, Function, FunctionSpace, NewtonSolver, NonlinearProblem, Constant, Form, log)
from dolfinx.fem import (locate_dofs_topological, assemble_matrix, assemble_vector, set_bc, apply_lifting)
from dolfinx.mesh import locate_entities_boundary
from ufl import (Identity, VectorElement, FiniteElement, MixedElement, TestFunction, TrialFunction,
derivative, dx, grad, div, dot, inner, sym, split)
from math import sqrt
log.set_log_level(log.LogLevel.INFO)
# Parameters
U = 0.1
beta = 1e-3
delta = 1e-3
# Mesh
mesh = RectangleMesh(
MPI.COMM_WORLD,
[np.array([0, -0.5, 0]), np.array([10, 0.5, 0])], [200, 20],
CellType.quadrilateral, dolfinx.cpp.mesh.GhostMode.none)
def inlet_boundary(x):
return np.isclose(x[0], 0.0)
def outlet_boundary(x):
return np.isclose(x[0], 10.0)
def wall_boundary(x):
return np.logical_or(np.isclose(x[1], -0.5), np.isclose(x[1], 0.5))
# Define Function Spaces
V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
N = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([V, Q, N]))
# Boundary conditions
inlet_facets = locate_entities_boundary(mesh, 1, inlet_boundary)
outlet_facets = locate_entities_boundary(mesh, 1, outlet_boundary)
wall_facets = locate_entities_boundary(mesh, 1, wall_boundary)
def u_bc_inlet(x):
return np.stack((U*np.ones(x.shape[1]), np.zeros(x.shape[1])))
def u_bc_wall(x):
return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1])))
u_bc1 = Function(W.sub(0).collapse())
u_bc1.interpolate(u_bc_inlet)
dofs1 = locate_dofs_topological((W.sub(0), FunctionSpace(mesh, V)), 1, inlet_facets)
bc1 = DirichletBC(u_bc1, dofs1, W.sub(0))
u_bc2 = Function(W.sub(0).collapse())
u_bc2.vector.set(0.0)
dofs2 = locate_dofs_topological((W.sub(0).sub(1), FunctionSpace(mesh, V).sub(1)), 1, outlet_facets)
bc2 = DirichletBC(u_bc2, dofs2, W.sub(0).sub(1))
u_bc3 = Function(W.sub(0).collapse())
u_bc3.interpolate(u_bc_wall)
dofs3 = locate_dofs_topological((W.sub(0), FunctionSpace(mesh, V)), 1, wall_facets)
bc3 = DirichletBC(u_bc3, dofs3, W.sub(0))
n_bc = Function(W.sub(2).collapse())
n_bc.vector.set(1.0)
dofs4 = locate_dofs_topological((W.sub(2), FunctionSpace(mesh, N)), 1, inlet_facets)
bc4 = DirichletBC(n_bc, dofs4, W.sub(2))
bcs = [bc1, bc2, bc3, bc4]
# Define test and trial functions
w_trial = TrialFunction(W)
w_test = TestFunction(W)
v, q, c = split(w_test)
# Define function
w = Function(W)
u, p, n = split(w)
w0 = Function(W)
u0, p0, n0 = split(w0)
dw = Function(W)
# Initial conditions
def w_init(x):
values = np.zeros((4, x.shape[1]))
values[0] = 0.0
values[1] = 0.0
values[2] = 0.0
values[3] = 1.0
return values
w.interpolate(w_init)
# Strain rate tensor
def D(u):
return sym(grad(u))
# Identity tensor
I = Identity(len(u))
# Total stress
T = -p*I + 2*beta*D(u)
# Body force
f = Constant(mesh, (0, 0))
# Stokes - Equation
F = inner(grad(v), T)*dx - inner(v, f)*dx + q*div(u)*dx
F += c*dot(grad(n), u)*dx + delta*inner(grad(c), grad(n))*dx + c*Constant(mesh, 0.01)*inner(D(u),D(u))*n*dx
# Jacobian
J = derivative(F, w, w_trial)
class NonlinearEquation(NonlinearProblem):
def __init__(self, a, L, bcs):
super().__init__()
self.L = Form(L)
self.a = Form(a)
self.bcs = bcs
self._F = None
self._J = None
def form(self, x):
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
def F(self, x):
if self._F is None:
self._F = assemble_vector(self.L)
else:
with self._F.localForm() as f_local:
f_local.set(0.0)
self._F = assemble_vector(self._F, self.L)
apply_lifting(self._F, [self.a], [self.bcs], [x], -1.0)
self._F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(self._F, self.bcs, x, -1.0)
return self._F
def J(self, x):
if self._J is None:
self._J = assemble_matrix(self.a, self.bcs)
else:
self._J.zeroEntries()
self._J = assemble_matrix(self._J, self.a, self.bcs)
self._J.assemble()
return self._J
problem = NonlinearEquation(J, F, bcs)
solver = NewtonSolver(MPI.COMM_WORLD)
w.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
# Files
u_file = XDMFFile(MPI.COMM_WORLD, "stokes/velocity.xdmf", "w")
u_file.write_mesh(mesh)
p_file = XDMFFile(MPI.COMM_WORLD, "stokes/pressure.xdmf", "w")
p_file.write_mesh(mesh)
n_file = XDMFFile(MPI.COMM_WORLD, "stokes/n.xdmf", "w")
n_file.write_mesh(mesh)
# Solver
solver.solve(problem, w.vector)
w.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
u_file.write_function(w.split()[0])
p_file.write_function(w.split()[1])
n_file.write_function(w.split()[2])
u_file.close()
p_file.close()
n_file.close()