Divergence when surface normal vector is flipped

Cases
1

I was running a 2D simulation that kept diverging over and over again (divergence occured after only a few time steps). I finally find out that it had to do with the mesh orientation in the z direction. simulation is in the xy plane but the surface normal vector was oriented in the negative z direction rather than the positive. flipping the orientation fixed the issue.

I think there should be some of warning when the orientation is wrong. or even better, an automatic fixing mechanism.

2

This is a common issue with meshes created with Gmsh. PyFR has code to detect this, here:

namely through the presence of negative Jacobian determinants.

Regards, Freddie.

3

Thank you for your answer Freddie,
You say that Pyfr has code to detect flipped mesh orientation, but in my simulations an error was not raised. do you have any idea why? or any hunch?

I am asking because I found this error very difficult to trace, as the simulation simply diverged.

thanks
Nathan

4

It is difficult to say without seeing the mesh/geo file.

Regards, Freddie.

5

Thank you freddie, you can find the mesh file here. If the .geo file is needed I can supply it as well

6

Which version of gmsh are you using?

Does your .geo file call the OpenCASCADE geometry kernel, i.e.:

SetFactory("OpenCASCADE");

PS: Yes, please add the .geo file.

7

Hi Nigel
I am using GMSH 4.8.0, with openCASCADE (I am using the GUI)
the text of the geo file is at the bottom of the message and the geometry file can be found here.
The process I went through to generate the mesh was:

  1. Set the “global model scaling” to 1e-3
  2. Merge the geo file
  3. Generate the mesh
  4. Recombine using “recombine 2d”

I am sure there is a better way to do this but this is what i did.
Geo file:

Merge "fluid_domain_scaled_supressuion-Surface-Plane1.STEP";

Field[1] = BoundaryLayer;
Field[1].CurvesList = {461,462,463,464,465,466};
Field[1].SizeFar = 0.004;
Field[1].Size = 0.00006;
Field[1].Thickness = 0.00025;
Field[1].Ratio = 1.1;
Field[1].Quads = 1; // no touch


BoundaryLayer Field = 1; // no touch

Mesh.MeshSizeMax = 0.004;
	
Field[2] = Box; // no touch
Field[2].VIn = 0.00015;
Field[2].VOut = 0.004;
Field[2].XMax = -0.335;
Field[2].XMin = -0.45;
Field[2].YMax = 0.030;
Field[2].YMin = 0;
Field[2].Thickness = 0.03; 
Background Field =2;  // no touch

Physical Curve("Stack") = {1:460};
Physical Curve("Top") = {463,464};
Physical Curve("Left") = {462};
Physical Curve("Bottom") = {461,466};
Physical Curve("Right") = {465};
Physical Surface("fluid") = {1};
1 Like
8

Looking at the .geo and .msh files (and the resulting extremely small quads) any chance the offending Jacobian determinants might have had negative values, but “less negative” than the hard-wired [-1.0e-5]… that is to say, did some have values in the range (-1e-5 : 0] ?

Which makes me wonder, should the code also trap zero-valued dets?

9

when looking at the surface orientation in GMSH, it seems the orientation is simply negative, not near zero. also as I said flipping the surface orientation fixed the issue

1 Like
10

I had a look at this and everything seems as is should in terms of the determinant, the lowest value is small but positive. And inspecting the smats, its seems to be performing as expected.

I think this issue is that in the code that finds the surface norms, _pnorm_fpts, I think there might be an assumption baked in that the mesh surface normal is in z +. And so I think, in 2D, it is possible that element have a positive jacobian, but the normals at the flux points are in the wrong direction. Does that seems possible @fdw? It should be possible to test the mesh surface normal by computing the cross product of two edge vectors and testing the sign.

11

Ok I realised I misled you all by sendind the wrong mesh. this is the mesh after the flipping. sorry!!!

this is a flipped mesh. the .geo file was correct, but the .msh file was after flipping.

sorry again and thanks for your help so far!

12

Ok, that is alot clearer, disregaurd all of what I said previously. The determinants are negative, but small enough to slip by the test.

It is probably worth looking at modifying this test, maybe scaling it by some value from the mesh. Although a completely robust way of doing this is not easy.

2 Likes
13

Picking these kinds of numbers can be tricky. I’m open to suggestions. Generally, I am suspicious of elements with small Jacobian’s as they often cause issues in terms of numerical precision.

Regards, Freddie.

14

I think that a warning couldn’t hurt, just a warning for any negative jacobian and an error above a certain value.