"Problems with the simulated results" - Navier-Stokes equations in 3D

General
1

Good day, community. During this time, I have tried to recreate the [Tutorial 02: Navier-Stokes problem] problem in 3D; however, the solutions I obtain are inconsistent. When running the program, no error occurs, but the visualization, whether in VISKEX or PARAVIEW, is not coherent, or the entire volume has the same values. The 2D problem has worked correctly, but in 3D, I believe there are missing parameters in the boundary conditions. Could you please help me find the error? I am attaching the code:

import typing
import basix.ufl
import dolfinx
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import numpy.typing
import petsc4py.PETSc
import ufl
import viskex
from petsc4py import PETSc
from dolfinx.io import VTXWriter


nu = 0.01

def u_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the flat velocity profile at the inlet."""
    values = np.zeros((3, x.shape[1]))
    values[0, :] = 1.0
    return values

def u_wall_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the zero velocity at the wall."""
    return np.zeros((3, x.shape[1]))


pre_step_length = 1.0
after_step_length = 2.0
pre_step_height = 1.0
after_step_height = 2.0
mesh_size = 0.1
height_z = 1.0  

gmsh.initialize()
gmsh.model.add("mesh")

p0 = gmsh.model.geo.addPoint(0.0, after_step_height - pre_step_height, 0.0, mesh_size)
p1 = gmsh.model.geo.addPoint(pre_step_length, after_step_height - pre_step_height, 0.0, mesh_size)
p2 = gmsh.model.geo.addPoint(pre_step_length, 0.0, 0.0, mesh_size)
p3 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, 0.0, 0.0, mesh_size)
p4 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, after_step_height, 0.0, mesh_size)
p5 = gmsh.model.geo.addPoint(0.0, after_step_height, 0.0, mesh_size)

p6 = gmsh.model.geo.addPoint(0.0, after_step_height - pre_step_height, height_z, mesh_size)
p7 = gmsh.model.geo.addPoint(pre_step_length, after_step_height - pre_step_height, height_z, mesh_size)
p8 = gmsh.model.geo.addPoint(pre_step_length, 0.0, height_z, mesh_size)
p9 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, 0.0, height_z, mesh_size)
p10 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, after_step_height, height_z, mesh_size)
p11 = gmsh.model.geo.addPoint(0.0, after_step_height, height_z, mesh_size)

l0 = gmsh.model.geo.addLine(p0, p1)
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p5)
l5 = gmsh.model.geo.addLine(p5, p0)

l6 = gmsh.model.geo.addLine(p6, p7)
l7 = gmsh.model.geo.addLine(p7, p8)
l8 = gmsh.model.geo.addLine(p8, p9)
l9 = gmsh.model.geo.addLine(p9, p10)
l10 = gmsh.model.geo.addLine(p10, p11)
l11 = gmsh.model.geo.addLine(p11, p6)

l12 = gmsh.model.geo.addLine(p0, p6)
l13 = gmsh.model.geo.addLine(p1, p7)
l14 = gmsh.model.geo.addLine(p2, p8)
l15 = gmsh.model.geo.addLine(p3, p9)
l16 = gmsh.model.geo.addLine(p4, p10)
l17 = gmsh.model.geo.addLine(p5, p11)

ll0 = gmsh.model.geo.addCurveLoop([l0, l13, -l6, -l12])
ll1 = gmsh.model.geo.addCurveLoop([l1, l14, -l7, -l13])
ll2 = gmsh.model.geo.addCurveLoop([l2, l15, -l8, -l14])
ll3 = gmsh.model.geo.addCurveLoop([l3, l16, -l9, -l15])
ll4 = gmsh.model.geo.addCurveLoop([l4, l17, -l10, -l16])
ll5 = gmsh.model.geo.addCurveLoop([l5, l12, -l11, -l17])

ll_bottom = gmsh.model.geo.addCurveLoop([l0, l1, l2, l3, l4, l5])
ll_top = gmsh.model.geo.addCurveLoop([l6, l7, l8, l9, l10, l11])

