Dirichlet boundary for pressure in Oasisx

dolfinx
1

Hi,

I am trying to use Oasisx for a problem where velocity is induced due to application of pressure at both inlet and outlet. I tried to modify this example from the forum. However, I always get zero velocity for the output. I would appreciate if someone could fix the problem with my attempt. Here is the MWE executed on 0.8.0 installed through conda:

from mpi4py import MPI
import dolfinx
import numpy as np
import oasisx
from oasisx import DirichletBC, LocatorMethod, FractionalStep_AB_CN
from typing import List
from tqdm.auto import tqdm


def inlet(x):
    return np.isclose(x[0], 0)


def wall(x):
    return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1))


def outlet(x):
    return np.isclose(x[0], 10)


domain = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([10, 1])], [
    20, 50], cell_type=dolfinx.mesh.CellType.triangle)

bcs_p: List[oasisx.PressureBC] = []
bc_inlet_x = DirichletBC(dolfinx.fem.Constant(domain, 1.),
                         method=LocatorMethod.GEOMETRICAL, marker=inlet)
bc_inlet_y = DirichletBC(dolfinx.fem.Constant(domain, 0.),
                         method=LocatorMethod.GEOMETRICAL, marker=inlet)
bc_wall = DirichletBC(dolfinx.fem.Constant(domain, 0.), method=LocatorMethod.GEOMETRICAL,
                      marker=wall)

inlet_facets = dolfinx.mesh.locate_entities_boundary(domain, domain.topology.dim - 1, inlet)
outlet_facets = dolfinx.mesh.locate_entities_boundary(domain, domain.topology.dim - 1, outlet)

ft = dolfinx.mesh.meshtags(domain, domain.topology.dim-1, np.hstack((inlet_facets, outlet_facets)),
                           np.hstack((np.full_like(inlet_facets, 1, dtype=np.int32), np.full_like(outlet_facets, 2, dtype=np.int32))))

bcs_p: List[oasisx.PressureBC] = [oasisx.PressureBC(8.0, marker=(ft, 1)), oasisx.PressureBC(0., marker=(ft, 2))]
bcs_u = [[bc_wall], [bc_wall]]


# 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, 0.1, 0.001
num_steps = int((T_end - T_start) // dt)
with dolfinx.io.VTXWriter(domain.comm, "poisseuille.bp", [solver.u], engine="BP4") as writer:
    for step in tqdm(range(num_steps), total=num_steps):
        solver.solve(dt, nu=1, max_iter=10)
        writer.write(dt*step)
2

You have not done the correct marking of the meshtags (I think).
Here is an adapted code:

from mpi4py import MPI
import dolfinx
import numpy as np
import oasisx
from oasisx import DirichletBC, LocatorMethod, FractionalStep_AB_CN
from typing import List
from tqdm.auto import tqdm


L = 10
H = 2

def inlet(x):
    return np.isclose(x[0], 0)


def wall(x):
    return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], H))


def outlet(x):
    return np.isclose(x[0], L)




domain = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([L, H])], [
    20, 20], cell_type=dolfinx.mesh.CellType.triangle)

bc_wall = DirichletBC(dolfinx.fem.Constant(domain, 0.), method=LocatorMethod.GEOMETRICAL,
                      marker=wall)

inlet_facets = dolfinx.mesh.locate_entities_boundary(domain, domain.topology.dim - 1, inlet)
outlet_facets = dolfinx.mesh.locate_entities_boundary(domain, domain.topology.dim - 1, outlet)


fdim = domain.topology.dim - 1
domain.topology.create_entities(fdim)
facet_map = domain.topology.index_map(fdim)
num_facets_local = facet_map.size_local + facet_map.num_ghosts
markers = np.full(num_facets_local, 0, dtype=np.int32)
markers[inlet_facets] = 1
markers[outlet_facets] = 2


ft = dolfinx.mesh.meshtags(domain, domain.topology.dim-1, np.arange(num_facets_local,dtype=np.int32), markers)


class InletPressure():
    def __init__(self, t:float):
        self.t = t
    def eval(self, x):
        scale = self.t if self.t < 1 else 1 
        return np.full_like(x[0], scale*8)
                

p_inlet = InletPressure(0.)
bcs_p: List[oasisx.PressureBC] = [oasisx.PressureBC(p_inlet.eval, marker=(ft, 1)), oasisx.PressureBC(0.0, marker=(ft, 2))]
bcs_u = [[bc_wall], [bc_wall]]


# 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, 10, 1e-3
num_steps = int((T_end - T_start) // dt)
p_inlet.t = T_start
with dolfinx.io.VTXWriter(domain.comm, "poisseuille.bp", [solver.u], engine="BP4") as writer, dolfinx.io.VTXWriter(domain.comm, "phi.bp", [solver._p], engine="BP4") as writer_phi:
    for step in tqdm(range(num_steps), total=num_steps):
        p_inlet.t += dt
        solver.solve(dt, nu=0.1, max_iter=10)
        writer.write(dt*step)
        writer_phi.write(dt*step)

However, the following does not give a good result, which I do believe is from the fact that having pressure inlet and outlet conditions in a splitting scheme is not a good idea.

1 Like
3

Thank you for the quick fix!

I am not familiar with CFD but to my eye you used “a” splitting method in your own tutorial here. Is that substantially different?

Could you point out whether Oasisx or another fenics based solver has an approach that could work efficiently for my case?

4

That splitting method is a historical artifact of the original FEniCS tutorial from 2016.
One should never set pressure conditions strong on inlets or outlets of a boundary when one consider Navier-Stokes.
See for instance: Pressure boundary conditions for blood flows | Chinese Annals of Mathematics, Series B

They are especially not well suited for splitting schemes, as the first step, the tentative velocity step, in absence of any non-zero Dirichlet conditions in pressure, is solely driven by the tentative pressure gradient.

1 Like