Fluid flow simulations, non-physical solutions

[p1rx]

Hi,

i have a troublesome problem with simulation of fluid flow for nontrivial geometries. I tried two approaches for solving transient NavierStokes eq: first IPCs + SUPG according with dolfinx tutorial (Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial), and second a direct approach with multiphenicsx blocked formulation and snesKSP solver. Both methods work fine for the examples from tutorials or for simple Poisseuile like channel flow, but when setting more complex geometry the simulations fail to produce any reasonable results. While Krylov subspace might not be the best choice for solver I cant manage to set any other solver without errors. Nether the less I thought that the stabilised incremental method from the tutorial would work without problem, but it reproduces wrong solutions as well.

I had tried multiple geometries, where i added obstacles or additional shunts to the pipe but I am getting rather incorrect results.

Bellow some images from tests

without the obstacle

with the obstacle

and the incremental scheme

As one of the tests i runned the geometry straight from dolfinx tutorial on my implementation in multiphenics but also got bizzare outcome

[p1rx]

dolfinx ver 0.7.1

Exemplary code:

import typing
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
from mpi4py import MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex
from dolfinx.mesh import locate_entities
from dolfinx.io import gmshio, XDMFFile, VTXWriter
from dolfinx import plot
import dolfinx
import numpy as np
from dolfinx import plot, default_scalar_type
import basix.ufl

import multiphenicsx
from petsc4py import PETSc 
from dolfinx import fem

import multiphenicsx.fem
import multiphenicsx.fem.petsc


g = 0.2
R = 1
gdim = 2
c_x, c_y = 1, 4
r = 0.5

# GMSH class
def create_mesh_with_gmsh(gdim=gdim):
    gmsh.initialize()
    gmsh.model.add("three_squares")

    # Define points for the squares and lines
    # Outer loop
    p1 = gmsh.model.occ.addPoint(0, 0, 0) 
    p2 = gmsh.model.occ.addPoint(1, 0, 0)
    p3 = gmsh.model.occ.addPoint(1, 3, 0)
    p4 = gmsh.model.occ.addPoint(3, 3, 0)
    p5 = gmsh.model.occ.addPoint(3, 3 + g, 0)
    p6 = gmsh.model.occ.addPoint(2 + np.sqrt(g*(2*R-g)), 3 + g, 0)
    p7 = gmsh.model.occ.addPoint(2, 3 + 2*R, 0) 
    p8 = gmsh.model.occ.addPoint(1, 3 + 2*R, 0) 
    p9 = gmsh.model.occ.addPoint(1, 7, 0) 
    p10 = gmsh.model.occ.addPoint(0, 7, 0) 
    # Inner loop
    p11 = gmsh.model.occ.addPoint(1, 3 + g, 0) 
    p12 = gmsh.model.occ.addPoint(2, 3 + g, 0) 
    p13 = gmsh.model.occ.addPoint(2, 3 + 2*R - g, 0) 
    p14 = gmsh.model.occ.addPoint(1, 3 + 2*R - g, 0) 
    
    c = gmsh.model.occ.addPoint(2, 3 + R, 0) 

    obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)

    # Create the squares
    l1 = gmsh.model.occ.addLine(p1, p2) # --> Inlet
    l2 = gmsh.model.occ.addLine(p2, p3) 
    l3 = gmsh.model.occ.addLine(p3, p4) 
    l4 = gmsh.model.occ.addLine(p4, p5) 
    l5 = gmsh.model.occ.addLine(p5, p6) 
    l6_arc1 = gmsh.model.occ.addCircleArc(p6, c, p7)
    l7 = gmsh.model.occ.addLine(p7, p8) 
    l8 = gmsh.model.occ.addLine(p8, p9) 
    l9 = gmsh.model.occ.addLine(p9, p10) # --> Outlet
    l10 = gmsh.model.occ.addLine(p10, p1) 
    

    loop1 = gmsh.model.occ.addCurveLoop([l1, l2, l3, l4, l5, l6_arc1, l7, l8, l9, l10])

    shape1 = gmsh.model.occ.addPlaneSurface([loop1])
    
    Omega_fluid = gmsh.model.occ.cut([(gdim, shape1)], [(gdim, obstacle)])
    gmsh.model.occ.synchronize()
    
    

    # Add physical groups
    
    fluid_marker = 1
    volumes = gmsh.model.getEntities(dim=gdim)
    assert (len(volumes) == 1)
    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
    gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")

    gmsh.model.addPhysicalGroup(1, [l1], 101) # --> Inlet
    gmsh.model.addPhysicalGroup(1, [l9], 102) # --> Outlet
    gmsh.model.addPhysicalGroup(1, [l2, l3, l4, l5, l6_arc1, l7, l8, l10], 201) # --> Walls
    # # gmsh.model.addPhysicalGroup(1, [l11, l12_arc2, l13, l14], 202) # --> InnerWalls

    obstacle = []
    boundaries = gmsh.model.getBoundary(volumes, oriented=False)
    for boundary in boundaries:
        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
        if np.allclose(center_of_mass, [c_x,  c_y, 0]):
            obstacle.append(boundary[1])
    obstacle_marker = 202
    gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
    gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")

    # Meshing options
    gmsh.option.setNumber("Mesh.ElementOrder", 2)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.2)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.2)
    
    gmsh.model.mesh.generate(2)

    # Convert Gmsh model to DOLFINx mesh
    gdim = 2
    gmsh_model_rank = 0
    mesh_comm = MPI.COMM_WORLD
    mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)
    mesh.topology.create_entities(mesh.topology.dim - 1)
    mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
    gmsh.finalize()

    return mesh, cell_tags, facet_tags

# Create and get the dolfinx mesh
mesh, subdomains, boundaries = create_mesh_with_gmsh()
cells_Omega1 = subdomains.indices[subdomains.values == 1]
dx = ufl.Measure("dx")(subdomain_data=subdomains)





# %% ############################################################
#  < Define function spaces > 

Vns_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim, ))
Qns_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)
Vns = dolfinx.fem.functionspace(mesh, Vns_element)
Qns = dolfinx.fem.functionspace(mesh, Qns_element)

