Integral within ufl form

General
1

Hi Dolfinx community!

I’m trying to implement a simple 0D resistance model downstream of my computational domain. The idea is to impose a Neumann-type boundary condition based on the mean pressure and flow rate:

P = QR; ~~~~Q = \int_{\Gamma} \underline{u} \cdot \hat{n} ~ d\Gamma

I’m testing this on a simple Stokes formulation:

sol = dolfinx.fem.Function(V)
(u, p) = ufl.split(sol) 
(w, q) = ufl.TestFunctions(V)

# parameters
mu = dolfinx.fem.Constant(msh, dolfinx.default_scalar_type(1.0))

# build weak form 
weak_form = mu * ufl.inner(ufl.grad(w), ufl.grad(u)) * ufl.dx \
            - p * ufl.div(w) * ufl.dx \
            - q * ufl.div(u) * ufl.dx

I need to add the following to the form for the 0D model:

-\int_{\Gamma_{out}} \underline{w} \cdot \underline{\hat{n}} \left( R \int_{\Gamma_{out}} \underline{u} \cdot \underline{\hat{n}} d\Gamma\right) d\Gamma

The problem is that the flow rate needs to be computed from \underline{u}. Precomputing (below) uses the default value velocity (\underline{0}) and negates the resistance term.

# RESISTANCE BC
right_facets = dolfinx.mesh.locate_entities_boundary(msh, fdim, lambda x: np.isclose(x[0], 1.0))
facet_values = np.ones_like(right_facets, dtype=np.int32)
facet_tags = dolfinx.mesh.meshtags(msh, fdim, right_facets, facet_values)
ds_out = ufl.Measure('ds', domain=msh, subdomain_data=facet_tags, subdomain_id=1)
n = ufl.FacetNormal(msh)

R = dolfinx.fem.Constant(msh, dolfinx.default_scalar_type(1.0))
Q_form = dolfinx.fem.form(ufl.inner(u, n) * ds_out)
Q = dolfinx.fem.assemble_scalar(Q_form)
print(f"Flow rate: {Q}")
res_form = -R * Q * ufl.inner(w, n) * ds_out
weak_form += res_form

Additionally, two UFL forms can’t be multiplied:

res_form = -R * ufl.inner(w, n) * ds_out
flow_form = ufl.inner(u, n) * ds_out
weak_form += (res_form * flow_form)

-> TypeError: unsupported operand type(s) for *: 'Form' and 'Form'

Question: What is the recommended approach in Dolfinx to implement a boundary condition that depends on the integral of the solution itself (i.e., an implicit 0D model)?

Full MWE:

import dolfinx 
import ufl
import numpy as np
import basix
from mpi4py import MPI
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from petsc4py import PETSc

def pressureBC(x):
    return np.zeros(x.shape[1])

def no_slip(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1])))

def inflow(x):
    return np.stack((4.0 * x[1] * (1.0 - x[1]), np.zeros(x.shape[1])))

# Create mesh 
msh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 20, 20, dolfinx.mesh.CellType.quadrilateral)

# Function space
PE = basix.ufl.element('CG', msh.basix_cell(), degree=1)
QE = basix.ufl.element('CG', msh.basix_cell(), degree=1, shape=(msh.topology.dim,))
ME = basix.ufl.mixed_element([QE, PE])
V = dolfinx.fem.functionspace(msh, ME)
u_sub, u_submap = V.sub(0).collapse()
p_sub, p_submap = V.sub(1).collapse()

# Create connectivity
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)

# Mark boundaries
facet_tags = dolfinx.mesh.meshtags(msh, fdim, np.array([0, 1, 2, 3]), np.array([1, 2, 3, 4]))

# Define boundary conditions
x = ufl.SpatialCoordinate(msh)

