optimize Fluid.BaseClasses.FlowModels.basicFlowFunction_dp

[mwetter]

The function Fluid.BaseClasses.FlowModels.basicFlowFunction_dp and basicFlowFunction_m_flow lead to many operations outside the turbulent region.
This issue should look at whether a simpler implementation leads to faster computation. One could for example use a linear function for |m_flow| < m_flow_turbulent/2 and use regularization between these two regions. This regularization could be done with a 3rd order polynomial. For fixed resistances, all polynomial coefficients are parameters, so there should be little overhead, whereas for valves, these coefficients need to be computed during the time stepping. Hence, maybe extending the linear function to
|m_flow| < 2/3*m_flow_turbulent may be more efficient.

[Mathadon]

I'm in favour of this change since the combined functions are indeed very complex and often they require many evaluations because they are used in non-linear systems. Moreover a linear regime would better correspond to the laminar flow regime that occurs in practice. I suggest to look for a solution in this direction. For valves things may be a bit more complicated, although a pragmatic approximation may be sufficiently accurate?

[Mathadon]

I created an (arguably dumb) implementation for Annex60.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp:

m_flow :=smooth(1, if dp>0
                    then Annex60.Utilities.Math.Functions.spliceFunction(x=dp-dp_turbulent,pos=if dp>dp_turbulent/4 then k*sqrt(dp) else 0, neg=dp/dp_turbulent*k*sqrt(dp_turbulent)/2, deltax=dp_turbulent/3)
  else 
     Annex60.Utilities.Math.Functions.spliceFunction(x=-dp_turbulent-dp,pos=if dp<-dp_turbulent/4 then -k*sqrt(abs(dp)) else 0, neg=dp/dp_turbulent*k*sqrt(dp_turbulent)/2, deltax=dp_turbulent/3));

Unit test result difference:

I then evaluated Annex60.Fluid.Examples.Performance.Example2(from_dp.k=true), which does pretty much nothing else than evaluate that function. Results are

before

Nonlin sys(5)  :    10,       20.357 s,       33.63 (      24.00 to      567.00),   605328

after

Nonlin sys(5)  :    10,        8.987 s,       14.85 (       5.00 to      354.00),   605328

so this does indicate that a lot of performance can be gained. We should however look for a better formulation. Can't we construct some non-piecewise function that approaches sqrt(x) but that does not have the division by zero around zero flow?

edit:

This is somewhat more elegant:

      m_flow := smooth(1, if dp>dp_turbulent then k*sqrt(dp) else 
                          if dp<-dp_turbulent then -k*sqrt(-dp) else 
                            k*sqrt(dp_turbulent)*sin(dp/dp_turbulent*Modelica.Constants.pi/2));

and results in

 Nonlin sys(5)  :    10,       6.5 s,       26.74 (       2.00 to     6583.00),   605328

we could tweak the switching point between the sinusoid and the square root to have a continuous first derivative. This implementation should work for valves too.

[Mathadon]

Note that these functions are not inlined even though LateInline=true. This may imply a performance penalty.

[Mathadon]

The functions are inlined when dp_turbulent and kSqu are eliminated:

  input Modelica.SIunits.Pressure dp(displayUnit="Pa") 
    "Pressure difference between port_a and port_b (= port_a.p - port_b.p)";
  input Real k(min=0, unit="") 
    "Flow coefficient, k=m_flow/sqrt(dp), with unit=(kg.m)^(1/2)";
  input Modelica.SIunits.MassFlowRate m_flow_turbulent(min=0) "Mass flow rate";
  output Modelica.SIunits.MassFlowRate m_flow 
    "Mass flow rate in design flow direction";

algorithm 
      m_flow := smooth(1, if dp>m_flow_turbulent^2/k/k then k*sqrt(dp) else 
                          if dp<-m_flow_turbulent^2/k/k then -k*sqrt(-dp) else 
                             k*sqrt(m_flow_turbulent^2/k/k)*sin(dp/(m_flow_turbulent^2/k/k)*Modelica.Constants.pi/2));

The speed is then:

Nonlin sys(5)  :    10,        3.607 s,        5.96 (       0.00 to      236.00),   605328

[mwetter]

These are nice speed-ups. It would be good to investigate how expensive evaluating sin() is compared to using a polynomial in this region.
As these functions are used so often, it may make sense to provide a different implementation for actuators and for fixed resistances if pre-computing polynomial coefficients is possible for fixed resistances only.
I think @marcusfuchs will also be looking into this.

You also need to use noEvent to guard against evaluating sqrt for negative numbers. See Modelica Specification 3.3, example 3.3.1

[Mathadon]

Based on the following example

model test
  Real y;
  Integer case = 1
    annotation(Evaluate=true);
algorithm 

  for i in 1:10000 loop
    if case ==1 then
      y:=3.44*time + 5.44*time^2+5.44*time^3;
    elseif case == 2  then
      y:=3.44*time + 5.44*time^2+5.44*time^3 +5.44*time^4;
    elseif case == 3 then
      y:=sin(1.34*time);
    else
      y:=time;
    end if;
  end for;
  annotation (
    uses(Modelica(version="3.2.1")),
    experiment(
      StopTime=100,
      __Dymola_fixedstepsize=0.1,
      __Dymola_Algorithm="Euler"),
    __Dymola_experimentSetupOutput(events=false));
end test;

I find:
Third degree polynomial with zero offset: 0.26 s
Fourth degree polynomial with zero offset: 0.44 s
sine: 0.21 s
overhead for loop: 0.02 s

[Mathadon]

The way I understood it the noEvent does not guard against division by zero. It guards against events being generated by the if clause that guards against division by zero. Is this correct?