s0 = gmsh.model.geo.addPlaneSurface([ll_bottom])
s1 = gmsh.model.geo.addPlaneSurface([ll_top])
s2 = gmsh.model.geo.addPlaneSurface([ll0])
s3 = gmsh.model.geo.addPlaneSurface([ll1])
s4 = gmsh.model.geo.addPlaneSurface([ll2])
s5 = gmsh.model.geo.addPlaneSurface([ll3])
s6 = gmsh.model.geo.addPlaneSurface([ll4])
s7 = gmsh.model.geo.addPlaneSurface([ll5])

surface_loop = gmsh.model.geo.addSurfaceLoop([s0, s1, s2, s3, s4, s5, s6, s7])
volume = gmsh.model.geo.addVolume([surface_loop])

gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(2, [s0, s1, s2, s3, s4, s5, s6, s7], 1)
gmsh.model.addPhysicalGroup(3, [volume], 2)

gmsh.model.mesh.generate(3)

mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=3)
gmsh.finalize()

mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)


boundaries_1 = boundaries.indices[boundaries.values == 1]
boundaries_2 = boundaries.indices[boundaries.values == 10]


viskex.dolfinx.plot_mesh(mesh)
#Function Space
V_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim, ))
Q_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)

#Formulation problem
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)

    # 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)

    # Boundary conditions
    W0 = W.sub(0)
    V, _ = W0.collapse()
    u_in = dolfinx.fem.Function(V)
    u_in.interpolate(u_in_eval)
    u_wall = dolfinx.fem.Function(V)
    u_wall.interpolate(u_wall_eval)
    bdofs_V_1 = dolfinx.fem.locate_dofs_topological((W0, V), mesh.topology.dim - 1, boundaries_1)
    bdofs_V_2 = dolfinx.fem.locate_dofs_topological((W0, V), mesh.topology.dim - 1, boundaries_2)
    inlet_bc = dolfinx.fem.dirichletbc(u_in, bdofs_V_1, W0)
    wall_bc = dolfinx.fem.dirichletbc(u_wall, bdofs_V_2, W0)
    bc = [inlet_bc, wall_bc]

    # 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

        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 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
            return x

        def update_solution(self, x: petsc4py.PETSc.Vec) -> None:  # type: ignore[no-any-unimported]
            """Update `self._solution` with data in `x`."""
            x.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
            with x.localForm() as _x, self._solution.x.petsc_vec.localForm() as _solution:
                _solution[:] = _x

        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_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=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
            dolfinx.fem.petsc.set_bc(F_vec, self._bcs, x, -1.0)

        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()
            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()

    # 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)

    # Solve
    snes = petsc4py.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()
    return up
    
up_m = run_monolithic()

(u_m, p_m) = (up_m.sub(0).collapse(), up_m.sub(1).collapse())

from pathlib import Path
folder = Path("resultsnew")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, folder / "poiseuille_u.bp", u_m, engine="BP4")
vtx_p = VTXWriter(mesh.comm, folder / "poiseuille_p.bp", p_m, engine="BP4")
vtx_u.write(0)
vtx_p.write(0)
vtx_u.close()
vtx_p.close()

viskex.dolfinx.plot_vector_field(u_m, "u")
viskex.dolfinx.plot_vector_field(u_m, "u", glyph_factor=0.25)
viskex.dolfinx.plot_scalar_field(p_m, "p")
2

This is definitely wrong. The boundary marker 10 doesn’t actually existing in your tags, hence:

  • you are not imposing the wall dirichlet BC anywhere
  • you are imposing Dirichlet BCs everywhere, which means that the pressure is defined up to a constant. A simple fix would be to drop the Dirichlet BC at the outlet.

For future posts, it would be appreciated if you could provide a minimal description of your problem (geometry, saying that it was about steady NS, boundary conditions) before showing us the code. I could immediately see that it was a 3D version of the backward facing step because I authored the original tutorial, but for everyone else it would have been helpful to have that short description.

3

Thank you very much for pointing out the error in the code. I have already corrected it, and now the solutions are consistent. However, I want to compare the solutions with those obtained in some literature, but the problem is that I don’t know the units or dimensions. In FEniCSX, is everything dimensionless? Although the kinematic viscosity gives me the impression that it is not dimensionless, could you tell me how dimensions and units work in FEniCSX?

I am attaching the code in case someone else needs it.

import typing
import basix.ufl
import dolfinx
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import numpy.typing
import petsc4py.PETSc
import ufl
import viskex
from petsc4py import PETSc
from dolfinx.io import VTXWriter