# Define degree of freedom locations
left_dofs = dolfinx.fem.locate_dofs_geometrical((V.sub(0), u_sub), lambda x: np.isclose(x[0], 0.0))
right_dofs = dolfinx.fem.locate_dofs_geometrical((V.sub(0), u_sub), lambda x: np.isclose(x[0], 1.0))
right_dofs_p = dolfinx.fem.locate_dofs_geometrical((V.sub(1), p_sub), lambda x: np.isclose(x[0], 1.0))
bottom_dofs = dolfinx.fem.locate_dofs_geometrical((V.sub(0), u_sub), lambda x: np.isclose(x[1], 0.0))
top_dofs = dolfinx.fem.locate_dofs_geometrical((V.sub(0), u_sub), lambda x: np.isclose(x[1], 1.0))

# Interpolate boundary functions to fenics function
u_bc_walls = dolfinx.fem.Function(u_sub)
u_bc_walls.interpolate(no_slip)
u_in = dolfinx.fem.Function(u_sub)
u_in.interpolate(inflow)
p_bc_right = dolfinx.fem.Function(p_sub)
p_bc_right.interpolate(pressureBC)

# Define Dirichlet boundary conditions
bc_u_l = dolfinx.fem.dirichletbc(u_in, left_dofs, V.sub(0))
bc_u_b = dolfinx.fem.dirichletbc(u_bc_walls, bottom_dofs, V.sub(0))
bc_u_t = dolfinx.fem.dirichletbc(u_bc_walls, top_dofs, V.sub(0))
bc_p_r = dolfinx.fem.dirichletbc(p_bc_right, right_dofs_p, V.sub(1))

bcs = [bc_u_l, bc_u_b, bc_u_t] #, bc_p_r]


# Define variational problem
sol = dolfinx.fem.Function(V)
(u, p) = ufl.split(sol) 
(w, q) = ufl.TestFunctions(V)

# parameters
mu = dolfinx.fem.Constant(msh, dolfinx.default_scalar_type(1.0))

# build weak form 
weak_form = mu * ufl.inner(ufl.grad(w), ufl.grad(u)) * ufl.dx \
            - p * ufl.div(w) * ufl.dx \
            - q * ufl.div(u) * ufl.dx 


# RESISTANCE BC
right_facets = dolfinx.mesh.locate_entities_boundary(msh, fdim, lambda x: np.isclose(x[0], 1.0))
facet_values = np.ones_like(right_facets, dtype=np.int32)
facet_tags = dolfinx.mesh.meshtags(msh, fdim, right_facets, facet_values)
ds_out = ufl.Measure('ds', domain=msh, subdomain_data=facet_tags, subdomain_id=1)
n = ufl.FacetNormal(msh)

R = dolfinx.fem.Constant(msh, dolfinx.default_scalar_type(1.0))
Q_form = dolfinx.fem.form(ufl.inner(u, n) * ds_out)
Q = dolfinx.fem.assemble_scalar(Q_form)
print(f"Flow rate: {Q}")
res_form = -R * ufl.inner(w, n) * ds_out
flow_form = ufl.inner(u, n) * ds_out
weak_form += (res_form * flow_form)

# stab
he = ufl.CellDiameter(msh)
beta = 1.0/3.0
tau = beta * ((he * he)/(4.0 * mu))
stab = -tau * ufl.inner(ufl.grad(q), ufl.grad(p)) * ufl.dx
weak_form += stab

problem = NonlinearProblem(F=weak_form, u=sol, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-7
solver.report = True

ksp = solver.krylov_solver
opts = PETSc.Options()
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"] = "umfpack"
ksp.setFromOptions()

n, converged = solver.solve(sol)

u = sol.sub(0).collapse()
u.name = 'u'
p = sol.sub(1).collapse()
p.name = 'p'

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "output/u.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh)
    xdmf.write_function(u, 0.0)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "output/p.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh)
    xdmf.write_function(p, 0.0)
    
Q = dolfinx.fem.assemble_scalar(Q_form)
print(f"Flow rate: {Q}")