Navier stokes equation with a spatially varing vector body force

General
1

Hi,

I’m trying to simulate a flow problem with no slip boundary conditions, free boundary conditions, and setting a pressure of 0 to one node.

I was able to set up the boundary conditions, but I am also trying to include a constant body force. I included this in the variational formulation:

Force_per_unit_vol = fem.Function(V)
Force_per_unit_vol.interpolate(lambda x: body_force_fun(x))

F = rho_air*inner(grad(u)*u, v)*dx + \
    mu_air*inner(grad(u), grad(v))*dx + \
    inner(Force_per_unit_vol, v) - \
    inner(p, div(v))*dx - \
    div(u)*q*dx

Is this correct? If so, what is the output of body_force_fun(x)? Is it a numpy array with the same shape as x? I’m looking to compute something like this:

def body_force_fun(x):
    
    F_max = 1.0e-6
    z_offset = 0.01

    r_xy = np.sqrt(x[0]**2 + x[1]**2)

    if r_xy < 1.0e-8:
        return [0, 0, 0]
    else:
        a = 0.005
        max_val = np.sqrt(a/2) * np.exp(-0.5)
        r_val   = r_xy*np.exp(-r_xy**2/(2*a^2)) / max_val
        
        b = 0.005
        z_val = np.exp(-(x[2]-z_offset)**2/(2*b^2))
    
        F_mag = F_max * r_val * z_val
        
        return F_mag / r_xy * [x[1], -x[0], 0]

The full code is below.

Thank you,
Alex

from dolfinx import mesh, fem, io
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

from mpi4py import MPI

import numpy as np

import ufl
from ufl import inner, div, grad, dx

import basix

import meshio





vol_mesh = meshio.read('mesh.vtk')

points = vol_mesh.points
cells  = vol_mesh.cells[1].data


c_el = ufl.Mesh(basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(points.shape[1],)))
msh = mesh.create_mesh(MPI.COMM_WORLD, cells, points, c_el)

fdim = msh.topology.dim - 1

msh.topology.create_connectivity(fdim, fdim+1)






# Array geometry variables
h_array_1 = 0.02     # Part of height of array's cylindrical section
h_array_2 = 0.035    # Height of array's truncated conical section 
h_gap     = 0.03     # Gap between array and horizontal surface where powder is falling 




rho_air = 1.225
mu_air  = 1.8e-2


mark_eps = 1.0e-8



# Function to mark array noslip
def noslip_boundary1(x):
    return np.logical_and(x[2] < h_array_2, x[2] > -h_array_1)


# Function to mark table surface noslip
def noslip_boundary2(x):
    return np.isclose(x[2], (h_array_2 + h_gap))




# the same as old ufl.VectorElement
# the same as old ufl.FiniteElement
P2 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 2, shape=(msh.geometry.dim, ))
P1 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 1)

V, Q = fem.functionspace(msh, P2), fem.functionspace(msh, P1)

# Create the function space
TH = basix.ufl.mixed_element([P2, P1])
W = fem.functionspace(msh, TH)
W0, _ = W.sub(0).collapse()



# No slip boundary condition 1
noslip  = fem.Function(V)
facets1 = mesh.locate_entities_boundary(msh, fdim, noslip_boundary1)
dofs1   = fem.locate_dofs_topological((W.sub(0), V), fdim, facets1)
bc0     = fem.dirichletbc(noslip, dofs1, W.sub(0))


# No slip boundary condition 2
facets2 = mesh.locate_entities_boundary(msh, fdim, noslip_boundary2)
dofs2   = fem.locate_dofs_topological((W.sub(0), V), fdim, facets2)
bc1     = fem.dirichletbc(noslip, dofs2, W.sub(0))



# Since for this problem the pressure is only determined up to a constant, we pin it at a point
zero = fem.Function(Q)
zero.x.array[:] = 0
dofs = fem.locate_dofs_geometrical(
    (W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, -h_array_1]).all(axis=1))
bc2 = fem.dirichletbc(zero, dofs, W.sub(1))




def body_force_fun(x):
    
    F_max = 1.0e-6
        
    z_offset = 0.01

    r_xy = np.sqrt(x[0]**2 + x[1]**2)


    if r_xy < 1.0e-8:
        return [0, 0, 0]
    else:
        a = 0.005
        max_val = np.sqrt(a/2) * np.exp(-0.5)
        r_val   = r_xy*np.exp(-r_xy**2/(2*a^2)) / max_val
        
        b = 0.005
        z_val = np.exp(-(x[2]-z_offset)**2/(2*b^2))
    
        F_mag = F_max * r_val * z_val
        
        return F_mag / r_xy * [x[1], -x[0], 0]