nu = 0.01

def u_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the flat velocity profile at the inlet."""
    values = np.zeros((3, x.shape[1]))
    values[0, :] = 1.0
    return values

def u_wall_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the zero velocity at the wall."""
    return np.zeros((3, x.shape[1]))


pre_step_length = 1.0
after_step_length = 2.0
pre_step_height = 1.0
after_step_height = 2.0
mesh_size = 0.25
height_z = 1.0  

gmsh.initialize()
gmsh.model.add("mesh")

p0 = gmsh.model.geo.addPoint(0.0, after_step_height - pre_step_height, 0.0, mesh_size)
p1 = gmsh.model.geo.addPoint(pre_step_length, after_step_height - pre_step_height, 0.0, mesh_size)
p2 = gmsh.model.geo.addPoint(pre_step_length, 0.0, 0.0, mesh_size)
p3 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, 0.0, 0.0, mesh_size)
p4 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, after_step_height, 0.0, mesh_size)
p5 = gmsh.model.geo.addPoint(0.0, after_step_height, 0.0, mesh_size)

p6 = gmsh.model.geo.addPoint(0.0, after_step_height - pre_step_height, height_z, mesh_size)
p7 = gmsh.model.geo.addPoint(pre_step_length, after_step_height - pre_step_height, height_z, mesh_size)
p8 = gmsh.model.geo.addPoint(pre_step_length, 0.0, height_z, mesh_size)
p9 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, 0.0, height_z, mesh_size)
p10 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, after_step_height, height_z, mesh_size)
p11 = gmsh.model.geo.addPoint(0.0, after_step_height, height_z, mesh_size)

l0 = gmsh.model.geo.addLine(p0, p1)
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p5)
l5 = gmsh.model.geo.addLine(p5, p0)

l6 = gmsh.model.geo.addLine(p6, p7)
l7 = gmsh.model.geo.addLine(p7, p8)
l8 = gmsh.model.geo.addLine(p8, p9)
l9 = gmsh.model.geo.addLine(p9, p10)
l10 = gmsh.model.geo.addLine(p10, p11)
l11 = gmsh.model.geo.addLine(p11, p6)

l12 = gmsh.model.geo.addLine(p0, p6)
l13 = gmsh.model.geo.addLine(p1, p7)
l14 = gmsh.model.geo.addLine(p2, p8)
l15 = gmsh.model.geo.addLine(p3, p9)
l16 = gmsh.model.geo.addLine(p4, p10)
l17 = gmsh.model.geo.addLine(p5, p11)

ll0 = gmsh.model.geo.addCurveLoop([l0, l13, -l6, -l12])
ll1 = gmsh.model.geo.addCurveLoop([l1, l14, -l7, -l13])
ll2 = gmsh.model.geo.addCurveLoop([l2, l15, -l8, -l14])
ll3 = gmsh.model.geo.addCurveLoop([l3, l16, -l9, -l15])
ll4 = gmsh.model.geo.addCurveLoop([l4, l17, -l10, -l16])
ll5 = gmsh.model.geo.addCurveLoop([l5, l12, -l11, -l17])

ll_bottom = gmsh.model.geo.addCurveLoop([l0, l1, l2, l3, l4, l5])
ll_top = gmsh.model.geo.addCurveLoop([l6, l7, l8, l9, l10, l11])

s0 = gmsh.model.geo.addPlaneSurface([ll_bottom])
s1 = gmsh.model.geo.addPlaneSurface([ll_top])
s2 = gmsh.model.geo.addPlaneSurface([ll0])
s3 = gmsh.model.geo.addPlaneSurface([ll1])
s4 = gmsh.model.geo.addPlaneSurface([ll2])
s5 = gmsh.model.geo.addPlaneSurface([ll3])
s6 = gmsh.model.geo.addPlaneSurface([ll4])
s7 = gmsh.model.geo.addPlaneSurface([ll5])

surface_loop = gmsh.model.geo.addSurfaceLoop([s0, s1, s2, s3, s4, s5, s6, s7])
volume = gmsh.model.geo.addVolume([surface_loop])

gmsh.model.geo.synchronize()


gmsh.model.addPhysicalGroup(2, [s7], 1)
gmsh.model.addPhysicalGroup(2, [s0, s1, s2, s3, s4, s6], 2)
gmsh.model.addPhysicalGroup(3, [volume], 0)

