Models using the IF97 water medium are slow due to lack of common subexpression elimination

[casella]

The Modelica.Media IF97 water/steam model was designed in order to provide efficient performance (i.e., no duplicate calls to core computations) assuming that the LateInline annotation and common subexpression elimination were supported by the compiler.

The former requirement has been satisfied in OMC for a while, the latter is currently not yet fulfilled. As a result, the function Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph (which is very computationally expensive) is called many times with exactly the same argument, leading to an outrageous slow-down of the simulation process, possibly leading to time-outs (see, e.g., ThermoPower.Examples.CISE.CISESim180503 or ThermoPower.Test.DistributedParameterComponents.TestWaterFlow1DFEM_A).

The analysis is confirmed by turning on +d=dumpSimCode and by running the profiler on test cases.

This is also the root cause of the very slow simulation of simulation/libraries/3rdParty/ThermoPower/Bug2537.mos, which was reported by Martin Sjoelund on Jan 12 as being the critical path preventing the testsuite speed-up by parallelization.

Since the IF97 model is essential for the modelling of most power production systems, and usually gobbles up most of the CPU time, I would recommend to implement at least a very rudimentary CSE scheme that is able to handle the case of calls to the same function with the same arguments, e.g., Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(hex.p, hex.h[1], 0.0, 0).

[lochel]

I contributed a prototype of a common subexpression elimination module last week. The results are very promising and it should be activated very soon.

[casella]

The CSE module can be activated with the omc flag +cseCall. It dramatically improves the performance of the slowest ThermoPower test cases, but unfortunately it still breaks a lot of MSL test cases, so it cannot yet be activated by default.

Work will resume after the Annual Meeting.

[casella]

Update: with +cseCall the situation has improved a lot. However, there are still some issues with duplicate calls not being detected. For example, in the code generated for ThermoPower.Test.CISE.CISESim120501, most waterBaseProp_ph() function calls are duplicated, e.g.:

1006: $CSE75 := Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(
                Plant.HeaderLower.p, Plant.Downcomer.h[1], 0.0, 0);
1009: $CSE69 := Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(
                Plant.HeaderLower.p, Plant.Downcomer.h[1], 0, 0);

Of course there is a difference in the two expressions: in one case the phase input is 0.0 (Real), while in the other case it is 0 (Integer).

Unfortunately the transformational debugger is not working on these CSE variables (please make sure it does!), so I cannot track from where the Integer input phase gets a Real value from, I guess there's something wrong in the CSE handling code.

In any case, if there is some casting involved when the function is called, CSE should obviously be applied downstream the casting. Incidentally, does OMC silently cast Reals to Integers? That sounds a bit weird for a strongly typed language as Modelica.

[sjoelund.se]

Replying to casella:

In any case, if there is some casting involved when the function is called, CSE should obviously be applied downstream the casting. Incidentally, does OMC silently cast Reals to Integers? That sounds a bit weird for a strongly typed language as Modelica.

This is often something in ExpressionSimplify going wrong, adding a bad zero variable somewhere.

[vitalij]

it seems the function calls came from

    fluidState = Medium.setState_ph(p, h); // $CSE69

    dM_dt = V*(Medium.density_derp_h(fluidState)*der(p) + Medium.density_derh_p(
      fluidState)*der(h)); // $CSE75

perhaps we handling function calls like

   fluidState = Medium.ThermodynamicState(d = Medium.density_ph(p, h), 
                                          T = Medium.temperature_ph(p, h),
                                          phase = phase, 
                                          h = h, 
                                          p = p); 
    //density_ph has waterBaseProp(p, h, phase, region)

   y = Medium.density_derp_h(fluidState); // waterBaseProp_ph(p, h, phase, region) inside

with extends object different?

[sjoelund.se]

Well. It's a bit more odd than I thought it would be...

input Integer phase = 0 "Phase: 2 for two-phase, 1 for one phase, 0 if unknown";

So... The integer 0 is correct and the 0.0 something weird.

[sjoelund.se]

