Investigating regfun3

[TManikantan]

          OK, I made some more enquiries using the very nice OMC [algorithmic debugger](https://openmodelica.org/doc/OpenModelicaUsersGuide/latest/debugger.html#the-algorithmic-debugger).

This MWE replicates the problem with regFun3() at the point were the mass flow rate function jumps:

model TestRegFun3
  Real x = 1.5023328;
  Real x0 = 0;
  Real x1 = 1.5023328811860801;
  Real y0 = 0;
  Real y1 = 0.0015270174250288026;
  Real yd0 = 0.0010172719108540143;
  Real yd1 = 0.0010164308084791662;
  Real y = Modelica.Fluid.Utilities.regFun3(x, x0, x1, y0, y1, yd0, yd1);
end TestRegFun3;

Since x is very close to x1, one would expect y to be very close to y1 = 0.0015270174250288026. Instead one gets y = 0.00137275.

What I understand is that the derivative yd1 is very, very close to (y1-y0)/(x1-x0), while the derivative yd0 is slightly larger. I can replicate the issue with round figures:

model TestRegFun3_round
  Real x = time;
  Real x0 = 0;
  Real x1 = 1;
  Real y0 = 0;
  Real y1 = 0.001;
  Real yd0 = 0.0010008;
  Real yd1 = 0.001;
  Real y = Modelica.Fluid.Utilities.regFun3(x, x0, x1, y0, y1, yd0, yd1);
  annotation(experiment(StopTime = 1, StartTime = 0, Tolerance = 1e-06, Interval = 0.002));
end TestRegFun3_round;

Variable y should go from y0 to y1, and quite obviously it doesn't do that.

Bottom line: regFun3() is numerically ill-conditioned when the two supplied derivatives are very close to the slope of the curve connecting the two points. I'm not sure whether this is inherent to the algorithm by Gasparo and Morandi or if it is just @sielemann's implementation that has some glitch in those conditions.

At this point this is no longer a Fluid issue, but a plain old computer science problem, so I gladly hand it over to someone else who's more competent than I am. @sielemann, @gkurzbach, @beutlich any suggestion? @hubertus65 do you think Michael would be interested at plunging back into this problem?

Originally posted by @casella in https://github.com/modelica/ModelicaStandardLibrary/issues/3758#issuecomment-1513468578

[TManikantan]

Also read https://github.com/modelica/ModelicaStandardLibrary/issues/3758#issuecomment-1513765997 by @beutlich

[HansOlsson]

Also read #3758 (comment) by @beutlich

I agree that there is something there, but I also realized a simple fix, that both seems to work and make sense. After computing xi2 add:

      if  xi1 < x0 or xi2>x1 then
        // The limits don't make sense, just use linear interpolation
        y := (y1-y0)*(x-x0)/(x1-x0) + y0;
      else
        ...
      end if;

[TManikantan]

Do we looker deeper into this issue for this milestone or different milestone ?as we have already found a working solution.

[casella]

I guess this has lower priority now. I would fix other pending issues first.