Backend issues with ThermalSeparation library

[Per Östlund]

Most of the ThermalSeparation library now flattens using the new frontend (-d=newInst), but the majority of the models fail in the backend due to various issues.

Most of the failing models gives errors such as (this one is from Examples.Absorption_IdealGases):

BackendDAEOptimize.makeEquationToResidualExp failed to transform equation: column1.balanceEquations.bool_eps[1] = false to residual form!

and seem to finally fail completely with the error:

Internal error function createNonlinearResidualEquations failed

The library can be found at: ​https://github.com/thom-marx/ThermalSeparation
The current status of the library can be found at: ​https://libraries.openmodelica.org/branches/newInst/ThermalSeparation/ThermalSeparation.html

[Karim Adbdelhak]

The problem is, that there are if equations which get processed as residual equations and some equations which cant get processed as residual equations (boolean equations), in the same condition. A minimal example is the following:

model minimal_bool_if
Boolean b;
Real x;
equation
if time>0.5 then
 b = false;
 x = cos(time);
else
 b = true;
 time = x;
end if;
end minimal_bool_if;

the second equation of the else part, time = x, will be reformed to a residual equation, because x is not on the lhs. Maybe we generally want to make a better approach here? The problem here is, that now all other equations in this if equation will get reformed to residual form, which is not possible. I changed that in my local branch such that each equation pair will get processed individually, but that requires a uniform order inside the if equation to always work as intended.

For the given model i am not quite sure if this still is the problem since it is quite complex, the new error message is as follows:

Error: Internal error function createNonlinearResidualEquations failed
[/home/karim/workspace/OpenModelica/OMCompiler/Compiler/SimCode/SimCodeUtil.mo:3335:7-3335:48:writable] Error: Internal error createOdeSystem2 failed for Generic Jacobian via directional derivatives

Next i will try to generally skip the simplification and see if that works.

[Francesco Casella]

Replying to Karim.Abdelhak:

What about asking the library developers to rewrite

model minimal_bool_if
  Boolean b;
  Real x;
equation
  if time>0.5 then
    b = false;
    x = cos(time);
  else
    b = true;
    time = x;
  end if;
end minimal_bool_if;

as

model minimal_bool_if
  Boolean b;
  Real x;
equation
  b = if time > 0.5 then false else true;
  x = if time > 0.5 then time else cos(time)
end minimal_bool_if;

?

[Karim Adbdelhak]

Replying to casella:

Willi had the same idea, maybe it is a good idea to give out the general advice to avoid mixing types in if cases, or to avoid multiple equations in if cases in general.

I will try it with one small model, since this is only a sub-problem. Mainly the residuals for the jacobian of non-linear systems cannot be generated, quite possibly there is more to it.

[Francesco Casella]

I mean, if a variable-structure model has some structure, one should try to help the tool exploit it by making it explicit in the Modelica source code.

Ideally, the back-end should be able to do that automatically, but maybe that's asking too much in general.

[Francesco Casella]

Added Thomas Marx-Schubach and Andreas Heuermann in cc:

[Francesco Casella]

@thomas.marx, in your e-mail you mention that ThermalSeparation.Components.HeatExchanger.CcounterFlowHeatExchanger is a single component, so it cannot be simulated because it's singular. In this case, please remove the experiment annotation from the source code, because it doesn't make sense.

[Willi Braun]

The model ThermalSeparation.Examples.Absorption_IdealGases starts to simulate with the following script and the changes in the the ​PR2.

loadFile("../ThermalSeparation/package.mo");
getErrorString();

setCommandLineOptions("-d=newInst,-NLSanalyticJacobian");
getErrorString();

simulate(ThermalSeparation.Examples.Absorption_IdealGases);
getErrorString();

Also we have to disable the analytical jacobian generation -d=-NLSanalyticJacobian for some reason. I will create a separate ticket.

[Francesco Casella]

@Karim, please make sure with @adrpo or @sjoelund.se that Hudson points to the right repo URL.

[Adrian Pop]

It does: ​https://github.com/OpenModelica/OMLibraries/blob/master/repos.json#L960

[thomas.marx]

Rewriting the if-equations solved the backend problems for many models.

Could someone please add the flag setCommandLineOptions("-NLSanalyticJacobian"); (as stated above) for the library tests in Hudson?

[Francesco Casella]

Replying to thomas.marx:

Rewriting the if-equations solved the backend problems for many models.

