Enable Tearing by default

[Christian Schubert]

Before the 1.9.0 release, we should enable Tearing for all systems by default, as we gain a good speedup during simulation from it (see presentation from OpenModelica Workshop)

Last time I checked, it has been disabled for linear systems (and can be enabled with +d=doLinearTearing). This was because a linear torn system is solved with a non-linear solver which sometimes has problems doing so (try EngineV6, for example).
The idea was to use, either a linear solver or a newton-type non-linear solver.

What is the current situation?

[Adrian Pop]

Could we enable tearing if the linear system is above a certain size? 100 or so?

[Willi Braun]

Yes, we can. I'll do it.

[Willi Braun]

If we activate tearing of linear system for system >50 following tests are failing:

./simulation/libraries/3rdParty/PlanarMechanics/PlanarMechanicsForTesting.Examples.SimpleCarWithDifferentialGear.mos
./simulation/libraries/msl31/Modelica.Mechanics.MultiBody.Examples.Loops.Fourbar1.mos
./simulation/libraries/msl32/Modelica.Mechanics.MultiBody.Examples.Loops.Engine1b.mos
./simulation/libraries/msl32/Modelica.Mechanics.MultiBody.Examples.Loops.Fourbar1.mos
./simulation/libraries/msl32/Modelica.Mechanics.MultiBody.Examples.Systems.RobotR3.Components.MechanicalStructure.mos

and for system >100 following tests.

./simulation/libraries/msl31/Modelica.Mechanics.MultiBody.Examples.Loops.Fourbar1.mos
./simulation/libraries/msl32/Modelica.Mechanics.MultiBody.Examples.Loops.Engine1b.mos
./simulation/libraries/msl32/Modelica.Mechanics.MultiBody.Examples.Loops.Fourbar1.mos

Results are than different, as far as I see there are some missing start value for the iteration variables. So we would than need to change the models.

[Lennart Ochel]

Why not have other tools these problems? I guess in Dymola it is working with tearing and without any model changes, or not?

[Willi Braun]

Replying to lochel:

Why not have other tools these problems? I guess in Dymola it is working with tearing and without any model changes, or not?

I guess they select other tearing variable and perhaps the models are adapted to dymola tearing variables.

[Lennart Ochel]

Has anything happened in the last weeks? Because this is a blocker, it should be fixed soon or the priority should be changed (but probably this should be fixed ;-)).

[Martin Sjölund]

Postponed until 1.9.1

[Christian Schubert]

I ran the whole msl32 testsuite with +d=doLinearTearing this night.
We got plenty mismatches due to the warning

Warning: There are iteration variables with default zero start attribute. Use +d=initialization for more information.

since linear torn systems are handled like nonlinear systems internally.

The only real differences are:

• Modelica.Magnetic.FundamentalWave.Examples.BasicMachines.AIMS_Start which gave "Files Equal" on Windows isntead of "Files not Equal"
• Modelica.Electrical.Machines.Examples.Transformers.Rectifier12pulse.mos went to "Files not Equal" on Windows

I had a look at the results of the last test and we have an oscillating difference of 0.07 max while the absolute value is 250 max. Thus it should be ok.

My suggestion would be to change the tolerances a little on those two tests, which would also help with differences between Windows and Linux.

Should we enable Tearing for linear systems by default then?

[Lennart Ochel]

Yes, we should enable tearing for linear systems as soon as possible.

[Christian Schubert]

What I did not look at are the time required for simulation. On windows the tests have no time limit as far as I know.
Volker pointed out to me, that Modelica.Electrical.Analog.Examples.CauerLowPassSC takes longer with tearing enabled.
Are there any other?

[Adrian Pop]

Maybe this should be enabled only for systems above a certain size?

[Christian Schubert]

In theory tearing always improves simulation speed, so there is no need to disable it for small systems.
When I did the tearing investigations this January with Willi, we noted that the nonlinear solver sometimes had troubles solving the torn linear system.

Firstly, we should identify the models that take longer to simulate.
Secondly, we should look into it why this is. First question would be if the torn linear system is solved in exactly one iteration.

Could someone do the first step? I don't know how this can be done efficiently except enabling linearTearing and waiting for the MSL Coverage to run this night.

[ptaeuber wbraun lochel]

@ptaeuber: Could you take a look at the Modelica.Electrical.Analog.Examples.CauerLowPassSC model and find out why it takes so much longer to simulate with +d=doLinearTearing enabled?

It seems that there is still a bug in the way we handle torn systems.

[Patrick Täuber]

Yes, I will have a look at it.

[Patrick Täuber]

