Perform euler case: subsonic flow around a sphere, with pyFR

General
1

Dear all,

My name is Antonio Garcia-Uceda, and I’m a researcher on Flux Reconstruction methods. I’m using pyFR as part of my project.

I’m trying to perform the euler case: subsonic flow around a sphere. You can find attached the mesh in .gmsh format, as well as a first attempt of .ini files.

Could you please help me out with the following? Thanks.

  • When running pyFR, an error pops-up asking for the following: [soln-bcs-sup_in_fa], [soln-bcs-sub_in_char], [soln-bcs-slip_wall]. The mesh attached only has two boundary faces: Slip-Wall, Subsonic Characteristic. What is the other BCs required? How pyFR detects the correct type? Perhaps I have to input the correct “physical name” in the .gmsh file? If so, which is this one? I can’t find anything related in the user guide section of pyFR website.

  • The bcs of type: [soln-bcs-sub_in_char] is not in the user guide. Could it be that I’m running a too old version of the code?

  • When runninf the .ini file attached, an error pops-up with the following: “ValueError: Invalid characters in expression”. What is wrong with it?

Many thanks in advance for your help.

Best regards,
Antonio

sphere_lev2o1_o1.ini (969 Bytes)

Sphere_lev2_o1.msh (384 KB)

2

Hi Antonio,

I’ve looked over both your PyFR configuration file and your mesh. With regards to your boundary conditions you have specified a supersonic inflow farfield boundary condition, a subsonic inflow characteristic boundary condition, and a slip wall boundary condition for your mesh. Note you have specified the type (with underscores instead of hyphens) instead of the actual name of the boundary. With PyFR you need to provide the boundary names as given in the mesh file. Then the boundary condition type is specified separately. As an example for a no-slip adiabatic surface using the name you have given your surface, “Slip_Wall”, an appropriate boundary condition setting might be:

[soln-bcs-Slip_Wall]

type = no-slp-adia-wall

I’m not sure whether a more appropriate slp-adia-wall or symmetry boundary condition has been implemented in PyFR for your specific situation yet. There is the char-riem-inv boundary option, but I believe that is intended for use at freestream outer boundaries. If there is a more suitable inviscid wall boundary, then I’m not sure it has reached the documentation. Give the types of flows PyFR is aimed at solving, it wouldn’t surprise me if this isn’t an option. I’m also not sure whether the boundary condition names can include underscores or recognize lower and upper case characters. Something that might be worthwhile stating explicitly somewhere in the documentation. I tend to stick with lower case names.

It appears you’re using SI units wherein u_inf has been set to 100 m/s, p_inf = 101325.0 Pa, and rho_inf = 1.2 kg/m3. This gives a Mach number ~ 0.3 which really isn’t supersonic flow, so you are much better off creating separate subsonic inflow and subsonic outflow boundaries in your model. In Gmsh, if you labeled the inflow boundary inflow, then an example for specifying this boundary in PyFR’s configuration file would be:

[soln-bcs-inflow]

type = sub-in-frv

rho = 1.0

u = 0.355

v = 0.0

w = 0.0

Here I’ve normalized/non-dimensionalized rho, u, v, and w by some reference values (e.g. u = 100m/s/[mach*sqrt(gam)]).

I also noticed in your solver-interfaces-quad and solver-elements-hex declarations that you have included a parameter vcjh-eta (probably related to a Vincent, Castonguay, Jameson, Huynh term), but this doesn’t appear in the documentation either. Lastly, I noticed that you set the order for the solver in your configuration file for the simulation to 1. You might consider changing this value to 3rd or 4th order since this is one of the advantages of Flux Reconstruction after all…

These are just some of the things I noticed straight away, but I’m sure the experts may have other comments as well.

Best Regards,

Zach Davis

3

Hi Antonio,

- When runninf the .ini file attached, an error pops-up with the
following: "ValueError: Invalid characters in expression". What is wrong
with it?

Underscores are not supported as variable names inside the .ini file.
When these are used inside of the initial condition expressions they
give rise to the invalid characters in expression error.

For those that are interested: we disallow underscores as it greatly
simplifies the logic required to ensure that initial condition
expressions are valid. If we did allow underscores then more work would
be required to ensure that the expression does not contain something
naughty in an attempt to execute arbitrary code:

  <http://nedbatchelder.com/blog/201206/eval_really_is_dangerous.html&gt;

Disallowing underscores and other 'fun' characters permits us to safely
use eval(). This is an order of magnitude simpler than pulling in a
proper expression evaluation library (or, worse still, rolling our own
with a custom parser/syntax tree).

Regards, Freddie.

4

Zach, Freddie - thanks for the input!

Antonio - please let us know if it helps.

Cheers

Peter

Dr Peter Vincent MSci ARCS DIC PhD
Senior Lecturer and EPSRC Early Career Fellow
Department of Aeronautics
Imperial College London
South Kensington
London
SW7 2AZ
UK

web: www.imperial.ac.uk/aeronautics/research/vincentlab
twitter: @Vincent_Lab

5

Hi Zach, Freddie, Peter,

