Get velocity field at a boundary

dolfinx
1

Dear all,

I am solving the convection-difussion equation on a rectangular domain.
The boundary conditions are:

  • Dirichlet BC at Inlet (parabolic velocity profile).
  • Homogeneous Dirichlet BC at Walls.
  • Zero stress BC at outlet.

The working code is shown below:

import numpy as np
from mpi4py import MPI
import dolfinx
from dolfinx.fem import (Constant, Function, FunctionSpace, apply_lifting, assemble_matrix,
                         assemble_scalar, assemble_vector, locate_dofs_topological, set_bc,
                         dirichletbc, form)
from dolfinx import fem
from ufl import (inner, div, dx, nabla_grad, sym, Identity, as_tensor, Dn, grad, ds, rhs, lhs, Measure, transpose, as_tensor,
                dot, FacetNormal, indices)
from dolfinx.io import gmshio
from petsc4py import PETSc
import ufl
import gmsh

rank = MPI.COMM_WORLD.rank
mesh_comm = MPI.COMM_WORLD
model_rank = 0

# Problem parameters
L = 1
H = L/10
nu = 0.01
C = 1
Re = C*L/nu

# Numerical parameters
T = 2
dt = 1e-3
num_steps = int(T/dt)

# Mesh generation
gdim = 2

gmsh.initialize()
if rank == 0:
    rectangle = gmsh.model.occ.addRectangle(0,0,0, L, H, tag=1)
    gmsh.model.occ.synchronize()

fluid_marker = 1
if rank == 0:
    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")

inlet_marker, outlet_marker, wall_marker = 2, 3, 4
inflow, outflow, walls = [], [], []
if rank == 0:
    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, [0, H/2, 0]):
            inflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L, H/2, 0]):
            outflow.append(boundary[1])
        else:
            walls.append(boundary[1])
    gmsh.model.addPhysicalGroup(1, walls, wall_marker)
    gmsh.model.setPhysicalName(1, wall_marker, "Walls")
    gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
    gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
    gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
    gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")

# Refinement
lc = 1
if mesh_comm.rank == model_rank:
    gmsh.model.mesh.field.add("Box", 1)
    gmsh.model.mesh.field.setNumber(1, "VIn", lc/50)
    gmsh.model.mesh.field.setNumber(1, "VOut", lc)
    gmsh.model.mesh.field.setNumber(1, "XMin", 0)
    gmsh.model.mesh.field.setNumber(1, "XMax", L)
    gmsh.model.mesh.field.setNumber(1, "YMin", 0)
    gmsh.model.mesh.field.setNumber(1, "YMax", H)
    gmsh.model.mesh.field.add("Min", 2)
    gmsh.model.mesh.field.setNumbers(2, "FieldsList", [1])
    gmsh.model.mesh.field.setAsBackgroundMesh(2)
    gmsh.model.occ.synchronize()

if mesh_comm.rank == model_rank:
    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.model.mesh.generate(gdim)
    gmsh.model.mesh.setOrder(2)
    gmsh.model.mesh.optimize("Netgen")

# Generate mesh and facets
mesh, _, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
ft.name = "Facet markers"

# Define surface differential at boundaries
ds_Outlet = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=outlet_marker)
ds_Inlet = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=inlet_marker)
ds_Wall = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=wall_marker)

# Define normal vector
n = FacetNormal(mesh)

# Define function spaces
from ufl import VectorElement, TrialFunction, TestFunction
# Velocity element space
v_cg2 = VectorElement('Lagrange', mesh.ufl_cell(), 2) # Order 2
V = FunctionSpace(mesh, v_cg2)
# Trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

# Create initial condition
u_n = Function(V)
u_n.name = "u_n"

# Inlet BC
class Inlet():
    def __init__(self, C, H):
        self.C = C
        self.H = H

    def __call__(self, x):
        values = np.zeros((gdim, x.shape[1]),dtype=PETSc.ScalarType)
        values[0] = (1.0/self.C)*x[1]*(self.H - x[1])
        values[1] = 0
        return values

# Walls BC
class Walls():
    def __init__(self, t):
        self.t = T

    def __call__(self, x):
        values = np.zeros((gdim, x.shape[1]),dtype=PETSc.ScalarType)
        values[0] = 0
        values[1] = 0
        return values


fdim = gdim - 1

outlet_facets = ft.indices[ft.values == outlet_marker]
inlet_facets = ft.indices[ft.values == inlet_marker]
wall_facets = ft.indices[ft.values == wall_marker]

t = 0
iterations = 0

# Dirichlet boundaries
inlet_velocity = Inlet(C, H)
noslip = Walls(t)

u_D_I = fem.Function(V)
u_D_W = fem.Function(V)
u_D_I.interpolate(inlet_velocity)
u_D_W.interpolate(noslip)

