How to calculate vorticity in 2D flow simulation

I am simulating flow past a cylinder with FEniCS 2019.1.0 and trying to calculate the vorticity as a scalar and plot it. I would be really grateful if anyone can tell me how can I modify the code to fix the error.

This is a MWE of my problem.

``````from fenics import *

class Env2DCylinder():
def __init__(self):

self.mesh = UnitSquareMesh(10,10)
self.V = VectorFunctionSpace(self.mesh,'P',2)
self.Q = FunctionSpace(self.mesh,'P',1)
self.W = FunctionSpace(self.mesh,'CG',1)

self.u = TrialFunction(self.V)
self.v = TestFunction(self.V)
self.p = TrialFunction(self.Q)
self.q = TestFunction(self.Q)
self.u_ = Function(self.V)
self.p_ = Function(self.Q)

self.w_ = Function(self.W)

self.L1 = dot(Constant((1,0)),self.v)*dx

self.L2 = Constant(0)*self.q*dx

def evolve(self):
solve(self.A1 == self.L1, self.u_)
solve(self.A2 == self.L2, self.p_)

self.w_ = self.compute_vorticity(self.u_)
plot(self.w_,cmap='viridis')
plt.show()

self.xdmffile_u.write(self.u_,0)
self.xdmffile_p.write(self.p_,0)
self.timeseries_u.store(self.u_.vector(),0)
self.timeseries_p.store(self.p_.vector(),0)

def compute_vorticity(self,u):
mesh = self.mesh
class VorticityExpression(UserExpression):
def __init__(self,ux_value,uy_value,degree=1,mesh=mesh):
self.ux_value = ux_value
self.uy_value = uy_value
self.YY = FunctionSpace(mesh,'CG',1)
self.ux_trial = Function(self.YY)
self.ux_trial.vector().set_local(self.ux_value)
self.uy_test = Function(self.YY)
self.uy_test.vector().set_local(self.uy_value)
super().__init__(degree=degree)
def eval(self,value,x):
value[0] = self.ux_trial.dx(1)(x)-self.uy_test.dx(0)(x)

def value_shape(self):
return ()

ux,uy = u.split(deepcopy = True)
ux_values = ux.compute_vertex_values(self.mesh)
uy_values = uy.compute_vertex_values(self.mesh)
#         vortex_expr = VorticityExpression('ux.dx[1]-uy.dx[0]',degree=1,ux=ux_,uy=uy_)
vortex_expr = VorticityExpression(degree=1,ux_value=ux_values,uy_value=uy_values)
vortex = interpolate(vortex_expr,self.W)
return vortex

env = Env2DCylinder()
env.evolve()
``````

And this is the reported error:

``````AssertionError                            Traceback (most recent call last)
Cell In[5], line 66
63         return vortex
65 env = Env2DCylinder()
---> 66 env.evolve()

Cell In[5], line 30, in Env2DCylinder.evolve(self)
27 solve(self.A1 == self.L1, self.u_)
28 solve(self.A2 == self.L2, self.p_)
---> 30 self.w_ = self.compute_vorticity(self.u_)
31 plot(self.w_,cmap='viridis')
32 plt.show()

Cell In[5], line 62, in Env2DCylinder.compute_vorticity(self, u)
60 #         vortex_expr = VorticityExpression('ux.dx[1]-uy.dx[0]',degree=1,ux=ux_,uy=uy_)
61         vortex_expr = VorticityExpression(degree=1,ux_value=ux_values,uy_value=uy_values)
---> 62         vortex = interpolate(vortex_expr,self.W)
63         return vortex

File ~/anaconda3/envs/fenics1/lib/python3.8/site-packages/dolfin/fem/interpolation.py:71, in interpolate(v, V)
69 # Compute interpolation
70 if hasattr(v, "_cpp_object"):
---> 71     Pv.interpolate(v._cpp_object)
72 else:
73     Pv.interpolate(v)

File ~/anaconda3/envs/fenics1/lib/python3.8/site-packages/dolfin/function/function.py:365, in Function.interpolate(self, u)
363     self._cpp_object.interpolate(u._cpp_object)
364 else:
--> 365     self._cpp_object.interpolate(u)

File ~/anaconda3/envs/fenics1/lib/python3.8/site-packages/dolfin/function/expression.py:53, in _InterfaceExpression.__init__.<locals>.wrapped_eval(self, values, x)
52 def wrapped_eval(self, values, x):
---> 53     self.user_expression.eval(values, x)

