# The application of ds() and assemble() to calculate the drag and lift coefficients of the flow around a cylinder

I am trying to compute the drag and lift coefficients for the flow around a cylinder, but unfortunately get stranded by this ArityMismatch error, can anyone tell me how should I modify the code? ( The body part of the code is almost the same as the example given in the fenics tutorial)

Thanks in advance for any generous help!

``````def compute_drag_lift_coefficients(u,p):
# Define normal vector along the cylinder surface
n = FacetNormal(mesh)
stress_tensor=sigma(U,p_n)

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)

class CylinderBoundary(SubDomain):
def inside(self,x,on_boundary):
tol=1E-14
return on_boundary and x[0]>0.1 and x[0]<0.3 and x[1] >0.1 and x[1] <0.3
Cy=MeshFunction("size_t",mesh,mesh.topology().dim()-1)
Cy.set_all(0)

Gamma_1=CylinderBoundary()
Gamma_1.mark(boundary_parts,1)
ds=Measure('ds',domain=mesh,subdomain_id=1)

force = dot(stress_tensor,n)
drag_force=assemble(force[0]*ds)
lift_force=assemble(force[1]*ds)
# Compute drag and lift coefficients
drag_coefficient = 2 * drag_force / (rho * 1.0 * 0.41 * 2.2)
lift_coefficient = 2 * lift_force / (rho * 1.0 * 0.41 * 2.2)

return drag_coefficient, lift_coefficient

drag_coefficient, lift_coefficient = compute_drag_lift_coefficients(u_,p_)
print("drag coefficient:",drag_coefficient,"lift coefficient",lift_coefficient)
``````

The error is like:

``````ArityMismatch                             Traceback (most recent call last)
Cell In[18], line 31
27     lift_coefficient = 2 * lift_force / (rho * 1.0 * 0.41 * 2.2)
29     return drag_coefficient, lift_coefficient
---> 31 drag_coefficient, lift_coefficient = compute_drag_lift_coefficients(u_,p_)
32 print("drag coefficient:",drag_coefficient,"lift coefficient",lift_coefficient)

