Previously solves the Navier-Stokes equations unsteady 3D [1]. Now, I need to solve the equations in steady state, I have implemented the following code, helping me with the following tutorials: [3D problem (sphere) Navier-Stokes], [Tutorial 02: Navier-Stokes problem], [Stokes equations using Taylor-Hood elements]
##NVSTK-SS-3D
import typing
import dolfinx
import numpy as np
from mpi4py import MPI
from dolfinx.cpp.mesh import CellType
from dolfinx.io import (XDMFFile)
from dolfinx.fem import (Constant, Function, )
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, set_bc, assemble_vector, create_vector)
from petsc4py import PETSc
import ufl
from ufl import (FacetNormal, Identity, Measure, TestFunction, TrialFunction,
as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
import basix
print(f"DOLFINx version: {dolfinx.__version__} is installed")
#-----------------------------------------------------------------------------------------------
#Mesh
filename="sphere-4.xdmf"
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "sphere-4.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
print(f"DOLFINx version: {dolfinx.__version__} is installed")
#------------------------------------------------------------------------------------------------------
#Parameters
eps=1e-14
nu=0.01
#------------------------------------------------------------------------------------------------------------
#Boundaries
def Inflow(x):
return np.abs(x[0] - 0.) < eps
def Outflow(x):
return np.abs(x[0] - 10.0) < eps
def Walls(x):
result = np.logical_or(np.abs(x[1] + 0.) < eps, np.abs(x[1] - 5.) < eps)
result2 = np.logical_or(np.abs(x[2] - 0.) < eps, np.abs(x[2] - 5.) < eps)
return np.logical_or(result, result2)
def Sphere(x):
result = np.logical_and(x[0] > 2.0 - eps, x[0] < 4.0 + eps)
result2 = np.logical_and(x[1] > 2.0 - eps, x[1] < 4.0 + eps)
result3 = np.logical_and(x[2] > 2.0 - eps, x[2] < 4.0 + eps)
return np.logical_and(np.logical_and(result, result2), result3)
def inflow_profile(x):
values = np.zeros((3, x.shape[1]), dtype=PETSc.ScalarType)
values[0, :] = 1.0
return values
print(f"DOLFINx version: {dolfinx.__version__} is installed")
#------------------------------------------------------------------------------------------------------------------
import dolfinx.fem as fem
V_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim, ))
Q_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)
#-------------------------------------------------------------------------------
def run_monolithic() -> dolfinx.fem.Function:
"""Run standard formulation using a mixed function space."""
# Function spaces
W_element = basix.ufl.mixed_element([V_element, Q_element])
W = dolfinx.fem.functionspace(mesh, W_element)
# Test and trial functions: monolithic
vq = ufl.TestFunction(W)
(v, q) = ufl.split(vq)
dup = ufl.TrialFunction(W)
up = dolfinx.fem.Function(W)
(u, p) = ufl.split(up)
print('After functions test and trial')
# Variational forms
F = (nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
+ ufl.inner(ufl.grad(u) * u, v) * ufl.dx
- ufl.inner(p, ufl.div(v)) * ufl.dx
+ ufl.inner(ufl.div(u), q) * ufl.dx)
J = ufl.derivative(F, up, dup)
print('After Variational forms')
# Boundary conditions
W0 = W.sub(0)
V, _ = W0.collapse()
u_in = fem.Function(V)
u_in.interpolate(inflow_profile)
u_zero = fem.Function(V)
u_zero.x.array[:] = 0
#p_zero = fem.Function(Q)
#p_zero.x.array[:] = 0
inflow_boundary_dofs = fem.locate_dofs_geometrical(V, Inflow)
#outflow_boundary_dofs = fem.locate_dofs_geometrical(Q, Outflow)
walls_boundary_dofs = fem.locate_dofs_geometrical(V, Walls)
sphere_boundary_dofs = fem.locate_dofs_geometrical(V, Sphere)
bcu_inflow = fem.dirichletbc(u_in, inflow_boundary_dofs)
#bcp_outflow = fem.dirichletbc(p_zero, outflow_boundary_dofs)
bcu_walls = fem.dirichletbc(u_zero, walls_boundary_dofs)
bcu_sphere = fem.dirichletbc(u_zero, sphere_boundary_dofs)
bc = [bcu_inflow, bcu_walls, bcu_sphere]
#bcp = [bcp_outflow]
print('After boundaries')
# Class for interfacing with SNES
class NavierStokesProblem:
"""Define a nonlinear problem, interfacing with SNES."""
def __init__( # type: ignore[no-any-unimported]
self, F: ufl.Form, J: ufl.Form, solution: dolfinx.fem.Function,
bcs: list[dolfinx.fem.DirichletBC], P: typing.Optional[ufl.Form] = None
) -> None:
self._F = dolfinx.fem.form(F)
self._J = dolfinx.fem.form(J)
self._obj_vec = dolfinx.fem.petsc.create_vector(self._F)
self._solution = solution
self._bcs = bcs
self._P = P
print('ending def init')
def create_snes_solution(self) ->PETSc.Vec: # type: ignore[no-any-unimported]
"""
Create a petsc4py.PETSc.Vec to be passed to petsc4py.PETSc.SNES.solve.
The returned vector will be initialized with the initial guess provided in `self._solution`.
"""
x = self._solution.x.petsc_vec.copy()
with x.localForm() as _x, self._solution.x.petsc_vec.localForm() as _solution:
_x[:] = _solution
print('ending def create_snes_solution')
return x
def update_solution(self, x: PETSc.Vec) -> None: # type: ignore[no-any-unimported]
"""Update `self._solution` with data in `x`."""
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
with x.localForm() as _x, self._solution.x.petsc_vec.localForm() as _solution:
_solution[:] = _x
print('ending def update_solution')
def obj( # type: ignore[no-any-unimported]
self, snes: PETSc.SNES, x: PETSc.Vec
) -> np.float64:
"""Compute the norm of the residual."""
self.F(snes, x, self._obj_vec)
print('ending obj')
return self._obj_vec.norm() # type: ignore[no-any-return]
def F( # type: ignore[no-any-unimported]
self, snes: PETSc.SNES, x: PETSc.Vec, F_vec: PETSc.Vec
) -> None:
"""Assemble the residual."""
self.update_solution(x)
with F_vec.localForm() as F_vec_local:
F_vec_local.set(0.0)
dolfinx.fem.petsc.assemble_vector(F_vec, self._F)
dolfinx.fem.petsc.apply_lifting(F_vec, [self._J], [self._bcs], x0=[x], scale=-1.0)
F_vec.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(F_vec, self._bcs, x, -1.0)
print('ending def F')
def J( # type: ignore[no-any-unimported]
self, snes: PETSc.SNES, x: PETSc.Vec, J_mat: PETSc.Mat,
P_mat: PETSc.Mat
) -> None:
"""Assemble the jacobian."""
J_mat.zeroEntries()
dolfinx.fem.petsc.assemble_matrix( # type: ignore[misc]
J_mat, self._J, self._bcs, diagonal=1.0) # type: ignore[arg-type]
J_mat.assemble()
if self._P is not None:
P_mat.zeroEntries()
dolfinx.fem.petsc.assemble_matrix( # type: ignore[misc]
P_mat, self._P, self._bcs, diagonal=1.0) # type: ignore[arg-type]
P_mat.assemble()
print('ending def J')
# Create problem
problem = NavierStokesProblem(F, J, up, bc)
F_vec = dolfinx.fem.petsc.create_vector(problem._F)
J_mat = dolfinx.fem.petsc.create_matrix(problem._J)
print(f"DOLFINx version: {dolfinx.__version__} is installed")
# Solve
snes = PETSc.SNES().create(mesh.comm)
snes.setTolerances(max_it=20)
snes.getKSP().setType("preonly")
snes.getKSP().getPC().setType("lu")
snes.getKSP().getPC().setFactorSolverType("mumps")
snes.setObjective(problem.obj)
snes.setFunction(problem.F, F_vec)
snes.setJacobian(problem.J, J=J_mat, P=None)
snes.setMonitor(lambda _, it, residual: print(it, residual))
up_copy = problem.create_snes_solution()
snes.solve(None, up_copy)
problem.update_solution(up_copy) # TODO can this be safely removed?
up_copy.destroy()
snes.destroy()
print(f"DOLFINx version: {dolfinx.__version__} is installed")
return up
up_m = run_monolithic()
(u_m, p_m) = (up_m.sub(0).collapse(), up_m.sub(1).collapse())
print(f"DOLFINx version: {dolfinx.__version__} is installed")
with XDMFFile(MPI.COMM_WORLD, "out_stokes/velocity.xdmf", "w") as ufile_xdmf:
u_m.x.scatter_forward()
P1 = element("Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,))
u1 = Function(functionspace(mesh, P1))
u1.interpolate(u_m)
ufile_xdmf.write_mesh(mesh)
ufile_xdmf.write_function(u1)
with XDMFFile(MPI.COMM_WORLD, "out_stokes/pressure.xdmf", "w") as pfile_xdmf:
p_m.x.scatter_forward()
pfile_xdmf.write_mesh(mesh)
pfile_xdmf.write_function(p_m)
However, I don’t get any results as seen in the images:
In the following link you can obtain all the files that were used for the simulation [Drive]
Could you please advise me on what went wrong and I can’t see any results. Please and thank you very much