gmsh.model.mesh.generate(3)

mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=3)
gmsh.finalize()

mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)


boundaries_1 = boundaries.indices[boundaries.values == 1]
boundaries_2 = boundaries.indices[boundaries.values == 2]


viskex.dolfinx.plot_mesh(mesh)
viskex.dolfinx.plot_mesh_tags(mesh, boundaries, "boundaries")

#Function Space
V_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim, ))
Q_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)

#Formulation problem
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)

    # 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)

    # Boundary conditions
    W0 = W.sub(0)
    V, _ = W0.collapse()
    u_in = dolfinx.fem.Function(V)
    u_in.interpolate(u_in_eval)
    u_wall = dolfinx.fem.Function(V)
    u_wall.interpolate(u_wall_eval)
    bdofs_V_1 = dolfinx.fem.locate_dofs_topological((W0, V), mesh.topology.dim - 1, boundaries_1)
    bdofs_V_2 = dolfinx.fem.locate_dofs_topological((W0, V), mesh.topology.dim - 1, boundaries_2)
    inlet_bc = dolfinx.fem.dirichletbc(u_in, bdofs_V_1, W0)
    wall_bc = dolfinx.fem.dirichletbc(u_wall, bdofs_V_2, W0)
    bc = [inlet_bc, wall_bc]

    # 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

        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 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
            return x

        def update_solution(self, x: petsc4py.PETSc.Vec) -> None:  # type: ignore[no-any-unimported]
            """Update self._solution with data in x."""
            x.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
            with x.localForm() as _x, self._solution.x.petsc_vec.localForm() as _solution:
                _solution[:] = _x

        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_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=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
            dolfinx.fem.petsc.set_bc(F_vec, self._bcs, x, -1.0)

        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()
            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()

    # 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)

    # Solve
    snes = petsc4py.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()
    return up
    
up_m = run_monolithic()

(u_m, p_m) = (up_m.sub(0).collapse(), up_m.sub(1).collapse())

from pathlib import Path
folder = Path("resultsnew")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, folder / "poiseuille_u.bp", u_m, engine="BP4")
vtx_p = VTXWriter(mesh.comm, folder / "poiseuille_p.bp", p_m, engine="BP4")
vtx_u.write(0)
vtx_p.write(0)
vtx_u.close()
vtx_p.close()

viskex.dolfinx.plot_vector_field(u_m, "u")
viskex.dolfinx.plot_vector_field(u_m, "u", glyph_factor=1)
viskex.dolfinx.plot_scalar_field(p_m, "p")
4

You decide the dimensions when you prepare the mesh. For instance, you could decide that one unit of length of the mesh is 1 mm, 1cm, 1 km, …, and then determine the constants in the equations accordingly.

5

Thank you very much, I understand the issue with the dimensions, but the kinematic viscosity is causing me problems. In the variational form, it seems to have units. I have used them in the international system, but when the viscosity is low, the solution becomes inconsistent or cannot be visualized. So, I don’t understand why this is happening. Could you help me, please?

...
nu = 0.01
...
F = (nu* ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
         + ufl.inner(ufl.grad(u)...
6

If the viscosity is small, then the Reynolds number is large. In such cases you may need to take a finer mesh and/or stabilize the FEM scheme. Furthermore, searching for a steady solution in those cases may not be physically meaningful.

7

Because I have been thinking that instead of using the kinematic viscosity, it is better to use the density and dynamic viscosity in their respective places within the variational form. This is because when I use the kinematic viscosity of air, ν = 1.5e-5 (with units of m²/s [L²*T⁝š]), the solutions become inconsistent. So, I wanted to know your opinion on changing the variational form, as well as being able to create a more structured mesh. The simulation aims to model a ventilation system to control the parameters.

8

I would scale the equation (non-dimensionalize it), as discussed in Test problem 1: Channel flow (Poiseuille flow) — FEniCSx tutorial
and throughly explained in Scaling of Differential Equations | SpringerLink

9

Thank you very much for your answers. One last question, regarding what Francesco said, a stable solution for this problem wouldn’t make physical sense. So, would it be recommended to switch to a transient solution?

10

Yes, you should change formulation

1 Like