The reason $CSE75 is not easily visible in the debugging file is:

1. It is a record, and the main record variable is not a real variable
2. The generated CSE call is an algorithm section instead of an equation section. This strips all transformations performed on the equation (even if it is added to the algorithm). I think it would be better to make CSE create equations. Or possibly to treat algorithm sections with a single assignment differently for printing debugging information.

[sjoelund.se]

OK. We already treat the first statement in an algorithm special. The CSE module simply strips out all equation information from the source. It should keep the old source information and add a CSE operation.

[sjoelund.se]

Inline:

-Plant.Downcomer.drdp[1] = ThermoPower.Water.Flow1DFV2ph$Plant$Downcomer.Medium.density_derp_h(Plant.Downcomer.fluidState[1])
+Plant.Downcomer.drdp[1] = Modelica.Media.Water.IF97_Utilities.ddph_props(Plant.Downcomer.fluidState[1].p, Plant.Downcomer.fluidState[1].h, Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(Plant.Downcomer.fluidState[1].p, Plant.Downcomer.fluidState[1].h, Plant.Downcomer.fluidState[1].phase, 0))

Substitute:

-Modelica.Media.Water.IF97_Utilities.ddph_props(Plant.Downcomer.fluidState[1].p, Plant.Downcomer.fluidState[1].h, Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(Plant.Downcomer.fluidState[1].p, Plant.Downcomer.fluidState[1].h, Plant.Downcomer.fluidState[1].phase, 0))
+Modelica.Media.Water.IF97_Utilities.ddph_props(Plant.HeaderLower.p, Plant.Downcomer.h[1], Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(Plant.HeaderLower.p, Plant.Downcomer.h[1], 0.0, 0))

Plant.Downcomer.fluidState[1].phase is substituted by 0.0. Which seems wrong since that is declared as Integer.
Plant.Downcomer.fluidState[1].phase also has fixed=false useStart=false...

[vwaurich]

The wrong substitution should be fixed with [24615].

[casella]

Replying to sjoelund.se:

The reason $CSE75 is not easily visible in the debugging file is:

1. It is a record, and the main record variable is not a real variable
2. The generated CSE call is an algorithm section instead of an equation section. This strips all transformations performed on the equation (even if it is added to the algorithm). I think it would be better to make CSE create equations. Or possibly to treat algorithm sections with a single assignment differently for printing debugging information.

Please make sure this is handled correctly. It is too bad it one cannot analyze a given model fully just because common subexpressions have been found. BTW, why putting CSEs in the algorithm section?

[casella]

Replying to vwaurich:

The wrong substitution should be fixed with [24615].

