Issue with the derivative annotation in the NF produces grossly incorrect results

[Francesco Casella]

Please check ​Modelica.Fluid.Examples.TraceSubstances.RoomCO2. The simulation fails after about 20 seconds. Analysis of the transient reveals that the temperature in the volume unexpectedly increases above 300 degC, while it is supposed to remain roughly constant throughout the transient, which mainly involves the concentration of CO2.

Further analysis of the -d=optdaedump output reveals that in the process of changing the state variables from volume.m, volume.U to volume.medium.p, volume.medium.T, carried out via the Pantelides algorithm, when the flattening is carried out by NF, the backend computes

boundary1.ports[2].h_outflow = 
  Modelica.Fluid.Examples.TraceSubstances.RoomCO2.Medium.specificEnthalpy_pTX(
  volume.medium.p, 293.15, {0.01})
der(volume.ports[2].h_outflow) = 0.0

This is plain wrong and causes the derivatives of pressure and temperature to be computed in a grossly incorrect way, leading to a totally unphysical transient.

When the OF is used, the Pantelides algorithm instead correctly computes

boundary1.ports[2].h_outflow = 
  Modelica.Fluid.Sources.MassFlowSource_T$boundary1.Medium.specificEnthalpy_pTX(
    volume.medium.p, 293.15, {0.01})
der(volume.ports[2].h_outflow) = 
  Modelica.Fluid.Vessels.ClosedVolume$volume.Medium.h_pTX_der(
  volume.medium.p, volume.medium.T, {volume.medium.Xi[1], 
  1.0 - volume.medium.Xi[1]}, der(volume.medium.p), 
  der(volume.medium.T), {der(volume.medium.Xi[1]), 
  -der(volume.medium.Xi[1])})

and eventually produces the correct simulation results.

It is apparent that

• the enthalpy function is flattened differently in the two cases
• its derivative is correctly computed when starting from the OF, while it is computed to zero when starting from the NF.

Unfortunately the flat output is stripped of annotations, so I can't see if the Derivative annotation was flattened correctly, but I guess there's some problem with it.

@perost, can you please investigate and report?

[Adrian Pop]

Run omc with --showAnnotations to display them.

[Per Östlund]

The issue turned out to be because with the NF we get (abbreviated for readability):

volume.ports[2].h_outflow = 
  specificEnthalpy_pTX(volume.medium.p, volume.medium.T, volume.medium.Xi)

while the OF gives:

volume.ports[2].h_outflow =
  specificEnthalpy_pTX(volume.medium.p, volume.medium.T, {volume.medium.Xi[1]})

That is, the OF expands the last argument in the call while the NF doesn't. This was then inlined and the volume.medium.Xi argument ended up in a cat-call, which the old ExpressionSimplify-module used by the backend failed to simplify when using the NF since it expected the arguments to be expanded. I didn't follow the cascade of failure any further than this, but somehow this made the derivative become zero in the end.

I made a fix in 5262ff52 that tries to expand one-dimensional crefs in ExpressionSimplify.simplifyCat. The reason why I do it only for one-dimensional crefs is because the expansion functions in the OF are very primitive compared to the NF, and I'd have to port a lot of NFExpandExp to the OF to be able to handle more complex expressions.

The fix seems to work for RoomCO2 when using the NF at least, but we'll see if messing with the old expression simplification breaks something or not.

[Per Östlund]

Btw, the backend doesn't use the derivative annotations directly. They're parsed by the frontend and the relationships between functions and their derivatives are encoded in the DAE, which is what the backend uses. This information doesn't show up in the flat model.

They will still show up in the flat model when using --showAnnotations though, since the frontend forwards all annotations to the backend. But at that point they've already served their purpose.

[Francesco Casella]

Replying to adrpo:

Run omc with --showAnnotations to display them.

Aha, I missed that. It would be nice to have a flag to only show annotations that are relevant for structural analysis and numerical solution. As far as I remember, they would be: Evaluate, Inline, LateInline, derivative, inverse, but that's kind of extra luxury. See #5299

[Francesco Casella]

Replying to perost:

Btw, the backend doesn't use the derivative annotations directly.

Sure

They're parsed by the frontend and the relationships between functions and their derivatives are encoded in the DAE, which is what the backend uses. This information doesn't show up in the flat model.

They will still show up in the flat model when using --showAnnotations though, since the frontend forwards all annotations to the backend. But at that point they've already served their purpose.

In fact, I'm always afraid that some annotation is lost, but I understand that's not the point, because there's no reason they should. It would be very well possible that the information in the annotation is not encoded in the DAE structure correctly, and I understand there is currently no way to have explicit debug info about that. In fact, I only saw the end result (after some searching), i.e. a derivative erroneously computed as zero.

Let's see what happens on the testsuite. I guess other Media/Fluid models will benefit, I hope this doesn't break other stuff.

2D array expansion shouldn't be necessary, I'm not aware of any model that would need that.

[Per Östlund]

There were some issues with my first attempt at the fix, because the old frontend thought a cref with type Real[:] should be expanded into {}. But the issue should hopefully be fixed now in 2c3fa1d.

[Francesco Casella]

Indeed it is, see ​report. Some other models have benefitted as well.