I think the fact that Modelica.Electrical.Analog.Examples.CauerLowPassSC takes longer with tearing enabled has something to do with the dassl solver. I figured out that the Jacobian for the dassl is at many points in time all over zero if tearing is enabled.

e.g.:
cj: 97.6563 and the states
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 analytical jacobian
at point in time : 0.17047
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

With tearing disabled:
cj: 97.6563 and the states
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 analytical jacobian
at point in time : 0.17047
-500000 0 0 0 0 0 0 0 0 0 0 0 -38272.4 0 11726.5 2315.82
0 -500000 0 0 0 0 0 0 0 0 0 0 38272.4 0 -11726.5 -2315.82
0 0 -500000 0 0 0 0 0 0 0 0 0 -11726.5 0 21906.4 4326.19
0 0 0 -500000 0 0 0 0 0 0 0 0 -11726.5 0 21906.4 4326.19
0 0 0 0 -500000 0 0 0 0 0 0 0 -2315.82 0 4326.19 38285.1
0 0 0 0 0 -500000 0 0 0 0 0 0 2315.82 0 -4326.19 -38285.1
0 0 0 0 0 0 -500000 0 0 0 0 0 2315.82 0 -4326.19 -38285.1
0 0 0 0 0 0 0 -500000 0 0 0 -58234.3 0 0 0 0
0 0 0 0 0 0 0 0 -500000 0 0 -58234.3 0 0 0 0
0 0 0 0 0 0 0 0 0 -500000 0 0 0 -38343.6 0 0
0 0 0 0 0 0 0 0 0 0 -500000 0 0 -38343.6 0 0
0 -500000 0 -500000 0 0 0 0 0 0 0 0 26545.8 0 10179.8 2010.37
-500000 0 0 0 0 0 0 0 5e-05 0 0 5.82343e-06 -38272.4 0 11726.5 2315.82
0 0 -500000 0 0 -500000 0 0 0 0 0 0 -9410.71 0 17580.2 -33958.9
0 0 0 0 0 0 0 5e-05 0 5e-05 0 5.82343e-06 0 3.83436e-06 0 0
0 0 0 0 500000 0 0 0 0 0 -5e-05 0 2315.82 -3.83436e-06 -4326.19 -38285.1

[Christian Schubert]

Hi Patrick,

thank you very much for your problem description!
I was wondering whether you could try to find out why this happens and maybe even fix it?

Then, finally, we could enable tearing and end this long development project.

Best wishes,

Christian

[Lennart Ochel]

Replying to ptaeuber:

I think the fact that Modelica.Electrical.Analog.Examples.CauerLowPassSC takes longer with tearing enabled has something to do with the dassl solver. I figured out that the Jacobian for the dassl is at many points in time all over zero if tearing is enabled.

@Willi: any comments?

[Patrick Täuber]

I doubt that I can solve the problem because I don't know the dassl code. But I found out some hints so maybe someone else (Willi?) could easier find out the reason for this problem...

The thing is that the simulation problem only occurs after time=1. At earlier points of time the simulation process is the same as the one with tearing disabled.

The second clue is something like an error message from the LOG_DDASRT-log. From time=0 to time=1 it is:
total number of error test failures: 0

From then on errors occur:
| | | | total number of error test failures: 725
| | | | total number of error test failures: 890
| | | | total number of error test failures: 7020
| | | | total number of error test failures: 12249
| | | | total number of error test failures: 3714
and so on...

Without tearing the number is all over 0 and so it is in all the other models I tested, too. So I guess this output has something to do with the problem.
I don't know the meaning of this errors but maybe this will help somebody to find out the problem.

[Willi Braun]

I had a look on the issue and will continue,
but I haven't found the reason now.

Patrick, thanks for the hint with the jacobians.

[Christian Schubert]

Are there any news on this topic?

[Christian Schubert]

I had a look at this issue and would like to document what I've found out so far.
The jacobian matrix which is zero as described by ptaeuber is being calculated using finite differences, i.e. column J_i of Jacobian J is calculated as

J_i = (f(x+eps*e_i)-f(x))/eps

with e_i being a vector full of zeros and a one at colunm i.

Now, when tearing is enabled

f(x+eps*e_i) == f(x)

and thus J_i is exactly zero. eps has the magnitude of about 1e-13.

My interpretation would be the following: Linear torn systems are currently handled by a non-linear Newton solver. I assume that the current/old/extrapolated solution is good enough for the non-linear solver and thus no iteration is done when the state vector x is being disturbed by eps*e_i. Consequently, the solution of the non-linear system would be the same and so would be f(x+eps*e_i).

@Willi: Can you confirm my theory?

The big question is: What can we do about it?

[Willi Braun]

Enabled linear tearing by default in r21558.