search

termolivre / Psychro.jl

Thermodynamic properties of moist air in Julia
Other
5 stars 6 forks source link

`calcz` doesn't handle many numeric types due to overly-restrictive type-signature in `cardan` #3

Open wnoise opened 1 year ago

wnoise commented 1 year ago

Example 1: bigfloats

spechum(MoistAir, big"300", WetBulb, 291.15, 101325.0)

gives

ERROR: MethodError: no method matching cardan(::Tuple{BigFloat, BigFloat, Float64, Float64})
Closest candidates are:
  cardan(::NTuple{4, T}) where T<:AbstractFloat at ~/external/Psychro.jl/src/utilities.jl:80

ForwardDiff also fails:

using ForwardDiff
spechum_t(t) = spechum(MoistAir, t, WetBulb, 291.0, 101325.0)
ForwardDiff.derivative(spechum_t, 293.0)

with

ERROR: MethodError: no method matching cardan(::Tuple{ForwardDiff.Dual{ForwardDiff.Tag{typeof(spechum_t), Float64}, Float64, 1}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(spechum_t), Float64}, Float64, 1}, Float64, Float64})
Closest candidates are:
  cardan(::NTuple{4, T}) where T<:AbstractFloat at ~/external/Psychro.jl/src/utilities.jl:80

The obvious thing to do is just remove the constraint (and the uses of T() inside the function.)

This works fine, (including still working with Unitful quantities, as those are dealt with higher up). However examining the function closer, cardan is passed information it cannot use, as the last element must always be 1. In addition, as cardan is not exported, and only called calcz, the second-to-last element is always -1.0, but isn't reflected in cardan.

I would recommend:

  1. rename cardan to make clear it isn't quite the normal cardano algorithm and only handles cases with 1.0x^3 + ...,
  2. specialize it to 1.0 x^3 - 1.0 x^2 + ..., and actually work out the consequences, though julia/llvm might be smart enough to propagate the constants effectively enough not to matter.
  3. make it take two arguments, rather than a tuple, with each argument separately descending from AbstractFloat.
  4. keep the constants as Float64 (as almost anything will be able to multiply by Float64s, though 1/3 as 0.33333... loses some precision with e.g. BigFloat. 1//3 might specialize well enough though?)
pjabardo commented 1 year ago

I will give it a look this weekend.

For a long time, the finding the root was done using an iterative procedure. A couple of years ago, someone made a pull request and implemented an analytical solution implementing Cardano's algorithm. I really didn't take a deeper look but it worked pretty well and is probably much faster, more accurate and robust than the old iterative procedure.

But I think this uses complex numbers (I haven't studied the implementation) and specifies the type (Float64). I think that both of these points is in conflict with algorithmic differentiation and more specifically ForwardDiff.jl. For my day to day use this is not a issue but it would be really nice to have the package working flawlessly with AD and SciML stuff.

I think we are at a point where we could have something as nice but more flexible and extensible than EES (engineering equation solver)

I will try to study this as soon as possible.

Paulo