Support of HelmholzMedia by OpenModelica

[casella]

@thorade, we've made good progress with the support of your library in OpenModelica. @perost did an excellent job with the new frontend, so we passed from two models only running with the old front end, to all except one of the 42 runnable models of the library getting successfully compiled into simulation code, see report. Unfortunately, about a quarter of them still fail during simulation because of numerical issues.

I haven't had the time so far to look into that, so I'm not sure whether this is because of some shortcomings of the OpenModelica backend/runtime. Could you have a quick look and report here?

[thorade]

Library report, Library coverage, last 7 runs

1. Examples.ConvergenceTest.Ancillary_Saturation
   This goes into a branch of an if-else that I thought is unreachable (this should also not happen)??? After updating the assert message, it is seen that h=nan

2. Examples.ConvergenceTest.setSat
   An iteration written by me does not converge, close to the critical point, but it works in Dymola. There is a warning about abs value being negative.

3. Examples.MediaTestModels.ButaneTestModel_dT_component_ph
   There are 2 more Butane examples with identical error message: log(delta - g[i,9]), but when searching the code for this expression, I cannot find it? Also not when searching for delta - g[i,9] or just delta. exp(g[i,6]*(delta - g[i,9]) is called several times from the Interfaces\PartialHelmholtzMedium\EoS\f_r...mo functions.
   Same error also for Helium, Isobutane and Propane, in total 6 cases.

4. Examples.MediaTestModels.CarbondioxideTestModel error message during initialization: functions.c:2865: Invalid root: (-0.996066)^(1.33333). Maybe this error message can be enhanced to give more information?

5. Examples.Validation.Derivatives_SaturationBoundary
   fails during initialization, with liq.s and vap.s having same value

6. Examples.Validation.idealGasLimit
   this should converge to a certain value, but starts to oscillate because of numerical noise, also with other tools, but at a later point. Might be solved by changing some tolerance or other setting?

[thorade]

Thanks @perost and @casella for the hard work. I already had a quick look some days ago, and wrote down some notes now. Not all of it makes sense to me right away. One error causes 50% of the failures, log(negative) during initialization.

[casella]

We'll look into this ASAP, given the growing interest in sCO2 cycles your library can be indeed very useful.

[branch171]

I have found the error associated with the log(delta - g[i,9]). This is assocaiated with the code f_r

• sum(g[i,1]*tau^g[i,2]delta^g[i,3]exp(g[i,6](delta - g[i,9])^g[i,4] + g[i,7](tau - g[i,8])^g[i,5]) for i in 1:nGauss)

(delta - g[i,9])^g[i,4] gets calculated as exp(g[i,4] *log(delta - g[i,9])) howver it is always a square so the negative argument it not an issue

it needs to be replaced with:

• sum(g[i,1]*tau^g[i,2]delta^g[i,3]exp(g[i,6](delta - g[i,9])^2 + g[i,7](tau - g[i,8])^2) for i in 1:nGauss)

i.e. remove g[i,4] and g[i,5] and replace with 2's. Note in other parts, i.e derivatives) this is already done g[i,4] and g[i,5] are not used.

The is also an error with the residualNonAnalytical terms which calculate a (-ve)^1.3333. I haven't found that yet but it can be fixed by setting all the NonAnalytical coefficients to zero. This is OK as residualNonAnalytical is only needed very close to the critical point.

[casella]

Excellent! @branch171, are you able to provide a PR to fix it?

[branch171]

what is a PR? I can show you which line to change.

[branch171]

The code to change is in function f_r

I have commented out the old line and below is the replacement:

algorithm
  f_residual :=
    sum(p[i,1]*tau^p[i,2]*delta^p[i,3] for i in 1:nPoly)
  + sum(b[i,1]*tau^b[i,2]*delta^b[i,3]*exp(-delta^b[i,4]) for i in 1:nBwr)
//+ sum(g[i,1]*tau^g[i,2]*delta^g[i,3]*exp(g[i,6]*(delta - g[i,9])^g[i,4] + g[i,7]*(tau - g[i,8])^g[i,5]) for i in 1:nGauss)
  + sum(g[i,1]*tau^g[i,2]*delta^g[i,3]*exp(g[i,6]*(delta - g[i,9])^2      + g[i,7]*(tau - g[i,8])^2)      for i in 1:nGauss)
  + sum(a[i,1]*Disb[i]*delta*Psi[i] for i in 1:nNonAna);

To fix the nNonAna (-ve)^1.3333 error in CarbonDioxide: set the residualNonAnalytical all to zero and comment out the previous values.

     residualNonAnalytical=fill(0.0, 0, 12)) 
     //residualNonAnalytical=[
     // -0.666422765408E+00,  0.000,   1.00,    2, 2,  0.875,  0.300, 0.70, 10.0, 275., 0.3, 3.5;
     //  0.726086323499E+00,  0.000,   1.00,    2, 2,  0.925,  0.300, 0.70, 10.0, 275., 0.3, 3.5;
     //  0.550686686128E-01,  0.000,   1.00,    2, 2,  0.875,  0.300, 0.70, 12.5, 275., 1.0, 3.])