Thanks a lot for your reply! In particular to Zach, it did help a lot your reply to start the code running. I have some other technical questions:

  • “order 1” corresponds to p=1 or p=0 (scheme order 1)?

  • Is it possible to monitor residuals? I see it’s required to input dt… Isn’t there a Global or Local time stepping scheme available?

  • How to specify the maximum time steps to run / convergence criteria?

  • If trying to run in parallel… Is it required to partition the .gmsh in as many blocks as procs to be used?

  • If I want a specific output (let’s say the L2 norm of entropy error on Euler cases). Is it possible to input in the .ini file? If not, Where in the code should I implement this output?

The first attempts to run this case with pyFR do not look like making any evolution of the solution in time. Perhaps I’m not setting correcltly the time integration scheme/ BCs… Could you please have a look and tell me whether it’s OK?

Thanks a lot for your help.

On the other hand, if I want to set periodic boundary conditions (not for this case though…) How should I specify the Boundaries in the .ini and .gmsh files? I see in the “euler_vortex” case uploaded in the website that no “bcs” tags are specified. Is it on the other hand necessary to name the boundaries in any way inside the .gmsh file?

Thanks once more!

Best regards,
Antonio

sphere_lev2o1_o1.ini (821 Bytes)

6

Hi Antonio,

I believe order 1 corresponds to p=1.

I think the ability to monitor residuals is on the developers to-do list. It isn’t straightforward, but they should be able to implement something in a similar vein to the way that they provide progress updates while the simulation runs.

You can specify the maximum number of time steps using the times parameter in the [soln-output] “namelist”. For example:

[soln-output]

times = range(0, 20, 21)

This will start the simulation at time t=0, and end at time t=20 using the dt value specified in the [solver-time-integrator] “namelist” while outputting 21 solution files over the course of that interval.

In order to run the solver in parallel, you first need to convert your Gmsh formatted mesh (i.e. your *.msh file) into PyFR’s native mesh format via:

pyfr mesh convert my_mesh.gmsh my_mesh.pyfrm

Then you can partition the mesh into however many parts you desire using:

pyfr mesh partition np my_mesh.pyfrm .

where np would be the number of partitions or a colon delimited list of weighting factors to decompose the graph. Note that you need to include either the absolute or relative path for the output directory in the command—here I’ve used . to represent the relative path for the current working directory.

I’m not sure setting up monitors for specific metrics is provided within the solver. Again, experts correct me if I’ve again stuck my foot in my mouth. You can however use Paraview on a solution file to calculate the specific scalar or vector field function that you’re interested in using the solution data PyFR provides. Others might have some insights on where in the code this sort of feature could be added.

The surface boundaries in your Gmsh *.msh file need to correspond to the names supplied in the your PyFR configuration file. There are no boundary conditions that need to be prescribed in the 2D Euler vortex example case. I’m not sure that a periodic boundary condition is implemented in PyFR, but I may be mis-stating the facts yet again :stuck_out_tongue: I haven’t seen the documentation if there is.

Best Regards,

Zach Davis

7

Hi Antonio,

I have seen your ini file. Regarding time step and residual, PyFR uses explicit unsteady time marching method, such as Runge-Kutta method. Thus, every cell marches the same time step dt which you imposed in your INI file.

CFL condition may roughly guide the allowable time step. From your ini file, freestream velocity and speed of sounds are 100 and 340 repectively. If the minimum characteristic length of cell is about 0.05, your time step is O(0.0001). As you know well, If you increase order, the allowable time step is decreased around O(0.0001/n). You may get NaN solution if your dt is too large. I recommend to use NaN value check option in pyfr-sim, such as -n 500 (frequency for NaN check), and find allowable time step for Runge-Kutta scheme.

Another option is to use adaptive time stepping, which find the maximum time step during computation. You need to modify [solver-time-integrator] section like below.

[solver-time-integrator]
scheme = rk45
controller = pi
t0 = 0.0
dt = 5e-5
atol = 1e-3
rtol = 1e-3

You need to use PI controller and rk45 or rk34 scheme and impose the proper value for atol and rtol. When I solve the non-dimensionalized flow condition, I use (atol, rtol) as 1e-6, but I think larger values are needed for your condition. I have run your simulation quickly and 1e-3 seems fine.

To compute and/or integrate some values, I think using Paraview is easy option, as Zach said. Actually, Paraview itself provides many features, so you can readily obtain these values. But sometimes, you may need to make a seperate code to read pyfrs solution and compute the values.

PyFR can also handle periodic boundary condition, as shown in Euler Vortex example. When you make a mesh using gmsh, you can just impose physical group name as ‘periodic_0_l’, ‘periodic_0_r’, instead of your Boundary name (slip wall and so on). You can give many periodic condition when change 0 in that name to other number. pyfr-mesh can make data structure for periodic boundary condition and you don’t need to change config file.

In addition, PyFR can also handle curved elements, thus I think that it would be better to use higher-order gmsh elements, especially for sphere surface. There is option to set order for your mesh in gmsh. In command line, you can siimply do it.

gmsh -3 -order 3 your.geo

Regards,

Jin Seok

8

Hi Zack, Jin,

Thanks a lot once again for your replies. I really appreciate your support. I’ll let you now if I run successfully this case.

Best regards,
Antonio