Cell In[18], line 23, in compute_drag_lift_coefficients(u, p)
18 ds=Measure('ds',domain=mesh,subdomain_id=1)
22 force = dot(stress_tensor,n)
---> 23 drag_force=assemble((force[0]*ds))
24 lift_force=assemble(force[1]*ds)
25 # Compute drag and lift coefficients

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/fem/assembling.py:198, in assemble(form, tensor, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, backend)
196     dolfin_form = form
197 else:
--> 198     dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
200 # Create tensor
201 comm = dolfin_form.mesh().mpi_comm()

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/fem/assembling.py:56, in _create_dolfin_form(form, form_compiler_parameters, function_spaces)
54     return form
55 elif isinstance(form, ufl.Form):
---> 56     return Form(form,
57                 form_compiler_parameters=form_compiler_parameters,
58                 function_spaces=function_spaces)
59 else:
60     raise TypeError("Invalid form type %s" % (type(form),))

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/fem/form.py:43, in Form.__init__(self, form, **kwargs)
40 if form_compiler_parameters is None:
41     form_compiler_parameters = {"external_include_dirs": dolfin_pc["include_dirs"]}
---> 43 ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
44                    mpi_comm=mesh.mpi_comm())
45 ufc_form = cpp.fem.make_ufc_form(ufc_form[0])
47 function_spaces = kwargs.get("function_spaces")

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/jit/jit.py:47, in mpi_jit_decorator.<locals>.mpi_jit(*args, **kwargs)
45 # Just call JIT compiler when running in serial
46 if MPI.size(mpi_comm) == 1:
---> 47     return local_jit(*args, **kwargs)
49 # Default status (0 == ok, 1 == fail)
50 status = 0

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/jit/jit.py:97, in ffc_jit(ufl_form, form_compiler_parameters)
95 p.update(dict(parameters["form_compiler"]))
96 p.update(form_compiler_parameters or {})
---> 97 return ffc.jit(ufl_form, parameters=p)

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ffc/jitcompiler.py:217, in jit(ufl_object, parameters, indirect)
214 kind, module_name = compute_jit_prefix(ufl_object, parameters)
216 # Inspect cache and generate+build if necessary
--> 217 module = jit_build(ufl_object, module_name, parameters)
219 # Raise exception on failure to build or import module
220 if module is None:
221     # TODO: To communicate directory name here, need dijitso params to call
222     #fail_dir = dijitso.cache.create_fail_dir_path(signature, dijitso_cache_params)

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ffc/jitcompiler.py:130, in jit_build(ufl_object, module_name, parameters)
123 params = dijitso.validate_params({
124     "cache": cache_params,
125     "build": build_params,
126     "generator": parameters,  # ffc parameters, just passed on to jit_generate
127 })
129 # Carry out jit compilation, calling jit_generate only if needed
--> 130 module, signature = dijitso.jit(jitable=ufl_object,
131                                 name=module_name,
132                                 params=params,
133                                 generate=jit_generate)
134 return module

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dijitso/jit.py:165, in jit(jitable, name, params, generate, send, receive, wait)
163 elif generate:
164     # 1) Generate source code
--> 165     header, source, dependencies = generate(jitable, name, signature, params["generator"])
166     # Ensure we got unicode from generate

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ffc/jitcompiler.py:65, in jit_generate(ufl_object, module_name, signature, parameters)
62 elif isinstance(ufl_object, ufl.Mesh):
63     compile_object = compile_coordinate_mapping
---> 65 code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
66         prefix=module_name, parameters=parameters, jit=True)
68 # Jit compile dependent objects separately,
69 # but pass indirect=True to skip instantiating objects.
70 # (this is done in here such that it's only triggered
71 # if parent jit module is missing, and it's done after
72 # compile_object because some misformed ufl objects may
73 # require analysis to determine (looking at you Expression...))
74 dependencies = []

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ffc/compiler.py:142, in compile_form(forms, object_names, prefix, parameters, jit)
139 def compile_form(forms, object_names=None,
140                  prefix="Form", parameters=None, jit=False):
141     """This function generates UFC code for a given UFL form or list of UFL forms."""
--> 142     return compile_ufl_objects(forms, "form", object_names,
143                                prefix, parameters, jit)

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ffc/compiler.py:185, in compile_ufl_objects(ufl_objects, kind, object_names, prefix, parameters, jit)
183 # Stage 1: analysis
184 cpu_time = time()
--> 185 analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
186 _print_timing(1, time() - cpu_time)
188 # Stage 2: intermediate representation

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ffc/analysis.py:89, in analyze_ufl_objects(ufl_objects, kind, parameters)
86 forms = ufl_objects
88 # Analyze forms
---> 89 form_datas = tuple(_analyze_form(form, parameters)
90                    for form in forms)
92 # Extract unique elements accross all forms
93 for form_data in form_datas:

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ffc/analysis.py:89, in <genexpr>(.0)
86 forms = ufl_objects
88 # Analyze forms
---> 89 form_datas = tuple(_analyze_form(form, parameters)
90                    for form in forms)
92 # Extract unique elements accross all forms
93 for form_data in form_datas:

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ffc/analysis.py:169, in _analyze_form(form, parameters)
162 if r == "uflacs":
163     # Temporary workaround to let uflacs have a different
164     # preprocessing pipeline than the legacy quadrature
165     # representation. This approach imposes a limitation that,
166     # e.g. uflacs and qudrature, representations cannot be mixed
167     # in the same form.
168     from ufl.classes import Jacobian
--> 169     form_data = compute_form_data(form,
170                                   do_apply_function_pullbacks=True,
171                                   do_apply_integral_scaling=True,
172                                   do_apply_geometry_lowering=True,
173                                   preserve_geometry_types=(Jacobian,),
174                                   do_apply_restrictions=True)
175 elif r == "tsfc":
176     try:
177         # TSFC provides compute_form_data wrapper using correct
178         # kwargs

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/algorithms/compute_form_data.py:418, in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
415 if not complex_mode:
416     preprocessed_form = remove_complex_nodes(preprocessed_form)
--> 418 check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)  # Currently testing how fast this is
420 # TODO: This member is used by unit tests, change the tests to
421 # remove this!
422 self.preprocessed_form = preprocessed_form

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/algorithms/check_arities.py:177, in check_form_arity(form, arguments, complex_mode)
175 def check_form_arity(form, arguments, complex_mode=False):
176     for itg in form.integrals():
--> 177         check_integrand_arity(itg.integrand(), arguments, complex_mode)

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/algorithms/check_arities.py:159, in check_integrand_arity(expr, arguments, complex_mode)
156 arguments = tuple(sorted(set(arguments),
157                          key=lambda x: (x.number(), x.part())))
158 rules = ArityChecker(arguments)
--> 159 arg_tuples = map_expr_dag(rules, expr, compress=False)
160 args = tuple(a[0] for a in arg_tuples)
161 if args != arguments:

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/corealg/map_dag.py:37, in map_expr_dag(function, expression, compress)
28 def map_expr_dag(function, expression, compress=True):
29     """Apply a function to each subexpression node in an expression DAG.
30
31     If *compress* is ``True`` (default) the output object from
(...)
35     Return the result of the final function call.
36     """
---> 37     result, = map_expr_dags(function, [expression], compress=compress)
38     return result

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/corealg/map_dag.py:86, in map_expr_dags(function, expressions, compress)
84     r = handlers[v._ufl_typecode_](v)
85 else:
---> 86     r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
88 # Optionally check if r is in rcache, a memory optimization
89 # to be able to keep representation of result compact
90 if compress:

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/algorithms/check_arities.py:48, in ArityChecker.sum(self, o, a, b)
46 def sum(self, o, a, b):
47     if a != b:
---> 48         raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
49     return a

