TwoPhase_setState fails on Openmodelica 1.18

[bilderbuchi]

We've encountered a bug that is also present in HelmholtzMedia.Examples.ConvergenceTest.TwoPhase_setState. See also #39, but I didn't want to pollute that one with details.

TwoPhase_setState fails in Openmodelica 1.18 (NF) on Windows 10 using either current HelmHoltzMedia master (18ff552) and MSL 3.2.3 or current msl4update (20bd884) and MSL 4.0.0

The error also happens on the OM 1.18 test server (also currently on the OM master):

Regular simulation: ./HelmholtzMedia_HelmholtzMedia.Examples.ConvergenceTest.TwoPhase_setState -abortSlowSimulation -alarm=480 -lv LOG_STATS
LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
assert            | debug   | HelmholtzMedia_HelmholtzMedia.Examples.ConvergenceTest.TwoPhase_setState_functions.c:11145: Invalid root: (nan)^(3)
stdout            | warning | Integrator attempt to handle a problem with a called assert.
assert            | debug   | HelmholtzMedia_HelmholtzMedia.Examples.ConvergenceTest.TwoPhase_setState_functions.c:11145: Invalid root: (nan)^(3)
stdout            | warning | Integrator attempt to handle a problem with a called assert.
assert            | debug   | HelmholtzMedia_HelmholtzMedia.Examples.ConvergenceTest.TwoPhase_setState_functions.c:11145: Invalid root: (nan)^(3)
stdout            | info    | model terminate | Simulation terminated by an assert at time: 10.0002

The test regressed at this point (also cited by @thorade in #39), and I suspect that this OM commit by @perost could be connected/the cause (old trac issue related to that).

The Modelica Specification issue modelica/ModelicaSpecification#2757 is probably connected but does not seem close to resolution.

Now, why all this exposition?
The error seems to occur, as far as we can tell, in HelmholtzMedia.Interfaces.PartialHelmholtzMedium.setState_pd, because

• here f.delta is NaN, with
• delta becoming NaN here during record constructions because it seems to see a d_crit of 0.
• d_crit maybe ends up being 0 a little above when
• Maybe MM is 0 here.

I'm not sure if those 0s occur because something is off in the ordering of binding equations when the function is evaluated, so a default/uninitialized value of 0.0 is used (hence the suspicion of OpenModelica/OpenModelica@14d191b423 above), or because something wonky happens when looking up fluidConstants[1].molarMass.

In any case, the minimal workaround that makes this test (and our internal model) work is to delay the assignment of f.delta:

protected
  //...
  EoS.HelmholtzDerivs f(d=d);
  //...
algorithm
  f.delta := d/d_crit;

This strengthens my suspicion that the ordering of operations could be off in OM.

Questions/remarks:

• Considering the spec issue above, maybe the desired behaviour of this corner of Modelica is not yet fully clear, maybe we should work around/avoid it in the meantime?
• Should I rather report this to OM to fix on their side (e.g. a tweak of OpenModelica/OpenModelica@14d191b423)?
• I'm not competent to assess the wider implications of this workaround, and if this has to be done elsewhere in the library, too. can anybody else contribute?

[perost]

I suspect that it shouldn't really have worked before my change either, it looks like we generate code that's using uninitialized variables so the behaviour is undefined. The actual issue looks to be that OM drops the modifiers on EoS.HelmholtzDerivs f(d=d, delta = d/d_crit);. If I flatten a toy example like this:

record R
  Real x;
  Real y;
end R;

function f
  input Real x;
  output Real y = r.x;
protected
  R r(x = x);
end f;

model M
  Real x = f(time);
end M;

I get:

function R "Automatically generated record constructor for R"
  input Real x;
  input Real y;
  output R res;
end R;

function f
  input Real x;
  output Real y = r.x;
  protected R r;
end f;

class M
  Real x = f(time);
end M;

As you can see the modifiers on f.r are gone. The issue is that modifiers are not a thing in the flat model, so we need to handle them in some other way. Our old frontend flattens the function to:

function f
  input Real x;
  output Real y = x;
  protected R r = R(x, <EMPTY(scope: <global scope>, name: y, ty: Real)>);
end f;

I.e. it generates a record constructor call to initialise z, but has to add some special empty expressions since not all record fields have values. I think having uninitialized variables is a bit dubious even if they're not used, but it seems the concensus is that it should be allowed (and it's used by other libraries too).

