reformulation of basicFlowFunctions

[Mathadon]

Similar to #723 IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp can be reformulated from

   m_flow := if noEvent(dp>dp_turbulent) then k*sqrt(abs(dp))
           elseif noEvent(dp<-dp_turbulent) then -k*sqrt(abs(-dp))
         else (k^2*5/4/m_flow_turbulent)*dp-k/4/(m_flow_turbulent/k)^5*dp^3;

into

     m_flow := if noEvent(abs(dp)>dp_turbulent)
               then sign(dp)*k*sqrt(abs(dp))
               else dp*m_flow_turbulent/dp_turbulent/4*(5-(dp/dp_turbulent)^2);

In a (best case) example model this yields a significant speedup. Consider the example:

model test
  parameter Integer i = 10000000;
  Real y;
algorithm 
  y:=0;
  for k in 1:i loop
    y :=y*IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp(
      0,
      0.5,
      2+time/1e10);
  end for;
end test;

which evaluates the function only in the regularisation region. I simulated this using

simulateModel("test", numberOfIntervals=1, method="Euler", fixedstepsize=1, resultFile="test");

The old formulation leads to:

Integration terminated successfully at T = 1
   CPU-time for integration      : 2.23 seconds
   CPU-time for one GRID interval: 2.23e+03 milli-seconds
...

Profiling information for the blocks.
Estimated overhead per call        0.00[us] total        0.000[s]
the estimated overhead has been subtracted below.
Name of block   , Block, Total CPU[s], Mean[us]    ( Min[us]    to Max[us]    ),   Called
Entire model    :     0,        0.000,        1.00 (       1.00 to        1.00),        3
Event handling  :     1,        2.520,  1260041.50 ( 1215051.00 to  1305032.00),        2
Empty timer     :     2,        0.000,        0.12 (       0.00 to        1.00),        8
Outside of model:     3,        0.000,       75.20 (       1.00 to      248.00),        5
AcceptedSection2:     4,        5.921,  1184230.00 ( 1086584.00 to  1305031.00),        5
 Auxiliary2 code:     5,        5.921,  1184229.00 ( 1086583.00 to  1305030.00),

The new formulations leads to:

(One-step solver with fixed stepsize and order 1)
Integration terminated successfully at T = 1
   CPU-time for integration      : 0.318 seconds
   CPU-time for one GRID interval: 318 milli-seconds
...
Profiling information for the blocks.
Estimated overhead per call        0.00[us] total        0.000[s]
the estimated overhead has been subtracted below.
Name of block   , Block, Total CPU[s], Mean[us]    ( Min[us]    to Max[us]    ),   Called
Entire model    :     0,        0.000,        1.00 (       1.00 to        1.00),        3
Event handling  :     1,        0.288,   143859.00 (  142819.00 to   144899.00),        2
Empty timer     :     2,        0.000,        0.12 (       0.00 to        1.00),        8
Outside of model:     3,        0.000,       34.00 (       1.00 to      100.00),        3
AcceptedSection2:     4,        0.743,   148631.20 (  137263.00 to   180118.00),        5
 Auxiliary2 code:     5,        0.743,   148630.00 (  137262.00 to   180117.00),        5

the cpu times seem to be self-inconsistent but they do show a clear advantage of the new formulation. Re-adding the else if construct does not affect computation time so the gains in computation time seem to be entirely related to the formulation of the else branch.

When changing the else branch from
dp*m_flow_turbulent/dp_turbulent/4*(5-(dp/dp_turbulent)^2) to
dp*m_flow_turbulent^5/m_flow_turbulent^4/dp_turbulent/4*(5-(dp/dp_turbulent)^2)
computation time rises by more than a factor of 10.

(m_flow_turbulent/k)^2 is a lot faster than (m_flow_turbulent/k)*(m_flow_turbulent/k).

Adding and eliminating
Modelica.SIunits.PressureDifference dp_div = dp/dp_turbulent "Help variable for computational efficiency";

Leads to CPU-time for integration : 0.255 seconds

Replacing /4 by *0.25 does not have a noticeable effect.

I'm not sure what's going on, but I suggest to go with this new formulation:

      m_flow := if noEvent(abs(dp)>dp_turbulent)
                then sign(dp)*k*sqrt(abs(dp))
                else 0.25*dp_div*m_flow_turbulent*(5-dp_div^2);

[Mathadon]

The same test for IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_m_flow

Old:

Old:

 dp :=if (m_flow>m_flow_turbulent) then (m_flow/k)^2
      elseif (m_flow<-m_flow_turbulent) then -(m_flow/k)^2
      else (m_flow_turbulent*m_flow+m_flow^3/m_flow_turbulent)/2/k^2;

CPU-time for integration : 0.622 seconds

New:

 dp :=if (abs(m_flow)>m_flow_turbulent)
      then sign(m_flow)*(m_flow/k)^2
      else 0.5*m_flow_turbulent*m_flow/k^2*(1+(m_flow/m_flow_turbulent)^2);

CPU-time for integration : 0.272 seconds

[Mathadon]

Note, inlining does not seem to work in dymola 2017.

[Mathadon]

IBPSA.Fluid.Examples.FlowSystem.Basic is about 15% faster

Old:

DynamicsSection :    11,        4.487,      639.97 (     304.00 to     3811.00),     7011
...
 Nonlin sys(29) :    15,        3.538,      504.69 (     211.00 to     3538.00),     7011

New:

DynamicsSection :    11,        3.854,      549.69 (     301.00 to     2541.00),     7011
...
 Nonlin sys(29) :    15,        3.005,      428.54 (     207.00 to     2239.00),     7011

[Mathadon]

I found another formulation that's a bit slower, but it is C2 continuous:

protected 
  Modelica.SIunits.PressureDifference dp_turbulent = m_flow_turbulent^2/k/k
    "Pressure where flow changes to turbulent";
  Real dpNorm=dp/dp_turbulent;
  Real dpNormSq=dpNorm*dpNorm;
algorithm 
   m_flow :=  if noEvent(abs(dp)>dp_turbulent) then sign(dp)*k*sqrt(abs(dp))
 else 
  (1.40625  + (0.15625*dpNormSq - 0.5625)*dpNormSq)*m_flow_turbulent*dpNorm;

this one takes 0.36 seconds in the test (instead of 0.255), but it's C2 continuous and for such an important function I think it makes sense to use it instead of the current C1 continuous implementation.

If you agree then I'll re-implement the derivative functions and basicFlowFunction_m_flow too and make a new pull request.

[mwetter]

I think that's fine. Maybe
dp_turbulent = (m_flow_turbulent/k)^2 would be a bit faster.

[Mathadon]

For future reference:
the total CPU time should not be used to compare the performance of individual functions, even in such simple tests. The simplification that @mwetter suggested leads to a CPU time increase in the total CPU time and to a CPU time decrease in the Advanced.GenerateBlockTimers=true output.

Based on this timer the C2 implementation is actually faster than 0.25*dpNorm*m_flow_turbulent*(5-dpNorm^2) (0.57s vs 0.59s).

edit: for IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_m_flow 0.61s vs 0.57s(old).