Errors on mesh with pyramid elements

General
1

Hi,

I am using Pointwise to generate fully 3D mesh. The mesh can be divided into two blocks, one is boundary layer mesh which is extruded from 2D, the other one, which in order to avoid the high aspect ratio in the farfield, is filled with 3D Tet mesh within surfaces. So inevitably, the pyramid type of mesh will appear on the boundary region of two blocks. I exported it of polynomial order of four to have better surface mesh with gmsh kernel. However when transferring to pyfr mesh, it gave me error:

ValueError: Pyramids with non-parallelogram bases are currently unsupported

I have talked with Pointwise support to decide make those pyramids a bit farfield to avoid being curved. But still not working. In the end, I changed the source code tolerance a bit:

    def _check_pyr_parallelogram(self, foeles):
        # Find PyFR node map for the quad face
        fnmap = self._petype_fnmap['pyr']['quad'][0]
        pfnmap = [self._nodemaps['quad', 4][i] for i in fnmap]

        # Face nodes
        fpts = self._nodepts[foeles[:, pfnmap]].swapaxes(0, 1)

        # Check if parallelogram or not
        if np.any(np.abs(fpts[0] - fpts[1] - fpts[2] + fpts[3]) > 1e-10):   # changed here from 1e-10 to 1e-8
            raise ValueError('Pyramids with non-parallelogram bases are '
                             'currently unsupported')

I changed tolerance to 10e-8 which can be compiled (POINTWISE out come is double precision). And I run the case, with terrible outcome:

z_cut
z_cut1920×983 42.7 KB

Those pyramids (the surfaces between pyramid and prism mesh to be exactly) will generate some fluctuation will blow up the code or cause very small dt and eventually leads to terrible flow field. how can I solve this problem properly. Is that because I changed the tolerance leads to this problem?

It is a LES case, dx+ dz+ is around 2 times of a decent DNS case.

Thanks in ahead,
Zhenyang

2

At the moment we only support pyramids whose bases are affine, ie the base has to some kind of parallelogram. This is due to a problem in the polynomial representation of a pyramid. It is well in the FEM literature that in order build a stable approximation space for a arbitrary square based pyramid, the approximation space has to be a rational polynomial space. As these kinds of basis are not supported in PyFR currently, we can’t stably handle anything other that affine based pyramids.

That being said, by changing the tolerance you are effectively ignoring this limitation. Regardless of the precision of Pointwise, if PyFR is telling you they are not sufficiently affine with the built in tolerance, then I’m afraid your pyramids aren’t sufficiently affine. It could be that you are exporting the mesh at high order and so the pyramids get slightly curved in the smoothing process?

3

Hi,

Okay, thank you. I have change the tolerance back to original. I have made a new set of mesh, magically this time no error appeared (I am pretty sure I have re-compiled the code and checked the parallelogram of those Pyramids). But the result is still very awful as showed in the picture. This will blow up the code very quickly after applying initial condition. (Those points with higher value than initial condition are in Pyramid meshes). Any suggestions about this situation?

z_cut_1
z_cut_12252×1153 85.9 KB

Best wishes,
Zhenyang

4

Can you try exporting the mesh with Q1 rather than Q4 and see if you have the same problem?

5

Yes, it is still there. Could it be some connectivity problem?

z_cut_q1
z_cut_q12252×1153 92.4 KB

6

No, if PyFR was unable to find a connecting point it would throw an error. Can you confirm for the Q1 mesh how large that affine metric is for your mesh?

7

The maximum of Pyramid of Q1 mesh should be 0.06*chord_length corresponding to the middle of aerofoil section, the minimum should be 0.02*chord_length at tailing edge. On the extruded direction. it is uniform 0.025*chord_length.

8

What I meant was for the Q1 mesh how big is this term in the above function for your mesh?

9

Oh sorry, it is the order of 10e-18.

10

@Zhenyang do you think you could send me your mesh? If possible a smaller resolution mesh that exhibits the same problem would be excellent. This will just make testing a bit quicker and easier, I can take a look and see if there is a problem. I’ll DM you with my email address.

1 Like
11

One thing that would be useful to check is if you can hold a free-stream. So set your initial conditions to be a free-stream of your choosing, set the system to Euler and take your boundary conditions to be sup-in-fa and have these boundaries also enforce the same free-stream. Then run your simulation for a few time steps and report back on your findings.

Regards, Freddie.

12

Hi Freddie,

The similar result appears with the free-stream setup

