How to extract pressure solution from inlet

General
1

Hello everybody,

I’ve successfully solved Stokes’ problem in the context of Poiseuille flow. Now, I’m looking to calculate the disparity between the input pressure and the outlet pressure. My question is how to extract the solution pressure from the inlet. While it’s straightforward to obtain the pressure solution at the outlet, I’m encountering difficulties when attempting to retrieve it from the inlet.

Could you please help me?
Thank you,

2

Then please post the code you have already obtained, so that people don’t have to start form scratch while trying to reply to you. It stands to reason that the code for computing it at the inlet and at the outlet will be similar.

1 Like
3

Below is my code for computing numerical solutions, i don’t know how to extract the pressure solution from only inlet and outlet boundaries. Specifically, I aim to calculate the difference in the mean pressure obtained from the inlet and outlet.
Thank you.

#======================================================
# Importing modules and packages
#======================================================


from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
import dolfinx.fem.petsc
import dolfinx.io
from basix.ufl import element, mixed_element
from dolfinx import mesh, fem
from dolfinx.fem import (Function, dirichletbc, form, functionspace,
                         locate_dofs_topological)
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner

import matplotlib.pyplot as plt
from dolfinx import plot
import pyvista
from dolfinx.plot import vtk_mesh
import warnings
warnings.filterwarnings("ignore")


#======================================================
# Numerical Resolution
#======================================================

def wall_boundary(x):
    return np.logical_or(np.isclose(x[1], -0.7), np.isclose(x[1], 0.7))

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

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


def u_expression(x):
    k=1
    h=0.7
    mu=1
    u_max=k*(h**2)/(2*mu)
    return np.stack((u_max*(1-(x[1]/h)**2), np.zeros(x.shape[1])))

def Dirichlet_inlet(msh,W,W0):
    u_in=Function(W0)
    u_in.interpolate(u_expression)
    facets = locate_entities_boundary(msh, 1, in_boundary)
    dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
    bc0 = dirichletbc(u_in, dofs, W.sub(0))
    return bc0

def Dirichlet_walls(msh,W,W0):
    u_wall= Function(W0)
    facets1 = locate_entities_boundary(msh, 1, wall_boundary)
    dofs1 = locate_dofs_topological((W.sub(0), W0), 1, facets1)
    bc1 = dirichletbc(u_wall, dofs1, W.sub(0))
    return bc1


def mixed_direct(TH, msh):

    # Define function spaces
    W = functionspace(msh, TH)
    W0, _ = W.sub(0).collapse()
    W1, W1_to_W = W.sub(1).collapse()

    # Dirichlet conditions
    bc0=Dirichlet_inlet(msh,W,W0)
    bc1=Dirichlet_walls(msh,W,W0)
    

    bcs = [bc0,bc1]

    # Define the variational form

    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    
    f = Function(W0)
    a = form((inner(grad(u), grad(v)) - inner(p, div(v)) + inner(div(u), q)) * dx)
    L = form(inner(f, v) * dx)

    # Assemble the bilinear form 'a' and linear form 'L'
    A = fem.petsc.assemble_matrix(a, bcs=bcs)
    A.assemble()
    b = fem.petsc.assemble_vector(L)

    fem.petsc.apply_lifting(b, [a], bcs=[bcs])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

    # Set Dirichlet boundary values in the right-hand side
    fem.petsc.set_bc(b, bcs)

    # Create and configure solver

    ksp = PETSc.KSP().create(msh.comm)
    ksp.setOperators(A)
    ksp.setType("preonly")
  

    # LU solver

    pc = ksp.getPC()
    pc.setType("lu")
    pc.setFactorSolverType("mumps")

    # Solve the system AU=b

    U = Function(W)
    try:
        ksp.solve(b, U.vector)
    except PETSc.Error as e:
        if e.ierr == 92:
            print("The required PETSc solver/preconditioner is not available. Exiting.")
            print(e)
            exit(0)
        else:
            raise e
    ksp.solve(b, U.vector)

    # Retrieve velocity (u) and pressure (p) solution

    u, p = U.sub(0).collapse(), U.sub(1).collapse()
    u.x.scatter_forward()
    p.x.scatter_forward()

    return u, p
4

If you look at the function you use to assign Dirichlet BCs

you already have some hints on how you could proceed:
locate the inlet boundary first to get facets, and then get the dofs on W0 with

dofs = locate_dofs_topological(W0, 1, facets)

Note that you don’t need the (W.sub(0), W0) pair in this case, because I am ussing you’ll want to use those dofs on the output u of the mixed_direct function, which is a dolfinx.fem.Function on W0. You can then get the values of the velocity on those dofs with u.x[dofs].

If your goal was just to get the average on the inlet and the average on the outlet, it may even be simpler to use a ds measure, and assemble u * ds(INLET_MARKER) with dolfinx.fem.assemble_scalar. Please look into the existing demos with these keywords (ds, assemble_scalar) and, if you need further help, post an updated example of what you managed to do.

1 Like
5

Ok, thank you so much, Sir. I will try to do what you’ve suggested.