Could someone please add the flag setCommandLineOptions("-NLSanalyticJacobian"); (as stated above) for the library tests in Hudson?

Done, see the last ​test report and ​coverage changes

Currently, 17/23 models compile and 6/23 simulate successfully.

[Karim Adbdelhak]

Replying to thomas.marx:

Most models that can't be simulated seem to fail with the same problem during initialisation:

messages = "Simulation execution failed for model: ThermalSeparation.Examples.Absorption_CO2_MEA
LOG_LS            | info    | initialize linear system solvers
|                 | |       | | 63 linear systems
LOG_LS            | info    | Start solving Linear System 965 (size 1) at time 0 with Lapack Solver
LOG_LS            | info    | Start solving Linear System 2150 (size 3) at time 0 with Lapack Solver
LOG_LS            | warning | Failed to solve linear system of equations (no. 2150) at time 0.000000, system is singular for U[2, 2].
LOG_LS            | info    | Matrix U
|                 | |       | |          0          0          0
|                 | |       | |          0          0          0
|                 | |       | |          0          0          0
LOG_LS            | info    | Output vector x
|                 | |       | | [ 1]    1.79765205516e-05
|                 | |       | | [ 2]    2.27687785313e-05
|                 | |       | | [ 3]     1.7983662067e-05
stdout            | warning | The default linear solver fails, the fallback solver with total pivoting is started at time 0.000000. That might raise performance issues, for more information use -lv LOG_LS.
LOG_LS            | info    | Start solving Linear System 2150 (size 3) at time 0 with Total Pivot Solver
LOG_LS            | warning | Matrix singular!
LOG_LS            | info    | rank =  0
LOG_LS            | warning | under-determined linear system not solvable!
stdout            | warning | Error solving linear system of equations (no. 2150) at time 0.000000.
stdout            | warning | Solving linear system 2150 fails at time 0. For more information use -lv LOG_LS.
LOG_LS            | warning | [1] Real column.balanceEquations.filmModel.diffCoeffVap[10].diffCoeff.D[5](start=0, nominal=1)
LOG_LS            | warning | [2] Real column.balanceEquations.filmModel.diffCoeffVap[10].diffCoeff.D[6](start=0, nominal=1)
LOG_LS            | warning | [3] Real column.balanceEquations.filmModel.diffCoeffVap[10].diffCoeff.D[4](start=0, nominal=1)

There are many matrices initialized with every element equal to 0, is this correct? In your library fill(0,...) is used often, so i guess this is no coincidence.

Furthermore there are many conflicting start values of alias variables, you can dump them by using "-d=aliasConflicts" while compiling. Here is one example of an alias set from ThermalSeparation.Examples.Absorption_CO2_MEA:

Warning: Alias set with conflicting start values
 * Candidate: column.p_v_in(start = 150000.0, confidence number = 7)
 * Candidate: column.mediumVapourIn.p(start = 100000.0, confidence number = 8)
 * Candidate: column.mediumVapourIn.state.p(start = 100000.0, confidence number = 9)
 * Candidate: column.mediumVapourIn.p_start(start = 160000.0, confidence number = 8)
 * Candidate: sourceGas_Vdot.medium.p(start = 100000.0, confidence number = 8)
 * Candidate: sourceGas_Vdot.medium.state.p(start = 100000.0, confidence number = 9)
 * Candidate: sourceGas_Vdot.medium.p_start(start = 160000.0, confidence number = 8)
=> Select value from column.p_v_in(start = 150000.0) for variable: column.p_v_in

These start values should all be equal, maybe the structural singularity of the initialization problem comes from these inconsistencies since it is mostly random which start values actually get chosen. I am not quite sure why Dymola seems to ignore this, but this should definitely be fixed in the models.

[Francesco Casella]

Regarding the alias conflicts, I think OMC is not behaving properly, see #5252.

[thomas.marx]

A few simulations failed due to a problem during the calculation of the diffusion coefficients.

The simulation of this test model below fails for nS= 4 as the resulting equation system is singular, but it works for nS=3.

model minimal_Diffcoeff

parameter Integer nS=4;
parameter Integer counter1[nS-1]={nS-i+1 for i in 2:nS};//für nS=4: {0,3,2,1};
parameter Integer aux[:] = {1,3,6,10,15, 21, 28, 36, 45};
parameter Integer a = aux[nS-1];

Real D[a];
Integer counter[nS];

equation

counter[1]=0;
counter[2:nS] = counter1;

