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

General
1

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)
    161     error("Please provide only one of generate or receive.")
    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
    167     header = as_unicode(header)

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 ().
2

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

3

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)

# Symmetric gradient
def epsilon(u):
    return sym(nabla_grad(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\
  +rho*dot(dot(u_n,nabla_grad(u_n)),v)*dx\
  +inner(sigma(U,p_n),epsilon(v))*dx\
  +dot(p_n*n,v)*ds-dot(mu*nabla_grad(U)*n,v)*ds\
  -dot(f,v)*dx
a1=lhs(F1)
L1=rhs(F1)
# Variational problem for step 2
a2=dot(nabla_grad(p),nabla_grad(q))*dx
L2=dot(nabla_grad(p_n),nabla_grad(q))*dx-(1/k)*div(u_)*q*dx
# Variational problem for step 3
a3=dot(u,v)*dx
L3=dot(u_,v)*dx-k*dot(nabla_grad(p_-p_n),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_)
4

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

5

Thank you! Now it’s working properly.

6

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!!