(v, q) = (ufl.TestFunction(Vns), ufl.TestFunction(Qns))
(du, dp) = (ufl.TrialFunction(Vns), ufl.TrialFunction(Qns))
(u, p) = (dolfinx.fem.Function(Vns), dolfinx.fem.Function(Qns))
u_n = dolfinx.fem.Function(Vns)  # Previous time step velocity

# Define restrictions
dofs_Vns_Omega1 = dolfinx.fem.locate_dofs_topological(Vns, subdomains.dim, cells_Omega1)
restriction_Vns_Omega1 = multiphenicsx.fem.DofMapRestriction(Vns.dofmap, dofs_Vns_Omega1)
dofs_Qns_Omega1 = dolfinx.fem.locate_dofs_topological(Qns, subdomains.dim, cells_Omega1)
restriction_Qns_Omega1 = multiphenicsx.fem.DofMapRestriction(Qns.dofmap, dofs_Qns_Omega1)
restriction = [restriction_Vns_Omega1, restriction_Qns_Omega1] #, restriction_M_Gamma


# %% ############################################################
# Parameters
t = 0
T = 8                     # Final time
dt=0.1
num_steps = int(np.floor(T-t) / dt)
nu = 0.01


# Blocked form
F = [(1/dt) * ufl.inner(u - u_n, v) * dx(1) + nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx(1) +
        ufl.inner(ufl.grad(u) * u, v) * dx(1) - ufl.inner(p, ufl.div(v)) * dx(1),
        ufl.inner(ufl.div(u), q) * dx(1)]
# Jacobian
J = [[ufl.derivative(F[0], u, du), ufl.derivative(F[0], p, dp)],
    [ufl.derivative(F[1], u, du), ufl.derivative(F[1], p, dp)]]



# %% ################################################
# < Define boundary conditions > 
gdim=2
class InletVelocity():
    def __init__(self, t):
        self.t = t

    def __call__(self, x):
        values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
        values[1] = 25*x[0]*(1-x[0])/(1**2) #, x[1]-x[1]
        return values


# Inlet
u_inlet = dolfinx.fem.Function(Vns)
inlet_velocity = InletVelocity(t)
u_inlet.interpolate(inlet_velocity)
# u_inlet.interpolate(u_in_eval)
selected_tags_inlet = [101]
boundaries_inlet = boundaries.indices[np.isin(boundaries.values, selected_tags_inlet)]
bcu_inflow = dolfinx.fem.dirichletbc(u_inlet, dolfinx.fem.locate_dofs_topological(Vns, mesh.topology.dim - 1, boundaries_inlet))


# OuterWalls
u_nonslip = np.array((0,) * gdim, dtype=PETSc.ScalarType)
selected_tags_walls_outer = [201]
boundaries_walls_outer = boundaries.indices[np.isin(boundaries.values, selected_tags_walls_outer)]
bcu_walls_outer = dolfinx.fem.dirichletbc(u_nonslip, dolfinx.fem.locate_dofs_topological(Vns, mesh.topology.dim - 1, boundaries_walls_outer), Vns)

# InnerWalls
selected_tags_walls = [202]
boundaries_walls = boundaries.indices[np.isin(boundaries.values, selected_tags_walls)]
bcu_walls_inner = dolfinx.fem.dirichletbc(u_nonslip, dolfinx.fem.locate_dofs_topological(Vns, mesh.topology.dim - 1, boundaries_walls), Vns)


# # velocity outlet
selected_tags_pressureoutlet = [102]
boundaries_veloutlet = boundaries.indices[np.isin(boundaries.values, selected_tags_pressureoutlet)]
bdofs_V_out = dolfinx.fem.locate_dofs_topological(Vns, mesh.topology.dim - 1, boundaries_veloutlet)
bcu_outlet = dolfinx.fem.dirichletbc(u_inlet, bdofs_V_out)

bcu = [bcu_inflow, bcu_walls_outer, bcu_walls_inner, bcu_outlet]


