I am running PyFR with very high order p = 6 or even higher with ILES method. But simulation blowed up but it is ok for lower polynomial order like 2 or 3. So I am thinking if it is the problem with quadrature rule and quadrature degree I chose for Pyramids elements. If I am right, I used gauss-legendre for solution points and witherden-vincent for quadrature points (I guess witherden-vincent is the only option for quadrature points right?). But could you tell me the formula to calculate the number of points that witherden-vincent need based on polynomial order? It seems it is not the same as what one can find from shape.py.

@WillT’s suggestion about making sure you reduce dt as you increase p (on a fixed mesh) is important. RE: the quadrature rules, an important point is to ensure that you have more quadrature points than solution points. What quad-deg are you using for the witherden-vincent points on pyramids? How many quadrature points does this give you? And how does this compare to the number of solution points you have at your given p?

Thanks for rapid replies. I did reduced time step and run it with adaptive time integration method. Maybe I will try even smaller time step later. But I am thinking, as @p.vincent said, one has to use more quadrature points than solution points. I am running p=5, which gives 91 gauss-lagendre solution points. But the maximum quadrature degree in the quadrature file is d10 with 83 points. So I am thinking if this is the problem.

I think that is likely the problem (or at least part of the problem). Could you try generating some higher strength rules with more points using Polyquad, and then use them?

Sorry for the late response. After increasing quadrature degree from files that Freddie sent to me, our simulation can run pretty well. Thanks again for helping.