In the demo “demo_stokes.py”, I want to change the original inlet condition to the average velocity inlet condition, i.e. -\frac{1}{|\Gamma_{inlet}|}\int_{\Gamma_{inlet}}\mathbf{u}\cdot\mathbf{n}ds=u_0. What should I do?
Is u_0 a known quantity or unknown quantity in your simulation?
Additionally, is your inlet going to be a straight inlet, or will it curve?
The u_0 is a given quantity and the inlet is a straight inlet.
Then you can easily compute what the magnitude of u should be, as you can compute the area of the inlet, and then compute the magnitude in the normal direction. In turn, you can either explicitly build the expression (if the normal direction is known analytically), or use
Maybe this is a bit a silly question, but does this code apply to taking the normal to a certain facet?
I saw a previous code based on approximating the normal at every degree of freedom and represent it with a projection; Obtain velocity normal to boundary - #4 by charithjeewantha
Right now I try to take such a velocity inlet condition and try to multiply the solution by the normal vector. For my case, I have a 3D geometry were the normal should be taken to the inlet facet. Is the code in the link an example of how to solve such a problem in dolfinx (rather than the old dolfin version)?
Yes, it uses ufl.FacetNormal to do an interpolation into the space whose degrees of freedom are located at a certain boundary.
One can of course also use projection routines, similar to:
to compute the normal vector as well.
Ok, thanks for your message. I tried to set up just a simple code to represent this;
from dolfinx import fem, mesh, io
from mpi4py import MPI
from petsc4py import PETSc
from ufl import (FacetNormal, VectorElement, FunctionSpace, TrialFunction,
TestFunction, ds, inner, dot, SpatialCoordinate, div, grad, dx)
from dolfinx.fem import Function
import numpy as np
def project_facet_normals(mesh):
# Create a suitable vector function space for normals
V_vector = FunctionSpace(mesh, VectorElement("CG", mesh.ufl_cell(), 1, dim=2))
# Define facet normals
n = FacetNormal(mesh)
# Define trial and test functions for projection
u_normal = TrialFunction(V_vector)
v_normal = TestFunction(V_vector)
# Define forms for projection
a_normal = inner(u_normal, v_normal) * n
L_normal = inner(n, v_normal) * n
# Solve projection problem
nh = Function(V_vector)
problem_normal = fem.petsc.LinearProblem(a_normal, L_normal)
nh = problem_normal.solve()
return nh
# Example usage
comm = MPI.COMM_WORLD
mesh = mesh.create_unit_square(comm, 30, 30)
# Project normals
nh = project_facet_normals(mesh)
# Output to check or use `nh` for further calculations
with io.XDMFFile(comm, "facet_normals.xdmf", "w") as file:
file.write_mesh(mesh)
file.write_function(nh)
This gave me the error;
Traceback (most recent call last):
File "/home/username/normal.py", line 37, in <module>
nh = project_facet_normals(mesh)
File "/home/username/normal.py", line 26, in project_facet_normals
nh = fem.Function(V_vector)
File "/home/username/miniconda3/envs/env_name/lib/python3.10/site-packages/dolfinx/fem/function.py", line 276, in __init__
self._cpp_object = functiontype(dtype)(V._cpp_object)
AttributeError: 'FunctionSpace' object has no attribute '_cpp_object'
I don’t know whether it was actually correct to take the multiplication by the normal n in those lines, but replacing it by ds gave me the following error;
# Define forms for projection
a_normal = inner(u_normal, v_normal) * ds
L_normal = inner(n, v_normal) * ds
Traceback (most recent call last):
File "/home/username/normal.py", line 37, in <module>
nh = project_facet_normals(mesh)
File "/home/username/normal.py", line 22, in project_facet_normals
a_normal = inner(u_normal, v_normal) * ds
File "/home/username/miniconda3/envs/env_name/lib/python3.10/site-packages/ufl/measure.py", line 408, in __rmul__
domains = extract_domains(integrand)
File "/home/username/miniconda3/envs/env_name/lib/python3.10/site-packages/ufl/domain.py", line 340, in extract_domains
return sorted(join_domains(domainlist))
File "/home/username/miniconda3/envs/env_name/lib/python3.10/site-packages/ufl/domain.py", line 304, in join_domains
gdims.add(domain.geometric_dimension())
AttributeError: 'Mesh' object has no attribute 'geometric_dimension'
In what way should this be adjusted to correctly store the numerical normal vectors?
Maybe I am a bit a silly, I am new to dolfinx and I don’t know how to use your code.
Here is my complete question. The origin boundary condition are as follow:
domain: \Omega=[0, 1]\times[0, 1],
inlet condition: \mathbf{u}= [0, -1]^T at y=1,
sysmetry condition: \mathbf{u}\cdot\mathbf{n}=\mathbf{u}.x=0 at x=0,
wall condition: \mathbf{u}=[0, 0]^T at y=0,
outlet condition: p=0 at x=1.
And now, I want to change the inlet condition to the average velocity inlet conditon as follow:
\frac{1}{|\Gamma_{inlet}|}\int_{y=1}\mathbf{u}\cdot\mathbf{n}ds=-1, \mathbf{u}\cdot\mathbf{t}=0
In my case, this can be reformula as follow:
\int_{y=1}\mathbf{u}.y ds=-1, \mathbf{u}.x=0 at y=1.
How should I do? Here is my code
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, la, plot
from dolfinx.fem import Constant, Function, dirichletbc, extract_function_spaces, form, functionspace, locate_dofs_topological
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary, meshtags, locate_entities
from ufl import div, dx, grad, inner, dot, FacetNormal, Measure, ds
# Create mesh
msh = create_rectangle(
MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [32, 32], CellType.triangle
)
def noslip_boundary(x):
return np.isclose(x[1], 0.0)
# Function to mark the lid (y = 1)
def inlet_boundary(x):
return np.isclose(x[1], 1.0)
# Lid velocity
def inlet_velocity_expression(x):
return np.stack((np.zeros(x.shape[1]), -np.ones(x.shape[1])))
def symmetry_boundary(x):
return np.isclose(x[0], 0.0)
def outlet_boundary(x):
return np.isclose(x[0], 1.0)
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)
def mixed_direct():
# Create the Taylot-Hood function space
TH = mixed_element([P2, P1])
W = functionspace(msh, TH)
W0 = W.sub(0)
V, _ = W0.collapse()
W1 = W.sub(1)
Q, _ = W1.collapse()
noslip = Function(V)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
dofs = locate_dofs_topological((W0, V), 1, facets)
bc0 = dirichletbc(noslip, dofs, W0)
inlet_velocity = Function(V)
inlet_velocity.interpolate(inlet_velocity_expression)
facets = locate_entities_boundary(msh, 1, inlet_boundary)
dofs = locate_dofs_topological((W0, V), 1, facets)
bc1 = dirichletbc(inlet_velocity, dofs, W0)
W00 = W0.sub(0)
V0, _ = W00.collapse()
symmetry = Function(V0)
facets = locate_entities_boundary(msh, 1, symmetry_boundary)
dofs = locate_dofs_topological((W00, V0), 1, facets)
bc2 = dirichletbc(symmetry, dofs, W00)
outlet_pressure = Function(Q)
facets = locate_entities_boundary(msh, 1, outlet_boundary)
dofs = locate_dofs_topological((W1, Q), 1, facets)
bc3 = dirichletbc(outlet_pressure, dofs, W1)
# Collect Dirichlet boundary conditions
bcs = [bc0, bc1, bc2, bc3]
# Define variational problem
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = Function(V)
a = form((inner(grad(u), grad(v)) - inner(p, div(v)) + inner(div(u), q)) * dx)
L = form(inner(f, v) * dx)
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)
# Set Dirichlet boundary condition values in the RHS
fem.petsc.set_bc(b, bcs)
# Create and configure solver
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
# Configure MUMPS to handle pressure nullspace
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")
pc.setFactorSetUpSolverType()
pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)
# Compute the solution
U = Function(W)
ksp.solve(b, U.x.petsc_vec)
mixed_direct()
Please state what version of DOLFINx you are using. This can be obtained with python3 -c "import dolfinx; print(dolfinx.__version__)".
Consider the following MWE:
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element, mixed_element
import dolfinx.fem.petsc
def noslip_boundary(x):
return np.isclose(x[1], 0.0)
# Function to mark the lid (y = 1)
def inlet_boundary(x):
return np.isclose(x[1], 1.0)
# Lid velocity
def inlet_velocity_expression(x):
return np.stack((np.zeros(x.shape[1]), -np.ones(x.shape[1])))
def symmetry_boundary(x):
return np.isclose(x[0], 0.0)
def outlet_boundary(x):
return np.isclose(x[0], 1.0)
def tangential_projection(u: ufl.Coefficient, n: ufl.FacetNormal) -> ufl.Coefficient:
"""
See for instance:
https://link.springer.com/content/pdf/10.1023/A:1022235512626.pdf
"""
return (ufl.Identity(u.ufl_shape[0]) - ufl.outer(n, n)) * u
def facet_vector_approximation(V: dolfinx.fem.FunctionSpace,
mt: dolfinx.mesh.MeshTags | None = None,
mt_id: int | None = None,
tangent: bool = False,
interior: bool = False,
jit_options: dict | None = None,
form_compiler_options: dict | None = None) -> dolfinx.fem.Function:
"""
Based on: https://gist.github.com/hherlyng/cb3ab37dc58205bcecbc9a6e15b1267f
by @hherlyng with a modification to avoid division by zero in norm computation
Approximate the facet normal or tangent on a surface by projecting it into the function space for a set of facets.
Args:
V: The function space to project into. Needs to be defined on discontinuous elements.
mt: The `dolfinx.mesh.MeshTags` containing facet markers.
mt_id: The id for the facets in `mt` we want to represent the normal or tangent on.
tangent: To approximate the tangent of the facets set this flag to `True`.
jit_options: Parameters used in CFFI JIT compilation of C code generated by FFCx.
See https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/jit.py#L22-L37
for all available parameters. Takes priority over all other parameter values.
form_compiler_options: Parameters used in FFCx compilation of this form. Run `ffcx - -help` at
the commandline to see all available options. Takes priority over all
other parameter values, except for `scalar_type` which is determined by
DOLFINx.
"""
jit_options = jit_options if jit_options is not None else {}
form_compiler_options = form_compiler_options if form_compiler_options is not None else {}
comm = V.mesh.comm # MPI Communicator
n = ufl.FacetNormal(V.mesh) # UFL representation of mesh facet normal
u, v = ufl.TrialFunction(V), ufl.TestFunction(V) # Trial and test functions
# Define integral measure
if interior:
# Create interior facet integral measure
dS = ufl.dS(domain=V.mesh) if mt is None else ufl.dS(domain=V.mesh, subdomain_data=mt, subdomain_id=mt_id)
else:
# Create exterior facet integral measure
ds = ufl.ds(domain=V.mesh) if mt is None else ufl.ds(domain=V.mesh, subdomain_data=mt, subdomain_id=mt_id)
# If tangent==True, the right-hand side of the problem should be a tangential projection of the facet normal vector.
if tangent:
if V.mesh.geometry.dim == 1:
raise ValueError("Tangent not defined for 1D problem")
elif V.mesh.geometry.dim == 2:
if interior:
a = (ufl.inner(u('+'), v('+')) + ufl.inner(u('-'), v('-'))) * dS
L = ufl.inner(ufl.as_vector([-n('+')[1], n('+')[0]]), v('+')) * dS \
+ ufl.inner(ufl.as_vector([-n('-')[1], n('-')[0]]), v('-')) * dS
else:
a = ufl.inner(u, v) * ds
L = ufl.inner(ufl.as_vector([-n[1], n[0]]), v) * ds
else:
c = dolfinx.fem.Constant(V.mesh, (1.0, 1.0, 1.0)) # Vector to tangentially project the facet normal vectors on
if interior:
a = (ufl.inner(u('+'), v('+')) + ufl.inner(u('-'), v('-'))) * dS
L = ufl.inner(tangential_projection(c, n('+')), v('+')) * dS \
+ ufl.inner(tangential_projection(c, n('-')), v('-')) * dS
else:
a = ufl.inner(u, v) * ds
L = ufl.inner(tangential_projection(c, n), v) * ds
# If tangent==false the right-hand side is simply the facet normal vector.
else:
if interior:
a = (ufl.inner(u('+'), v('+')) + ufl.inner(u('-'), v('-'))) * dS
L = (ufl.inner(n('+'), v('+')) + ufl.inner(n('-'), v('-'))) * dS
else:
a = ufl.inner(u, v) * ds
L = ufl.inner(n, v) * ds
# Find all boundary dofs, which are the dofs where we want to solve for the facet vector approximation.
# Start by assembling test functions integrated over the boundary integral measure.
ones = dolfinx.fem.Constant(V.mesh, dolfinx.default_scalar_type((1,) * V.mesh.geometry.dim)) # A vector of ones
if interior:
local_val = dolfinx.fem.form((ufl.dot(ones, v('+')) + ufl.dot(ones, v('-')))*dS)
else:
local_val = dolfinx.fem.form(ufl.dot(ones, ufl.TestFunction(V))*ds)
local_vec = dolfinx.fem.assemble_vector(local_val)
# For the dofs that do not lie on the boundary of the mesh the assembled vector has value zero.
# Extract these dofs and use them to deactivate the corresponding block in the linear system we will solve.
bdry_dofs_zero_val = np.flatnonzero(np.isclose(local_vec.array, 0))
deac_blocks = np.unique(bdry_dofs_zero_val // V.dofmap.bs).astype(np.int32)
# Create sparsity pattern by manipulating the blocks to be deactivated and set
# a zero Dirichlet boundary condition for these dofs.
bilinear_form = dolfinx.fem.form(a, jit_options=jit_options,
form_compiler_options=form_compiler_options)
pattern = dolfinx.fem.create_sparsity_pattern(bilinear_form)
pattern.insert_diagonal(deac_blocks)
pattern.finalize()
u_0 = dolfinx.fem.Function(V)
u_0.vector.set(0)
bc_deac = dolfinx.fem.dirichletbc(u_0, deac_blocks)
# Create the matrix
A = dolfinx.cpp.la.petsc.create_matrix(comm, pattern)
A.zeroEntries()
# Assemble the matrix with all entries
form_coeffs = dolfinx.cpp.fem.pack_coefficients(bilinear_form._cpp_object)
form_consts = dolfinx.cpp.fem.pack_constants(bilinear_form._cpp_object)
dolfinx.fem.petsc.assemble_matrix(A, bilinear_form, constants=form_consts, coeffs=form_coeffs, bcs=[bc_deac])
# Insert the diagonal with the deactivated blocks.
if bilinear_form.function_spaces[0] is bilinear_form.function_spaces[1]:
A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)
A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH)
dolfinx.cpp.fem.petsc.insert_diagonal(A=A, V=bilinear_form.function_spaces[0], bcs=[bc_deac._cpp_object], diagonal=1.0)
A.assemble()
# Assemble the linear form and the right-hand side vector.
linear_form = dolfinx.fem.form(L, jit_options=jit_options,
form_compiler_options=form_compiler_options)
b = dolfinx.fem.petsc.assemble_vector(linear_form)
# Apply lifting to the right-hand side vector and set boundary conditions.
dolfinx.fem.petsc.apply_lifting(b, [bilinear_form], [[bc_deac]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, [bc_deac])
# Setup a linear solver using the Conjugate Gradient method.
solver = PETSc.KSP().create(MPI.COMM_WORLD)
solver.setType("cg")
solver.rtol = 1e-8
solver.setOperators(A)
# Solve the linear system and perform ghost update.
nh = dolfinx.fem.Function(V) # Function for the facet vector approximation
solver.solve(b, nh.vector)
nh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
# Normalize the vectors to get the unit facet normal/tangent vector.
nh_norm = ufl.sqrt(ufl.inner(nh, nh)) # Norm of facet vector
cond_norm = ufl.conditional(ufl.gt(nh_norm, 1e-10), nh_norm, 1.0) # Avoid division by zero
if V.mesh.geometry.dim == 1:
nh_norm_vec = ufl.as_vector((nh[0]/cond_norm))
elif V.mesh.geometry.dim == 2:
nh_norm_vec = ufl.as_vector((nh[0]/cond_norm, nh[1]/cond_norm))
elif V.mesh.geometry.dim == 3:
nh_norm_vec = ufl.as_vector((nh[0]/cond_norm, nh[1]/cond_norm, nh[2]/cond_norm))
nh_normalized = dolfinx.fem.Expression(nh_norm_vec, V.element.interpolation_points())
n_out = dolfinx.fem.Function(V)
n_out.interpolate(nh_normalized)
return n_out
def mixed_direct():
# Create mesh
msh = dolfinx.mesh.create_rectangle(
MPI.COMM_WORLD, [np.array([0, 0]), np.array([3, 1])], [32, 32], dolfinx.mesh.CellType.triangle
)
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)
TH = mixed_element([P2, P1])
W = dolfinx.fem.functionspace(msh, TH)
W0 = W.sub(0)
V, _ = W0.collapse()
W1 = W.sub(1)
Q, _ = W1.collapse()
facets = dolfinx.mesh.locate_entities_boundary(msh, 1, inlet_boundary)
ft = dolfinx.mesh.meshtags(msh, 1, facets, np.full_like(facets, 7))
area_form = dolfinx.fem.form(dolfinx.fem.Constant(msh, dolfinx.default_scalar_type(1))*ufl.ds(domain=msh, subdomain_data=ft, subdomain_id=7))
inlet_area = msh.comm.allreduce(dolfinx.fem.assemble_scalar(area_form), op=MPI.SUM)
nh = facet_vector_approximation(V, mt=ft, mt_id=7, tangent=False, interior=False)
inlet_velocity = dolfinx.fem.Function(V)
u_val = dolfinx.fem.Constant(msh, dolfinx.default_scalar_type(1))
area = dolfinx.fem.Constant(msh, inlet_area)
facet_expr = dolfinx.fem.Expression(-u_val/area*nh, V.element.interpolation_points())
boundary_cells = dolfinx.mesh.compute_incident_entities(msh.topology, facets, 1, msh.topology.dim)
inlet_velocity.interpolate(facet_expr, cells=boundary_cells)
inlet_velocity.x.scatter_forward()
with dolfinx.io.VTXWriter(msh.comm, "u.bp", [inlet_velocity], engine="BP5") as bp:
bp.write(0.0)
if __name__ == "__main__":
mixed_direct()
yielding:
Hi @dokken, I have tested it on both 0.7.0 and 0.7.3
Thank you for your help!
If I understand it correctly, your code seems to impose the condition \mathbf{u}\cdot\mathbf{n} = -1 at the inlet boundary.
But I don’t think the condition is equivalent to \frac{1}{|\Gamma_{inlet}}|\int_{\Gamma_{inlet}}\mathbf{u}\cdot\mathbf{n}ds=-1. The inlet condition in my case also called “full developed flow”.
Let the domain \Omega is [0, 1]\times[0, 1], if the no-slip condition (\mathbf{u}=[0,0]^T) is imposed at \{(x, y):x=0\}\cup\{(x, y): x=1\}, and the inlet boundary at \{(x, y): y=1\}.
If impose the full developed flow condition, the distribution of inlet velocity is unimodal, being high in the middle and low on two sides. But if using your code, the boundary condition is equivalent to \mathbf{u}=[0, -1]^T, which means the distribution of the inlet velocity is uniform.
Then you haven’t properly specified your bc in your original post (there is no way for me to interpet that you wanted a developed flow).
The bc you are describing would go together with a uniform pressure BC (where i think you would have to use a Lagrange multiplier formulation).
Alternatively, if you know what kind of quadratic profile you would like to enforce, you can easily do a similar calculation as I have done above with the exact profile.
I.e determine C such that
1/\vert\Gamma\vert\int_{\Gamma} C (1-x)(x)~\mathrm{d}s=1