How to retrieve numeric solution from boundary

[quinso]

I’m currently investigating the Stokes problem (poiseuille case) to ensure that the flow from the inlet to the outlet conserves mass. To achieve this, I’ve used the finite element method to calculate the velocity solution. Now, I aim to retrieve the outlet solution (were Neumann are applied) and use it as the inlet boundary condition for resolving the problem in the opposite direction.

My question is: How can I obtain the solution on the outlet?"

The problem i’m considering (from inlet to outlet) is :

-\Delta u+ \nabla p=f
\nabla.u=0
u=0 in \Gamma_{walls}
u= (0.245(1-\frac{y}{0.7^2}),0) in \Gamma_{in}
\nabla u.n-pn=0 in \Gamma_{out}

Here is my code to solve the problem but only from inlet to outlet:

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


# Boundary conditions

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

    # Définir les espaces de fonctions
    W = functionspace(msh, TH)
    W0, _ = W.sub(0).collapse()
    W1, W1_to_W = W.sub(1).collapse()

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

    bcs = [bc0,bc1]

    # Define the variational formulation

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

    # Définir les valeurs de condition de bord de Dirichlet dans le second membre
    fem.petsc.set_bc(b, bcs)

    # Create and configure solver

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


    # Create a null space

    null_vec = dolfinx.fem.Function(W)
    null_vec.x.array[W1_to_W] = 1.0
    dolfinx.la.orthonormalize([null_vec.vector])
    nsp = PETSc.NullSpace().create(vectors=[null_vec.vector])
    A.setNearNullSpace(nsp)
    A.setNullSpace(nsp)

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

    # compute velocity (u) and pressure (p) solutions

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


    p_avg = W.mesh.comm.allreduce(
        fem.assemble_scalar(fem.form(p*dx)), op=MPI.SUM)
    p.x.array[:] -= p_avg

    return u, p

the numeric solution is shown below :

I want to get the solution like this :

[francesco-ballarin]

You can access the degrees of freedom with u.x.array. Flip the sign by multiplying the content of the array by minus 1, and then use u exactly as you do with u_in in Dirichlet_inlet or u_wall in Dirichlet_walls.

[quinso]

Thank you Sir for your answer. Do you mean like this:

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

def u_expression_inv(x):
   u_inv=-u_h.x.array
   return np.stack((u_inv, np.zeros(x.shape[1])))

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

[francesco-ballarin]

I can’t try it myself because the example you posted is not reproducible (for instance, it doesn’t build the mesh). The example could even be more minimal that what you write: create a mesh, create a P2-P1 function space, interpolate a function on that space, and then extract the u component.

Still, I think that if you write

u_h.x.array[:] *= -1.0

it will flip the sign of all DOFs. Note that both x and y component will be flipped, but that doesn’t matter because the y component of your velocity is zero anyways.

In a more general case, you would have to extract only the x component out of u_h.
To help you with that, you should use the second output of

W00, map0x = W0.sub(0).collapse()

This gives you the list of DOFs which correspond to the x component of the velocity, so if you want to modify just them (and not the y component) you can use

u_h.x.array[map0x] *= -1.0

Note that you are doing something similar when setting up the null space (just for the W space instead of W0)