The old frontend doesn't work at all with HelmholtzMedia so I can't check if it simulates with it, but it seems very likely that this is the issue. So l'll try to implement something similar for our new frontend and see if that fixes it.

[perost]

After wasting some time implementing my proposal I realized it wasn't actually needed and we already do handle the modifiers, and the issue is actually my previous changes like you said.

The issue is that d_crit is used to initialize f, but my sorting of the local elements doesn't realize that because it doesn't look at the modifiers. So d_crit gets sorted below f, which is wrong. I didn't see that at first since we don't print the modifiers in the flat model, but looking at the generated code it's obvious(ish). I will improve the sorting to also take the modifiers into account, that should fix it for real.

[perost]

The issue should be fixed once my PR has been merged, the model now simulates successfully at least.

[thorade]

Thanks @bilderbuchi and @perost for the analysis. If any changes in the library are needed afterwards, I am also happy to try that. I think I tried to change this part before, and then it started to fail in another tool, but I can try again if needed.

[bilderbuchi]

Awesome, thank you for the (very!) quick fix, @perost !
I'll check the next CI library job when it becomes available to see how much this has improved things (and update #39).

[perost]

Awesome, thank you for the (very!) quick fix, @perost ! I'll check the next CI library job when it becomes available to see how much this has improved things (and update #39).

The report looks good, a few more models simulate and nothing broke (just some unstable models). The remaining models that don't simulate all fail with similar errors to this issue, so there might still be some issue with how we initialise variables. All the models that broke due to my earlier commit all work now though.

[bilderbuchi]

Thanks! I think this can be closed now.

[thorade]

Very nice, thank you!

[bilderbuchi]

All the models that broke due to my earlier commit all work now though.

Sorry to disappoint, but I iiuc the earlier commit (or at least one in that changeset) broke HelmholtzMedia.Examples.ConvergenceTest.SinglePhase_setState_psine_Tramp and HelmholtzMedia.Examples.ConvergenceTest.TwoPhase_setState which are still failing.

[perost]

All the models that broke due to my earlier commit all work now though.

Sorry to disappoint, but I iiuc the earlier commit (or at least one in that changeset) broke HelmholtzMedia.Examples.ConvergenceTest.SinglePhase_setState_psine_Tramp and HelmholtzMedia.Examples.ConvergenceTest.TwoPhase_setState which are still failing.

They show up as working for me, but I had to manually refresh the page with F5 for some reason. Maybe you're having the same issue. There should be 5 models still failing during simulation.

[bilderbuchi]

There should be 5 models still failing during simulation.

exactly, and those two are 2 of these 5 from what I see (after even Ctrl-F5). They dropped from simulate->compile in the earlier commit. I'm confused now?

[perost]

There should be 5 models still failing during simulation.

exactly, and those two are 2 of these 5 from what I see (after even Ctrl-F5). They dropped from simulate->compile in the earlier commit.

The ones I see that are still failing are:

HelmholtzMedia.Examples.ConvergenceTest.Ancillary_Saturation 
HelmholtzMedia.Examples.ConvergenceTest.SinglePhase_setState_pramp_Tsine 
HelmholtzMedia.Examples.ConvergenceTest.setSat
HelmholtzMedia.Examples.Validation.Derivatives_SaturationBoundary
HelmholtzMedia.Examples.Validation.idealGasLimit

The simulation log for SinglePhase_setState_psine_Tramp shows some limits being violated (though it seems to just be due to floating point precision issues) but the simulation is successful in the end. TwoPhase_setState on the other hand simulates without any warnings.

[bilderbuchi]

pramp_Tsine != psine_Tramp, setSat != setState .. my apologies for wasting your time, all good! Clearly, I should stop working for today 😿