Nasty bug with package constants in Modelica.Media gas mixture models

[Francesco Casella]

Consider this test case:

model TestMixture
  package Medium = Modelica.Media.IdealGases.MixtureGases.FlueGasSixComponents;
  Medium.BaseProperties medium;
  Medium.MolarMass MM[:] = Medium.data.MM;
equation 
  medium.p = 1e5;
  medium.T = 300;
  medium.X[1] = time;
  medium.X[2] = 1 - time;
  medium.X[3:6] = zeros(4);
end TestMixture;

The record Medium.data.MM should contain the molar masses of the six components of the gas. Unfortunately, all six components of the array turn out to be equal to the molar mass of the last component (CO2 in this case).

This result is obtained using r22544 of OMC. I have been reported that if older versions of OMC with the previous Inst module were used, the value of the first component turned up, though I haven't checked this myself.

Apparently something goes wrong when the modifier

     data={Common.SingleGasesData.N2,
    Common.SingleGasesData.H2,
    Common.SingleGasesData.CO,
    Common.SingleGasesData.O2,
    Common.SingleGasesData.H2O,
    Common.SingleGasesData.CO2}

is expanded in the medium model.

This is a particularly nasty bug, because the compiler does not produce any error, but simple gives wrong numerical results. Basically, all models using Modelica.Media with ideal gas mixtures are invalid.

I would strongly urge to fix this issue before releasing 1.9.1. - if that is not possible, then at least the compiler should give an error when trying to compile these models, instead of producing wrong simulation code.

Thanks to Andrea Brambilla of ETH Zurich for reporting this problem.

[Martin Sjölund]

I get the following in the model: MM = {0.0440095, 0.0440095, 0.0440095, 0.0440095, 0.0440095, 0.0440095};. What is the expected results?

[Francesco Casella]

import Data = Modelica.Media.IdealGases.Common.SingleGasesData;
Real MM[:] = 
  {Data.N2.MM, Data.H2.MM, Data.CO.MM, Data.O2.MM, Data.H2O.MM, Data.CO2.MM};

[Adrian Pop]

Yes, it should be:

MM = {0.00201588, 0.0280101, 0.0319988, 0.01801528, 0.0440095};

It seems that accessing components of record constants put in an array is a bit off.

[Adrian Pop]

Easiest model to debug with:

package Constants
  record Data
    Real MM;
  end Data;
  
  constant Data C1 = Data(MM=1);
  constant Data C2 = Data(MM=2);
  constant Data C3 = Data(MM=3);
  constant Data C4 = Data(MM=4);
  
  constant Data data[:] = {C1, C2, C3, C4};
 
  constant Real x[:] = data.MM;
end Constants;

model Test
  parameter Real x[:] = Constants.x;
end Test;

Flattening wrongly gives:

adrpo@ida-liu050 ~/dev/OpenModelica/build/bin/media
$ ../omc ConstantRecordArray.mo 
function Constants.Data "Automatically generated record constructor for Constants.Data"
  input Real MM;
  output Data res;
end Constants.Data;

class Test
  parameter Real x[1] = 4.0;
  parameter Real x[2] = 4.0;
  parameter Real x[3] = 4.0;
  parameter Real x[4] = 4.0;
end Test;

[Adrian Pop]

Found where the problem is (Lookup.lookupVarF) and looking for a way to fix it.

[Adrian Pop]

Ouch, we treat constant arrays or records so badly is incredible we haven't detected this before.
If you use parameters then is fine as the bindings in the DAE are correct but the bindings we have
in the environment are totally off and when we handle constants we use those.

I partially fixed it so the constant in the model

Medium.MolarMass MM[:] = Medium.data.MM;

gets the correct value:

MM = {0.0280134, 0.00201588, 0.0280101, 0.0319988, 0.01801528, 0.0440095};

However the problem still remains in functions (or when indexing data with a non constant index):

When we expand this:

MM := 1/sum(state.X[j]/data[j].MM for j in 1:size(state.X, 1));

in function molarMass, is still wrong:

function TestMixture.Medium.molarMass "Inline before index reduction" "Return molar mass of mixture"
  input TestMixture.Medium.ThermodynamicState state "Thermodynamic state record";
  output Real MM(quantity = "MolarMass", unit = "kg/mol", min = 0.001, max = 0.25, nominal = 0.032) "Mixture molar mass";