[thorade]

Thanks for the input!
PR stands for Pull Request, you can learn how to do it here:
https://docs.github.com/en/github/collaborating-with-issues-and-pull-requests/proposing-changes-to-your-work-with-pull-requests
Each page should have an edit button also, for small changes this would be the easiest way to propose a change. For the first time it seems to be difficult, but it is worth learning it!

I am not 100% sure what would be the best solution here. RefProp used the coefficients g4 and g5 instead of hardcoded 2 because they see the possibility for other vaues. Before changing this, we should probably check no fluid (in RefProp) ever uses this. Maybe OpenModelica can also do some additional transformations/evaluations with constants!? Also, the values for density and also delta could be limited to a certain range.

For the non-analytical terms, this change would have to be moved into some if-else probably. These terms are only used for water and carbondioxide, where the critical region can be the operating region.

[branch171]

Just as a point the code for derivatives does NOT use g4 and g5 and has hardcoded 2 for these already. So if we use g4 and g5 then the fix needs to propagate out to these as well. It would seem strange to change from hardcode 2 as by definition a Gaussian is exp(-x^2).

Yes I agree the analytical term should have if-else. Its not really important until you are very close to the critical point.

I will look at the Pull Request and help on this.

[thorade]

I am on vacation until end of July, and will have a look at the code when I am back.
Noting two more ideas, so I do not forget about it:

• OpenModelica treats Reals and Integers differently during translation, I believe, so one probably easy to implement (in OpenModelica) trick could be to check the values of Real constants: If they have an Integer value, OpenModelica could apply the same transformations also here.
• On library side, we could split up the large g[i,9] matrix of Reals, and instead have vectors of Reals or Integers, as in CoolProp, see these json files: https://github.com/CoolProp/CoolProp/tree/master/dev/fluids

[casella]

* OpenModelica treats Reals and Integers differently during translation, I believe, so one probably easy to implement (in OpenModelica) trick could be to check the values of Real constants: If they have an Integer value, OpenModelica could apply the same transformations also here.

@thorade, I'm not sure what you mean here. I'm on holiday in August, but I should be online from Aug 15, we can continue the discussion then.

[thorade]

According to the analysis by branch171:

g[i,9]^g[i,4]  // does NOT work
g[i,9]^2  // works

Even if we have constant Real g[i,4] = 2.0. So I was wondering whether OM could be smarter here, because for the human eye both cases are more or less identical.

Of course this could also be fixed on library side, by using hardcoded 2, but before doing this change, I would like to check whether RefProp ever uses a value that is different from 2. The normal Gaussian bell uses a 2, and the derivatives already use hardcoded 2, as pointed out by branch171.

Another thing that might be worth trying is to just split up the matrix, and use Integers where possible. Before doing this, I would like to check whether it even solves the issue, and then check with RefProp or CoolProp which exponents can be restricted to Integers instead of Reals. For checking in OM, just use

