dt = Constant(mesh, PETSc.ScalarType(0.1))
nu = Constant(mesh, PETSc.ScalarType(0.001/1))
T = 10
theta = 1
## Inlet velocity
U_0 = 0.5
# Inlet velocity profile
def u_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[
PETSc.ScalarType]:
"""Return the flat velocity profile at the inlet."""
values = np.zeros((2, x.shape[1]))
values[0, :] = 1.5
return values
def u_wall_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[
PETSc.ScalarType]:
"""Return the zero velocity at the wall."""
return np.zeros((2, x.shape[1]))
# Variational formulation
V_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim, ))
Q_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)
W_element = basix.ufl.mixed_element([V_element, Q_element])
W = dolfinx.fem.functionspace(mesh, W_element)
V, _ = W.sub(0).collapse()
Q, _ = W.sub(1).collapse()
fdim = mesh.topology.dim - 1
def Inflow(x):
# return np.abs(x[0] - 0.) < DOLFIN_EPS
return np.isclose(x[0], 0)
def Walls(x):
return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], H))
def Obstacle(x):
return np.isclose(np.sqrt(x[0]**2 + x[1]**2), r)
u_inlet = Function(V)
u_inlet.interpolate(u_in_eval)
inflow_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim,
Inflow)
inflow_dofs = dolfinx.fem.locate_dofs_topological((W.sub(0), V), fdim,
ft.indices[ft.values == 2])
bcu_inflow = dolfinx.fem.dirichletbc(u_inlet, inflow_dofs, W.sub(0))
u_nonslip = dolfinx.fem.Function(V)
u_nonslip.interpolate(u_wall_eval)
wall_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim,
Walls)
wall_dofs = dolfinx.fem.locate_dofs_topological((W.sub(0), V), fdim,
ft.indices[ft.values == 4])
bcu_walls = dolfinx.fem.dirichletbc(u_nonslip, wall_dofs, W.sub(0))
Obstacle_facets = dolfinx.mesh.locate_entities_boundary(mesh,fdim,
Obstacle)
Obstacle_dofs = dolfinx.fem.locate_dofs_topological((W.sub(0), V),fdim,
ft.indices[ft.values == 5])
bcu_obstacle = dolfinx.fem.dirichletbc(u_nonslip, Obstacle_dofs, W.sub(0))
# Boundary conditions
bcu = [bcu_inflow, bcu_walls, bcu_obstacle]
p_outlet = PETSc.ScalarType(0)
outflow_dofs = dolfinx.fem.locate_dofs_topological(Q, fdim,
ft.indices[ft.values == 3])
bcp_outlet = dirichletbc(p_outlet, outflow_dofs, W.sub(1))
bcp = [bcp_outlet]
up = Function(W)
u, p = split(up)
up_old = Function(W)
u_old, p_old = split(up_old)
(v, q) = TestFunctions(W)
F = ((1/dt)*dot(u - u_old, v)
+ (theta)*nu * inner(grad(u), grad(v))
+ (theta)*dot(dot(grad(u),u),v)
+ (1-theta)*nu * inner(grad(u_old), grad(v))
+ (1-theta)*dot(dot(grad(u_old),u_old),v)
- p * div(v)
- q * div(u)
) * dx
# Solver parameters
problem = NonlinearProblem(F, up, bcu)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.rtol = 1e-12 # Tolerancia relativa
solver.atol = 1e-15 # Tolerancia absoluta
solver.max_it = 100 # Iteraciones máximas permitidas
solver.convergence_criterion = "residual" # También puede ser "increment"
solver.report = True # Muestra en consola
solver.error_on_nonconvergence = False # No lanzar error si no converge
log.set_log_level(log.LogLevel.INFO)
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()=97
# Write the solution
xdmf_u = XDMFFile(mesh.comm, "velocity.xdmf", "w")
xdmf_u.write_mesh(mesh)
t = 0.0
up_old.x.array[:] = 0.0
ii = 0
while t < T:
ii += 1
t += float(dt)
r = solver.solve(up)
# Residual calculation
residual = fem.petsc.assemble_vector(fem.form(F))
fem.petsc.apply_lifting(residual, [fem.form(F)], bcs=[bc])
residual.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(residual, [bc])
up_old.x.array[:] = up.x.array
xdmf_u.close()
vtkfile_p.close()
I’ve tried in various ways but I can’t manage to compute the residual of the momentum equation of the Navier–Stokes system. Is there any alternative? I must be doing something wrong — I’ve been stuck for quite a while trying to get the residual of this equation, and the main issue is that when I apply the boundary condition, it doesn’t allow me to proceed. Thank you very much.