I’m trying to do LESes of flows past a large number (~100) of different bluff bodies using the acnavierstokes
solver at around Re=400. My BCs are periodic in the zdirection (the bluff bodies are extruded along the zdirection), slip BCs above and below the object, acinfv
at the inlet with (u, v, w) = (1.0, 0.0, 0.0) and acoutfp
at the outlet with p=1.0. My meshes have ~70k elements and I use the 3rd order methods (tried 4th as well but those NaN even quicker)
However, I cannot seem to find a good set of solver settings, as the settings I explored so far either lead to NaNs about half way into the simulation, or greatly baloon the runtime. Here’s what I tried so far:
 A ‘standard’ setup with the RK45 integrator. This takes a while to ‘settle’ in the beginning, and goes quicker once the residuals die down but still takes ~24hrs on really beefy hardware like two A100 GPUs
[solvertimeintegrator]
formulation = dual
scheme = sdirk33
pseudoscheme = rk45
controller = none
pseudocontroller = localpi
tstart = 0.0
tend = 60.0
dt = 5e3
pseudodt = 2.5e4
pseudonitersmin = 1
pseudonitersmax = 250
pseudoresidnorm = l2
pseudoresidtol = 5e4
pseudoresidtolp = 2.5e2
atol = 1e1
pseudodtmaxmult = 2.5
[solverdualtimeintegratormultip]
pseudodtfact = 1.75
cycle = [(3, 1), (2, 2), (1, 4), (0, 8), (1, 4), (2, 2), (3, 1)]

I noticed that things go a lot quicker if I raise the dt/pseudodt ratio to ~40, especially with the
none
pseudocontroller. Predicted runtime was about 34 hours, especially if I up the # of iterations in the lower p levels of the multigrid method. However, this would cause NaNs about 2030 seconds of physical time into the simulation, which corresponds to the onset of instability in the zdirection I believe. In the snapshots saved just before the NaN occurred, I noticed that the pressure in some cells spiked up to extreme levels. (I also experimented with using other pseudoschemes likevermeire
andtvdrk3
but these did not really change anything) 
If I lowered the dt/pseudodt ratio with the
none
pseudocontroller, the simulation would go very quickly until, again, 2030 seconds, but after that it’d slow down to a crawl, requiring the max number of iterations in every time step to progress. Going by the progress bar, this’d mean the simulation would take days even at this relatively modest Re. 
I went back to the
localpi
pseudo controller with RK45, increased the ratio of dt/pseudodt to 40, and also massively increased the # of allowed iterations per timestep to 10000. This of course was glacial in the beginning, though it seemed to go very quickly later on. However, like the previous case, the runtime ballooned to over a day at around the 20 second mark.
Finally, these trends seem to persist even when I try to change the value of aczeta
; for instance, the slowdown with attempt #3 above at around t=20 occurs with both large (812) and small (24) values of aczeta
.
So, at least within the space of settings I tried, there seems to be only two possibilities: things go fast but eventually NaN, or they take a very long time for the Reynolds #, and it all seems to be connected to the onset of instability along the zdirection at around the 20 second mark. Is there a suggested way to make things go faster while not compromising numerical stability? I know PyFR’s numerical scheme can be sensitive to solver settings, so I’d appreciate any help.
Here are some images of my mesh. I can share my mesh file if needed.
(As a final note, these observations are all from runs using single precision  I tried double precision too, but in the cases of settings getting NaNs it did not fix the issue, and in the cases of settings that took too long it simply doubled the runtime)