Hi, I am trying to extract vorticity from the velocity field for fluid simulation. This is my code and extracted field. I am wondering why the vorticity field is so large that it ranges from -160 to 160. Am I doing the computation wrongly?
// We have a known mesh in the file fname
mesh = gen_mesh_from_file(fnames)
nodes = mesh.coordinates()
cell = mesh.cells()
// define function space
V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)
u_n = Function(V)
// We have a vector field u_n in second order VectorFunctionSpace, and we compute vorticity from that
ux, uy = u_n.split(deepcopy=True)
vortex = ux.dx(1)-uy.dx(0)
vor = Function(Q)
vor = project(vortex,Q)
// We plot the vorticity
np_solution = vor.compute_vertex_values()
plot_fig = plt.tricontourf(nodes[:,0], nodes[:,1],cell,np_solution,levels=100,cmap=cmocean.cm.thermal)
If u_n is in a second order continuous Lagrange space, then the vorticity is in a Discontinuous Lagrange space of order 1.
Thus your projection to CG-1 is introducing an error between the actual and shown vorticity.
Please use a DG-1 space for the projection, and XDMFFile.write_checkpoint to being able to visualize the DG-1 space in Paraview, as compute_vertex_values only does a CG-1 interpolation of any function.
I am trying to simulate the fluid dynamics and get the velocity field out of it. Now I am assuming sensors are placed at the node of the mesh. To get the v-field, I am doing this:
V = VectorFunctionSpace(mesh, "Lagrange", 2)
u_1 = Function(V)
// some process to set value for u_1
ux, uy = u_1.split()
u_solution[0, :] = ux.compute_vertex_values()
u_solution[1, :] = uy.compute_vertex_values()
Do you think it is a correct way of doing this? Here I have u_1 in CG-2 and try to compute_vertex_values() to get the field value.
Another question is that with u_solution, we won’t be able to get u_1 right? As there are many dof are missing. Is it possible to get a u_solution, so that:
Why don’t you use the V.tabulate_dof_coordinates() to get the coordinates of all the degrees of freedom?
As to your second question, you are correct that you cannot recover u_1 by this method, as you have discarded all information about the higher order space when computing the vertex values. You can easily verify this by comparing the underlying arrays of u and u_1 by doing print(u.vector().get_local()-u_1.vector().get_local())
I would still suggest you use XDMFFile.write_checkpoint and Paraview to look at your vorticity.
It’s a great idea to get V.tabulate_dof_coordinates() to get the coordinates for all the degrees of freedom. Can I also get the value for those nodes? Do you know how should I get this?
One thing I can come up with is to value the field at those coordiantes. Though I am not sure if this is the correct way of doing this.
