Calculation of Enstrophy/Vorticity

Question

Hi,

I'm very new to FEniCS and was running the 2D NS flow around a cylinder program given in section 3.4 of the FEniCS tutorial vol1. In addition to the velocity and pressure want to plot the some quantities like vorticity or enstrophy ( 0.5*(u_x^2 + v_y^2) where the solution is U = U(u,v,t) ). I tried adding

 #calculate enstrophy
n_y = Expression(('0.0','1.0'),degree = 0)
uy = dot(grad(u_),n_y)
n_x = Expression(('1.0','0.0'),degree = 0)
ux = dot(grad(u_),n_x)
ens = dolfin.project(dot(ux,ux) + dot(uy,uy))
# Save to file (.pvd)
vfile << ens

to the time-stepping loop.

This does not cause an error and I was almost able to plot it in Paraview but not quite (only the first frame had color, values on the scale looked too large).

Any advice on how to properly calculate such quantities would be greatly appreciated.

Thanks,

Eric

Answer 1

The 'problem' you observe is probably caused by the way the data is structured in vtk rather than a mistake in your computation. To each .vtk field you produce, a label is assigned. These labels should be the same if you store multiple .vtk files of a time-dependent variable. This exactly, is where it goes wrong if you use a copy of your code in a time loop, since 'ens' is considered to be a different variable every time instance and hence, the .vtk fields you produce are labeled differently every time step.
To clarify this somewhat abstract story, consider the following minimal working example (including a work-around --> Option 2!):

from dolfin import *
nx = 20; ny = 20;
mesh = RectangleMesh(Point(-1,-1),Point(1,1),nx, ny,"crossed")
V = VectorFunctionSpace(mesh,'CG',1)
# Define a velocity field
U_exact=('-U*exp(-2*nu*pow(pi,2)*t)*(cos(pi*(x[0]))*sin(pi*(x[1])))',
         ' U*exp(-2*nu*pow(pi,2)*t)*(cos(pi*(x[1]))*sin(pi*(x[0])))')

curlvtk = Function(V)
vtkplot = File('vorticity.pvd')

# Then do some arbitrary time stepping
for t in range(10):
    u00   = Expression(U_exact, degree=3, U=float(1.), \
                              nu=float(1E-2),t=float(t))
    uh    = interpolate(u00, V)
    # Compute curl
    curlu = curl(uh)
    # Project and assign to func
    curlup = project(curlu)
    curlvtk.assign(curlup)
    # Option 1: does produce .vtk files with different labels, 
    #      reproducing your issue
    # vtkplot << curlup

    # Option 2: each .vtk file has the same label!
    vtkplot << curlvtk

Comments

This worked, thanks!!