Hi everyone,
I have noticed an issue with the function setState(rho,rhoe,0) using the EquilStateModel for conditions that have zero (specific) energy.
Example test
For a mixture of air5p ≈ 4127.07 [Pa]T ≈ 416.4035 [K]
one gets
ρ ≈ 0.0343909 [kg/m^3]e ≈ 0.0737387 [J/kg]ρe ≈ 0.0025359 [J/m^3]
Trying to setState with these values (or similarly, with any combination with energy density abs(rhoe)<1e-2 in my testing) will cause either a division by zero (if you are very very unlucky), or an endless while loop with no convergence.
The cause
In file src/thermo/EquilStateModel.cpp at line 110 there is a division f1/rhoe while the value rhoe, being the energy density, takes both positive and negative values, hence it may be zero or infinitely close to zero.
Bad solution
Most of the times I (or someone else) had this problem (never looked for its cause before) the solution was to change the tolerance tol to a larger value (eg. 1e-6 from 1e-12). I had done that while working with COOLFluiD two years ago in my code, and I know that this is what they did at Ames where they use Mutation++ with the ARCHeS code. This does not solve the problem (just narrows the problematic regime around zero).
Solution 1
The easiest quick-fix is to replace line 110 of the file
while (std::max(std::abs(f1/rhoe),std::abs(f2/rho)) > tol) {
with
while ((std::max(std::abs(f1/rhoe),std::abs(f2/rho)) > tol) && (std::max(std::abs(f1),std::abs(f2)) > tolAbs)) {
This way, for values of energy
away from zero the relative convergence limit will cause the while loop to end as before, but
near zero the absolute convergence limit will take over (since the relative one would never be satisfied).
This does not fix the rare issue with exactly zero energy, which could be avoided by instead dividing f1/(rhoe+tolAbs) and f2/(rho+tolAbs) to avoid the division by zero.
So line 110 has to become
while ((std::max(std::abs(f1/(rhoe+tol)),std::abs(f2/(rho+tol))) > tol) && (std::max(std::abs(f1),std::abs(f2)) > tol)) {
For me, tolAbs = 1e-10; tol = 1e-12 has worked, but it might depend on the conditions and the system.
Solution 2
A possibly better idea would be to use the iterative residual as a convergence criterion, making line 110 like
while ((std::max(std::abs((f1-f1_old)/(rhoe+tol)),std::abs((f2-f2_old)/(rho+tol))) > tol) && (std::max(std::abs(f1-f1_old),std::abs(f2-f2_old)) > tol)) {
and of course defining f1_old, f2_old accordingly.
Notes:
There could still be the issue with division by zero or near-zero for relative tolerance whether you implement either Solution 1 or Solution 2, unless you add something to the denominator.
It is not a good idea to check only the absolute convergence, since it will be unnecessarily long in the case of values away from zero.
Hope this summary makes it easy to fix it. I talked to @grgbellas and he said he'll take a look soon.
jbscoggi
commented
4 years ago
Hi everyone, I have noticed an issue with the function
setState(rho,rhoe,0)using theEquilStateModelfor conditions that have zero (specific) energy.Example test
For a mixture of
air5p ≈ 4127.07 [Pa]T ≈ 416.4035 [K]one getsρ ≈ 0.0343909 [kg/m^3]e ≈ 0.0737387 [J/kg]ρe ≈ 0.0025359 [J/m^3]Trying tosetStatewith these values (or similarly, with any combination with energy densityabs(rhoe)<1e-2in my testing) will cause either a division by zero (if you are very very unlucky), or an endlesswhileloop with no convergence.The cause
In file
src/thermo/EquilStateModel.cppatline 110there is a divisionf1/rhoewhile the valuerhoe, being the energy density, takes both positive and negative values, hence it may be zero or infinitely close to zero.Bad solution
Most of the times I (or someone else) had this problem (never looked for its cause before) the solution was to change the tolerance
tolto a larger value (eg.1e-6from1e-12). I had done that while working with COOLFluiD two years ago in my code, and I know that this is what they did at Ames where they use Mutation++ with the ARCHeS code. This does not solve the problem (just narrows the problematic regime around zero).Solution 1
The easiest quick-fix is to replace
line 110of the filewhile (std::max(std::abs(f1/rhoe),std::abs(f2/rho)) > tol) {withwhile ((std::max(std::abs(f1/rhoe),std::abs(f2/rho)) > tol) && (std::max(std::abs(f1),std::abs(f2)) > tolAbs)) {This way, for values of energy
whileloop to end as before, butThis does not fix the rare issue with exactly zero energy, which could be avoided by instead dividing
f1/(rhoe+tolAbs)andf2/(rho+tolAbs)to avoid the division by zero.So
line 110has to becomewhile ((std::max(std::abs(f1/(rhoe+tol)),std::abs(f2/(rho+tol))) > tol) && (std::max(std::abs(f1),std::abs(f2)) > tol)) {For me,
tolAbs = 1e-10; tol = 1e-12has worked, but it might depend on the conditions and the system.Solution 2
A possibly better idea would be to use the iterative residual as a convergence criterion, making
line 110likewhile ((std::max(std::abs((f1-f1_old)/(rhoe+tol)),std::abs((f2-f2_old)/(rho+tol))) > tol) && (std::max(std::abs(f1-f1_old),std::abs(f2-f2_old)) > tol)) {and of course definingf1_old, f2_oldaccordingly.Notes:
Hope this summary makes it easy to fix it. I talked to @grgbellas and he said he'll take a look soon.
Cheers, Thanos