Illegitimate transformation during translation of HelmholtzMedia function

[Francesco Casella]

As discussed ​here, we have some issues in the HelmholtzMedia library that can be ultimately traced to a symbolic transformation carried out in a function, namely

a^b -> exp(log(a)*b)

The Modelica Specification, ​Section 10.6.7 states:

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.

So, from what I understand, it is perfectly legitimate to compute a^b, where a = -2.5 and b = 2.0 are two Real numbers, as the ANSI C function can recognize the even integer exponent and handle this correctly. And in fact, even the case of Real a = -2.5 and Integer b = 2 is first turned into the previous case, at least conceptually, before further processing.

Hence, the above-mentioned symbolic transformation is in general invalid and should be avoided, unless min(a) >= 0, that is.

This issue affects models of HelmholtzMedia that are currently failing for this reason, see e.g. ​HelmholtzMedia.Examples.MediaTestModels.ButaneTestModel_dT_component_ph

[matthis.thorade@…]

subscribing

[Francesco Casella]

In the meantime, a ​workaround was discussed, namely to add integer to the exponent to make sure that this transformation is not applied. But it would be good to fix it anyway.

[Karim Adbdelhak]

I dug deeper into this and found out that the actual call for this log that causes the error is in

$DER.HelmholtzMedia.Examples.MediaTestModels.HeliumTestModel.volume.Medium.EoS.f_r`

which is just the time derivative of

HelmholtzMedia.Examples.MediaTestModels.HeliumTestModel.ambient.Medium.EoS.f_r

The actual function name in the end should be

$DER$HelmholtzMedia$PExamples$PMediaTestModels$PHeliumTestModel$Pvolume$PMedium$PEoS$Pf_r

due to c-compliant character replacement.

It is not a simplification but just a differentiation process. The original functions stay the same. The rules for this differentiation are quite simple.

(x^y)' = x^(y-1) * ( x*ln(x)*y'+(y*x'))

This produces the logarithmic function. If the exponent is constant on the other hand, we can use a simplified version:

(x^r)' = r*(x^(r-1))

The new frontend does not evaluate the exponents (which are parameters) in this case, so the first rule has to be applied. The old frontend did evaluate them so it worked. I suppose we could use this rule for constants also for parameters? Is there an argument against it for parameters that could change during simulation?

[matthis.thorade@…]

Thanks for looking into this!
Here, the coefficients and exponents are even declared as constant, see
​https://github.com/thorade/HelmholtzMedia/blob/master/HelmholtzMedia/Interfaces/PartialHelmholtzMedium/EoS/f_r.mo

Whether the simplification should be applied for parameters also I am not so sure. They do not change during simulation by definition, to my understanding, but can change between simulation runs without recompilation.

[Karim Adbdelhak]

I implemented two new rules for exponential differentiation with parameters involved in ​PR6685.

The changes fix the problem described in this ticket, but jenkins reports 6 failing tests related to jacobian differentiation. I suppose for partial differentiation the rules have to be adapted. For now i will exclude these rules from partial differentiation and try to push the changes. Later on we can think about the jacobians and if the new rules are applicable in that case.

[Karim Adbdelhak]

I had to adapt some structures such that the known variables are correctly passed in the previous mentioned pull request, but everything works fine now. I also did not have to exclude jacobians, the rule is now always applied.

I added a minimal model to the testsuite that represents the problem:

model FuncLocalConstantsDer
  Real x,y(stateSelect = StateSelect.prefer);
  function f
    input Real x;
    output Real y;  
  protected
    final constant Real[2,2] r = fill(2.0, 2, 2);
  algorithm
    y := x^r[1,2];
  end f;
  equation
  x = sin(time);
  y = f(x);
end FuncLocalConstantsDer;

(The StateSelect attribute forces the compiler to try and differentiate the function).

This does not work with the current OMC, but it runs with my changes (also the constant could be a parameter).

The changes are pushed. I tested it with HelmholtzMedia.Examples.MediaTestModels.HeliumTestModel, but since they are all structurally equal i suppose all of them work now.

[Karim Adbdelhak]

Replying to matthis.thorade@…:

Whether the simplification should be applied for parameters also I am not so sure. They do not change during simulation by definition, to my understanding, but can change between simulation runs without recompilation.

I fiddled around with it and i don't see a problem. I actually do not replace any constants or parameters by their values, i just assume that they are not time dependent and can therefore apply the better differentiation rules without using the natural logarithm. It would still not work if the base becomes negative and the exponent is non integer (not by type, just by value), but if that is a realistic behavior of the model, it has to be modeled with complex numbers. Otherwise it just does not have a solution.

For all that i can see it should also work for discrete variables since they are (locally) independent of time. But that is another story.

[Francesco Casella]

Thank @Karim for the analysis and the quick fix!

I guess this issue was quite pervasive in HelmholtzMedia, so I hope there will be a significant improvement in the ​test report after the next run. I'm leaving for some holidays tomorrow, so I won't follow this for a while.

Please close this ticket if the test report has improved in a couple of days and this kind of issues no longer take place.

[matthis.thorade@…]

Next test has been run, and all 6 cases with this error message now simulate successfully.

[Francesco Casella]

Hooray! :)

There is still a related issue in ​HelmholtzMedia.Examples.MediaTestModels.CarbondioxideTestModel, which fails with

CarbondioxideTestModel_functions.c:2829: Invalid root: (-0.996066)^(1.33333)

@Karim, would you mind having a look?

[matthis.thorade@…]

Some first investigation for non-analytic terms is here
​https://github.com/thorade/HelmholtzMedia/issues/34#issuecomment-667541330

[Karim Adbdelhak]

Replying to matthis.thorade@…:

Some first investigation for non-analytic terms is here
​https://github.com/thorade/HelmholtzMedia/issues/34#issuecomment-667541330

These investigations were done by myself, but i will try to revisit what i did and try to shortly summarize what i think is going on. :)

The error once again happens in ​HelmholtzMedia.Interfaces.PartialHelmholtzMedium.EoS.f_r and revolves around the binding from Theta:

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

The crucial part is ((delta-1)^2)^(0.5/a[i,7]) which can be simplified to (delta-1)^(1/a[i,7]). Looking up the value for a[i,7] yields 0.3 for all i. As far as i can see the actual error reporting value (-0.996066)^(1.33333) could be wrong, since 1/0.3 = 3.3333 <> 1.3333, but i might be missing some simplifications here. Also the actual error is thrown from a partially differentiated version of this function, so that might be the origin for this discrepancy.

Nevertheless the same error should occur in the original function since delta has a minimum value of 0 defined in the function itself. That allows the base to become negative and since the exponent is not of integer value, this error occurs.

I only see four options for the origin of this error:

1. a[i,7] = 0.3 is an incorrect value (defined in ​HelmholtzMedia.HelmholtzFluids.Carbondioxide as residualNonAnalytical)
2. delta(min = 0.0) is an incorrect lower boundary
3. delta is never supposed to actually be lower than 1.0 but due to issues during solving of algebraic loops it temporarily gets assigned a value lower than 1.0.
4. The result of this is actually allowed to become complex. In that case you should use the Complex record from the MSL.

I guess 2. and 3. are pretty much connected.

Please check if any of these could be the case, otherwise i don't know how the simulation could work since we do not have complex numbers as a primitive type in any lower target language.

[Karim Adbdelhak]

I actually found the error i made and presumably the compiler also makes the same mistake.

It is not

(a^2)^(0.5) = a

but rather

(a^2)^(0.5) = |a|

The base never gets negative, seems to be a simplification error. I will try to fix this!

[Karim Adbdelhak]

Indeed it is that very simple error. Easy to reproduce with following example:

model M
  Real x;
equation
  x = (sin(time)^2)^0.5;
end M;

That simulates a perfect sine curve, but it should simulate the absolute value of a sine curve. I am a little stunned that this error exists.

I fixed it in ​PR6714, that fix also fixes your carbon dioxide model. If the tests run fine (which i assume) i will merge it and close this ticket again.

[Karim Adbdelhak]

​PR6714 is pushed and the second issue is fixed!

If you have new issues with the Helmholtz library please create a new ticket. We do not want to get this one crammed, might be confusing later on! :)

[Francesco Casella]

Thanks @Karim for the quick fix!

[Francesco Casella]

Replying to Karim.Abdelhak:

​PR6714 is pushed and the second issue is fixed!

If you have new issues with the Helmholtz library please create a new ticket. We do not want to get this one crammed, might be confusing later on! :)

#6088 collects remaining all remaining issues with the library. We may open specific tickets for each of them, though it is probably not necessary to do so a priori.