Force_per_unit_vol = fem.Function(V)
Force_per_unit_vol.interpolate(lambda x: body_force_fun(x))





# Collect Dirichlet boundary conditions
bcs = [bc0, bc1, bc2]

# Define variational problem
w = fem.Function(W)

(u, p) = ufl.split(w)
(v, q) = ufl.TestFunctions(W)

# Tentative velocity step
# https://github.com/bragostin/Fenics/blob/main/HTC_SquareDuct_Steady_Neumann_Minimal_Pub.py
F = rho_air*inner(grad(u)*u, v)*dx + \
    mu_air*inner(grad(u), grad(v))*dx + \
    inner(Force_per_unit_vol, v) - \
    inner(p, div(v))*dx - \
    div(u)*q*dx

dw = ufl.TrialFunction(W)
dF = ufl.derivative(F, w, dw)

problem = NonlinearProblem(F, w, bcs=bcs, J=dF)

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol   = 1e-6
solver.report = True
solver.error_on_nonconvergence = True
#solver.max_it = 10


# Compute the solution
solver.solve(w)


# Split the mixed solution and collapse
u = w.sub(0).collapse()
p = w.sub(1).collapse()


# Write the solution to file
with io.XDMFFile(MPI.COMM_WORLD, "pressure.xdmf", "w") as pfile_xdmf:
    p.x.scatter_forward()
    pfile_xdmf.write_mesh(msh)
    pfile_xdmf.write_function(p)

with io.XDMFFile(MPI.COMM_WORLD, "velocity.xdmf", "w") as ufile_xdmf:
    u.x.scatter_forward()
    P1 = basix.ufl.element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
    u1 = fem.Function(fem.functionspace(msh, P1))
    u1.interpolate(u)
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(u1)
2

The function you send in to interpolate, should take in x as a (3, num_points) numpy array,
and return a (num_components, num_points) numpy array. This is to ensure that the code is efficient when working with many points (i.e. all interpolation points in the function space). This is for instance explained in:

3

Hi,

Thank you for the fody force info I have working code. Right now I’m able to run it in 1.2 mins, but I will need to use higher body forces which will require a finder mesh. The scaling I’m getting is proportional to the number of cells^6 which I think comes from the size of the stiffness matrix.

Are there any things you could recommend to speed up the code?

Thank you,
Alex

The latest running code is below:

from dolfinx import mesh, fem, io, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

from mpi4py import MPI

import numpy as np

import ufl
from ufl import inner, div, grad, dx

import basix

import time


height = 0.06
width  = 0.08

n_h = 10
n_w = 10


msh = mesh.create_box(MPI.COMM_WORLD, ((-width/2, -width/2, 0.0), (width/2,width/2,height)), (n_h,n_h,n_w), mesh.CellType.tetrahedron)

fdim = msh.topology.dim - 1


rho_air = fem.Constant(msh, 1.225)
mu_air  = fem.Constant(msh, 1.8e-5)


# Function to mark noslip
def noslip_boundary1(x):
    return np.isclose(x[2], height)



P2 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 2, shape=(msh.geometry.dim, ))
P1 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 1)

V, Q = fem.functionspace(msh, P2), fem.functionspace(msh, P1)

# Create the function space
TH = basix.ufl.mixed_element([P2, P1])
W = fem.functionspace(msh, TH)
W0, _ = W.sub(0).collapse()


# No slip boundary condition 1
noslip  = fem.Function(V)
facets1 = mesh.locate_entities_boundary(msh, fdim, noslip_boundary1)
dofs1   = fem.locate_dofs_topological((W.sub(0), V), fdim, facets1)
bc0     = fem.dirichletbc(noslip, dofs1, W.sub(0))


