Good morning,
I would like to solve the linearized version of the steady Navier-Stokes system without using Multiphenics or the implemented Newton solver because then I would like to move to a more involved formulation with the Level Set. I have adapted the tutorial on the Stokes problem demo_stokes. I don’t get any errors, but what I don’t understand is that all the vectors become zeros if I solve inside the loop, and nothing is updated. Do you have any suggestions on why this is happening? My code is the following:
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, fem, la
from dolfinx.fem import (
Constant,
Function,
dirichletbc,
extract_function_spaces,
form,
functionspace,
locate_dofs_topological,
)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner, sym, nabla_grad, dot
def analitical_velocity(x):
values = np.zeros((2, x.shape[1]))
values[0] = np.sin(np.pi*x[1])*np.cos(np.pi*x[1])*np.sin(np.pi*x[0])**2
values[1] = -np.sin(np.pi*x[0])*np.cos(np.pi*x[0])*np.sin(np.pi*x[1])**2
return values
def analitical_pressure(x):
return x[0]*x[1]*(1-x[0]-x[1])
def analitical_force(x, nu):
values = np.zeros((2, x.shape[1]))
values[0] = ((1-2*x[0])*x[1]*(1-x[1])
- nu*2*(np.pi**2) * ((np.cos(np.pi*x[0])**2)-3*np.sin(np.pi*x[0])**2)*np.sin(np.pi*x[1])*np.cos(np.pi*x[1])
+ np.pi*(np.sin(np.pi*x[0])**3) * (np.sin(np.pi*x[1])**2) * np.cos(np.pi*x[0]))
values[1] = ((1-2*x[1])*x[0]*(1-x[0])
+ nu*2*(np.pi**2) * ((np.cos(np.pi*x[1])**2)-3*np.sin(np.pi*x[1])**2)*np.sin(np.pi*x[0])*np.cos(np.pi*x[0])
+ np.pi*(np.sin(np.pi*x[1])**3) * (np.sin(np.pi*x[0])**2) * np.cos(np.pi*x[1]))
return values
# Function to mark x = 0, x = 1, y = 0 and y = 1
def noslip_boundary(x):
return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) | np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)
# Function to mark the lid (y = 1)
def lid(x):
return np.isclose(x[1], 1.0)
# Create mesh
msh = create_rectangle(
MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [10, 10], CellType.triangle
)
P2 = element(
"Lagrange", msh.basix_cell(), degree=2, shape=(msh.geometry.dim,), dtype=default_real_type
)
P1 = element("Lagrange", msh.basix_cell(), degree=1, dtype=default_real_type)
V, Q = functionspace(msh, P2), functionspace(msh, P1)
# No-slip condition on boundaries where x = 0, x = 1, and y = 0
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType) # type: ignore
u_fixed = Function(V)
u_fixed.interpolate(analitical_velocity)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc0 = dirichletbc(u_fixed, locate_dofs_topological(V, 1, facets))
# In case you want to try without null-space
p_fixed = Function(Q)
p_fixed.interpolate(analitical_pressure)
facets = locate_entities_boundary(msh, 1, lid)
bc1 = dirichletbc(p_fixed, locate_dofs_topological(Q, 1, facets))
# Collect Dirichlet boundary conditions
bcs = [bc0]
# Define variational problem
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
u_n = fem.Function(V)
u_n.x.petsc_vec.set(0)
nu = Constant(msh, PETSc.ScalarType(1e-1))
a = form([
[(nu*inner(sym(nabla_grad(u)), sym(nabla_grad(v)))
+ dot(dot(u_n, nabla_grad(u)), v)
+ dot(dot(u, nabla_grad(u_n)), v)) * dx, - inner(p, div(v)) * dx],
[inner(q, div(u)) * dx, Constant(msh, PETSc.ScalarType(0))*p*q * dx(999999)],
])
f = fem.Function(V)
f.interpolate(lambda x: analitical_force(x, nu))
L = form([(dot(dot(u_n, nabla_grad(u_n)), v)
+ inner(f, v)) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx])
a_p11 = form(inner(p, q) * dx)
a_p = [[a[0][0], None], [None, a_p11]]
def block_operators():
"""Return block operators and block RHS vector for the Stokes
problem"""
# Assembler matrix operator, preconditioner and RHS vector into
# single objects but preserving block structure
A = assemble_matrix_block(a, bcs=bcs)
A.assemble()
P = assemble_matrix_block(a_p, bcs=bcs)
P.assemble()
b = assemble_vector_block(L, a, bcs=bcs)
########## If you set the pressure somewhere you don't need the null space ##########
# Set the nullspace for pressure (since pressure is determined only
# up to a constant)
null_vec = A.createVecLeft()
offset = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
null_vec.array[offset:] = 1.0
null_vec.normalize()
nsp = PETSc.NullSpace().create(vectors=[null_vec])
assert nsp.test(A)
A.setNullSpace(nsp)
return A, P, b
def block_iterative_solver():
"""Solve the Navier-Stokes problem using blocked matrices and an iterative
solver."""
# Assembler the operators and RHS vector
A, P, b = block_operators()
# Build PETSc index sets for each field (global dof indices for each
# field)
V_map = V.dofmap.index_map
Q_map = Q.dofmap.index_map
offset_u = V_map.local_range[0] * V.dofmap.index_map_bs + Q_map.local_range[0]
offset_p = offset_u + V_map.size_local * V.dofmap.index_map_bs
is_u = PETSc.IS().createStride(
V_map.size_local * V.dofmap.index_map_bs, offset_u, 1, comm=PETSc.COMM_SELF
)
is_p = PETSc.IS().createStride(Q_map.size_local, offset_p, 1, comm=PETSc.COMM_SELF)
# Create a GMRES Krylov solver and a block-diagonal preconditioner
# using PETSc's additive fieldsplit preconditioner
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A, P)
ksp.setTolerances(rtol=1e-9)
ksp.setType("gmres")
ksp.getPC().setType("fieldsplit")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
ksp.getPC().setFieldSplitIS(("u", is_u), ("p", is_p))
# Configure velocity and pressure sub-solvers
ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.MULTIPLICATIVE) # Better for nonsymmetric systems
ksp_u.setType("preonly")
ksp_u.getPC().setType("hypre") # Use HYPRE BoomerAMG for nonsymmetric systems
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi") # Or "hypre" if pressure block is stiff
ksp.getPC().setUp()
Pu, _ = ksp_u.getPC().getOperators()
Pu.setBlockSize(msh.topology.dim)
# Create a block vector (x) to store the full solution and solve Stokes
# we have an initial guess
x = A.createVecRight()
ksp.solve(b, x)
# Create Functions to split u and p
u_h, p_h = Function(V), Function(Q)
offset = V_map.size_local * V.dofmap.index_map_bs
u_h.x.array[:offset] = x.array_r[:offset]
p_h.x.array[: (len(x.array_r) - offset)] = x.array_r[offset:]
print("u_h: ", u_h.x.array[:offset])
u_n.x.array[:] = u_h.x.array
print("u_n: ", u_n.x.array[:])
# Save solution to file in XDMF format for visualization, e.g. with
# ParaView. Before writing to file, ghost values are updated using
# `scatter_forward`.
with XDMFFile(MPI.COMM_WORLD, "out_stokes/velocity.xdmf", "w") as ufile_xdmf:
u_h.x.scatter_forward()
P1 = element(
"Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,), dtype=default_real_type
)
u1 = Function(functionspace(msh, P1))
u1.interpolate(u_h)
ufile_xdmf.write_mesh(msh)
ufile_xdmf.write_function(u1)
with XDMFFile(MPI.COMM_WORLD, "out_stokes/pressure.xdmf", "w") as pfile_xdmf:
p_h.x.scatter_forward()
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(p_h)
max_iterations = 2
i = 1
uD = fem.Function(V)
uD.interpolate(analitical_velocity)
L2_error = fem.form(dot(u_n - uD, u_n - uD) * dx)
while i < max_iterations:
A.zeroEntries()
with b.localForm() as b_loc:
b_loc.set(0)
A, P, b = block_operators()
x = A.createVecRight()
ksp.solve(b, x)
print("x: ", x[:])
u_h.x.array[:offset] = x.array_r[:offset]
print("u_h: ", u_h.x.array[:offset])
p_h.x.array[: (len(x.array_r) - offset)] = x.array_r[offset:]
u_n.x.array[:] = u_h.x.array
print("u_n: ", u_n.x.array[:])
i += 1
# Compute norm of update
L2_norm = np.sqrt(msh.comm.allreduce(fem.assemble_scalar(L2_error), op=MPI.SUM))
print(f"Iteration {i}: L2 error {L2_norm}")
if L2_norm < 1e-10:
break
u1.interpolate(u_n)
u1.x.scatter_forward()
ufile_xdmf.write_function(u1, i)
# Compute the $L^2$ norms of the solution vectors
norm_u, norm_p = la.norm(u_h.x), la.norm(p_h.x)
if MPI.COMM_WORLD.rank == 0:
print(f"(B) Norm of velocity coefficient vector (blocked, iterative): {norm_u}")
print(f"(B) Norm of pressure coefficient vector (blocked, iterative): {norm_p}")
return norm_u, norm_p
# Solve using PETSc block matrices and an iterative solver
norm_u_1, norm_p_1 = block_iterative_solver()