algorithm
  MM := 1.0 / (state.X[1] / 0.0440095 + state.X[2] / 0.0440095 + state.X[3] / 0.0440095 + state.X[4] / 0.0440095 + state.X[5] / 0.0440095 + state.X[6] / 0.0440095);
end TestMixture.Medium.molarMass;

If you look on how Dym handles this function (or constant arrays of records) you go crazy.

Because in Modelica the type of the array elements is the same except for the binding
we don't actually keep all the elements of the array after we build the DAE, we keep just
one, that's why we get the same binding whatever index we use in data[NN].component.
This problem appears because constants do not generate a visible DAE and the binding
we have internally is just for one element of the array.

The entire array data[:] does have the correct binding but we do not use it when
we go and find components inside it and then we find the same component with the same
binding. Is not very clear to me how to properly fix this in an easy way. The only way
I can think of is to never access anything below data with dot notation but by
using some sort of record indexing:

data[:].MM -> {RECORD_INDEX(data[j], "MM") for j in 1:size(data, 1)}

This way we never go below data and only use the correct binding.

Martin and Per, any suggestion about this?

[Martin Sjölund]

Wouldn't it make more sense to use the index and dot operators for the record indexing? And just re-write the bad parts or remove the bad information altogether? So we always start from the data[:] during the lookup/ceval or whatever?

[Per Östlund]

Things like this still does not work either (testsuite/flattening/modelica/scodeinst/const4.mo):

class P
  constant Integer i = 3;
end P;

model A
  P p[2](i = {1, 2});
  Real x[2] = p.i; // Flattening gives x = {1.0, 1.0}.
end A;

I don't remember the details anymore, but I looked into this a long time ago and couldn't come up with a solution for the current instantiation. In scodeInst it's not an issue since the modifier isn't split until after instantiation.

[Adrian Pop]

Per, I thought about not splitting the modifier as you do in scodeInst
which is quite nice but it doesn't work for constant arrays of records
or any arrays of structures.
Consider:

constant Data data[:] = {C1, C2};

How would you handle this? push the binding inside Data or just keep
it at this level and then the components of data have no binding?

You don't split the modifier and you go with BINDING( {C1, C2} propagated dim 1 ) inside Data.
What do you do with it as suddenly you have inside Data component M and record binding C1
and C2 for it.

So you need to do record indexing when you push a structure modifier into
its class. So binding should have

BINDING( {C1.MM, C2.MM}, propagated dim 1).

But how should the binding look if you have more than one component in the
structure? You would need to unfold C1 and C2 to their components:

BINDING( { MM = C1.MM, OTHER = C1.OTHER, .... } propagated dim 1 ).

and if you have arrays of records of records of arrays this gets quite complicated.

[Adrian Pop]

I think the easiest way to fix this is to not do constant evaluation (well except for structural parameters). This will actually help a lot with other things we try to evaluate during flattening and we shouldn't.
We should generate DAE for constants and use the component references everywhere and not their bindings. I.e. data.MM stays the same (maybe vectorized and prefixed with the proper prefix) but not constant evaluated.

Because this works fine:

class P
  parameter Integer i = 3; // not a constant
end P;

model A
  P p[2](i = {1, 2});
  Real x[2] = p.i;
end A;

flattening gives:

class A
  parameter Integer p[1].i = 1;
  parameter Integer p[2].i = 2;
  Real x[1];
  Real x[2];
equation
  x = {/*Real*/(p[1].i), /*Real*/(p[2].i)};
end A;

We could constant evaluate things on the DAE between front-end and back-end but not in the environment.

[Martin Sjölund]

The "except for structural parameters" part means you need to consider those cases anyway, so it does not help.
You need to use the constants for example in functions that you need to interpret for structural parameters.

[Adrian Pop]

Right :) Even thou' I can't think of a structural parameter depending on a binding of an structured array element that doesn't mean it cannot exist.

Then the only other way I can think of is to use the binding of data whenever possible and when not possible generate DAE and leave the component reference not constant evaluated.

package Constants
  record Data
    Real MM;
  end Data;
  
  constant Data C1 = Data(MM=1);
  constant Data C2 = Data(MM=2);
  
  constant Data data[:] = {C1, C2};
 
  constant Real x[:] = data.MM;
end Constants;

model Test
  parameter Real x[:] = Constants.x;
end Test;

