search

CICE-Consortium / Icepack

Development repository for sea-ice column physics
Other
25 stars 131 forks source link

Undefined behaviour in `calculate_Tin_from_qin` #482

Open kieranricardo opened 7 months ago

kieranricardo commented 7 months ago

Hi there,

I've come across some undefined behaviour in calculate_Tin_from_qin. Depending on the compiler and compilation options, when the sqrt argument is negative the following will either yield a NaN or Tmltk. Replacing min with ieee_min_num fixes this for my use case, and will consistently return Tmltk when the sqrt argument is negative, independent of the compiler.

https://github.com/CICE-Consortium/Icepack/blob/f6ff8f7c4d4cb6feabe3651b13204cf43fc948e3/columnphysics/icepack_therm_shared.F90#L80C1-L81

Tin =  min((-bb1 - sqrt(bb1*bb1 - c4*aa1*cc1)) /  &
                         (c2*aa1),Tmltk)

Happy to submit a PR from my fork https://github.com/kieranricardo/Icepack/tree/bugfix-calculate_Tin_from_qin-undefined-behaviour :)

eclare108213 commented 7 months ago

The quantity under the square root probably should not ever be negative. What are the values of aa1, bb1, cc1 when this fails? Is it a roundoff problem or a bigger issue that needs to be looked into?

dabail10 commented 2 months ago

Any updates on this one?

anton-seaice commented 1 month ago

Kieran is on leave this week, we resolved why the quantity was negative.

We could just ignore this, or if we think returning Tmltk is better than returning Nan when the first argument is Nan, then make the change?

dabail10 commented 1 month ago

I am happy to put in a fix for this. Maybe we should just put in a check for the argument in the square root and set it to zero when negative. Then if the result is greater than Tmltk, then you just get that. I suspect this is a roundoff issue as highlighted by @eclare108213 .

anton-seaice commented 1 month ago

@apcraig - Do you have any experience with ieee_min_num and if it is better than min? I am not sure about compiler support for it?

dabail10 commented 1 month ago

I did some looking into this and ieee_min_num is a fortran 2018 feature, and so we need to be careful here. I think the typical approach here is to do the following instead:

Tin =  min((-bb1 - sqrt(max(bb1*bb1 - c4*aa1*cc1),puny)) /  &
                         (c2*aa1),Tmltk)
anton-seaice commented 1 month ago

@kieranricardo resolved why the quantity was negative some time ago, we think it was an issue with how the coupling was getting set-up.

My view is we shouldn't change anything, because we don't think its valid for the quantity under the square root to be negative, and giving an error helps with debugging.

I guess a more graceful handling would be to check if bb1*bb1 - c4*aa1*cc1 is negative, and abort with a meaningful error if it is.