Dear all
I am simulating an adverse pressure gradient problem using the navier-stokes solver. I started the computation at p=2 and progressively increased to p=5. At p=5, I encounter an error associated with the time step calculation.
The original setup with scheme=rk45, controller=pi, errest-norm=uniform and dt-max gives a ZeroDivisionError: 0.0 cannot be raised to a negative power. This error is at fac = err**-expa * self._errprev**expb.
When I try with a fixed dt the solution diverge to NaNs. Any changes made to the errest-norm option result in Minimum time step rejected error.
How can I resolve this issue?
Although I think there is a fix for your problem, I think this might be a symptom of an issue with your case.
The error message that you’re getting indicates that you are err goes to zero in your simulation. This is slightly strange, have you looked at the output from the dtstats plugin? Is it behaving as expected? Given what you have told me about the case you are running, it might be that your case is unstable. What or how much anti-aliasing are you using? Have you looked a the p=2 flow? Does it look reasonable?
A fix would be to add an epsilon on to err to ensure it can’t go to zero, but the error shouldn’t really be going to zero so I am reluctant to suggest this.
I do not notice anything abnormal in p=2 and p=4 solutions. I use anti-alias = flux,surf-flux.
The dtstats output is oscillating. I managed to control that to a certain degree by setting dtmax.
What is the degree of anti-aliasing you are using at p=4 and p=5? Also do you have a plot of the error and dt at p=4 and p=5?
So quad-deg for p5 hexes is set to 10? And I assume you are using GL points? If so then I think you will have less quadrature points than solution points (or maybe the same number). Which would explain what you are seeing.
If you are running p5, then you will have 6 solution points in each direction in a hex. Thus you probably want 7 quadrature points in each direction in a hex. Hence you want to set quad-deg for a hex to 13. You should also change quad-deg for the prisms, and the quad and tri faces, accordingly.
Yep, that makes sense. If you want to check out for more complicated element types what quadrature degree you might need as a minimum you can check the data files in pyfr.quadrules.*
For example here: PyFR/pyfr/quadrules/hex at develop · PyFR/PyFR · GitHub
n in the file name is the total number of points and d is the degree.
But yes try upping the anti-aliasing order. Given that your plots seem to show for p=5 the dt plummeting and the error sky rocketing, this indicates to me that the adaptive time stepper isn’t necessarily the problem. Instead the current configuration is unstable.