bc_I = fem.dirichletbc(u_D_I, fem.locate_dofs_topological(V, fdim, inlet_facets))
bc_W = fem.dirichletbc(u_D_W, fem.locate_dofs_topological(V, fdim, wall_facets))
bc = [bc_I, bc_W]

unit_vector_x = dolfinx.fem.Constant(mesh, np.array([1.0, 0.0]))
k = dolfinx.fem.Constant(mesh, dt)

''' Variational Form '''
# Unknown terms of the variational form
# Temporal component
T_comp = ufl.inner(u,v)*ufl.dx
# Spatial component
S_comp = ufl.inner(ufl.dot(u,unit_vector_x),ufl.div(v))*ufl.dx
S_comp -= ufl.inner((1/Re)*ufl.grad(u),ufl.grad(v))*ufl.dx
# Boundary component
B_comp = ufl.inner(ufl.dot(ufl.grad(u),n),v)*ds_Outlet
B_comp = ufl.inner(ufl.dot(ufl.grad(u),n),v)*ds_Wall

# Unknown terms of the variational form
# Temporal component
T_comp_old = ufl.inner(u_n,v)*ufl.dx
# Spatial component
S_comp_old = ufl.inner(ufl.dot(u_n,unit_vector_x),ufl.div(v))*ufl.dx
S_comp_old -= ufl.inner((1/Re)*ufl.grad(u_n),ufl.grad(v))*ufl.dx
# Boundary component
B_comp_old = ufl.inner(ufl.dot(ufl.grad(u_n),n),v)*ds_Outlet
B_comp_old = ufl.inner(ufl.dot(ufl.grad(u_n),n),v)*ds_Wall

# Bilinear and linear form
Form = (1/k)*(T_comp - T_comp_old) - (1/2)*(S_comp + S_comp_old) - (1/2)*(B_comp + B_comp_old)
a = fem.form(ufl.lhs(Form))
L = fem.form(ufl.rhs(Form))

# Assemble left-hand side matrix A and right-hand side vector b
A = fem.petsc.assemble_matrix(a, bcs=bc)
A.assemble()
b = fem.petsc.create_vector(L)

# Linear system solver
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.GMRES)

import tqdm.notebook
if rank == 0:
    progress = tqdm.tqdm(desc="Solving PDE", total=num_steps)

for iterations in range(num_steps):
    if rank == 0:
        progress.update(1)

    t += dt

    # Update the right hand side reusing the initial vector
    with b.localForm() as loc_b:
        loc_b.set(0)
    fem.petsc.assemble_vector(b, L)

    # Apply Dirichlet boundary condition to the vector
    apply_lifting(b, [a], [bc])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b, bc)

    # Solve linear problem
    solver.solve(b, u_n.vector)
    u_n.x.scatter_forward()

I need to get the velocity field at the right boundary, i.e., the outlet boundary.
I know the dofs of the boundary using fem.locate_dofs_topological(V, fdim, outlet_facets), but if I just do u_n.x.array[fem.locate_dofs_topological(V, fdim, outlet_facets)], the values are not in the correct order.

Does anybody know how to solve this issue?

Many thanks in advance.

2

You are using a function space with a block size (this is due to the usage of vectorelement). This means that the dofs are not unrolled.
Consider this MWE:

import dolfinx
import numpy as np
from dolfinx import fem
from dolfinx.fem import Function, FunctionSpace, locate_dofs_topological
from mpi4py import MPI
from ufl import VectorElement

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

def boundary(x):
    return x[0]<0.5

boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim-1, boundary)

v_cg2 = VectorElement('Lagrange', mesh.ufl_cell(), 2) # Order 2
V = FunctionSpace(mesh, v_cg2)
u = fem.Function(V)

dof_coords = V.tabulate_dof_coordinates()
bs = V.dofmap.bs
boundary_dofs = locate_dofs_topological(V,mesh.topology.dim-1, boundary_facets)


def velocity(x):
    return (np.sin(x[0]), x[1])
# Create initial condition
u_n = Function(V)
u_n.name = "u_n"
u_n.interpolate(velocity)

for dof in boundary_dofs:
    coord = dof_coords[dof]
    dofs = np.array([dof*bs+i for i in range(bs)], dtype=np.int32)
    exact = [np.sin(coord[0]), coord[1]]
    print("Coordinate", coord, "Function", u_n.x.array[dofs], "Exact", exact[0], exact[1])
    for i in range(bs):
        assert np.isclose(u_n.x.array[dofs[i]], exact[i])

Note that it is often useful to try to reduce your problem to something on a built-in mesh, and with a known function (usually something found through interpolation).

3

Thank you very much, @dokken . I’ve adapted my code based on your response and it works!