constant Integer g4 = 2;
constant Integer g5 = 2;
algorithm
  f_residual :=
    sum(p[i,1]*tau^p[i,2]*delta^p[i,3] for i in 1:nPoly)
  + sum(b[i,1]*tau^b[i,2]*delta^b[i,3]*exp(-delta^b[i,4]) for i in 1:nBwr)
//+ sum(g[i,1]*tau^g[i,2]*delta^g[i,3]*exp(g[i,6]*(delta - g[i,9])^g[i,4] + g[i,7]*(tau - g[i,8])^g[i,5]) for i in 1:nGauss)
  + sum(g[i,1]*tau^g[i,2]*delta^g[i,3]*exp(g[i,6]*(delta - g[i,9])^g4     + g[i,7]*(tau - g[i,8])^g5)     for i in 1:nGauss)
  + sum(a[i,1]*Disb[i]*delta*Psi[i] for i in 1:nNonAna);

[branch171]

The issue here is how to deal with x^y. Where x^y the y can be real it maps to exp(yln(x)) for a c compiler. This is a problem when x is negative as ln(-ve) is undefined. However, x^2 can be interpreted as xx and this of cause works for any x even -ve. Within the helmholtz function the Gaussian can apply to -ve and +ve x as it a function of distance of delta and tau normally either side of the critical point. You could fix this with exp(yln(abs(x))). I would be surprised if this was ever anything other than y = 2. According to R Span's book Multiparameter Equations of State the Guassian terms are always exp(-a(delta- b)^2 - c(tau-d)^2) with 2 fixed. Its not a parameter. See page 151. If y was fractional it would give complex numbers and if odd it would give unsymmetrical values that changed sign this would be bad in the eos. Also if evaluated at the critical point exp(y*ln(0)) is undefined for any y.

[casella]

The definition of x^y in Modelica is is given in Section 10.6.7 of the language specification:

Exponentiation ”a^b” is defined as pow(double a,double b) in the ANSI C library if both ”a” and ”b” are Real scalars. A Real scalar value is returned. If ”a” or ”b” are Integer scalars, they are automatically promoted to ”Real”. Consequences of exceptional situations, such as (a==0.0 and b<=0.0, a<0 and b is not an integer) or overflow are undefined

As I understand, the definition of pow in ANSI C is the following

double pow(double x, double y);

The pow function returns x raised to the power of y. If x is negative and y is not an integer value, the pow function will return a domain error.

I guess this relies on the fact that the IEEE double precision standard allows to identify Real numbers that happen to be integer (e.g. 2.0) and hence handle them appropriately.

From this point of view, the transformation made by the OpenModelica backend is illegitimate, because it assumes x > 0, which is not necessarily true. I opened #6068 on this topic.

As a workaround while OMC is fixed, maybe writing g[i,9]^integer(g[i,4]) should be enough to avoid the (illegal) transformation, because it would force the exponent to be treated as an Integer equals to two, which apparently works fine. Would you mind trying this out?

[thorade]

Using g[i,9]^integer(g[i,4]) also works.
I tried creating a minimal example, but failed to do so.

[casella]

Using g[i,9]^integer(g[i,4]) also works.

That sounds good.

I tried creating a minimal example, but failed to do so.

Me too. Never mind.

[thorade]

I would be surprised if this was ever anything other than y = 2. According to R Span's book Multiparameter Equations of State the Guassian terms are always exp(-a(delta- b)^2 - c(tau-d)^2) with 2 fixed.

If I read the RefProp Fortran source code and the .fld files correctly, some fluids like R125 use exponents that are different from 2. The book by Span was published in 2000.

[branch171]

is the RefProp source code on GIT. I'll take a look and see what is happening.

[thorade]

RefProp source code is not on github or anywhere else online, but if you install RefProp, you should have received a copy of the code. I looked at lines 147-181 of the file "C:\Program Files (x86)\REFPROP\FORTRAN\CORE_FEQ.FOR".

[branch171]

is there a chance you could attach the file or just those lines. I don't have RefProp.

[kabdelhak]

I fixed the problem in PR6685. The workaround should not be necessary anymore.

For anyone interested, the explanation of what is happening can be found in #6068.

[thorade]

