search

OpenModelica / OpenModelica

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

NB simulations fail because of division by zero involving parameters #9414

Closed casella closed 1 year ago

casella commented 2 years ago

In the newInst-newBackend ScalableTestSuite report, there are many models failing because of division by zero involving parameters, e.g. CounterCurrentHeatExchangerEquations_N_10 or OneDHeatTransfer_TT_FD.

@kabdelhak, I guess that's because of some bug in code generation, I would suggest you to have a quick look.

casella commented 2 years ago

This seems to be related to the following error reported by the new backend in all those models:

Error: Internal error NBResolveSingularities.balanceInitialization failed because following
non-fixable variables could not be solved:
        [DER-] (10) Real[10] $DER.T
hkiel commented 1 year ago

A MWE for the above NewBackend error is the following model, which works fine for the old backend. It works again if you replace the second equation by y = der(sin(time));

model test
  Real x;
  Real y;
equation
  x = sin(time);
  y = der(x);
end test;
casella commented 1 year ago

@hkiel thanks for the really minimal MWE!

@kabdelhak I'd say the array connection algorithm has definitely nothing to do with this MWE. Can you please check and try to resolve this issue, so we can make some progress with testing?

Thanks!

casella commented 1 year ago

@kabdelhak, here is the regression report after your PRs #9953 and #9954. Note that most of the regressions are in Buildings_8 because Martin removed the reference files, that were not appropriate for that version of the library.

The results are encouraging, but we are still not there. The OneDHeatTransferTI_Modelica and OneDHeatTransferTT_Modelica models now simulate, but unfortunately the OneDHeatTransfer_TI_FD and OneDHeatTransfer_TT_FD still fail because of division by zero, and so still do all the heat exchanger models.

Furthermore, the simulation performance shows that there is some problem with wrong ODE Jacobians, similarly to #9843.

For example, ScalableTestSuite.Thermal.HeatConduction.ScaledExperiments.OneDHeatTransferTT_Modelica_N_320 has these simulation statistics with the new backend

LOG_STATS         | info    | ### STATISTICS ###
|                 | |       | | timer
|                 | |       | | |    0.0346164s          reading init.xml
|                 | |       | | |   0.00147384s          reading info.xml
|                 | |       | | |  0.000869611s [  0.0%] pre-initialization
|                 | |       | | |   0.00015479s [  0.0%] initialization
|                 | |       | | |    0.0440796s [  0.0%] steps
|                 | |       | | |      174.009s [ 90.2%] solver (excl. callbacks)
|                 | |       | | |    0.0939878s [  0.0%] creating output-file
|                 | |       | | |   0.00179789s [  0.0%] event-handling
|                 | |       | | |   0.00280491s [  0.0%] overhead
|                 | |       | | |      18.7898s [  9.7%] simulation
|                 | |       | | |      192.942s [100.0%] total
|                 | |       | | events
|                 | |       | | |     0 state events
|                 | |       | | |     0 time events
|                 | |       | | solver: dassl
|                 | |       | | | 467125 steps taken
|                 | |       | | | 935978 calls of functionODE
|                 | |       | | | 700610 evaluations of jacobian
|                 | |       | | |     0 error test failures
|                 | |       | | | 233533 convergence test failures

while it has these ones 

LOG_STATS         | info    | ### STATISTICS ###
|                 | |       | | timer
|                 | |       | | |    0.0433434s          reading init.xml
|                 | |       | | |   0.00263026s          reading info.xml
|                 | |       | | |   0.00109046s [  0.8%] pre-initialization
|                 | |       | | |  0.000323908s [  0.2%] initialization
|                 | |       | | |   1.9767e-05s [  0.0%] steps
|                 | |       | | |    0.0219177s [ 16.2%] solver (excl. callbacks)
|                 | |       | | |    0.0127131s [  9.4%] creating output-file
|                 | |       | | |  0.000248615s [  0.2%] event-handling
|                 | |       | | |   0.00105418s [  0.8%] overhead
|                 | |       | | |    0.0976644s [ 72.3%] simulation
|                 | |       | | |     0.135032s [100.0%] total
|                 | |       | | events
|                 | |       | | |     0 state events
|                 | |       | | |     0 time events
|                 | |       | | solver: dassl
|                 | |       | | |   193 steps taken
|                 | |       | | |   223 calls of functionODE
|                 | |       | | |    25 evaluations of jacobian
|                 | |       | | |     2 error test failures
|                 | |       | | |     0 convergence test failures