# %% ############################################################
class NonlinearProblemV1V2(object):
    """Define a nonlinear problem, interfacing with SNES."""

    def __init__(  # type: ignore[no-any-unimported]
        self, F: typing.List[ufl.Form], J: typing.List[typing.List[ufl.Form]],
        solutions: typing.Tuple[dolfinx.fem.Function, dolfinx.fem.Function],
        bcs: typing.List[dolfinx.fem.DirichletBC],
        V1: dolfinx.fem.FunctionSpace, V2: dolfinx.fem.FunctionSpace,
        restriction: typing.List[multiphenicsx.fem.DofMapRestriction],
        P: typing.Optional[typing.List[typing.List[ufl.Form]]] = None
        
    ) -> None:
        self._F = dolfinx.fem.form(F)
        self._J = dolfinx.fem.form(J)
        self.V1 = V1
        self.V2 = V2
        self.restriction = restriction
        self._obj_vec = multiphenicsx.fem.petsc.create_vector_block(self._F, self.restriction)
        self._solutions = solutions
        self._bcs = bcs
        
        self._P = P
        
        

    def create_snes_solution(self) -> petsc4py.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 guesses provided in `self._solutions`,
        properly stacked together and restricted in a single block vector.
        """
        x = multiphenicsx.fem.petsc.create_vector_block(self._F, restriction=self.restriction)
        with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(x, [self.V1.dofmap, self.V2.dofmap], restriction) as x_wrapper:
            for x_wrapper_local, sub_solution in zip(x_wrapper, self._solutions):
                with sub_solution.vector.localForm() as sub_solution_local:
                    x_wrapper_local[:] = sub_solution_local
        return x

    def update_solutions(self, x: petsc4py.PETSc.Vec) -> None:  # type: ignore[no-any-unimported]
        """Update `self._solutions` with data in `x`."""
        x.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
        with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(x, [self.V1.dofmap, self.V2.dofmap], restriction) as x_wrapper:
            for x_wrapper_local, sub_solution in zip(x_wrapper, self._solutions):
                with sub_solution.vector.localForm() as sub_solution_local:
                    sub_solution_local[:] = x_wrapper_local

    def obj(  # type: ignore[no-any-unimported]
            self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec
    ) -> np.float64:
        """Compute the norm of the residual."""
        self.F(snes, x, self._obj_vec)
        return self._obj_vec.norm()  # type: ignore[no-any-return]

    def F(  # type: ignore[no-any-unimported]
        self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec, F_vec: petsc4py.PETSc.Vec
    ) -> None:
        """Assemble the residual."""
        self.update_solutions(x)
        with F_vec.localForm() as F_vec_local:
            F_vec_local.set(0.0)
        multiphenicsx.fem.petsc.assemble_vector_block(  # type: ignore[misc]
            F_vec, self._F, self._J, self._bcs, x0=x, scale=-1.0,
            restriction=self.restriction, restriction_x0=self.restriction)

    def J(  # type: ignore[no-any-unimported]
        self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec, J_mat: petsc4py.PETSc.Mat,
        P_mat: petsc4py.PETSc.Mat
    ) -> None:
        """Assemble the jacobian."""
        J_mat.zeroEntries()
        multiphenicsx.fem.petsc.assemble_matrix_block(
            J_mat, self._J, self._bcs, diagonal=1.0,  # type: ignore[arg-type]
            restriction=(self.restriction, self.restriction))
        J_mat.assemble()
        if self._P is not None:
            P_mat.zeroEntries()
            multiphenicsx.fem.petsc.assemble_matrix_block(
                P_mat, self._P, self._bcs, diagonal=1.0,  # type: ignore[arg-type]
                restriction=(self.restriction, self.restriction))
            P_mat.assemble()


# %% ############################################################
# Create a problem

# problem = NonlinearProblemV1V2M(F, J, (u1, u2, Rb), bcs)
problem = NonlinearProblemV1V2(F, J, (u, p), bcu, Vns, Qns, restriction)
F_vec = multiphenicsx.fem.petsc.create_vector_block(problem._F, restriction=restriction)
J_mat = multiphenicsx.fem.petsc.create_matrix_block(problem._J, restriction=(restriction, restriction))
snes = petsc4py.PETSc.SNES().create(mesh.comm)
snes.setTolerances(max_it=100)
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)
solution = problem.create_snes_solution()



# < Solve system >


from tqdm import tqdm
for i in tqdm(range(num_steps)):
    t += dt
    snes.solve(None, solution)
    u_n.x.array[:] = u.x.array[:]

# Finalize
solution.destroy()
snes.destroy()

[francesco-ballarin]

Please format the code properly: use “```” instead of the quote command.

[p1rx]

sry, done
as for the problem, I get that the NS tends to be numerically unstable, but can that really be for the presented probelms? It seems that the instabilities concentrate around sharp edges, that for slow flows produces nonphysical solutions and for high fluid flows explode the solution.

My intent for the model is to use in more engineering project in which fluid flow is one of multipple physics in some pipe network. Initially i wanted to implement RANS but it is too complex to code for just one project (adding the fact that the flow is just one of the physics)
I guess there is also an option to switch to Stokes as long as the flow magnitude is small.

What are your thoughts?

[hherlyng]

I just skimmed through the code, and if I interpret your code correctly, you’re setting an inlet velocity profile Dirichlet BC, noslip BCs on the walls, and additionally an outlet velocity profile Dirichlet BC. This seems like an over-constrained problem to me, as you’re forcing the nature of both the inlet and the outlet of the flow.

I would try with not setting a BC on the velocity at the outlet. One possibility is just dropping the Dirichlet BC you are settingat the outlet, this would be the same as having a homogeneous Neumann BC at the outlet. Another alternative is setting a pressure BC at the outlet (e.g. p=0).

By the way, it would be very helpful if the code you provide is as minimal as possible. In this case, you could e.g. just use the built-in rectangle mesh of FEniCSx so that you don’t have to include the whole gmsh code. Also, in my personal opinion, it would help a lot if you provided a mathematical statement of your problem, to make it clear what the variational form is and which boundary conditions you are setting on which parts of the geometry.

Cheers

[p1rx]

hey, Thx for the answer
you are right that there is inlet and outlet set as velocity dirichletbc. The inlet is at the bottom horizontal edge and outlet at the top. Rest of the edges are set to be the noslip walls. This is a residual of many tests. Initially i tested a free flux and 0 pressure outlet (the last did crash in the multiphenicsx implementation - probably i just dont understand how to set the bc for 2 lin spaces).

As for the mathematical formulation of the problem - it is a straight forward transient nonlinear Navier Stokes equation for incompressible newtonian fluid:

F = [(1/dt) * ufl.inner(u - u_n, v) * dx(1) + nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx(1) +
        ufl.inner(ufl.grad(u) * u, v) * dx(1) - ufl.inner(p, ufl.div(v)) * dx(1),
        ufl.inner(ufl.div(u), q) * dx(1)]

Also the presented implementation utilises multiphenicsx “restricition”. This example has just one domain so there is no need to use restrictions, but also this is a residual of tests in my project.

[p1rx]

as for the geometry, my main complain is that for simple rectangle domains the solutions seems ok - as long as the flow is stable. When there is some complexity added to the geometry the solution explodes.

Below i post the implementation of the incremental scheme from the tutorial link with only the geometry changed and inlet/outlet adjusted to the geometry (inlet is on the bottom edge and pressure outlet on the top

# < Imports > 
import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import tqdm.autonotebook

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                               create_vector, create_matrix, set_bc)
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from dolfinx.mesh import create_mesh, meshtags_from_entities

from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)


# %% ################################################
# < G M S H > 

# Geometry Parameters
g = 0.3
R = 1
gdim = 2
def create_mesh_with_gmsh(gdim=gdim):
    gmsh.initialize()
    gmsh.model.add("three_squares")

    # Define points for the squares and lines
    # Outer loop
    p1 = gmsh.model.occ.addPoint(0, 0, 0) 
    p2 = gmsh.model.occ.addPoint(1, 0, 0)
    p3 = gmsh.model.occ.addPoint(1, 3, 0)
    p4 = gmsh.model.occ.addPoint(3, 3, 0)
    p5 = gmsh.model.occ.addPoint(3, 3 + g, 0)
    p7 = gmsh.model.occ.addPoint(2, 3 + 2*R, 0) 
    p8 = gmsh.model.occ.addPoint(1, 3 + 2*R, 0) 
    p9 = gmsh.model.occ.addPoint(1, 7, 0) 
    p10 = gmsh.model.occ.addPoint(0, 7, 0) 
    # Inner loop
    p11 = gmsh.model.occ.addPoint(1, 3 + g, 0) 
    p12 = gmsh.model.occ.addPoint(2, 3 + g, 0) 
    p13 = gmsh.model.occ.addPoint(1.5, 3 + 2*R - g, 0) 
    p14 = gmsh.model.occ.addPoint(1, 3 + 2*R - g, 0) 
    
    c = gmsh.model.occ.addPoint(2, 3 + R, 0) 


    # Create the squares
    l1 = gmsh.model.occ.addLine(p1, p2) # --> Inlet
    l2 = gmsh.model.occ.addLine(p2, p3) 
    l3 = gmsh.model.occ.addLine(p3, p4) 
    l4 = gmsh.model.occ.addLine(p4, p5) 
    l6_arc1 =  gmsh.model.occ.addLine(p5, p7)
    l7 = gmsh.model.occ.addLine(p7, p8) 
    l8 = gmsh.model.occ.addLine(p8, p9) 
    l9 = gmsh.model.occ.addLine(p9, p10) # --> Outlet
    l10 = gmsh.model.occ.addLine(p10, p1) 
    l11 = gmsh.model.occ.addLine(p11, p12) 
    l12_arc2 = gmsh.model.occ.addLine(p13, p12) 
    l13 = gmsh.model.occ.addLine(p13, p14) 
    l14 = gmsh.model.occ.addLine(p14, p11) 
    
    loop1 = gmsh.model.occ.addCurveLoop([l1, l2, l3, l4, l6_arc1, l7, l8, l9, l10])
    loop2 = gmsh.model.occ.addCurveLoop([l11, l12_arc2, l13, l14])

    shape1 = gmsh.model.occ.addPlaneSurface([loop1])
    shape2 = gmsh.model.occ.addPlaneSurface([loop2])
    gmsh.model.occ.synchronize()
    Omega_fluid = gmsh.model.occ.cut([(gdim, shape1)], [(gdim, shape2)])
    gmsh.model.occ.synchronize()

    # Add physical groups
    gmsh.model.addPhysicalGroup(2, [Omega_fluid[0][0][1]], 1)

    gmsh.model.addPhysicalGroup(1, [l1], 101) # --> Inlet
    gmsh.model.addPhysicalGroup(1, [l9], 102) # --> Outlet
    gmsh.model.addPhysicalGroup(1, [l2, l3, l4,  l6_arc1, l7, l8, l10], 201) # --> Walls
    gmsh.model.addPhysicalGroup(1, [l11, l12_arc2, l13, l14], 202) # --> InnerWalls


    # Meshing options
    gmsh.option.setNumber("Mesh.Algorithm", 8)
    gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
    gmsh.option.setNumber("Mesh.RecombineAll", 1)
    gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.2)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.5)
    
    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.setOrder(2)
    gmsh.model.mesh.optimize("Netgen")

    # Convert Gmsh model to DOLFINx mesh
    gdim = 2
    gmsh_model_rank = 0
    mesh_comm = MPI.COMM_WORLD
    mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)
    mesh.topology.create_entities(mesh.topology.dim - 1)
    mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
    gmsh.finalize()

    return mesh, cell_tags, facet_tags

# Create and get the dolfinx mesh
mesh, subdomains, ft = create_mesh_with_gmsh()
ft.name = "Facet markers"

# < Parameters and lin spaces > 
t = 0
T = 5                       # Final time
dt = 1 / 1000                 # Time step size
num_steps = int(T / dt)
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(0.001))  # Dynamic viscosity
rho = Constant(mesh, PETSc.ScalarType(1))     # Density

v_cg2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

fdim = mesh.topology.dim - 1

u = TrialFunction(V)
v = TestFunction(V)
u_ = Function(V)
u_.name = "u"
u_s = Function(V)
u_n = Function(V)
u_n1 = Function(V)
p = TrialFunction(Q)
q = TestFunction(Q)
p_ = Function(Q)
p_.name = "p"
phi = Function(Q)

# %% ################################################
# < Define boundary conditions> 
gdim=2
class InletVelocity():
    def __init__(self, t):
        self.t = t

    def __call__(self, x):
        values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
        values[1] = np.sin(self.t * np.pi / 5) * 1*x[0]*(1-x[0])/(1**2) 
        return values


# Inlet
u_inlet = Function(V)
inlet_velocity = InletVelocity(t)
u_inlet.interpolate(inlet_velocity)
# u_inlet.interpolate(u_in_eval)
selected_tags_inlet = [101]
boundaries_inlet = ft.indices[np.isin(ft.values, selected_tags_inlet)]
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, mesh.topology.dim - 1, boundaries_inlet))

# OuterWalls
u_nonslip = np.array((0,) * gdim, dtype=PETSc.ScalarType)
selected_tags_walls_outer = [201]
boundaries_walls_outer = ft.indices[np.isin(ft.values, selected_tags_walls_outer)]
bcu_walls_outer = dirichletbc(u_nonslip, locate_dofs_topological(V, mesh.topology.dim - 1, boundaries_walls_outer), V)

# InnerWalls
selected_tags_walls = [202]
boundaries_walls = ft.indices[np.isin(ft.values, selected_tags_walls)]
bcu_walls_inner = dirichletbc(u_nonslip, locate_dofs_topological(V, mesh.topology.dim - 1, boundaries_walls), V)

bcu = [bcu_inflow, bcu_walls_outer, bcu_walls_inner] 

# pressure Outlet
selected_tags_pressureoutlet = [102]
ft_pressureoutlet = ft.indices[np.isin(ft.values, selected_tags_pressureoutlet)]
bcp_outlet = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(Q, mesh.topology.dim - 1, ft_pressureoutlet), Q)

bcp = [bcp_outlet] 

# %% ################################################
# < V A R I A T I O N A L   F O R M >
# initial guess for u
f = Constant(mesh, PETSc.ScalarType((0, 0)))
F1 = rho / k * dot(u - u_n, v) * dx
F1 += inner(dot(1.5 * u_n - 0.5 * u_n1, 0.5 * nabla_grad(u + u_n)), v) * dx
F1 += 0.5 * mu * inner(grad(u + u_n), grad(v)) * dx - dot(p_, div(v)) * dx
F1 += dot(f, v) * dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
A1 = create_matrix(a1)
b1 = create_vector(L1) 

# pressure correction
a2 = form(dot(grad(p), grad(q)) * dx)
L2 = form(-rho / k * dot(div(u_s), q) * dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)

# u corection
a3 = form(rho * dot(u, v) * dx)
L3 = form(rho * dot(u_s, v) * dx - k * dot(nabla_grad(phi), v) * dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)

# %% ################################################
# < Set S O L V E R S >

# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.JACOBI)

# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.MINRES)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")

# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)

# %% ################################################
# < S O L V E >
t = 0
progress = tqdm.autonotebook.tqdm(desc="Solving PDE", total=num_steps)
for i in range(num_steps):
    progress.update(1)
    # Update current time step
    t += dt
    # Update inlet velocity
    inlet_velocity.t = t
    u_inlet.interpolate(inlet_velocity)

    # Step 1: Tentative velocity step
    A1.zeroEntries()
    assemble_matrix(A1, a1, bcs=bcu)
    A1.assemble()
    with b1.localForm() as loc:
        loc.set(0)
    assemble_vector(b1, L1)
    apply_lifting(b1, [a1], [bcu])
    b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b1, bcu)
    solver1.solve(b1, u_s.vector)
    u_s.x.scatter_forward()

    # Step 2: Pressure corrrection step
    with b2.localForm() as loc:
        loc.set(0)
    assemble_vector(b2, L2)
    apply_lifting(b2, [a2], [bcp])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b2, bcp)
    solver2.solve(b2, phi.vector)
    phi.x.scatter_forward()

    p_.vector.axpy(1, phi.vector)
    p_.x.scatter_forward()

    # Step 3: Velocity correction step
    with b3.localForm() as loc:
        loc.set(0)
    assemble_vector(b3, L3)
    b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    solver3.solve(b3, u_.vector)
    u_.x.scatter_forward()


    # Update variable with solution form this time step
    with u_.vector.localForm() as loc_, u_n.vector.localForm() as loc_n, u_n1.vector.localForm() as loc_n1:
        loc_n.copy(loc_n1)
        loc_.copy(loc_

n)

Outcome to the presented code is:

To bre frank, this is probably the numerical instability of the method, although i faced it faster than i expected

[dokken]

@p1rx

Initially i tested a free flux and 0 pressure outlet (the last did crash in the multiphenicsx implementation - probably i just dont understand how to set the bc for 2 lin spaces).

You should not prescribe the outlet velocity, as you now force your scheme to have a particular flow profile (as @hherlyng says). You should rather pursue having an open boundary (outlet) and figure out why that didn’t work in your case.

[p1rx]

i have retested the simulations using Oasisx package to rule out mistakes on my part in formulations (i have used two to this point). Same as previously, while simple geometry works fine the non-trivial geometry invites errors to the solution. Here is a solution to geometry with looping pipe network, where the velocity inlet (constant ux=1, uy=0) is at the bottom edge and pressure outlet at the bottom boundary:

here is simplified code for said example:

from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
from dolfinx import fem, mesh
from dolfinx.io import gmshio, XDMFFile, VTXWriter
import oasisx
from oasisx import DirichletBC, LocatorMethod, FractionalStep_AB_CN
import ufl
import gmsh
from tqdm import tqdm
from typing import List
import numpy as np

# < G M S H > 
gmsh.initialize()
g = 0.3
R = 1
gdim = 2

gmsh.initialize()

p1 = gmsh.model.occ.addPoint(0, 0, 0) 
p2 = gmsh.model.occ.addPoint(1, 0, 0)
p3 = gmsh.model.occ.addPoint(1, 3, 0)
p4 = gmsh.model.occ.addPoint(3, 3, 0)
p5 = gmsh.model.occ.addPoint(3, 3 + g, 0)
p7 = gmsh.model.occ.addPoint(2, 3 + 2*R, 0) 
p8 = gmsh.model.occ.addPoint(1, 3 + 2*R, 0) 
p9 = gmsh.model.occ.addPoint(1, 7, 0) 
p10 = gmsh.model.occ.addPoint(0, 7, 0) 
# Inner loop
p11 = gmsh.model.occ.addPoint(1, 3 + g, 0) 
p12 = gmsh.model.occ.addPoint(2, 3 + g, 0) 
p13 = gmsh.model.occ.addPoint(1.5, 3 + 2*R - g, 0) 
p14 = gmsh.model.occ.addPoint(1, 3 + 2*R - g, 0) 


l1 = gmsh.model.occ.addLine(p1, p2) # --> Inlet
l2 = gmsh.model.occ.addLine(p2, p3) 
l3 = gmsh.model.occ.addLine(p3, p4) 
l4 = gmsh.model.occ.addLine(p4, p5) 
l6_arc1 =  gmsh.model.occ.addLine(p5, p7)
l7 = gmsh.model.occ.addLine(p7, p8) 
l8 = gmsh.model.occ.addLine(p8, p9) 
l9 = gmsh.model.occ.addLine(p9, p10) # --> Outlet
l10 = gmsh.model.occ.addLine(p10, p1) 
l11 = gmsh.model.occ.addLine(p11, p12) 
l12_arc2 = gmsh.model.occ.addLine(p13, p12) 
l13 = gmsh.model.occ.addLine(p13, p14) 
l14 = gmsh.model.occ.addLine(p14, p11) 
loop1 = gmsh.model.occ.addCurveLoop([l1, l2, l3, l4, l6_arc1, l7, l8, l9, l10])
loop2 = gmsh.model.occ.addCurveLoop([l11, l12_arc2, l13, l14])
shape1 = gmsh.model.occ.addPlaneSurface([loop1])
shape2 = gmsh.model.occ.addPlaneSurface([loop2])
gmsh.model.occ.synchronize()

Omega_fluid = gmsh.model.occ.cut([(gdim, shape1)], [(gdim, shape2)])
gmsh.model.occ.synchronize()

# Add physical groups
gmsh.model.addPhysicalGroup(2, [Omega_fluid[0][0][1]], 1)
inlet_marker, outlet_marker, wall_marker, wall2_marker = 101, 102, 201, 202
gmsh.model.addPhysicalGroup(1, [l1], inlet_marker) # --> Inlet
gmsh.model.addPhysicalGroup(1, [l9], outlet_marker) # --> Outlet
gmsh.model.addPhysicalGroup(1, [l2, l3, l4,  l6_arc1, l7, l8, l10], wall_marker) # --> Walls
gmsh.model.addPhysicalGroup(1, [l11, l12_arc2, l13, l14], wall2_marker) # --> InnerWalls


# Meshing options
gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.2)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.5)

gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")

# Convert Gmsh model to DOLFINx mesh
gdim = 2
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
model_rank = 0
domain, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
domain.topology.create_entities(domain.topology.dim - 1)
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim - 1)
gmsh.finalize()
facet_tags.name = "Facet markers"


# < BCu >
bcs_p: List[oasisx.PressureBC] = []
bc_inlet_x = DirichletBC(dolfinx.fem.Constant(domain, 0.),
                          method=LocatorMethod.TOPOLOGICAL, marker = (facet_tags, inlet_marker)) # method= marker = 
bc_inlet_y = DirichletBC(dolfinx.fem.Constant(domain, 1.),
                         method=LocatorMethod.TOPOLOGICAL, marker = (facet_tags, inlet_marker))
bc_wall = DirichletBC(dolfinx.fem.Constant(domain, 0.), method=LocatorMethod.TOPOLOGICAL,
                      marker=(facet_tags, wall_marker))
bc_wall2 = DirichletBC(dolfinx.fem.Constant(domain, 0.), method=LocatorMethod.TOPOLOGICAL,
                      marker=(facet_tags, wall2_marker))
bcs_u = [[bc_inlet_x, bc_wall, bc_wall2], [bc_inlet_y, bc_wall, bc_wall2]]

# < BCp >
bcs_p: List[oasisx.PressureBC] = [oasisx.PressureBC(
    0.,  marker=(facet_tags, outlet_marker))]


# fractional step solver
solver = FractionalStep_AB_CN(
    mesh=domain,
    u_element=("Lagrange", 2),
    p_element=("Lagrange", 1),
    bcs_u=bcs_u,
    bcs_p=bcs_p,
    solver_options={"tentative": {"ksp_type": "preonly", "pc_type": "lu"}, "pressure": {
        "ksp_type": "preonly", "pc_type": "lu"}, "scalar": {"ksp_type": "preonly", "pc_type": "lu"}},
    body_force=None
)

# Time-stepping
T_start, T_end, dt = 0.0, 1, 0.001
num_steps = int((T_end - T_start) // dt)
with dolfinx.io.VTXWriter(domain.comm, "_OASISX_Test1-PipeNetwork_v1.bp", [solver.u], engine="BP4") as writer:
    for step in tqdm(range(num_steps)):
        solver.solve(dt, nu=0.01, max_iter=100)
        writer.write(dt*step)

here is an output from Ansys Fluent (yes i know that the Fluent uses Rans + SIMPLE + k/eps omega whichc is different from formulations of the discussed program implemented in dolfinx. Take it more as a reference of how the solution should look like)

As final note i would like to say that this problem is not prevlent only for this geometry. More complex geometries tend to have even more unphysical solutions, that usually forms around sharp edges

[dokken]

You have marked you boundaries wrong. You have set the outlet as the top wall in the right-channel. CHanging this yields

from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
from dolfinx import fem, mesh
from dolfinx.io import gmshio, XDMFFile, VTXWriter
import oasisx
from oasisx import DirichletBC, LocatorMethod, FractionalStep_AB_CN
import ufl
import gmsh
from tqdm import tqdm
from typing import List
import numpy as np

# < G M S H >
gmsh.initialize()
g = 0.3
R = 1
gdim = 2

gmsh.initialize()

p1 = gmsh.model.occ.addPoint(0, 0, 0)
p2 = gmsh.model.occ.addPoint(1, 0, 0)
p3 = gmsh.model.occ.addPoint(1, 3, 0)
p4 = gmsh.model.occ.addPoint(3, 3, 0)
p5 = gmsh.model.occ.addPoint(3, 3 + g, 0)
p7 = gmsh.model.occ.addPoint(2, 3 + 2*R, 0)
p8 = gmsh.model.occ.addPoint(1, 3 + 2*R, 0)
p9 = gmsh.model.occ.addPoint(1, 7, 0)
p10 = gmsh.model.occ.addPoint(0, 7, 0)
# Inner loop
p11 = gmsh.model.occ.addPoint(1, 3 + g, 0)
p12 = gmsh.model.occ.addPoint(2, 3 + g, 0)
p13 = gmsh.model.occ.addPoint(1.5, 3 + 2*R - g, 0)
p14 = gmsh.model.occ.addPoint(1, 3 + 2*R - g, 0)


l1 = gmsh.model.occ.addLine(p1, p2)  # --> Inlet
l2 = gmsh.model.occ.addLine(p2, p3)
l3 = gmsh.model.occ.addLine(p3, p4)
l4 = gmsh.model.occ.addLine(p4, p5)
l6_arc1 = gmsh.model.occ.addLine(p5, p7)
l7 = gmsh.model.occ.addLine(p7, p8)
l8 = gmsh.model.occ.addLine(p8, p9)
l9 = gmsh.model.occ.addLine(p9, p10)  # --> Outlet
l10 = gmsh.model.occ.addLine(p10, p1)
l11 = gmsh.model.occ.addLine(p11, p12)
l12_arc2 = gmsh.model.occ.addLine(p13, p12)
l13 = gmsh.model.occ.addLine(p13, p14)
l14 = gmsh.model.occ.addLine(p14, p11)
loop1 = gmsh.model.occ.addCurveLoop([l1, l2, l3, l4, l6_arc1, l7, l8, l9, l10])
loop2 = gmsh.model.occ.addCurveLoop([l11, l12_arc2, l13, l14])
shape1 = gmsh.model.occ.addPlaneSurface([loop1])
shape2 = gmsh.model.occ.addPlaneSurface([loop2])
gmsh.model.occ.synchronize()

Omega_fluid = gmsh.model.occ.cut([(gdim, shape1)], [(gdim, shape2)])
gmsh.model.occ.synchronize()

# Add physical groups
gmsh.model.addPhysicalGroup(2, [Omega_fluid[0][0][1]], 1)
inlet_marker, outlet_marker, wall_marker, wall2_marker = 101, 102, 201, 202
gmsh.model.addPhysicalGroup(1, [l1], inlet_marker)  # --> Inlet
gmsh.model.addPhysicalGroup(1, [l4], outlet_marker)  # --> Outlet
gmsh.model.addPhysicalGroup(
    1, [l2, l3,  l6_arc1, l8, l7, l9, l10], wall_marker)  # --> Walls
gmsh.model.addPhysicalGroup(
    1, [l11, l12_arc2, l13, l14], wall2_marker)  # --> InnerWalls


# Meshing options
gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.2)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.5)

gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")

# Convert Gmsh model to DOLFINx mesh
gdim = 2
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
model_rank = 0
domain, cell_tags, facet_tags = gmshio.model_to_mesh(
    gmsh.model, mesh_comm, model_rank, gdim=gdim)
domain.topology.create_entities(domain.topology.dim - 1)
domain.topology.create_connectivity(
    domain.topology.dim, domain.topology.dim - 1)
gmsh.finalize()
facet_tags.name = "Facet markers"
with dolfinx.io.XDMFFile(mesh_comm, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facet_tags, domain.geometry)

# < BCu >
bcs_p: List[oasisx.PressureBC] = []
bc_inlet_x = DirichletBC(dolfinx.fem.Constant(domain, 0.),
                         method=LocatorMethod.TOPOLOGICAL, marker=(facet_tags, inlet_marker))  # method= marker =
bc_inlet_y = DirichletBC(dolfinx.fem.Constant(domain, 1.),
                         method=LocatorMethod.TOPOLOGICAL, marker=(facet_tags, inlet_marker))
bc_wall = DirichletBC(dolfinx.fem.Constant(domain, 0.), method=LocatorMethod.TOPOLOGICAL,
                      marker=(facet_tags, wall_marker))
bc_wall2 = DirichletBC(dolfinx.fem.Constant(domain, 0.), method=LocatorMethod.TOPOLOGICAL,
                       marker=(facet_tags, wall2_marker))
bcs_u = [[bc_inlet_x, bc_wall, bc_wall2], [bc_inlet_y, bc_wall, bc_wall2]]

# < BCp >
bcs_p: List[oasisx.PressureBC] = [oasisx.PressureBC(
    0.,  marker=(facet_tags, outlet_marker))]


# fractional step solver
solver = FractionalStep_AB_CN(
    mesh=domain,
    u_element=("Lagrange", 2),
    p_element=("Lagrange", 1),
    bcs_u=bcs_u,
    bcs_p=bcs_p,
    solver_options={"tentative": {"ksp_type": "preonly", "pc_type": "lu"}, "pressure": {
        "ksp_type": "preonly", "pc_type": "lu"}, "scalar": {"ksp_type": "preonly", "pc_type": "lu"}},
    body_force=None
)

# Time-stepping
T_start, T_end, dt = 0.0, 1, 0.001
num_steps = int((T_end - T_start) // dt)
with dolfinx.io.VTXWriter(domain.comm, "_OASISX_Test1-PipeNetwork_v1.bp", [solver.u], engine="BP4") as writer:
    for step in tqdm(range(num_steps)):
        solver.solve(dt, nu=0.01, max_iter=100)
        writer.write(dt*step)

Please always inspect your boundary tags visually before making simulations.
Resulting velocity field

[p1rx]

wow, that is embarrassing. Idk why the numbering of the edges was messed in gmsh but indeed my final problem with running simulaitons was related to wrong selection of boundary facets.

Btw, would anyone wonder, method to check boundaries is available in viskex but also you can use pyvista (i couldnt find anywhere example to do so, so i am posting here a short script):

# < visualise dofs of facets with pvista >

V = dolfinx.fem.FunctionSpace(domain, ("CG", 1))
grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(domain))
tag = inlet_marker  # Example tag

# Locate degrees of freedom associated with the specific tag
facet_indices = facet_tags.indices[facet_tags.values == tag]
dofs = dolfinx.fem.locate_dofs_topological(V, domain.topology.dim - 1, facet_indices)

# return coords of highlighted nodes:
coords = V.tabulate_dof_coordinates()
dof_coords = coords[dofs]

# highlight dof coordinates:
points = pyvista.PolyData(dof_coords)
points["dofs"] = np.zeros(dof_coords.shape[0]) 

plotter = pyvista.Plotter()
plotter.add_mesh(grid, color="lightgray", show_edges=True)  # Show the mesh
plotter.add_points(points, color="red", point_size=10) # Show highlighted points
plotter.camera_position = [(7, 6, 1), (7, 6, 0), (0, 0, 0)] 
plotter.camera.zoom(0.18)  
plotter.show()

[QibangLiu]

Hi, I copy and paste the code to run, but the results of velocity is 0?

I am using dolfinx 0.9.0 and the oasisx verion is 1.2.0.

Thanks.

[dokken]

Hm, seems like the GMSH mesh generation has changed in the latest release, meaning that facets aren’t tagged in the correct way.
Here is an updated code that runs nicely with the 0.9 docker image of dolfinx with oasisx 1.2:

from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
from dolfinx import fem, mesh
from dolfinx.io import gmshio, XDMFFile, VTXWriter
import oasisx
from oasisx import DirichletBC, LocatorMethod, FractionalStep_AB_CN
import ufl
import gmsh
from tqdm import tqdm
from typing import List
import numpy as np

# < G M S H >
gmsh.initialize()
g = 0.3
R = 1
gdim = 2

gmsh.initialize()

p1 = gmsh.model.occ.addPoint(0, 0, 0)
p2 = gmsh.model.occ.addPoint(1, 0, 0)
p3 = gmsh.model.occ.addPoint(1, 3, 0)
p4 = gmsh.model.occ.addPoint(3, 3, 0)
p5 = gmsh.model.occ.addPoint(3, 3 + g, 0)
p7 = gmsh.model.occ.addPoint(2, 3 + 2 * R, 0)
p8 = gmsh.model.occ.addPoint(1, 3 + 2 * R, 0)
p9 = gmsh.model.occ.addPoint(1, 7, 0)
p10 = gmsh.model.occ.addPoint(0, 7, 0)
# Inner loop
p11 = gmsh.model.occ.addPoint(1, 3 + g, 0)
p12 = gmsh.model.occ.addPoint(2, 3 + g, 0)
p13 = gmsh.model.occ.addPoint(1.5, 3 + 2 * R - g, 0)
p14 = gmsh.model.occ.addPoint(1, 3 + 2 * R - g, 0)


l1 = gmsh.model.occ.addLine(p1, p2)  # --> Inlet
l2 = gmsh.model.occ.addLine(p2, p3)
l3 = gmsh.model.occ.addLine(p3, p4)
l4 = gmsh.model.occ.addLine(p4, p5)
l6_arc1 = gmsh.model.occ.addLine(p5, p7)
l7 = gmsh.model.occ.addLine(p7, p8)
l8 = gmsh.model.occ.addLine(p8, p9)
l9 = gmsh.model.occ.addLine(p9, p10)  # --> Outlet
l10 = gmsh.model.occ.addLine(p10, p1)
l11 = gmsh.model.occ.addLine(p11, p12)
l12_arc2 = gmsh.model.occ.addLine(p13, p12)
l13 = gmsh.model.occ.addLine(p13, p14)
l14 = gmsh.model.occ.addLine(p14, p11)
loop1 = gmsh.model.occ.addCurveLoop([l1, l2, l3, l4, l6_arc1, l7, l8, l9, l10])
loop2 = gmsh.model.occ.addCurveLoop([l11, l12_arc2, l13, l14])
shape1 = gmsh.model.occ.addPlaneSurface([loop1])
shape2 = gmsh.model.occ.addPlaneSurface([loop2])
gmsh.model.occ.synchronize()

Omega_fluid = gmsh.model.occ.cut([(gdim, shape1)], [(gdim, shape2)])
gmsh.model.occ.synchronize()

# Add physical groups
gmsh.model.addPhysicalGroup(2, [Omega_fluid[0][0][1]], 1)
inlet_marker, outlet_marker, wall_marker = 101, 102, 201
bndry = gmsh.model.getBoundary(Omega_fluid[0], combined=False, oriented=False)
inlet_curves = []
outlet_curves = []
wall_curves = []
for tag in bndry:
    cmm = gmsh.model.occ.getCenterOfMass(tag[0], tag[1])
    if np.allclose(cmm, [0.5, 0.0, 0.0]):
        inlet_curves.append(tag[1])
    elif np.allclose(cmm, [0.5, 7.0, 0.0]):
        outlet_curves.append(tag[1])
    else:
        wall_curves.append(tag[1])


gmsh.model.addPhysicalGroup(1, inlet_curves, inlet_marker)  # --> Inlet
gmsh.model.addPhysicalGroup(1, outlet_curves, outlet_marker)  # --> Outlet
gmsh.model.addPhysicalGroup(1, wall_curves, wall_marker)  # --> Walls


# Meshing options
gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.2)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.5)

gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")
gmsh.write("mesh_d.msh")
# Convert Gmsh model to DOLFINx mesh
gdim = 2
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
model_rank = 0
domain, cell_tags, facet_tags = gmshio.model_to_mesh(
    gmsh.model, mesh_comm, model_rank, gdim=gdim
)
domain.topology.create_entities(domain.topology.dim - 1)
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim - 1)
gmsh.finalize()
facet_tags.name = "Facet markers"
with dolfinx.io.XDMFFile(mesh_comm, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facet_tags, domain.geometry)

# < BCu >
bcs_p: List[oasisx.PressureBC] = []
bc_inlet_x = DirichletBC(
    dolfinx.fem.Constant(domain, 0.0),
    method=LocatorMethod.TOPOLOGICAL,
    marker=(facet_tags, inlet_marker),
)  # method= marker =
bc_inlet_y = DirichletBC(
    dolfinx.fem.Constant(domain, 1.0),
    method=LocatorMethod.TOPOLOGICAL,
    marker=(facet_tags, inlet_marker),
)
bc_wall = DirichletBC(
    dolfinx.fem.Constant(domain, 0.0),
    method=LocatorMethod.TOPOLOGICAL,
    marker=(facet_tags, wall_marker),
)
bcs_u = [[bc_inlet_x, bc_wall], [bc_inlet_y, bc_wall]]

# < BCp >
bcs_p: List[oasisx.PressureBC] = [
    oasisx.PressureBC(0.0, marker=(facet_tags, outlet_marker))
]


# fractional step solver
solver = FractionalStep_AB_CN(
    mesh=domain,
    u_element=("Lagrange", 2),
    p_element=("Lagrange", 1),
    bcs_u=bcs_u,
    bcs_p=bcs_p,
    solver_options={
        "tentative": {"ksp_type": "preonly", "pc_type": "lu"},
        "pressure": {
            "ksp_type": "preonly",
            "pc_type": "lu",
            "ksp_error_if_not_converged": True,
        },
        "scalar": {
            "ksp_type": "preonly",
            "pc_type": "lu",
        },
    },
    body_force=None,
)

# Time-stepping
T_start, T_end, dt = 0.0, 1, 0.001
num_steps = int((T_end - T_start) // dt)
with dolfinx.io.VTXWriter(
    domain.comm, "_OASISX_Test1-PipeNetwork_v1.bp", [solver.u], engine="BP4"
) as writer:
    for step in tqdm(range(num_steps)):
        solver.solve(dt, nu=0.01, max_iter=100)
        writer.write(dt * step)