The code contains the smooth function. Dymola automatically converts this to a noEvent (see dsmodel.mof) but we can explicitly add it too.
Code excerpt:

smooth(1, (if noEvent(res[2].dp > 0.2699999999999999)
       then 0.5773502691896258*sqrt(res[2].dp) else (if noEvent(res[2].dp < 
      -0.2699999999999999) then (-0.5773502691896258)*sqrt( -res[2].dp) else 0.3
      *sin(5.817764173314433*res[2].dp))));

[Mathadon]

Code example is on https://github.com/iea-annex60/modelica-annex60/tree/issue279_flowFunctions

Also, I have the impression that the old implementation does not comply with the Modelica specification:

The inverse requires that for all valid values of the input arguments of f2(...,y, ...) and uk being calculated as uk := f2(..., y, ...) implies the equality y = f1(..., uk, ...,) up to a certain precision.

This is currently not the case in the regularisation region as illustrated by Annex60.Fluid.BaseClasses.FlowModels.Examples.InverseFlowFunction (see variable deltaDp). The new version does have this property.

[Mathadon]

Here is an alternative that has a continuous first derivative:

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

computation time is marginally faster than when using the sinusoid function.

A possible problem is that this function's inverse is pretty expensive to evaluate:

(1/3)*((-54*m*dp_turbulent^(5/2)+3*dp_turb^2*sqrt(-3* dp_turbulent*(125* dp_turbulent*k^2-108*m^2)))*k^2)^(1/3)/k+5* dp_turbulent ^2*k/((-54*m* dp_turbulent ^(5/2)+3* dp_turbulent ^2*sqrt(-3* dp_turbulent*(125* dp_turbulent*k^2-108*m^2)))*k^2)^(1/3)

a solution may be to allow that f(dp) is not the exact inverse of f(m_flow) in the turbulent region as is currently the case.

[Mathadon]

Some more analyses:

implementation for m_flow

Annex60.Fluid.BaseClasses.FlowModels.basicFlowFunction_m_flow can be implemented using either:

 dp :=smooth(2, if noEvent(m_flow>m_flow_turbulent) then (m_flow/k)^2 else
                if noEvent(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);

or even:

 dp :=smooth(10, if noEvent(m_flow>0)
                then (m_flow/k)^2
                else -(m_flow/k)^2);

*** Warning: InverseFlowFunction.mat: deltaDp has absolute and relative error = 2.665e-02, 8.523e-01.

When running the first option together with the proposed implementation of Annex60.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp and LateInline=true unit tests give identical results. So it seems that this implementation is (almost?) identical to what MSL does, although it is more efficient.

Inlining

Setting LateInline=true is not sufficient for the function to be inlined. The annotation must contain Inline=true. However, when using Inline=true the function will always be inlined and therefore the function's inverse will not be used, where this may otherwise be possible (i.e. only for non-linear systems of 1 equation as far as I know).
I suggest to use Inline=true and to use from_dp where needed for avoiding a function inverse to be required.

This causes reference results to change since in some models the function inverse will not be used, where this was the case earlier. Since the normal function and its inverse are not identical, this will change the reference results.

I'm happy with this solution, so I'll already make a pull request.

Edit: My building model's speed increased by ~ 11% due to this measure.

[marcusfuchs]

Looks like a good improvement from @Mathadon , looking forward to discuss this on the call later today.

[marcusfuchs]

Also, I tested this new flow function for the pipes example in #264 and this implementation seems to solve the problem we were experiencing with different mass flow rates for 2 x 50 m pipes vs. 1 x 100 m pipe.

[mwetter]

@marcusfuchs will review, ensure that it is twice continuously differentiable and add a unit test for the 1st and 2nd derivative.

[marcusfuchs]

I added a test case for the 1st- and 2nd-order derivative each.

The 1st-order derivative check works for from_dp = true and from_dp = false.

The 2nd-order derivative check works for from_dp = true but fails for from_dp = false. So if the tests are implemented correctly, the flowFunction_dp would be twice differentiable, the flowFunction_mflow only once. I will double-check this tomorrow.

[marcusfuchs]

Also the 2nd-order derivative check works with from_dp = false on the dassl solver if the solver tolerance is decreased from 1e-4 to 1e-6.

[mwetter]

@marcusfuchs
I renamed the examples using

import buildingspy.development.refactor as r

r.move_class("Annex60.Fluid.BaseClasses.FlowModels.Examples.InverseFlowFunction",
             "Annex60.Fluid.BaseClasses.FlowModels.Validation.InverseFlowFunctions")
r.move_class("Annex60.Fluid.BaseClasses.FlowModels.Examples.TestFlowFunctions_1st_derivative",
             "Annex60.Fluid.BaseClasses.FlowModels.Validation.FirstDerivative")
r.move_class("Annex60.Fluid.BaseClasses.FlowModels.Examples.TestFlowFunctions_2nd_derivative",
             "Annex60.Fluid.BaseClasses.FlowModels.Validation.SecondDerivative")
r.move_class("Annex60.Fluid.BaseClasses.FlowModels.Examples.TestFlowFunctions",
             "Annex60.Fluid.BaseClasses.FlowModels.Validation.FlowFunctions")

as they are validation tests and not examples. I also added scripts for the regression tests.

[bramvdh91]

Currently we are again experiencing different mass flows during laminar/transient flow regimes for two 50 m pipes in series compared to a single pipe of 100 m, whereas they should be the same. Please refer to bramvdh91#24 for an example of the issue.

The problem only appears when use_dh is set to true. When use_dh = false and the nominal mass flow is adapted such that m_flow_turbulent is the same for both cases, the problem also does not appear.