with the old backend. As explained in #9843, the problem here is that the Jacobian is not computed correctly, so the convergence of the BDF algorithm, which depends on (I - h*df/dx), is hampered if h is not very small, forcing the solver to take many very short steps. Please check the values of the jacobian with the smallest test and compare it to the new backend case.

kabdelhak commented 1 year ago

I now understand why this does not work. It requires index reduction for balanced initialization, which is not yet implemented in the new backend.

  x = sin(time);
  y = cos(time); 
  der(x) = y;

or 

  x = sin(time);
  der(x) = cos(time); 
  y = der(x);

is how this should be solved. The second equation is the result of differentiating the first one due to no unknown being present. I do have a work in progress solution for index reduction, but i have to clean it up to make future work regarding arrays easier.

casella commented 1 year ago

@kabdelhak I'm not sure I understand.

This is the OneDHeatTransferTI_FD model, which is one of the models failing at runtime:

model OneDHeatTransferTI_FD
  "One end at a fixed temperature, one end is insulated; implemented by FD method"
  import SIunits = Modelica.Units.SI;
  parameter SIunits.Length L "Length";
  parameter Integer N = 2 "number of nodes";
  parameter SIunits.Temperature T0 "Initial temperature";
  parameter SIunits.Temperature TN "Temperature at the last node";
  parameter SIunits.SpecificHeatCapacity cp "Material Heat Capacity";
  parameter SIunits.ThermalConductivity lambda
    "Material thermal conductivity";
  parameter Modelica.Units.SI.Density rho "Material Density";
  final parameter Modelica.Units.SI.Length dx=L/(N - 1) "Element length";
  SIunits.Temperature T[N] "temperature of the nodes";
initial equation
  for i in 1:N - 1 loop
    T[i] = T0;
  end for;
equation
  T[N] = TN;
  for i in 2:N - 1 loop
    der(T[i]) = lambda * ((T[i + 1] - T[i]) / dx + ((-T[i]) + T[i - 1]) / dx) / cp / rho / dx;
  end for;
  der(T[1]) = lambda * ((T[2] - T[1]) / dx) / cp / rho / dx;
end OneDHeatTransferTI_FD;

This is an index-1 model, in fact, almost index-0, except for the equation T[N] = T0. The state vector is T[1:N-1], and there are explicit initial equations setting their initial values to a parameter, and explicit equations to compute its derivative.

I don't understand why index reduction should be required here at all, nor why your example with sin(time) and cos(time) should be relevant. The only tricky point could be that we have a variable array, but only a slice of it is actually a state. Is that a problem?

The runtime error is

division by zero in partial equation: dx
|                 | |       | at Time=0.000000

it seems to me that for some reason the binding equation dx=L/(N - 1) is not computed before dx is used.

Are we on the same page?

kabdelhak commented 1 year ago

Are we on the same page?

My comments referred to the MWE by hkiel, i thought it represents the same problem.

Regarding the proposed model ScalableTestSuite.Thermal.HeatConduction.ScaledExperiments.OneDHeatTransferTT_Modelica_N_10 (smaller version): It is an older problem which i mentioned to you last summer. The new backend currently cannot handle mixed state arrays. Only index 1-9 of T are states and index 10 is not. The following model works fine with the new backend making use of dummy variables to represent the state slice. 

model OneDHeatTransferTI_FD
  "One end at a fixed temperature, one end is insulated; implemented by FD method"
  import SIunits =
         Modelica.Units.SI;
  parameter SIunits.Length L "Length";
  parameter Integer N = 2 "number of nodes";
  parameter SIunits.Temperature T0 "Initial temperature";
  parameter SIunits.Temperature TN "Temperature at the last node";
  parameter SIunits.SpecificHeatCapacity cp "Material Heat Capacity";
  parameter SIunits.ThermalConductivity lambda
    "Material thermal conductivity";
  parameter Modelica.Units.SI.Density rho "Material Density";
  parameter Modelica.Units.SI.Length dx=L/(N - 1) "Element length";
  SIunits.Temperature T[N] "temperature of the nodes";
  SIunits.Temperature X[N-1] "temperature of the nodes";