for i in 1:nS-1 loop
  for k in i+1:nS loop
   D[k-i+ sum(counter[1:i])] = 1;
end for;
end for; 
end Fuller;

The problem can be solved by defining counter[nS] as a parameter and deleting the first two equations:

parameter Integer counter[nS]=cat(1,{0},counter1);

Of course, I think it is not desirable to use a variable for an array index. However, it seems odd that the model works for nS=3 and not for nS=4. Maybe this is a bug in OpenModelica or I am missing some point.

[Adrian Pop]

The NF sets all variables and parameters to DAE.POTENTIAL (even if they are not connectors), of course if they are not stream or flow.

Because of this BackendVariable.calcAliasKey will be different with the NF as it uses isVarConnector which checks for DAE.NON_CONNECTOR (which is never generated by the NF).

Opened #5336 with this issue.

[Karim Adbdelhak]

Sadly the fix of #5336 did not fix this one.

But i think i found the root of the problem, in NFEvalConstants.mo in function evaluateExpTraverser i found following code:

        // TODO: The runtime doesn't handle array literals well when used as
        //       arguments of external functions, since it sometimes tries to
        //       write to them (e.g. when trying to pack them). Until that's
        //       fixed we keep them as they are here.
        if ComponentRef.nodeVariability(cref) <= Variability.STRUCTURAL_PARAMETER and not 
           (isExternalArg and Type.isArray(ty))
            then
          // Evaluate all constants and structural parameters.
          outExp := Ceval.evalCref(cref, outExp, Ceval.EvalTarget.IGNORE_ERRORS(), evalSubscripts = false);
          outChanged := true;
        elseif outChanged then
          // If the cref's subscripts changed, recalculate its type.
          outExp := Expression.CREF(ComponentRef.getSubscriptedType(cref), cref);
        end if;

I replaced the condition with true:

        // TODO: The runtime doesn't handle array literals well when used as
        //       arguments of external functions, since it sometimes tries to
        //       write to them (e.g. when trying to pack them). Until that's
        //       fixed we keep them as they are here.
        if true then
          // Evaluate all constants and structural parameters.
          outExp := Ceval.evalCref(cref, outExp, Ceval.EvalTarget.IGNORE_ERRORS(), evalSubscripts = false);
          outChanged := true;
        elseif outChanged then
          // If the cref's subscripts changed, recalculate its type.
          outExp := Expression.CREF(ComponentRef.getSubscriptedType(cref), cref);
        end if;

This fixed it for the minimal example, this for sure isn't the right attempt, but when we fix this ToDo, i think we fix the ticket.

[Per Östlund]

Replying to Karim.Abdelhak:

But i think i found the root of the problem, in NFEvalConstants.mo in function evaluateExpTraverser i found following code:
...

With your change the new frontend will evaluate all variables, regardless of whether they are constants, parameter (structural or not) or time-varying. I.e. you removed not only the check for external arguments that the TODO is talking about, but also the check whether the cref refers to a constant or structural parameter.

Since there aren't any external functions in the minimal example I can only assume that it started working because all the parameters were evaluated, which can be achieved without any change to the compiler by just adding Evaluate = true annotations to the parameters. In fact, it seems like adding Evaluate = true to only the counter1 parameter is enough to get the model to compile.

[Karim Adbdelhak]

Replying to perost:

Since there aren't any external functions in the minimal example I can only assume that it started working because all the parameters were evaluated, which can be achieved without any change to the compiler by just adding Evaluate = true annotations to the parameters. In fact, it seems like adding Evaluate = true to only the counter1 parameter is enough to get the model to compile.

For the minimal case it works, thomas.marx could you try that with your libraries?

To clarify what went wrong:
The NF only evaluates structural parameters and constants, which is expected behavior but different from before. Since counter1 is used for indexing D in an equation all equations belonging to D form an algebraic loop since it is not clear which element gets solved in which equation. The jacobian for that loops resolves to be 0 in every element because all equations do not explicitly contain any of the unknown elements of D.

[thom.marx]

Adding the Evaluate=true annotation to the counter1 parameter works as well in the larger models.

[Francesco Casella]

@thom.marx, @perost's recent commit had a quite positive impact on the coverage of the ThermalSeparation library, see the ​report.

All models now pass successfully the frontend, backend and templates phases. Five more models successfully simulate, ten more pass the compilation phase successfully but then fail during simulation for some reason.

I would encourage you to check these new results with the latest 1.16.0 nightly build and report here if you have any indication regarding the reason of the simulation failures.