Hello All,
I am hacking around with PyFR trying to see if I can implement reacting flows into the solver.
My first step has been to just get species mass fractions as passive scalars in the euler solver and implement mixture thermodynamics (i.e. calculate mixture Cp, Rg, Cv, etc. to then calculate P, T, etc.)
This has all gone really well but I have hit a snag that I have narrowed down to simply catastrophic round-off/cancellation errors in floating point arithmetic in the C code templates (which I still am having trouble believing 100% myself but I have exhausted my debugging abilities and I have sat and watched the values slowly lose precision as the simulation runs)
Here has been my procedure (all of this is with the openMP backend, in serial, with 1 thread, NaNs checked at every step, and all optimizations turned off in the compilers):
In the euler flux.mako template I am replacing the line
p = ${c['gamma'] - 1}*(E - 0.5*invrho*${pyfr.dot('rhov[{i}]', i=ndims)});
with
//Calculate internal energy
U = (E - 0.5*invrho*${pyfr.dot('rhov[{i}]', i=ndims)});
From there it is just a straightforward calculation of thermodynamics from internal energy assuming CPG for now (I add to the template arguments the number of species, as well as the variable “species_data” which is a list of lists in which the first list is the first species’ [MW, Cp] and the second list is the second species’ [MW, Cp] and so on.
fpdtype_t mw_mix = 0.0; //initialize mixture molecular weight
fpdtype_t cp_mix = 0.0; //initialize mixture Cp
fpdtype_t Runi = 8.3144598; //initialize universal gas constant
fpdtype_t Rg_mix; //initialize specific gas constant
fpdtype_t cv_mix; //initialize mixture cv
// Compute mixture MW (MW units in kg/mol)
% for i in range(nspecies):
mw_mix += s[${ndims + 2 + i}]*${species_data[i][0]}; \\After the E index of 's' I have the species mixture fraction conservative vars in the array
% endfor
// Compute mixture cp
% for i in range(nspecies):
cp_mix += s[${ndims + 2 + i}]*${species_data[i][1]};
% endfor
// Compute mixture Rg
Rg_mix = Runi / mw_mix;
// Compute mixture cv
cv_mix = cp_mix - Rg_mix;
// Compute the pressure using calorically perfect gas assumption
p = U * Rg_mix / cv_mix;
I am running a simple case (same mesh as euler vortex example) with uniform pressure, uniform species mixture fraction (0.79 N2 and 0.21 O2). If I just manually calculate Rg_mix and cp_mix and put them in the constants section of the ini file and replace the calculations above with the constant values it runs perfectly.
** But when I implement the above calculations, specifically the calculations of mixture Cp and molecular weight. The values start off perfect (I can see from print outs to the terminal) and slowly lose precision until NaNs are detected.***
This is surprising on a few levels, first I am only adding 13 floating point operations, and second, none of them are violating the numerical NO-NO’s (that I know of) like subtracting two numbers of near equal value. For example in the calculation cv_mix = cp_mix - Rg_mix, the magnitude of cp_mix ~1,000 and Rg_mix ~288.
Anyway, I am just reaching out to the community to see if anyone else has run into a similar issue or anyone has an idea of what’s going on here.
Thanks to all who work on this project, it is a wonderful code base and has been a joy to utilize.
Kyle