initial equation
  for i in 1:N - 1 loop
    T[i] = T0;
  end for;
equation
  T[N] = TN;
  for i in 1:N - 1 loop
    X[i] = T[i];
  end for;
  for i in 2:N - 1 loop
    der(X[i]) = lambda * ((T[i + 1] - T[i]) / dx + ((-T[i]) + T[i - 1]) / dx) / cp / rho / dx;
  end for;
  der(X[1]) = lambda * ((T[2] - T[1]) / dx) / cp / rho / dx;
end OneDHeatTransferTI_FD;

I did not yet decide on how to handle sliced array states in the new backend and we agreed on postponing that to when we handle index reduction since that is closely related to one another.

@casella would it be fine for now to update the scalable testsuite accordingly? and might this be the issue for other real models as well?

casella commented 1 year ago

Regarding the proposed model ScalableTestSuite.Thermal.HeatConduction.ScaledExperiments.OneDHeatTransferTT_Modelica_N_10 (smaller version): It is an older problem which i mentioned to you last summer.

I apologize, but there are really too many pending issues, and often they differ in intricate and subtle ways. We need to keep track of them in writing in these tickets 😅

The new backend currently cannot handle mixed state arrays. Only index 1-9 of T are states and index 10 is not. The following model works fine with the new backend making use of dummy variables to represent the state slice.

OK, good.

model OneDHeatTransferTI_FD
  "One end at a fixed temperature, one end is insulated; implemented by FD method"
  import SIunits =
         Modelica.Units.SI;
  parameter SIunits.Length L "Length";
  parameter Integer N = 2 "number of nodes";
  parameter SIunits.Temperature T0 "Initial temperature";
  parameter SIunits.Temperature TN "Temperature at the last node";
  parameter SIunits.SpecificHeatCapacity cp "Material Heat Capacity";
  parameter SIunits.ThermalConductivity lambda
    "Material thermal conductivity";
  parameter Modelica.Units.SI.Density rho "Material Density";
  parameter Modelica.Units.SI.Length dx=L/(N - 1) "Element length";
  SIunits.Temperature T[N] "temperature of the nodes";
  SIunits.Temperature X[N-1] "temperature of the nodes";
initial equation
  for i in 1:N - 1 loop
    T[i] = T0;
  end for;
equation
  T[N] = TN;
  for i in 1:N - 1 loop
    X[i] = T[i];
  end for;
  for i in 2:N - 1 loop
    der(X[i]) = lambda * ((T[i + 1] - T[i]) / dx + ((-T[i]) + T[i - 1]) / dx) / cp / rho / dx;
  end for;
  der(X[1]) = lambda * ((T[2] - T[1]) / dx) / cp / rho / dx;
end OneDHeatTransferTI_FD;

Excellent.

I did not yet decide on how to handle sliced array states in the new backend and we agreed on postponing that to when we handle index reduction since that is closely related to one another.

Sure. Quaternions in the MultiBody library are a nice example of that. Though I guess the case in question here is a really trivial one.

@casella would it be fine for now to update the scalable testsuite accordingly?

Absolutely. The point of ScalableTestSuite is to stress compilers with increasingly large models, not necessarily to check how well array-preserving models can handle complicated sliced matching and index reduction. BTW:

It may be useful to also have models to test sliced variable arrays, but that's a different purpose than ScalableTestSuite's.

I understand array-preserving index reduction is probably the most challenging topic in array-preserving code generation, and there are lots of very interesting problems (e.g. power grid models, thermal chip models) where this is really not necessary and the choice of states is trivial. It definitely makes sense to start testing them first without worrying about this issue.

I will update ScalableTestSuite myself, making sure that all the tests that do not involve index reduction do not have sliced

and might this be the issue for other real models as well?

I will check.

Thanks!

casella commented 1 year ago

After updating ScalableTestSuite to avoid arrays containing both state and non-state variables, this problem is solved both in heat transfer and heat exchanger models, see test report.