OpenModelica / OpenModelica

OpenModelica is an open-source Modelica-based modeling and simulation environment intended for industrial and academic usage.
https://openmodelica.org
Other
820 stars 305 forks source link

Avoid equations in the form record = function(var1, var2, ...) to end up in the residual equations when applying tearing #11823

Open casella opened 8 months ago

casella commented 8 months ago

When using Modelica.Media, a common modelling pattern is to write something like

  Medium.ThermodynamicState state;
  Medium.Pressure p;
  Medium.Temperature t;
  ...
equation
  state = Medium.setState_pT(p, T)
  ...

When solving initial equations, sometimes OMC tears implicit equations in such a way that this equation ends up in the residual equations set. This is never a good idea, particularly when using ExternalMedia, which has many variables in the state record, including an Integer one, which makes it impossible to solve the residual (outer) equations with the Newton algorithm.

Hence, it would be nice if the tearing algorithms had some provision to avoid this situation. This may actually turn out to be a difficult problem, in which case an alternative heuristic would be to always select as tearing (iteration) variables the arguments of functions that return records. I guess this would more or less guarantee that this equation ends up in the torn (inner) part, not in the residual part.

See also #11564; the proposal here also deals with purely static models, e.g. for process design.

It would be nice to have this in the current backend. For sure, we should have it in the new one.

Keeping @FedericaAscione in the loop.

casella commented 8 months ago

Discussion with @phannebohm: in fact, selecting p and T as states is not by itself a guarantee that state := setState_pT(p, T) is picked as a torn/inner equation. If there is one equation that can be solved for an element of state and it gets solved before in the sequence of torn equations, then you can't put state := setState_pT(p, T) in the torn section anymore.

In fact, this is likely the case here. The ThermodynamicState state contains state.p and state.T and there are equations in the model stating p = state.p and T = state.T. If you pick p and T as tearing variables, state.p := p and state.T := T may end up being solved first, and then you can't have state := setState_pT(p, T) as torn equation because some variables of state have been solved for already.

phannebohm commented 8 months ago

That's correct! We can control the selection of iteration variables. But selecting inner equations is done by running a matching algorithm. We could try to prefer record equations of that form if given a choice. But I'm not sure how easily the matching can be adapted for that.

For the new backend we want to do Tearing differently, by explicitly selecting inner variables together with their equations according to some rules like the one in this issue, because often we can know which equation to choose and which not to.

casella commented 8 months ago

In any case, 0 = state - setState_pT(p,T) should be avoided at all costs, because it is really a bad tearing choice, and it doesn't work for media that have an Integer phase in the state record.

I'm not really sure how exactly records and parts thereof are represented in the backend, but one idea could be, before starting the tearing algorithm:

In this way, if the input arguments of the function are all marked as tearing variables via the annotation, the fact that state = setState_ph(p, h) lands in the torn/inner section would be guaranteed.

@phannebohm, @kabdelhak, what do you think? Could this be implemented in the old backend without too much effort?

casella commented 8 months ago

For the new backend we want to do Tearing differently, by explicitly selecting inner variables together with their equations according to some rules like the one in this issue, because often we can know which equation to choose and which not to.

As I wrote above, if some variables are coming from a record that appears in an equation record = function(...), you can avoid matching to subsets of it - either the whole thing or nothing.