Opened 4 years ago

Closed 4 years ago

Last modified 4 years ago

#6068 closed defect (fixed)

Illegitimate transformation during translation of HelmholtzMedia function

Reported by: casella Owned by: Karim.Abdelhak
Priority: blocker Milestone: 1.17.0
Component: Backend Version:
Keywords: Cc: vitalij, matthis.thorade@…

Description

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

Change History (18)

comment:1 Changed 4 years ago by matthis.thorade@…

subscribing

comment:2 Changed 4 years ago by casella

  • Cc matthis.thorade@… added

comment:3 Changed 4 years ago by 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.

comment:4 Changed 4 years ago by Karim.Abdelhak

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?

comment:5 follow-up: Changed 4 years ago by 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.

comment:6 Changed 4 years ago by Karim.Abdelhak

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.

comment:7 Changed 4 years ago by Karim.Abdelhak

  • Resolution set to fixed
  • Status changed from new to closed

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.

comment:8 in reply to: ↑ 5 Changed 4 years ago by Karim.Abdelhak

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.

comment:9 Changed 4 years ago by 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.

comment:10 Changed 4 years ago by matthis.thorade@…

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

comment:11 Changed 4 years ago by casella

  • Resolution fixed deleted
  • Status changed from closed to reopened

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?

comment:12 follow-up: Changed 4 years ago by matthis.thorade@…

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

comment:13 in reply to: ↑ 12 Changed 4 years ago by Karim.Abdelhak

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.

Last edited 4 years ago by Karim.Abdelhak (previous) (diff)

comment:14 Changed 4 years ago by Karim.Abdelhak

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!

Last edited 4 years ago by Karim.Abdelhak (previous) (diff)

comment:15 Changed 4 years ago by Karim.Abdelhak

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.

Last edited 4 years ago by Karim.Abdelhak (previous) (diff)

comment:16 follow-up: Changed 4 years ago by Karim.Abdelhak

  • Resolution set to fixed
  • Status changed from reopened to closed

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! :)

comment:17 Changed 4 years ago by casella

Thanks @Karim for the quick fix!

comment:18 in reply to: ↑ 16 Changed 4 years ago by 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.

Last edited 4 years ago by casella (previous) (diff)
Note: See TracTickets for help on using tickets.