I have rerun the ThermoPower_Experimental task on hudson (see ​https://test.openmodelica.org/libraries/ThermoPower_Experimental/BuildModelRecursive.html) using r24619 and +cseCall but I don't see any substantial improvement in simulation time. CISESim120501 is still at 17 seconds. Unfortunately I cannot see the code dump of those runs, nor I can download a nightly build and try it out myself. Have you checked that you no longer get all those duplicate CSE computations?

[vwaurich]

Replying to casella:

Replying to vwaurich:

The wrong substitution should be fixed with [24615].

I have rerun the ThermoPower_Experimental task on hudson (see ​https://test.openmodelica.org/libraries/ThermoPower_Experimental/BuildModelRecursive.html) using r24619 and +cseCall but I don't see any substantial improvement in simulation time. CISESim120501 is still at 17 seconds. Unfortunately I cannot see the code dump of those runs, nor I can download a nightly build and try it out myself. Have you checked that you no longer get all those duplicate CSE computations?

I just fixed the issue with the real expression that should be an integer.The bug was in RemoveSimpleEquations. The two function calls from comment 3 are now equal but there are still 2 different cses generated. Maybe some one from the cse-module developers could have a look again.

[sjoelund.se]

The following is not optimized the correct way either (+cseCall +cseEachCall):

model test
  Real x,y;
  final parameter Real c = 0;
equation
  x = sin(time) + 2;
  y = sin(time) + c;
end test;

It becomes:

y = sin(time);
CSE3 = sin(time);
x = CSE3 + 2;

[lochel]

Replying to sjoelund.se:

The following is not optimized the correct way either (+cseCall +cseEachCall):

model test
  Real x,y;
  final parameter Real c = 0;
equation
  x = sin(time) + 2;
  y = sin(time) + c;
end test;

It becomes:

y = sin(time);
CSE3 = sin(time);
x = CSE3 + 2;

I think this is working correctly now. Maybe it was broken a few days back.
If +cseCall is used then nothing should happen, since both equations are separated in different partitions.
If +cseEachCall is use then both calls are substituted with an own cse variable, due to the same reason:

y = $cse1
$cse1 = sin(time)
x = 2.0 + $cse2
$cse2 = sin(time)

[sjoelund.se]

Replying to lochel:

If +cseEachCall is use then both calls are substituted with an own cse variable, due to the same reason:

y = $cse1
$cse1 = sin(time)
x = 2.0 + $cse2
$cse2 = sin(time)

Then you should move these equations to the same partition, I guess. Or a special partition for cse without dependent vars? Or ... something. Because for slow functions, we do have this very issue that we execute them multiple times with the identical input.

[lochel]

Replying to sjoelund.se:

Then you should move these equations to the same partition, I guess. Or a special partition for cse without dependent vars? Or ... something. Because for slow functions, we do have this very issue that we execute them multiple times with the identical input.

Yes, that is true for independent function calls. This is one thing I will add later on since the first goal is just to get it working in general. I will set up a ticket to keep track of it.

BTW: Can one summarize what the actual issue of this ticket is?

[casella]

The original motivation of this ticket was that test cases using the Modelica.Media IF97 water model were unnecessarily slow because of repeated calls to the computationally-intensive waterBaseProp_ph function with the same arguments. If CSE is implemented correctly, this should never be the case.

Issue 1: in the test model by Martin (comment:15), the expression sin(time) shows up twice. I would expect it to be only computed once in an auxiliary variable, which is then used twice. What I see, instead, is two CSE auxiliaries each computing sin(time), which is not useful at all.

Issue 2: as of r24759, code generated for the ThermoPower.Examples.CISE.CISESim120501 model still contains duplicated calls, e.g.:

$CSE75 := Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(
            Plant.HeaderLower.p, Plant.Downcomer.h[1], 0, 0);
$CSE69 := Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(
            Plant.HeaderLower.p, Plant.Downcomer.h[1], 0, 0);

regardless whether +cseCall or +cseEachCall is activated. Only one CSE auxiliary variable should be created for each such function call.

[lochel]

I just added a flag to disable partitioning for now (+d=disablePartitioning). That is a workaround for issue 1:

########### Updated Equation List ########### (3, 3)
========================================
1/1 (1): x = 2.0 + $cse1
2/2 (1): y = $cse1
3/3 (1): $cse1 = sin(time)

I didn't tested it with the model from issue 2, but maybe the duplicated calls will also disappear in this case.

[lochel]

Replying to casella:

Issue 2: as of r24759, code generated for the ThermoPower.Examples.CISE.CISESim120501 model still contains duplicated calls, e.g.:

$CSE75 := Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(
            Plant.HeaderLower.p, Plant.Downcomer.h[1], 0, 0);
$CSE69 := Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(
            Plant.HeaderLower.p, Plant.Downcomer.h[1], 0, 0);

regardless whether +cseCall or +cseEachCall is activated. Only one CSE auxiliary variable should be created for each such function call.

I just run ThermoPower.Examples.CISE.CISESim120501 with +cseEachCall and I got the mentioned function call only once:

861/1351 (16): $cse38 = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(Plant.HeaderLower.p, Plant.Downcomer.h[1], 0, 0)

[sjoelund.se]

Milestone changed to 1.9.3 since 1.9.2 was released.

[sjoelund.se]

Moved to new milestone 1.9.4

[sjoelund.se]

Milestone pushed to 1.9.5

[sjoelund.se]

Milestone renamed

[casella]

Consider the model ThermoPower.Examples.CISE.CISESim120501. The simulation time with Dymola on my PC is 2.45 seconds, and it takes 248 DASSL steps and 89 Jacobian computations.

I have switched fromt +cseCall +cseEachCall to --postOptModules+=wrapFunctionCalls, as suggested by omc's warning message.

The simulation time on my pc is then 18 seconds. DASSL takes 311 steps and 87 Jacobian computations, roughly the same as Dymola. Colored jacobians are used in both tools, so the difference in performance is all in the rhs computation of the ODEs. Of these 18 seconds, 15.6 are spent computing waterBaseProp_ph 477888 times, and one more second is spent computing dewcurve_p and boilingcurve_p 72383 times each.

If I look at the odeEquations section of the dumpSimCode output, there are many cases where CSE is working as expected, e.g.:

2087: $cse29 := Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(Plant.Risers.p, Plant.HeaderUpper.h, 0, 0);
2088: $cse30=Modelica.Media.Water.IF97_Utilities.ddph_props(Plant.Risers.p, Plant.HeaderUpper.h, $cse29)[Real ]
2089: $cse31=Modelica.Media.Water.IF97_Utilities.ddhp_props(Plant.Risers.p, Plant.HeaderUpper.h, $cse29)[Real ]

However, there are still a lot of other cases where it doesn't and duplicate function calls show up, e.g.:

1764: Plant.Risers.T[1]   =Modelica.Media.Water.IF97_Utilities.T_props_ph  (Plant.Risers.p, Plant.HeaderLower.h, Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(Plant.Risers.p, Plant.HeaderLower.h, 0, 0))[Real ]
1765: Plant.Risers.rho[1] =Modelica.Media.Water.IF97_Utilities.rho_props_ph(Plant.Risers.p, Plant.HeaderLower.h, Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(Plant.Risers.p, Plant.HeaderLower.h, 0, 0))[Real ]
1766: Plant.Risers.drdp[1]=Modelica.Media.Water.IF97_Utilities.ddph_props  (Plant.Risers.p, Plant.HeaderLower.h, Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(Plant.Risers.p, Plant.HeaderLower.h, 0, 0))[Real ]
1767: Plant.Risers.drdh[1]=Modelica.Media.Water.IF97_Utilities.ddhp_props  (Plant.Risers.p, Plant.HeaderLower.h, Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(Plant.Risers.p, Plant.HeaderLower.h, 0, 0))[Real ]

This is quite weird, as the original Modelica code reads:

    fluidState[j] = Medium.setState_ph(p, h[j]);
    T[j] = Medium.temperature(fluidState[j]);
    rho[j] = Medium.density(fluidState[j]);
    drdp[j] = Medium.density_derp_h(fluidState[j]);
    drdh[j] = Medium.density_derh_p(fluidState[j]);

there is probably something wrong with function inlining. Another such case is:

2110: Plant.SH.rhol=Modelica.Media.Water.IF97_Utilities.BaseIF97.Regions.rhovl_p($cse14.psat, Modelica.Media.Water.IF97_Utilities.BaseIF97.Regions.boilingcurve_p($cse14.psat))[Real ]
2112: Plant.SH.hl=  Modelica.Media.Water.IF97_Utilities.BaseIF97.Regions.hvl_p  ($cse14.psat, Modelica.Media.Water.IF97_Utilities.BaseIF97.Regions.boilingcurve_p($cse14.psat))[Real ]

stemming from

  rhol = Medium.bubbleDensity(sat);
  hl = Medium.bubbleEnthalpy(sat);

I would say that, on average, each function call is repeated 3 times, which means we could reduce the simulation time of this model (and of most models using IF97) by a factor two or more by avoiding these unnecessarily repeated function calls.

The efficient support of the IF97 model is an essential feature for most people using Modelica.Media and for all power plant models. Therefore I would set this ticket as a blocker for 2.0.0. Please have a look at this issue and find out how hard it is to avoid these multiple function calls.

[rfranke]

Even worse: the command line option +cseCall does not seem to do anything anymore. It provided a speedup of a factor of 2 for the DrumBoiler optimization example used in the publication at Modelica 2015: ​http://www.ep.liu.se/ecp/118/036/ecp15118339.pdf

This is gone with the current nightly builds. I discovered this recently when working on the Newton solver and comparing the new results with the old ones.

[casella]

Replying to rfranke:

Even worse: the command line option +cseCall does not seem to do anything anymore.

Try +cseEachCall instead

[rfranke]

+cseEachCall produced similar results as +cseCall when writing the paper last year. It does not seem to work anymore either now.

[adrpo]

I find this really strange as we do have tests for it which they pass fine. Maybe we broke this particular model somehow.

[casella]

Replying to adrpo:

I find this really strange as we do have tests for it which they pass fine. Maybe we broke this particular model somehow.

The model is not broken, it simply takes much longer to simulate. I don't think our automatic testing procedures catch this.

[adrpo]

By broken I mean the modifications that should be applied by cse are not.

[rfranke]

It worked with OpenModelica 1.9.6. Just simulate Modelica.Fluid.Examples.DrumBoiler.DrumBoiler in OMEdit and browse for _cse variables.

[casella]

Adrian, can you point me out the tests? I understand there are some cases where CSE works and some where it doesnt', I'd like to see them.

[lochel]

Both flags, cseCall and cseEachCall, are deprecated and are processed as --postOptModules+=wrapFunctionCalls. Should we remove the old flags? They were just kept for backwards compatibility and to not break existing scripts.

Tests for module wrapFunctionCalls are mainly located in testsuite/simulation/modelica/commonSubExp.

Jan is currently working on this module and further improvements and bug fixes will be available soon.

[rfranke]

Good to know. It should become a default setting asap (hopefully without introducing these _cse variables in the plot browser).

I just tried the DrumBoiler example with

setCommandLineOptions("--postOptModules+=wrapFunctionCalls");

but didn't see an impact of this option either.

[lochel]

Replying to rfranke:

Good to know. It should become a default setting asap (hopefully without introducing these _cse variables in the plot browser).

Currently, there are issues with 4 MSL models (besides some known bugs). Once that is fixed, we plan to make this a default module.

[casella]

.mos file with IF97 CSE test cases

[casella]

I have added an attachment with three typical use patterns of the IF97 medium model. At the moment, no CSE at all is reported by OMC, while it should, in all three cases.

[jhagemann]

in all three test cases (for example in Test1: T = Medium.temperature(state)) is T a cse variable and is not necessary to create a "new" cse.

see also ticket: 3910

[rfranke]

OpenModelica with the new --postOptModules+=wrapFunctionCalls still does not seem to identify any common subexpressions in the example Modelica.Fluid.DrumBoiler.DrumBoiler. Please compare with the released version 1.9.6. Your new developments should not make things worse.

[jhagemann]

Patrick and I fixed one bug (​a1f6080). We hope that the module creates for every function call a cse variable again. The idea that for the case “CREF=CALL” no new cse variable is introduced but CREF is used instead, still exists and will be implemented in a correct way. Please test your examples and report if your problems are solved or not.

The handling of cse variables in different partitions is still not implemented in a good way. Please use +d=disablePartitioning for now. (see also: ticket 3812)

[casella]

After Jan's commit, the attached test case testIF97CSE.mos​ now produces the expected results:

Test1:
1/1 (1): cp = $cse2   [dynamic]
2/2 (1): rho = $cse3   [dynamic]
3/3 (1): T = $cse4   [dynamic]
4/4 (1): der(h) = 200000.0 - h   [dynamic]
5/5 (1): der(p) = 2000000.0 - p   [dynamic]
6/6 (1): $cse4 = Modelica.Media.Water.IF97_Utilities.T_props_ph(p, h, $cse1)   [binding]
7/7 (1): $cse3 = Modelica.Media.Water.IF97_Utilities.rho_props_ph(p, h, $cse1)   [binding]
8/8 (1): $cse2 = Modelica.Media.Water.IF97_Utilities.cp_props_ph(p, h, $cse1)   [binding]
9/9 (16): $cse1 = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(p, h, 0, 0)   [binding]

Test2:
1/1 (1): der(p) = 2000000.0 - p   [dynamic]
2/2 (1): der(h) = 200000.0 - h   [dynamic]
3/3 (1): T = $cse2   [dynamic]
4/4 (1): rho = $cse3   [dynamic]
5/5 (1): $cse3 = Modelica.Media.Water.IF97_Utilities.rho_props_ph(p, h, $cse1)   [binding]
6/6 (1): $cse2 = Modelica.Media.Water.IF97_Utilities.T_props_ph(p, h, $cse1)   [binding]
7/7 (16): $cse1 = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(p, h, 0, 0)   [binding]

Test 3:
1/1 (1): der(p) = 2000000.0 - p   [dynamic]
2/2 (1): der(h) = 200000.0 - h   [dynamic]
3/3 (1): T = $cse2   [dynamic]
4/4 (1): rho = $cse3   [dynamic]
5/5 (1): cp = $cse4   [dynamic]
6/6 (1): $cse4 = CSE.Test3.Medium.specificHeatCapacityCp($cse1)   [binding]
7/7 (1): $cse3 = CSE.Test3.Medium.density($cse1)   [binding]
8/8 (1): $cse2 = CSE.Test3.Medium.temperature($cse1)   [binding]
9/9 (5): $cse1 = CSE.Test3.Medium.setState_ph(p, h, 0)   [binding]

In each case, the computationally-intensive functions (setState_ph and waterBaseProp_ph) are only called once as expected.

In response to @rfranke's concerns, I would suggest to add the attached testIF97CSE.mos​ to the testsuite, to make sure some future developments do not break this again.

[casella]

I have re-run the ThermoPower_Experimental coverage test, and it is obvious that Jan's commit has had a very positive impact on the performance of all the models using IF97, which are the majority, just compare the simulation times in
​https://test.openmodelica.org/libraries/history/ThermoPower_Experimental-2016-05-10.html
​https://test.openmodelica.org/libraries/history/ThermoPower_Experimental-2016-05-19.html

Some tests, e.g. CISESim180503, had an almost 3X improvement, most had about 2X improvement. There are still some where the improvement was marginal, for example TestWaterFlow1DFEM_G went down from 30.84 to only 27.87 (-10%). I will investigate why with some profiling, but I suspect this is because the bottleneck in this case is the nonlinear solver, see #3111.

I guess also Ruediger's boiler models should have benefitted from this change.

[jhagemann]

Thank you for your detailed report.

[casella]

I have now analyzed the dumpSimCode logs of some ThermoPower models. For example, check the dumpSimCode log of ThermoPower.Test.DistributedParameterComponents.TestWaterFlow1DFEM_F looking for calls to waterBaseProp_ph(hexB.p, hexB.h[1], 0, 0).

Before Jan's commit, the ode equations contained four repeated calls to this function. After the commit, the function is only called once to set a CSE variable, so the performance in simulation is now as expected.

However, in the initial equation section there are still four calls to that function, and four more in the initial_lambda0 section. In the case of steady-state initialization with homotopy, also the initial equations could end up getting evaluated a large number of times. I have opened #3921 on this topic.

[casella]

One last comment. As reported in comment:44, the TestWaterFlow1DFEM_G only went down from 30 to 27 seconds after CSE was improved. However, we have tried running it on three different machines and two OSs here at Politecnico, turning on profiling, and the simulation time was around 5 seconds, of which only 1 second was spent by the steam property function calls, that need CSE, which is a good result if we look at CSE alone.

Unfortunately, we don't have access to detailed profiling of the coverage test results, so we can't say why the OSMC server takes so much more time to run the simulation, and to do what. My guess is that this is related to convergence problems with the nonlinear equation solver.

[sjoelund.se]

1.11 will be released before 2.0