Impose Velocity Profile on Inlet Face using Array

[bcu]

I am trying to impose a time-varying velocity profile to the inlet face of a rectangular domain as a Dirichlet Boundary Condition. The velocity profile is in the form of numpy array.

t2    = 0.0
inlet = 0.0
class InletVelocity(UserExpression):
    def __init__(self, t2, inlet, **kwargs):
        super().__init__(kwargs)
        self._t2 = t2
        self._vel = inlet
    def eval_cell(self, value, x, ufc_cell):        
        value[0,:] = inlet
        value[1,:] = 0
    def value_shape(self):
        return (2,)
    
v_in = InletVelocity(t2,inlet)

# Define boundaries
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 8)'
walls = 'near(x[1], 0) || near(x[1], 2)'

for n in range(num_steps):
    # Updating inlet and outlet boundary conditions
    print('t =', t)
    v_in.t2     = t2
    ui          = vel_prof(ai,omega,mu2,rho2,waves,dp_dz,t2) 
    v_in.inlet  = ui
    # Update current time
    t += dt
    # Define boundary conditions
    bcu_inflow = DirichletBC(W.sub(0), v_in, inflow)

How can I use the time-varying numpy array as a velocity profile for the inlet at each time step and imposing it on the inlet face?

[dokken]

Is the velocity profile going to vary both in space and time (is it constant per time step?). If it is constant per time step, I do not see your issue. If it is varying in time, is the data points ordered from lowest y-value to highest y-value, and is the number of data points in the numpy array one to one with the number of dofs in the function space?

If so, you could do the following:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(4,4)

V = VectorFunctionSpace(mesh, "CG", 2)
u_inlet = Function(V)

# As we have a four cells along each side,
# we have 9 degrees of freedom in a "CG-2" space along the left side
values_x = np.array([0.1,0.3,0.5,0.8,1.1,0.8,0.7,0.3,0.2])
values_y = np.zeros(9)
dof_coords = V.tabulate_dof_coordinates()

# Identify coordinates where x[0]=0
dofs_on_left_side = np.flatnonzero(np.isclose(dof_coords[:,0], 0))
left_dofs_coords = dof_coords[dofs_on_left_side]
#Sort these coord by increasing y-coordinate
sorted_indices = np.argsort(left_dofs_coords[:,1])

# From the sorted dofs, decide which is the x and y component of the space
x_dofs = V.sub(0).dofmap().dofs()
y_dofs = V.sub(1).dofmap().dofs()
x_c, y_c = 0, 0
for index in sorted_indices:
    dof = dofs_on_left_side[index]
    if dof in x_dofs:
        u_inlet.vector().vec().setValueLocal(dof, values_x[x_c])
        x_c +=1
    elif dof in y_dofs:
        u_inlet.vector().vec().setValueLocal(dof, values_y[y_c])
        y_c +=1
XDMFFile("u_inlet.xdmf").write_checkpoint(u_inlet, "u", 0)

And repeat this insertion into u_inlet once per time step (which you can have assigned to a boundary condition).

[bcu]

Thank you for this. The velocity profile varies in space and time but only in the normal direction x-direction. It’s a velocity profile based on pulsatile flow. For Example these are the velocity profiles at three different time frames.

# =====================================================================
# Create mesh
# =====================================================================

channel = Rectangle(Point(0, 0), Point(14, 0.6))
domain = channel
mesh = generate_mesh(domain, 64)

# =====================================================================
# Define function spaces
# =====================================================================

V_e = VectorElement("CG", mesh.ufl_cell(), 1)
Q_e = FiniteElement("CG", mesh.ufl_cell(), 1)
W_e = V_e * Q_e
W = FunctionSpace(mesh, W_e)
# =====================================================================
u_inlet = Function(W)

for n in range(num_steps):
    # Updating inlet and outlet boundary conditions
    print('t =', t)
    inlet_velocity = <<<ARRAY solution at Time t>>>>
    values_x = inlet_velocity
    values_y = np.zeros(14)
    dof_coords = W.tabulate_dof_coordinates()
    # Identify coordinates where x[0]=0
    dofs_on_left_side = np.flatnonzero(np.isclose(dof_coords[:,0], 0))
    left_dofs_coords = dof_coords[dofs_on_left_side]
    print(left_dofs_coords)
    #Sort these coord by increasing y-coordinate
    sorted_indices = np.argsort(left_dofs_coords[:,1])
    
    # From the sorted dofs, decide which is the x and y component of the space
    x_dofs = W.sub(0).dofmap().dofs()
    y_dofs = W.sub(1).dofmap().dofs()
    x_c, y_c = 0, 0
    for index in sorted_indices:
        dof = dofs_on_left_side[index]
        if dof in x_dofs:
            u_inlet.vector().vec().setValueLocal(dof, values_x[x_c])
            x_c +=1
        elif dof in y_dofs:
            u_inlet.vector().vec().setValueLocal(dof, values_y[y_c])
            y_c +=1

    # Update current time
    t += dt
    # Define boundary conditions
    bcu_inflow = DirichletBC(W.sub(0), u_inlet, inflow)
    # Update current time
    t += dt
    # Define boundary conditions
    bcu_inflow = DirichletBC(W.sub(0), u_inlet, inflow)
    # free to slip walls
    bcu_walls = DirichletBC(W.sub(0).sub(1), Constant(0), walls)
    #bcp_outflow = DirichletBC(W.sub(1), outflow_profile, outflow)
    bcs = [bcu_inflow, bcu_walls]
    problem         = NonlinearVariationalProblem(F, w, bcs, J=J)
    solver          = NonlinearVariationalSolver(problem)

I’m getting the following error.

