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.A1 = inner(grad(self.u),grad(self.v))*dx
self.L1 = dot(Constant((1,0)),self.v)*dx
self.A2 = inner(grad(self.p),grad(self.q))*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!