Strange issue: no vector display when plotting

gpr · June 28, 2021, 6:15pm

Hello everyone !

I haven’t been able to find a topic that deals with this weird issue, so here it is:
I’m trying to solve a 3D Navier-Stokes problem using IPCS ( Incremental Pressure Correction Scheme ). It is the same core principle as in Fenics tutorial ft07.

Here is my method:
-First, I create a STEP file with Onshape (here, a cylinder)
-I mesh the STEP file and tag all the boundaries in GMSH and then convert it to XDMF.
-Finally I import the XDMF file in Fenics, compute Navier-Stokes solution and plot the results with Matplotlib and Paraview.

Here is my issue:
I can see no arrows to represent the velocity field in matplotlib, only points. When plotting in Paraview arrows seems to be unit vectors on the x direction, which is very unrealistic.
However, the norm values seems quite relevant: zero at the cylinder’s wall and maximum velocity near the center.

Any idea on why this is happening ?

Scripts:
(dowload the github repository)
Minimal working example:

github.com

Fra26802/discourse-post-2/blob/main/flow_cylinder_3D.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jun 22 23:04:14 2021

@author: gpr
"""

#USEFUL LIBRARIES---------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
import meshio


#USEFUL FUNCTIONS---------------------------------------------------
#tensors needed for flow equation:
def epsilon(u):
    return sym(nabla_grad(u))

def sigma(u, p):

This file has been truncated. show original

Script named that is needed to import the create_mesh function:

github.com

Fra26802/discourse-post-2/blob/main/gmsh_interface.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 21 11:53:47 2021

@author: gpr

source J.Dokken tutorial https://jsdokken.com/converted_files/tutorial_pygmsh.html
"""
import meshio
import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

.msh file needed:

github.com

Fra26802/discourse-post-2/blob/main/cylinder.msh

$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
4
2 11 "walls"
2 22 "inflow"
2 33 "outflow"
3 1 "cylinder"
$EndPhysicalNames
$Entities
2 3 3 1
1 5 3.06161699786838e-16 60 0 
2 5 3.06161699786838e-16 0 0 
1 -5.0000001 -5.0000001 59.9999999 5.0000001 5.0000001 60.0000001 0 2 1 -1 
2 4.99999 -9.999999999081516e-06 -9.999999999621423e-06 5.00001 1.000000000091849e-05 60.00001 0 2 1 -2 
3 -5.0000001 -5.0000001 -1e-07 5.0000001 5.0000001 1e-07 0 2 2 -2 
1 -5.0000001 -5.0000001 -1.000000011686097e-07 5.0000001 5.0000001 60.0000001 1 11 4 1 2 3 -2 
2 -5.0000001 -5.0000001 59.9999999 5.0000001 5.0000001 60.0000001 1 22 1 1 
3 -5.0000001 -5.0000001 -1e-07 5.0000001 5.0000001 1e-07 1 33 1 3

This file has been truncated. show original

Kei_Yamamoto · June 29, 2021, 3:07pm

Since you are only applying constant pressure gradient in x-direction, I think it is not that strange that the velocity vector only has x component. One thing I am curious is that I don’t see any arrows inside the mesh but only on the boundary parts of the mesh. How do you plot in Paraview ?

gpr · July 2, 2021, 8:11am

Hi @Kei_Yamamoto !

Thanks for your answer.Tell me if i’m wrong, but the pressure gradient is only along the z-axis, with a pression of 18.6 at the inflow and 0 at the outflow. That’s why I don’t understand why the velocity field is along the x-axis.
For plotting in paraview I just open the file velocity_field.xdmf with Xdmf3-ReaderS then select 3D glyphs and coloring : f_28

here is a screenshot:

Kei_Yamamoto · July 2, 2021, 8:48am

Okay, I misunderstood your Paraview plot and it is indeed strange. I have done similar simulation and the variational formulation looks no problem to me. One difference is that I have never used keep_diagonal=True and A1.ident_zeros(). I don’t really know what they are doing but could you tell me why they are needed ??

Also, maybe you already know it but
https://fenicsproject.org/pub/tutorial/html/._ftut1009.html
this might be helpful.

gpr · July 2, 2021, 10:25am

before building this code I made a 2D one that wasn’t working without those two lines :
keep_diagonal=True and A1.ident_zeros()

Following what you said I just tried without A1.ident_zeros() and turning keep_diagonal to False. It makes no difference.

Actually i knew about this tutorial (Solving PDEs in Python - <br> The FEniCS Tutorial Volume I) and used it to implement the IPCS scheme.

Actually what is really weird is that the two-dimensionnal simulation and plotting works fine. I think i did something wrong when switching from 2D to 3D

Kei_Yamamoto · July 5, 2021, 7:38am

Hi gpr,

I think I figured out why it does not work.

When you create triangle_mesh on line 38,39 , you do not set prune_z=true. This needs to be done to create 2D mesh that consists of boundary part.

You should modify your code as follows.

triangle_mesh = create_mesh(mesh3D_from_msh, “triangle”, prune_z=True)
meshio.write("facet_mesh3D_from_msh.xdmf", triangle_mesh)

gpr · July 19, 2021, 3:31pm

Hi @Kei_Yamamoto ! I don’t think the z coordinate needs to be pruned. We still need it for the surfaces as they have a z-coordinate (they aren’t in the x,y plane). I tried setting prune_z to true and it doesn’t make any difference. Still the idea was relevant and i thank you for your suggestion

Eventually i found the solution. The problem lied in how i applied the filters in Paraview . It didn’t had anything to do with the Fenics code.

EngineerMama · July 19, 2021, 8:43pm

@gpr
Out of curiosity, what was the filter problem & resolution (in case someone modelling ft07 runs into the same issue)?

gpr · July 20, 2021, 8:56pm

At the beginning I only selected “glyph” in the top menu without applying any filters.
It indeed displayed glyphs, but the orientation of the vector field wasn’t taken into account.
In order to display properly a vector field, one must select a glyph filter and enable scaling by magnitude. Then you have to play a bit with the scale factor for the arrows to appear.

This might seem trivial for a regular paraview user, but i was a complete beginner so i got a bit confused

Topic		Replies	Views
Navier-Stokes Secondary Flows (Dean Drag Forces) in a curved channel variational formulation	4	1668	November 21, 2020
Incompressible flow Unable to solve linear system using PETSc Krylov solver Errors	2	427	April 21, 2023
Unable to solve linear system using PETSc Krylov solver for incompressible flow simulation Errors	0	247	April 21, 2023
How to plot solutions while running fenics code in parallel? General	9	389	November 13, 2023
XDMF file causes a crash but .pvd files work General	3	446	November 30, 2022

Strange issue: no vector display when plotting

Related topics