*** Error: Unable to create Dirichlet boundary condition.
*** Reason: Illegal value dimension (3), expecting (2).
*** Where: This error was encountered inside DirichletBC.cpp.

[dokken]

Please format your code without the > as there is no way of nicely copying your code in its current state.

Note that there are several things with your code that could be improved.
First of all, to set the boundary condition, make u_inlet a function of the subspace of velocity and not the mixed space (solution sketched below):

from mshr import *
channel = Rectangle(Point(0, 0), Point(14, 0.6))
domain = channel
mesh = generate_mesh(domain, 64)

V_e = VectorElement("CG", mesh.ufl_cell(), 1)
Q_e = FiniteElement("CG", mesh.ufl_cell(), 1)
W_e = V_e * Q_e
W = FunctionSpace(mesh, W_e)
u_inlet = Function(W.sub(0).collapse())
u_inlet.vector()[:] = 1
bc = DirichletBC(W.sub(0), u_inlet, "on_boundary")
q = Function(W)
bc.apply(q.vector())
File("u.pvd") << q.sub(0, deepcopy=True)

As other issues, I would create the problem, solver, and therefore the variational formulation and BCS outside the temporal loop. You can update the boundary condition by changing the vector of u_inlet (and not redefining the function). Similarly for other temporal expressions.

[bcu]

Thank you. I’m still a little lost as I’m fairly new to Fenics. I’ve made the updates you suggested but still having some trouble. It might help to see the full code.

# =====================================================================
# Create mesh
# =====================================================================
channel = Rectangle(Point(0, 0), Point(14, 0.6))
domain = channel
mesh = generate_mesh(domain, 64)

# =====================================================================
# Define function spaces
# =====================================================================
V_e = VectorElement("CG", mesh.ufl_cell(), 1)
Q_e = FiniteElement("CG", mesh.ufl_cell(), 1)
W_e = V_e * Q_e
W = FunctionSpace(mesh, W_e)

# =====================================================================
# COLLECT BOUNDARY CONDITIONS
# =====================================================================
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 14)'
walls = 'near(x[1], 0) || near(x[1], 0.6)'

u_inlet = Function(W.sub(0).collapse())
u_inlet.vector()[:] = 1

# Define boundary conditions
bcu_inflow = DirichletBC(W.sub(0), u_inlet, inflow)
# free to slip walls
bcu_walls = DirichletBC(W.sub(0).sub(1), Constant(0), walls)
#bcp_outflow = DirichletBC(W.sub(1), outflow_profile, outflow)
bcs = [bcu_inflow, bcu_walls]


# define Jacobian
dW              = TrialFunction(W)
J               = derivative(F,w,dW)
a = lhs(F)
L = rhs(F)

# =====================================================================

problem         = NonlinearVariationalProblem(F, w, bcs, J=J)
solver          = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['maximum_iterations'] = 20
solver.parameters['newton_solver']['relative_tolerance'] = 1.0e-6

# =====================================================================
# create files for storing
# =====================================================================
uFile         = File(outputFileV, 'compressed')
pFile         = File(outputFileP, 'compressed')
# =====================================================================
t       = dt
count   = 0

for n in range(num_steps):
    # Updating inlet and outlet boundary conditions
    # As we have a four cells along each side,
    # we have 9 degrees of freedom in a "CG-2" space along the left side
    values_x = vel_prof(ai,omega,mu2,rho2,waves,dp_dz,t)
    values_y = np.zeros(14)
    dof_coords = W.tabulate_dof_coordinates()
    
    # Identify coordinates where x[0]=0
    dofs_on_left_side = np.flatnonzero(np.isclose(dof_coords[:,0], 0))
    left_dofs_coords = dof_coords[dofs_on_left_side]
    #Sort these coord by increasing y-coordinate
    sorted_indices = np.argsort(left_dofs_coords[:,1])
    
    # From the sorted dofs, decide which is the x and y component of the space
    x_dofs = W.sub(0).dofmap().dofs()
    y_dofs = W.sub(1).dofmap().dofs()
    x_c, y_c = 0, 0
    for index in sorted_indices:
        dof = dofs_on_left_side[index]
        if dof in x_dofs:
            u_inlet.vector().vec().setValueLocal(dof, values_x[x_c])
            x_c +=1
        elif dof in y_dofs:
            u_inlet.vector().vec().setValueLocal(dof, values_y[y_c])
            y_c +=1

    # Solve for velocity and pressure fields
    print("Solving.....")
    solver.solve()
    
    (u,p) = w.split()
    N1 = 20 #save every N time steps
    if (count/N1).is_integer():
        uFile << u
        pFile << p
        
    print("Written Velocity and Pressure Data")
    w0.assign(w)
    t += dt

I’m still getting the following error.

u_inlet.vector().vec().setValueLocal(dof, values_x[x_c])
  File "PETSc/Vec.pyx", line 688, in petsc4py.PETSc.Vec.setValueLocal (src/petsc4py.PETSc.c:98842)
petsc4py.PETSc.Error: error code 63

[dokken]

@bcu

 x_dofs = W.sub(0).dofmap().dofs()
    y_dofs = W.sub(1).dofmap().dofs()

These should be

V = W.sub(0).collapse()
x_dofs = V.sub(0).dofmap().dofs()

etc.
Please reduce your example to a minimal problem, there are so many functions in this code that are not necessary to reproduce your error. See the guidelines in: Read before posting: How do I get my question answered? - #3

[bcu]

Thank you. I was able to fix it using your comments.

[RokyC]

x_dofs = V.sub(0).dofmap().dofs() , it doesn’t seem to work at Dolfinx
Traceback (most recent call last):
File “”, line 1, in
TypeError: ‘DofMap’ object is not callable