Standards for EOS return values

[adamjermyn]

I don't know precisely what this will look like, but I'd like to propose that we figure out some standards for what the EOS is allowed to return, and say that return values that violate those conditions must raise ierr /= 0. This would be the job of individual EOS's to enforce, not a blend-level enforcement.

A few standards that come to mind:

1. No EOS return quantity should be NaN.
2. No negative pressures.
3. No negative sound speeds.
4. No negative adiabatic indices.

Is there anything else we want to check? We can also make a module-level test that runs a few grids of (rho, T) at different compositions as an easy way of catching these...

[fxt44]

for a GIGO model, no negative or NAN input quantities - notionally T, Rho, composition. on the output side of GIGO, maybe add no negative entropies and specific heats.

[adamjermyn]

I've made a new branch (eos_better_tests) for this. It currently implements your suggestions on the input side. Working on the output side now.

[adamjermyn]

Implemented output checks in this commit: https://testhub.mesastar.org/eos_better_tests/commits/a249bd4

These checks all throw an ierr /= 0 and return immediately.

[adamjermyn]

Here are the current output checks:

         do k=1,num_eos_basic_results
            if (is_bad(res(k) .or. is_bad(d_dlnd(k)) .or. is_bad(d_dlnT(k))) then
               ierr = 1
               return
            end if
            do l=1,species
               if (is_bad(d_dxa(k,l))) then
                  ierr = 1
                  return
               end if
            end do
         end do

         if (res(i_grad_ad) < 0d0 .or. res(i_chiRho) < 0d0 .or. res(i_chiT) < 0d0 &
            .or. res(i_Cp) < 0d0 .or. res(i_Cv) < 0d0 .or. res(i_gamma1) < 0d0 &
            .or. res(i_gamma2) < 0d0 .or. res(i_gamma3) < 0d0) then
            ierr = 1
            return
         end if

[fxt44]

another option on the input side is to only consider values within a valid range - no temperatures less than 10 K or larger than 1e13 K, no densities larger than 1e12 g/cc, abundances greater than 1, and so on. might get false positives if the validity range is made too small. i also wonder how expensive all this checking will get. maybe its an option?

[adamjermyn]

Added.

I'd be shocked if these checks matter for runtime. EOS calls take microseconds, these checks can't possibly add more than ~hundreds of cycles...

Copying and pasting from slack, we have a problem in net/test where it was silently passing garbage inputs to the EOS (yet passing!):

I've added GIGO checks to the EOS inputs (throwing ierr). Now net fails its test because test_one_zone_burn_small_net tries to run an EOS call with negative abundances.

Here's the tail of the output:
  **************** approx21_cr60_plus_co56 **************** 
                                log temp         4.623300792
                                 log rho       -10.746410108
                               rate_raw rco56ec_to_fe56    1.0402929319524918D-07
                               rate_raw rni56ec_to_co56    1.3152697923338621D-06
 test_one_zone_burn_small_net
 logT=   8.3000000000000007     
 logRho=   7.0000000000000000     
 x=   1.1195903029095758E-030  -1.7186989648912257E-031   1.0000000020735527E-030   3.0005476214493425E-030   1.2802186851313698E-002  0.98714342252602971        5.3226249706257383E-026   5.4390620909983984E-005   2.3809678445156945E-012   1.6462501154694145E-018   6.4280022786219520E-029   3.2000000002366084E-029   3.6000000000001163E-029   4.0000000000000003E-029   4.4000000000000004E-029   4.8000000000000004E-029   5.6000000000000005E-029   5.2000000000000004E-029   5.4000000000000005E-029   7.9348354368784564E-029   5.5999992429863682E-029
(I added printing logT, logRho,and xa)
Could someone who knows what this test is doing tell me what to patch here?
It looks like it's just numerical noise that let the burn solver push to negative abundances, but then it doesn't clean that up before the EOS call.
So we could just clean that up. But I'm not sure where the best place is for that...
(But passing negative xa to the EOS is definitely not okay!)

[fxt44]

usually between steps i do things like x(i) = max(x(i), floor) where here floor is 1d-30 and optionally renormalize the positive definite mass fractions.

[adamjermyn]

@warrickball The composition being set in set_composition in atm/test/test_atm_support.f90 had a substantial negative entry (of order -1d-3). This appears to be because one component was set by subtracting the sum of the rest from 1. I've fixed this by putting a floor of 0d0 on all components and then normalizing-by-division afterwards. Let me know if you object to this.

[warrickball]

For the record, I think @rhdtownsend is really the ur-master of atm (and is, in fact, the one who wrote the offending composition back in r11869) but I'm happy with that change. To be specific, I presume the offending composition is

     X = 0.70d0
     Z = 1d-2
     XC = 3.2724592105263235d-03
     XN = 9.5023842105263292d-04
     XO = 8.8218000000000601d-03

where the XC, XN and XO are copied from the Z=0.02 cases. set_composition() computes Y from 1-X-Z and ultimately gets a ²⁴Mg abundance from Z-XC-XN-XO, which is the issue. The intention may have also been to halve those abundances, in line with the changed Z, which I think would also fix the issue.

[adamjermyn]

Ok, great. As long as someone more informed than me has looked at it I'm happy!

-Adam

On Mon, Jun 28, 2021 at 3:52 PM, Warrick Ball @.***> wrote:

For the record, I think @rhdtownsend https://github.com/rhdtownsend is really the ur-master of atm (and is, in fact, the one who wrote the offending composition back in r11869) but I'm happy with that change. To be specific, I presume the offending composition is

 X = 0.70d0

 Z = 1d-2

 XC = 3.2724592105263235d-03

 XN = 9.5023842105263292d-04

 XO = 8.8218000000000601d-03

where the XC, XN and XO are copied from the Z=0.02 cases. set_composition() computes Y from 1-X-Z and ultimately gets a ²⁴Mg abundance from Z-XC-XN-XO, which is the issue. The intention may have also been to halve those abundances, in line with the changed Z, which I think would also fix the issue.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/MESAHub/mesa/issues/285#issuecomment-869985933, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABPR5HY6HU2YQBXRJLKNAZDTVDHJRANCNFSM47ATBT5Q .