z_cut_euler
z_cut_euler2252×1153 148 KB

13

Are there any differences between p = 2 and p = 1? This will help us narrow down where the issue might be.

Regards, Freddie.

14

Do you mean the p degree when I output from Pointwise or in the solver setup?

Best wishes,
Zhenyang

15

The solver itself.

Regards, Freddie.

16

Actually I alway run the code with p = 4, I just changed to p = 1 and 2 and both of them can run for several time steps with freestream setup. And outcome is very normal freestream field. Then I tried the p = 1 with navier-stokes solver and previous setups and it is working! so, perhaps the problem is with high p order with this mesh type.

z_cut_p1
z_cut_p12150×1096 47.1 KB

17

Okay, there is a good chance I know where the issue is. Specifically, one of these expressions is likely wrong:

Regards, Freddie.

1 Like
18

Something does seem to be a bit off with the linear Jacobian expressions, at least to me. I am assuming x is a vector with the solution point location in reference space.

From the gmsh importer I think that the node locations should be

V = [-1,-1,-1; 
      1,-1,-1;
      1, 1,-1;
     -1, 1,-1;
      0, 0, 1];

If I then set x to be a point in the middle of this reference element, I get J=diag(0,1,1). This seems a bit off to me, can you confirm the node ordering?

Also I have used pyramids before on a number of cases, however the base has always been on a cardinal plane. This makes me think that there is maybe an issue with the off diagonal terms? Also I presume the lack of x(3) dependence is due to parallelogram base constraint?

Matlab test script 
d2r = pi/180;
thetax = 0*d2r;
thetay = 0*d2r;
thetaz = 0*d2r;
Rx = [1, 0, 0; 0, cos(thetax), -sin(thetax); 0, sin(thetax), cos(thetax)];
Ry = [cos(thetay), 0, sin(thetay); 0, 1, 0; -sin(thetay), 0, cos(thetay)];
Rz = [cos(thetax), -sin(thetax), 0; sin(thetax), cos(thetax), 0; 0, 0, 1];

V = [-1,-1,-1; 
      1,-1,-1;
      1, 1,-1;
     -1, 1,-1;
      0, 0, 1];
  
V = (Rx*Ry*Rz*V')';


scatter3(V(:,1),V(:,2),V(:,3),'.r')
xlabel('x')
ylabel('y')
zlabel('z')

x = [0, 0, 0];

J = [((1 - x(2))*V(2,1) + (-x(2) - 1)*V(3,1) + (x(2) - 1)*V(1,1) + (x(2) + 1)*V(4,1))/4,...
     ((1 - x(2))*V(2,2) + (-x(2) - 1)*V(3,2) + (x(2) - 1)*V(1,2) + (x(2) + 1)*V(4,2))/4,...
     ((1 - x(2))*V(2,3) + (-x(2) - 1)*V(3,3) + (x(2) - 1)*V(1,3) + (x(2) + 1)*V(4,3))/4;
     ((1 - x(1))*V(3,1) + (-x(1) - 1)*V(2,1) + (x(1) - 1)*V(1,1) + (x(1) + 1)*V(4,1))/4,...
     ((1 - x(1))*V(3,2) + (-x(1) - 1)*V(2,2) + (x(1) - 1)*V(1,2) + (x(1) + 1)*V(4,2))/4,...
     ((1 - x(1))*V(3,3) + (-x(1) - 1)*V(2,3) + (x(1) - 1)*V(1,3) + (x(1) + 1)*V(4,3))/4;
     (-V(1,1) - V(2,1) - V(3,1) - V(4,1) + 4*V(5,1))/8,...
     (-V(1,2) - V(2,2) - V(3,2) - V(4,2) + 4*V(5,2))/8,...
     (-V(1,3) - V(2,3) - V(3,3) - V(4,3) + 4*V(5,3))/8];
19

Standard element on a pyramid is:

In [1]: from pyfr.shapes import PyrShape

In [2]: PyrShape.std_ele(1)
Out[2]:
[(-1.0, -1.0, -1.0),
(1.0, -1.0, -1.0),
(-1.0, 1.0, -1.0),
(1.0, 1.0, -1.0),
(0.0, 0.0, 1.0)]
20

If the issue is related to the Jacobian expression of linear elements, it may be possible to force PyFR to import all elements as non-linear and then tests whether free stream is preserved or not. It could also be an issue of Pointwise when exporting pyramids. Is this mesh available publicly?

1 Like