Check == and <> usages in functions

[modelica-trac-importer]

Reported by otter on 19 Feb 2016 13:31 UTC
In https://trac.modelica.org/Modelica/ticket/1867#comment:41 the following was reported:

The LMS Amesim compiler found several "warnings" related to equality comparison of Real variables. For example, there are lines like if (d_vap <> d_liq) or if abs(m) <= tol or fb == 0.0 . I suspect most other tools just ignore the warnings, like we do. However, it might be nice to clean up as many of them as we can before the v3.2.2 release. This is most likely just a partial list…

Functions that use == comparisons on Real variables:

Modelica.Media.Water.IF97_Utilities.BaseIF97.Basic.f3
Modelica.Media.Common.OneNonLinearEquation.solve 

Functions that use <> comparisons on Real variables:

Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph
Modelica.Media.Water.IF97_Utilities.waterBaseProp_ps
Modelica.Media.Water.IF97_Utilities.waterBaseProp_pT
Modelica.Media.Water.IF97_Utilities.BaseIF97.Isentropic.hofps4
Modelica.Media.Water.IF97_Utilities.BaseIF97.TwoPhase.waterSat_ph
Modelica.Media.Water.IF97_Utilities.BaseIF97.TwoPhase.waterR4_ph
Modelica.Media.Water.IF97_Utilities.BaseIF97.TwoPhase.waterR4_dT 

Although == and <> are allowed in Modelica functions, there might be issues and the functions should be inspected.

I have analyzed the first two (Hubertus, please have a look at the other ones):

Modelica.Media.Common.OneNonLinearEquation.solve

Here the following code is present

 if abs(m) <= tol or fb == 0.0 then
    // root found (interval is small enough)

This code fragment checks whether the interval of the independent variable (abs(m) <= tol) is small enough or whether the residue is exactly zero. If this is the case, the root was found. The "fb == 0.0" is correct here: We cannot use an eps, abs(fb) <= eps, because the "size" of fb is not known (could be very large or very small). However, if the residue is exactly zero, it is impossible to further reduce the x-interval. Therefore, this check must be present.

Modelica.Media.Water.IF97_Utilities.BaseIF97.Basic.f3

Here the following code fragement is present:

f.delta := if (d == data.DCRIT and T == data.TCRIT) then 
                1 - Modelica.Constants.eps 
           else abs(d/data.DCRIT);

Looks a bit surprising to me and it might be better to have an eps here, such as (but I am unsure):

f.delta := if ( abs(d-data.DCRIT)<=eps and 
                abs(T-data.TCRIT)<=eps then 
                1 - Modelica.Constants.eps 
...

---

Migrated-From: https://trac.modelica.org/Modelica/ticket/1915

[modelica-trac-importer]

Comment by hansolsson on 19 Feb 2016 14:01 UTC
Just so that the previous comment isn't lost:
https://trac.modelica.org/Modelica/ticket/1867#comment:44

It is formally correct Modelica to use equality tests - but it might indicate a problem, and might cause problems for e.g. automatic differentiation.

Consider f3:

  f.delta := if (d == data.DCRIT and T == data.TCRIT) then 1 - Modelica.Constants.eps
     else abs(d/data.DCRIT);

Around the critical density f.delta goes from -1/data.DCRIT to 1/data.DCRIT.

IF we want f.delta to be continuously differentiable we will have to create a complicated expression. I don't know if it is worth it, and I don't see the point in any other change. (I don't even know if we want to differentiate this function in some cases.)

[modelica-trac-importer]

Comment by hansolsson on 19 Feb 2016 15:30 UTC
Replying to [comment:1 hansolsson]:

Just so that the previous comment isn't lost:
https://trac.modelica.org/Modelica/ticket/1867#comment:44

It is formally correct Modelica to use equality tests - but it might indicate a problem, and might cause problems for e.g. automatic differentiation.

Consider f3:

  f.delta := if (d == data.DCRIT and T == data.TCRIT) then 1 - Modelica.Constants.eps
     else abs(d/data.DCRIT);

Around the critical density f.delta goes from -1/data.DCRIT to 1/data.DCRIT.

Oh, sorry - this is wrong.

The f.delta derivative is 1/data.DCRIT both before and after the critical point.

I thought the code was more complicated.

That makes it easier to support automatic differentiation (assuming DCRIT is positive!):

   f.delta := if (d == data.DCRIT and T == data.TCRIT) then 1 - Modelica.Constants.eps
      + (d/DCRIT-1) else abs(d/data.DCRIT);

But that also makes me wonder what problem would f.delta=1 at the critical point cause?

[modelica-trac-importer]

Comment by hubertus on 19 Feb 2016 22:17 UTC
There are derivatives that jump from +1/x to -1/x when crossing the critical point (cp does that). This is a type of function where I would really not recommend to use AD to compute derivatives. In no symbolic math engine that I ever used (Mathematica, Maple, and an open source one) this ever worked properly without my massaging the code manually before and after taking derivatives. That aside, it is worthwhile checking that the current code is fully numerically robust. But it has been used for around 15 years now without me ever getting a bug report on this, so I would strongly prefer to see a failure-example before I change anything.