Cell In[5], line 52, in Env2DCylinder.compute_vorticity.<locals>.VorticityExpression.eval(self, value, x)
51 def eval(self,value,x):
---> 52     value[0] = self.ux_trial.dx(1)(x)-self.uy_test.dx(0)(x)

File ~/anaconda3/envs/fenics1/lib/python3.8/site-packages/ufl/exproperators.py:341, in _call(self, arg, mapping, component)
339     return _restrict(self, arg)
340 else:
--> 341     return _eval(self, arg, mapping, component)

File ~/anaconda3/envs/fenics1/lib/python3.8/site-packages/ufl/exproperators.py:331, in _eval(self, coord, mapping, component)
329     mapping = {}
330 index_values = StackDict()
--> 331 return f.evaluate(coord, mapping, component, index_values)

File ~/anaconda3/envs/fenics1/lib/python3.8/site-packages/ufl/indexed.py:117, in Indexed.evaluate(self, x, mapping, component, index_values, derivatives)
115     return A.evaluate(x, mapping, component, index_values, derivatives)
116 else:
--> 117     return A.evaluate(x, mapping, component, index_values)

File ~/anaconda3/envs/fenics1/lib/python3.8/site-packages/ufl/differentiation.py:176, in Grad.evaluate(self, x, mapping, component, index_values, derivatives)
174 component, i = component[:-1], component[-1]
175 derivatives = derivatives + (i,)
--> 176 result = self.ufl_operands[0].evaluate(x, mapping, component,
177                                        index_values,
178                                        derivatives=derivatives)
179 return result

File ~/anaconda3/envs/fenics1/lib/python3.8/site-packages/ufl/core/terminal.py:71, in Terminal.evaluate(self, x, mapping, component, index_values, derivatives)
69 # If it has an ufl_evaluate function, call it
70 if hasattr(self, 'ufl_evaluate'):
---> 71     return self.ufl_evaluate(x, component, derivatives)
72 # Take component if any
73 warning("Couldn't map '%s' to a float, returning ufl object without evaluation." % str(self))

File ~/anaconda3/envs/fenics1/lib/python3.8/site-packages/dolfin/function/function.py:257, in Function.ufl_evaluate(self, x, component, derivatives)
254 """Function used by ufl to evaluate the Expression"""
255 # FIXME: same as dolfin.expression.Expression version. Find
256 # way to re-use.
--> 257 assert derivatives == ()   # TODO: Handle derivatives
259 if component:
260     shape = self.ufl_shape

AssertionError:

``````

Many thanks in advance for the help!

You can simply reproduce the error with:

``````from fenics import *

from dolfin import *
mesh = UnitSquareMesh(10, 10)
Y = FunctionSpace(mesh, 'CG', 1)
ux = Function(Y)
ux([0.5, 0.5])
ux.dx(0)([0.5, 0.5])
``````

and trace

``````Traceback (most recent call last):
File "/root/shared/mwe3.py", line 8, in <module>
ux.dx(0)([0.5, 0.5])
File "/usr/lib/python3/dist-packages/ufl_legacy/exproperators.py", line 330, in _call
return _eval(self, arg, mapping, component)
File "/usr/lib/python3/dist-packages/ufl_legacy/exproperators.py", line 320, in _eval
return f.evaluate(coord, mapping, component, index_values)
File "/usr/lib/python3/dist-packages/ufl_legacy/indexed.py", line 107, in evaluate
return A.evaluate(x, mapping, component, index_values)
File "/usr/lib/python3/dist-packages/ufl_legacy/differentiation.py", line 164, in evaluate
result = self.ufl_operands[0].evaluate(x, mapping, component,
File "/usr/lib/python3/dist-packages/ufl_legacy/core/terminal.py", line 60, in evaluate
return self.ufl_evaluate(x, component, derivatives)
File "/usr/lib/python3/dist-packages/dolfin/function/function.py", line 275, in ufl_evaluate
assert derivatives == ()   # TODO: Handle derivatives
AssertionError
``````

Legacy DOLFIN does not have the capabilities of evaluating the derivative of a function at a point (out of the box). This functionality has been added in DOLFINx with the introduction of `dolfinx.fem.Expression`.