# Since for this problem the pressure is only determined up to a constant, we pin it at a point
zero = fem.Function(Q)
zero.x.array[:] = 0
dofs = fem.locate_dofs_geometrical(
    (W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = fem.dirichletbc(zero, dofs, W.sub(1))



def body_force_fun(x):
        
    F_max    = 1.0e-2
    z_offset = 0.02

    r_xy = np.sqrt(x[0]**2 + x[1]**2)

    a = 0.005
    max_val = a * np.exp(-0.5)
    Fr_val  = np.exp(-r_xy**2/(2*a**2)) / max_val
        
    b = 0.005
    Fz_val = np.exp(-(x[2]-z_offset)**2/(2*b**2))

    F_mag = F_max * Fr_val * Fz_val
      
    return -F_mag * np.vstack((x[1], -x[0], np.zeros_like(x[2])))
        

Force_per_unit_vol = fem.Function(V)
Force_per_unit_vol.interpolate(lambda x: body_force_fun(x))




# Collect Dirichlet boundary conditions
bcs = [bc0, bc2]

# Define variational problem
w = fem.Function(W)

(u, p) = ufl.split(w)
(v, q) = ufl.TestFunctions(W)


# Tentative velocity step

# https://github.com/bragostin/Fenics/blob/main/HTC_SquareDuct_Steady_Neumann_Minimal_Pub.py

# https://math.stackexchange.com/questions/1952588/
# weak-form-of-steady-navier-stokes-equations-with-special-boundary-condition

F = rho_air*inner(grad(u)*u, v)*dx + \
    mu_air*inner(grad(u), grad(v))*dx + \
    inner(Force_per_unit_vol, v)*dx - \
    inner(p, div(v))*dx - \
    div(u)*q*dx

dw = ufl.TrialFunction(W)
dF = ufl.derivative(F, w, dw)

problem = NonlinearProblem(F, w, bcs=bcs, J=dF)

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol   = 1e-6
#solver.max_it = 2
solver.report = True
solver.error_on_nonconvergence = False



# Compute the solution
log.set_log_level(log.LogLevel.INFO)
start = time.time()
solver.solve(w)
end   = time.time()

print('\n')
print(msh.topology.index_map(msh.topology.dim).size_local)
print((end - start) / 60)


# Split the mixed solution and collapse
u = w.sub(0).collapse()
p = w.sub(1).collapse()


# Write the solution to file
with io.XDMFFile(MPI.COMM_WORLD, "pressure.xdmf", "w") as pfile_xdmf:
    p.x.scatter_forward()
    pfile_xdmf.write_mesh(msh)
    pfile_xdmf.write_function(p)

with io.XDMFFile(MPI.COMM_WORLD, "velocity.xdmf", "w") as ufile_xdmf:
    u.x.scatter_forward()
    P1 = basix.ufl.element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
    u1 = fem.Function(fem.functionspace(msh, P1))
    u1.interpolate(u)
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(u1)
    

with io.XDMFFile(MPI.COMM_WORLD, "body_force.xdmf", "w") as ufile_xdmf:
    Force_per_unit_vol.x.scatter_forward()
    P1 = basix.ufl.element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
    u1 = fem.Function(fem.functionspace(msh, P1))
    u1.interpolate(Force_per_unit_vol)
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(u1)
4

You have not customized the solver options at all.
That is the first place to start.
I would start with using KSP type «preonly» and an LU preconditioner (mumps or superlu_dist for factorisation), and then check if you get any speedup when going from 1 to 2 to 4 processes.

See for instance Dolfinx seems much slower than dolfin in solving nonlinear mechanics - #3 by dokken

5

Thank you!

Adding those options made the code 9x faster:

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol   = 1e-6
#solver.max_it = 2
solver.report = True
solver.error_on_nonconvergence = False

opts = PETSc.Options()  # type: ignore
ksp = solver.krylov_solver
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

solver.solve(w)

What do you mean processes? Is that in parallel?

6

If you run your code with mpirun -n 4 python3 mycode.py the computational load is distributed across 4 processes using the message passing interface (mpi4py is what is user-facing in python).

This can potentially speed up your code even further, if you have more than 1 million dofs in your problem

7

Hi, I tried running this and I get the following error:

mesh = _cpp.mesh.create_mesh(comm, cells, cmap._cpp_object, x, partitioner)
RuntimeError: A facet is connected to more than two cells.

I guess its because one of the exterior (?) facets is connected to more than two cells, but I haven’t run this line yet:

fdim = msh.topology.dim - 1
msh.topology.create_connectivity(fdim, fdim+1)

What do you recommend, try to clean up the mesh before I import it? What software would be good for that?

8

You need to be a bit careful when reading in the mesh. See:

for details