ArityMismatch: Adding expressions with non-matching form arguments ('v_1',) vs ().

``````

Please make a minimal working/reproducible example, as your function depends on input that you haven’t supplied

Thanks for the reminder!
The input is like:

``````from fenics import *
from mshr import *
import numpy as np

T=5
num_steps=5001
dt=T/num_steps
mu=0.001
rho=1

#Create mesh
channel = Rectangle(Point(0,0),Point(2.2,0.41))
cylinder=Circle(Point(0.2,0.2),0.05)
domain=channel-cylinder
mesh=generate_mesh(domain,64)

# Function Space
V=VectorFunctionSpace(mesh,'P',2)
Q=FunctionSpace(mesh,'P',1)

#Define boundaries
inflow = 'near(x[0],0)'
outflow = 'near(x[0],2.2)'
walls = 'near(x[1],0)||near(x[1],0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1] >0.1 && x[1] <0.3'

# Inflow profile
inflow_profile = ('4.0*1.5*x[1]*(0.41-x[1])/pow(0.41,2)','0')

# BCs
bcu_inflow=DirichletBC(V,Expression(inflow_profile,degree=2),inflow)
bcu_walls=DirichletBC(V,Constant((0,0)),walls)
bcu_cylinder=DirichletBC(V,Constant((0,0)),cylinder)
bcp_outflow=DirichletBC(Q,Constant(0),outflow)
bcu=[bcu_inflow, bcu_walls, bcu_cylinder]
bcp=[bcp_outflow]

#Trial and Test functions
u=TrialFunction(V)
v=TestFunction(V)
p=TrialFunction(Q)
q=TestFunction(Q)

# Functions for solutions at previous and current time steps
u_n=Function(V)
u_=Function(V)
p_n=Function(Q)
p_=Function(Q)

# Expressions used in variational forms
U=0.5*(u_n+u)
n=FacetNormal(mesh)
f=Constant((0,0))
k=Constant(dt)
mu=Constant(mu)

def epsilon(u):

# Stress tensor
def sigma(u,p):
return  2*mu*epsilon(u)-p*Identity(len(u))

# Variational problem for step 1
F1= rho*dot((u-u_n)/k,v)*dx\
+inner(sigma(U,p_n),epsilon(v))*dx\
-dot(f,v)*dx
a1=lhs(F1)
L1=rhs(F1)
# Variational problem for step 2
# Variational problem for step 3
a3=dot(u,v)*dx

# Assemble matrices
A1=assemble(a1)
A2=assemble(a2)
A3=assemble(a3)

#Apply bcs to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for visualization output
xdmffile_u =XDMFFile('NS_cylinder/velocity.xdmf')
xdmffile_p =XDMFFile('NS_cylinder/pressure.xdmf')

# Create time series (for use in reaction_system.py)
timeseries_u=TimeSeries('NS_cylinder/velocity_series')
timeseries_p=TimeSeries('NS_cylinder/pressure_series')

# Save mesh to file
File('NS_cylinder/cylinder.xml.gz') << mesh

# Create progress bar
# progress=Progress('Time-stepping')
# set_log_level(PROGRESS)

#Time-stepping
t=0

for n in range(num_steps):
t+=dt
#1 Tentative velocity step
b1=assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

#2 pressure correction step
b2=assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg' )

#3 Velocity correction step
b3=assemble(L3)
solve(A3, u_.vector(), b3, 'cg','sor')

xdmffile_u.write(u_,t)
xdmffile_p.write(p_,t)

timeseries_u.store(u_.vector(),t)
timeseries_p.store(p_.vector(),t)

u_n.assign(u_)
p_n.assign(p_)
``````

Why is the input non-captialized, while the `U` in sigma is capitalized?

Thank you! Now it’s working properly.

Hi, Dr. Dokken! I am having another issue with my code. Initially I accidentally integrated over the whole domain and get some values of the drag and lift coefficients, but I later found out that when I integrate near the boundary of the cylingder, the drag and lift coefficients are always 0 no matter the timestep.
My code now can be found here: Getting 0 drag coefficient and lift coefficient when working on the flow over a cylinder - General - FEniCS Project.

I would be more than grateful if you could help me again!!