Thank you @kabdelhak !!
Large step forward towards full support, coverage went up from 29 to 35 (out of 42).

[kabdelhak]

Thank you @kabdelhak !!
Large step forward towards full support, coverage went up from 29 to 35 (out of 42).

No problem!

I also had a short look at 4. (Examples.MediaTestModels.CarbondioxideTestModel error message during initialization: functions.c:2865: Invalid root: (-0.996066)^(1.33333). Maybe this error message can be enhanced to give more information?)
since it seemed related to what i already did.

Unfortunately it would be hard to give a better error message here, although i totally agree it should. We would need a new way of reporting and writing these errors during simulation since they rely entirely on temporary variables.

But i was right, it is related. The function in question is once again HelmholtzMedia.Interfaces.PartialHelmholtzMedium.EoS.f_r. In this case it is the partial derivative w.r.t. delta (but the problem should also occur in the original function, the derivative is just triggered first). The part which is problematic is the binding of variable Theta:

  Real[nNonAna] Theta = {(1-tau) + a[i,8]*((delta-1)^2)^(0.5/a[i,7]) for i in 1:nNonAna} "Theta";

Here we have the exponent ^2^(0.5/a[i,7]) which can be simplified to ^1/a[i,7] (weird way of writing it in the first place, but i guess it is to mirror the other bindings which also have an exponent of 2).

First of all the minimum value for delta is only set to 0 in this model and not to 1, so the base delta-1 can become negative. This would be not problematic with integer valued exponents, but in this case they are not. a is bound by residualNonAnalytical, which is defined in HelmholtzMedia.HelmholtzFluids.Carbondioxide as

     residualNonAnalytical=[
      -0.666422765408E+00,  0.000,   1.00,    2, 2,  0.875,  0.300, 0.70, 10.0, 275., 0.3, 3.5;
       0.726086323499E+00,  0.000,   1.00,    2, 2,  0.925,  0.300, 0.70, 10.0, 275., 0.3, 3.5;
       0.550686686128E-01,  0.000,   1.00,    2, 2,  0.875,  0.300, 0.70, 12.5, 275., 1.0, 3.]

Here we can see: a[i,7] pointing to the 7th row is always 0.3 which results in a non integer exponent. As far as i can see there are two things going wrong here. The first thing is that either the minimum value for delta or this matrix have to be incorrect and have to be adapted. The second thing is that we do not seem to index shift correct when writing the code. We still point to the 7th element, but since have 0-based indexing in C that is a 0.7, but that could just be a debug message fault.

I hope this helps!

If the model still fails after you have updated these values i can have a look on why the index shift does not work, but as far as i can see it is just an error message fault and we actually use the correct values for computation.

[thorade]

I am double checking now if I typed the equation for non-analytcial terms correctly.
It is based on the book by Span 2000, table 3.8 (snippet below):

And on these two notebooks:
https://github.com/CoolProp/CoolProp/blob/master/doc/notebooks/Derivatives_of_NA_Term.ipynb
https://github.com/thorade/HelmholtzMedia/blob/master/docs/HelmholtzDerivatives/Helmholtz_energy_derivatives.ipynb

[kabdelhak]

I also had a short look at 4. (Examples.MediaTestModels.CarbondioxideTestModel error message during initialization: functions.c:2865: Invalid root: (-0.996066)^(1.33333). Maybe this error message can be enhanced to give more information?)

My analysis for issue 4. was a little off here, it wasn't a modeling error, but a simplification error during compliation!

• Fixed in: PR6714
• Further information: Ticket #6068

[thorade]

Thanks again Karim!
And thanks @casella for collecting all remaining issues here:
https://trac.openmodelica.org/OpenModelica/ticket/6088
OpenModelica/OpenModelica#6088

One more note for users interested in CO2: There are two equations of state available, the very accurate one, and a shorter equation of state without analytical terms.

[thorade]

This ticket has grown quite long, and some very nice improvements have been done.
These improvements are now available in the OpenModelica 1.16 release.
I am closing this ticket, and have created #39 instead to track the remaining issues.
Thanks all for the contributions!