Currently when we instantiate data we instantiate each array element with the correct modification (generate the DAE) and only keep one in the component environment (the first or
the last in the array with its binding).

Lookup.lookupVar on data.MM lookups data and ignores the correct binding and then lookups MM inside the component environment of data and returns its wrong binding.

I changed Lookup.lookupVar to use the binding from data and that works fine for constant
indexes or whole dimensions but not for iterators. When you get data[j].MM you cannot fail
and you cannot index in the binding of data with it and you cannot return the binding of MM as is incorrect. What I did is return itself data[j].MM as binding which is correct
but then you get references to constants that do not have any DAE in the flattened code.

So we actually need to *always* generate DAE from constants that we cannot constant evaluate.
Combined with my Lookup.lookupVar fix it should work fine.

Ok. I'll try this now and see where I get.

[Adrian Pop]

Partial fix for this in r22588.
Does not yet work with for iterators indexing the record array.
I will do the remaining part later.

[Francesco Casella]

Re-tried the test case reported initially with r22604. The array MM is now computed correctly, and so is the density medium.d. However, the mixture molar mass medium.MM is still wrong, as it is computed with the function Modelica.Media.IdealGases.Common.MixtureGasNasa.molarMass that uses an iterator on the array.

This function is inlined, I'm not sure whether this is relevant for the problem.

[Adrian Pop]

That the function is inlined is not relevant. I'm working of fully fixing this bug. Hopefully will be ready this week as it requires some changes on how bindings are handled in the compiler.

[Francesco Casella]

The TextMixtureLoop model above computes the molar mass in four different ways. Currently, two of them work fine, two don't.

[Martin Sjölund]

This ticket was not closed for 1.9.1, which has now been released. It was batch modified for milestone 1.9.2 (but maybe an empty milestone was more appropriate; feel free to change it).

[Martin Sjölund]

Milestone changed to 1.9.3 since 1.9.2 was released.

[Francesco Casella]

Model with four different ways of using molar masses

[Francesco Casella]

The problem with package constants in ideal gas mixtures being wrongly evaluated is still sitting there. I would urge to find a fix to this issue. Getting models to compile and run with plain wrong values is really bad, as it undermines the trust one can put on the compiler.

[Martin Sjölund]

Moved to new milestone 1.9.4

[Martin Sjölund]

Milestone pushed to 1.9.5

[Martin Sjölund]

Milestone renamed

[Francesco Casella]

It seems that we'll only get this with the new frontend, so changing the milestone

[Christoph <buchner@…>]

What new frontend? I couldn't find any discussion on this.
Looking at comment:19, isn't this too big an issue to leave until 2.0? When is that planned to arrive?

[Francesco Casella]

Adrian Pop and Per Östlund have been working for a while on a completely new front-end, which should be both faster and more comprehensive as far as coverage is concerned than the current one. There is a consensus among the developers on the fact that trying to fix this kind of bugs in the current front-end is like trying to patch an old rusted bucket with many holes and fissures, i.e., a waste of time.

According to the plans the new front-end should be available on the master branch before the end of this year. Significant core OSMC resources are currently being allocated to this effort.

The new front-end will live in parallel with the old one for some time, in order to conduct enough beta testing and iron out all the nasty details. As soon as the new front-end is available and well debugged, I expect that most blocker tickets for 2.0.0, which stem from erroneous handling of the code by the front-end, will be solved automatically. My understanding is that 2.0.0, for which the requirement is full coverage of industrially relevant libraries, will need it.

Version 1.10.0 will be released earlier, in order to include all the improvements made so far w.r.t. 1.9.6, in particular the 64-bit Windows version and the many improvements on scalability of the back-end.

Adrian will post a ticket on this topic shortly.

[Christoph <buchner@…>]

I see, thank you for the clarification/info, Francesco!

[Francesco Casella]

I am confident that the new front-end can handle this correctly. Unfortunately, because of #5010 it is still not possible to run this model with it.

I'll get back to this once #5010 is fixed.

[Francesco Casella]

#5010 is now fixed, but we still have a backend issue, see #5200. However, that is not related to the problem originally reported here.

[Francesco Casella]

For the record, after #5200 was fixed, the test model originally reported in this ticket now works flawlessly, as long as the new front-end -d=newInst is used. This will become the default from 2.0.0, but it should already be usable with 1.14.0, possibly even with 1.13.0.