Jump in Modelica.Fluid.Pipes.BaseClasses.WallFriction.LaminarAndQuadraticTurbulent.massFlowRate_dp_staticHead

[uschna]

As can be seen in the Funktion.txt model, there is a jump in Modelica.Fluid.Pipes.BaseClasses.WallFriction.LaminarAndQuadraticTurbulent.massFlowRate_dp_staticHead when the pressure difference exceeds dp_small. This causes problems in ModelicaTest.Fluid.TestPipesAndValves.BranchingPipes18 (version 4.0.0) with pipe2 in this case. Also for this case there should exist a continuous transition in this function between dp >= dp_a and dp_b < dp < dp_a. Funktion.txt

[beutlich]

Some observations:

• Also happens for Modelica.Fluid.Pipes.BaseClasses.WallFriction.Laminar.massFlowRate_dp_staticHead.
• The jump happens when dp is equal to dp_small (You can set dp_small=1.5 (instead of 1.5023328811860799) in the given example to observe the jump at time = 0.)
• The problem is in Modelica.Fluid when it first was added by 3487ffff6f5539b3e2608383cb1babeb51d797e2 for MSL 3.1.

[uschna]

Also for Modelica.Fluid.Pipes.BaseClasses.WallFriction.Detailed.massFlowRate_dp_staticHead?

[gkurzbach]

Is there anybody who cares about this?

[uschna]

There is also a jump in these functions at dp = 0. I don't know if there is an error in the interval 0 <= dp <= dp_small.

[uschna]

This problem is only for g_times_height_ab = 0 or small values (<= 1e-11), for g_times_height_ab >= 1e-10 this is continuous.

[uschna]

Reasons for these:

Incorrect values: In function massFlowRate_dp_staticHead: dp_zero < dp < dp_a: Computation of dm_flow_ddp_fric_zero using Modelica.Fluid.Utilities.regFun3: abs((y1d+y0d)-2 * Delta0)<100 * Modelica.Constants.eps Computation of m_flow using Modelica.Fluid.Utilities.regFun3: abs((y1d+y0d)-2 * Delta0)>=100 * Modelica.Constants.eps, abs(y0d-y1d)>100 * Modelica.Constants.eps, (mu>0 and eta<h0 and Delta0 * omega<=0 or abs(aux01)<abs(aux02)and aux02 * Delta0>=0 or abs(aux01)<abs(0.10000000000000001 * Delta0))and not useSingleCubicPolynomial, xi1 < x < xi2

If dp_b < dp < dp_zero or dp_zero < dp < dp_a and in the computation of dm_flow_ddp_fric_zero also abs((y1d+y0d)-2 * Delta0)>=100 * Modelica.Constants.eps: correct values

[TManikantan]

@casella kindly look into it

[casella]

I'll have a look at it ASAP. The issue is non-trivial, and I did not write the original code myself.

[casella]

This is kind of odd.

The proposed test case sets the 8th input of massFlowRate_dp_staticHead g_times_height_ab = 0, i.e., the pipe is horizontal. This simplifies the analysis of massFlowRate_dp_staticHead considerably, because

dp_grav_a = 0;
dp_grav_b = 0;
dp_zero = 0;
dp_a = dp_small;

with this choice of input.

Consider now the if-statement starting at this line. When dp_a = dp_small the first branch is active and you get

m_flow := Internal.m_flow_of_dp_fric(dp_small, rho_a, rho_b, mu_a, mu_b, length, diameter, crossArea, Re1, Re2, Delta);

if you take an infintesimally smaller value, then the third branch is active, and you get

(m_flow_a, dm_flow_ddp_fric_a) := Internal.m_flow_of_dp_fric(dp_a-dp_grav_a, rho_a, rho_b, mu_a, mu_b, length, diameter, crossArea, Re1, Re2, Delta);
(m_flow_b, dm_flow_ddp_fric_b) := Internal.m_flow_of_dp_fric(dp_b-dp_grav_b, rho_a, rho_b, mu_a, mu_b, length, diameter, crossArea, Re1, Re2, Delta);
(m_flow, dm_flow_ddp_fric_zero) := Utilities.regFun3(dp_zero, dp_b, dp_a, m_flow_b, m_flow_a, dm_flow_ddp_fric_b, dm_flow_ddp_fric_a);
if dp>dp_zero then
    m_flow := Utilities.regFun3(dp, dp_zero, dp_a, m_flow_zero, m_flow_a, dm_flow_ddp_fric_zero, dm_flow_ddp_fric_a);
else
    m_flow := Utilities.regFun3(dp, dp_b, dp_zero, m_flow_b, m_flow_zero, dm_flow_ddp_fric_b, dm_flow_ddp_fric_zero);
end if;

Now, considering that dp = dp_small, dp_grav_a = 0, dp_grav_b = 0, and dp_zero = 0, the third branch can be simplified to

(m_flow_a, dm_flow_ddp_fric_a) := Internal.m_flow_of_dp_fric(dp_small, rho_a, rho_b, mu_a, mu_b, length, diameter, crossArea, Re1, Re2, Delta);
// irrelevant (m_flow_b, dm_flow_ddp_fric_b) := Internal.m_flow_of_dp_fric(dp_b-dp_grav_b, rho_a, rho_b, mu_a, mu_b, length, diameter, crossArea, Re1, Re2, Delta);
// irrelevant (m_flow, dm_flow_ddp_fric_zero) := Utilities.regFun3(dp_zero, dp_b, dp_a, m_flow_b, m_flow_a, dm_flow_ddp_fric_b, dm_flow_ddp_fric_a);
m_flow := Utilities.regFun3(dp_small, dp_zero, dp_small, m_flow_zero, m_flow_a, dm_flow_ddp_fric_zero, dm_flow_ddp_fric_a);

Now, function regFun3 is such that when the first input is equal to the second input, the output is equal to the fourth input; when the first input is equal to the third input, the output is equal to the fifth input.

Hence, the third branch can be further simplified as:

(m_flow_a, dm_flow_ddp_fric_a) := Internal.m_flow_of_dp_fric(dp_small, rho_a, rho_b, mu_a, mu_b, length, diameter, crossArea, Re1, Re2, Delta);
m_flow := m_flow_a;

and then further on as

m_flow := Internal.m_flow_of_dp_fric(dp_small, rho_a, rho_b, mu_a, mu_b, length, diameter, crossArea, Re1, Re2, Delta);

So, basically, one should get the same result from the first and third branch when dp is very close to dp_small. Apparently this doesn't happen, which means the regFun3 function is not working as expected.

[casella]

OK, I made some more enquiries using the very nice OMC 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?

[beutlich]

I spent some time with the test model TestRegFun3_round provided by @casella. What I observed:

1. The abscissas, ordinates and slopes are actually well-posed, meaning that the original cubic interpolation polynomial already satisfies the required monotonicity. This can be seen by the following verfications: a. Set useSingleCubicPolynomial=true in https://github.com/modelica/ModelicaStandardLibrary/blob/7d2e8d7a7e2738f9eaa0fe45ca21af54eb1e6ac9/Modelica/Fluid/Utilities.mo#L603 and check TestRegFun3_round.y (and the derivative). b. Call Modelica.Fluid.Utilities.cubicHermite instead of Modelica.Fluid.Utilities.regFun3 and check TestRegFun3_round.y (and the derivative). c. Call my test function cubicMonotoneInterpolation instead of Modelica.Fluid.Utilities.regFun3 and check TestRegFun3_round.y (and the derivative).

2. Reducing the magic threshold 1e-6 in the test for correction of https://github.com/modelica/ModelicaStandardLibrary/blob/7d2e8d7a7e2738f9eaa0fe45ca21af54eb1e6ac9/Modelica/Fluid/Utilities.mo#L694-L697 resolves the issue. I changed it to 1e-12, such that the correction no longer applies. That correction does not look correct to me and might be the root-cause of the bug.

3. The term "co-montone" or "comonotone" as used in Modelica.Fluid.Utilities.regFun3 is not a common term.

I propose to further check the thresholds of Modelica.Fluid.Utilities.regFun3 and especially simplify the well-posed case where the cubic Hermite interpolant satisfies monotonicity.

[sielemann]

Thanks for all the findings and work! Per https://github.com/modelica/ModelicaStandardLibrary/issues/3758#issuecomment-854501677 and https://github.com/modelica/ModelicaStandardLibrary/issues/3758#issuecomment-859379868, I dont see a lot of practical relevance of fixing this. If g_times_height_ab <= 1e-11, i.e., a horizontal pipe, the existing regularization functions can be used (without static head). I will therefore not work on this myself short term. I understand that, for consistency and flexibility of all corner cases, it is very desirable to solve this nonetheless. Hopefully, https://github.com/modelica/ModelicaStandardLibrary/issues/3758#issuecomment-1513765997 and https://github.com/modelica/ModelicaStandardLibrary/issues/3758#issuecomment-1513323969 already identified a way forward.

[beutlich]

@sielemann If I read your comment correctly, you claim that it is still a Fluid issue by using the regularization functions inappropriately in corner cases. This is in contradiction to @casella's comment who locates the issue in the regularization function itself. Someone needs to moderate. @TManikantan

[TManikantan]

Issue started with a discontinuous transition in m_flow calculation from dp. This might have been caused by using zero value for static head (seen in Funktion.txt) as pointed out by @uschna, so we agree with @sielemann that this is an unintended use of this function and should be closed (maybe with clear documentation in Modelica.Fluid.Pipes.BaseClasses.WallFriction.LaminarAndQuadraticTurbulent.massFlowRate_dp_staticHead that this function should not be used for horizontal tubes), especially when another function exists without the use of static head (LaminarAndQuadraticTurbulent.massFlowRate_dp). @casella investigated it and found a tangential issue (in our opinion) that regfun3 is not working as intended. We agree that this is not creating a continuous transition. So, our suggestion is to close the existing ticket this ticket with improved documentation regarding the use of static head but open another issue to investigate regfun3.

[beutlich]

Since this issue originates from the test model ModelicaTest.Fluid.TestPipesAndValves.BranchingPipes18 this models needs to be adapated as well (either by avoiding horizontal tubes or by fixing the calls of rgularization functions or by fixing regfun3).

[TManikantan]

We agree that ModelicaTest.Fluid.TestPipesAndValves.BranchingPipes18 would have to be adapted to have a non-zero static-head value. Should we go ahead and make those changes and commit? @beutlich Should we also add an assert statement in `Laminar.massFlowRate_dp_staticHead' for static-head value as well? Additionally, we have created an issue to investigate regfun3.

[TManikantan]

@beutlich @casella @hubertus65 Do you agree with our proposed solution?

[hubertus65]

Hello Mani! Yes, I agree with your proposed solution. It's very good to now have met you actually in person!

[TManikantan]

Thankyou for your response @hubertus65,Likewise I've been looking forward to meeting you too. Thanks for all great interactions.

[beutlich]

@beutlich @casella @hubertus65 Do you agree with our proposed solution?

Delegating to the issue authors @uschna and @gkurzbach whos should be more suitable in answering.

[casella]

Since this issue originates from the test model ModelicaTest.Fluid.TestPipesAndValves.BranchingPipes18 this models needs to be adapated as well (either by avoiding horizontal tubes or by fixing the calls of rgularization functions or by fixing regfun3).

Avoiding horizontal tubes is IMHO not a good idea. In most cases of thermofluid models, tubes are indeed horizontal, models that fail in that case are not very good quality.

[casella]

We agree that ModelicaTest.Fluid.TestPipesAndValves.BranchingPipes18 would have to be adapted to have a non-zero static-head value.

Again, unless I am missing some point (this issue is quite intricate) I don't really agree with this. Horizontal pipe models are not singular by any means, if some singularity is introduced by regFun, this should be handled by improving (or avoiding) regFun, not by giving some inclination to the pipes to sweep the problem under the rug 😅

[casella]

If (as suggested in #3758 comment) changing some thresholds allows to handle the horizontal tube case successfully, I am definitely in favour of that.

[uschna]

Because the creation of this ticket was some time ago, I am not sure. But I suspect the problem does not only affect ModelicaTest.Fluid.TestPipesAndValves.BranchingPipes18 but also ModelicaTest.Fluid.TestComponents.Pipes.DynamicPipesWithTraceSubstances (pipe8.flowModel.WallFriction.massFlowRate_dp_staticHead -> pipe8.port_a.m_flow). As an interim workaround, both models could be reviewed and adjusted accordingly. The goal should be to model the MSL types in a way that ensures a steady transition. This also applies to functions (massFlowRate_dp_staticHead).

[casella]

IMHO the right approach here is to make sure that regFun3 always produces the expected results, as stated in its documentation. If it doesn't, there is a hole in the library, and all kind of weird things can potentially happen when that function is used, which will take forever to debug.

If regFun3 doesn't work well in some corner case, but it's actually known what its output should be in those cases, just add an if clause to its the algorithm that checks for those corner case conditions and gives the correct result. Everything else will work out automatically.

As I wrote earlier, I'm not really sure how to do that, but it seems to me a quintessential computer science / numerical mathematics problem, not a physical modelling problem.

I'd leave the physical models alone, as well as their parameters. Also implementing tweaks to the MSL examples to avoid generating this issue seems to me a bit like sweeping the dirt under the rug 😅

[uschna]

I didn't understand offhand the algorithm in Modelica.Fluid.Utilities.regFun3. Therefore, I don't know what the problem is with the discontinuity at the edges and whether this can be fixed in Modelica.Fluid.Utilities.regFun3 or needs to be corrected in the calls in the various massFlowRate_dp_staticHead functions. If I understand correctly, the problem is handled in https://github.com/modelica/ModelicaStandardLibrary/issues/4128. My question is: Until this is resolved, will the models be adjusted? Thanks.

[casella]

I didn't understand offhand the algorithm in Modelica.Fluid.Utilities.regFun3.

Me neither 😅

Therefore, I don't know what the problem is with the discontinuity at the edges and whether this can be fixed in Modelica.Fluid.Utilities.regFun3 or needs to be corrected in the calls in the various massFlowRate_dp_staticHead functions.

I am strongly in favour of the programming by contract approach. Every piece of software assumes certain preconditions, and given them, guarantees certain postconditions. If the preconditions are not met by other pieces of software, that's not the fault of this piece of software. This is the only way to manage responsibility in large software projects, otherwise you get spaghetti programming and enduring chaos. This is even more important for a modular software library such as MSL.

From this point of view, the problem is very clear: the pipe model carefully avoids discontinuities by using specialized smoothing functions. If they don't work as they claim, as shown here, it's them that should be fixed.

What is less clear is how to fix regFun3, and that's the topic of #4128.

My question is: Until this is resolved, will the models be adjusted? Thanks.

If we consider this as an urgent issue (I'm not sure), then we could patch this somehow, e.g. by adding some inclination to the pipes. But it's not really urgent, we are not planning a new release of MSL tomorrow (are we?).

I'm not really in favour of this approach, because then it's very likely that #4128 will never be addressed, and the underlying Modelica.Fluid models will remain numerically fragile forever, for a very complicated reason (took us a while to figure out) in a condition (being horizontal) which is true in the large majority of applications. That's the worst possible outcome.

My recommendation is to get through #4128 ASAP. This shouldn't be a big deal, the problem has been clearly identified with a simple MWE and can be solved locally without the need of understanding the intricacies of the pipe model, which is quite involved. I would leave those models alone: they are perfectly fine, provided regFun3 does what it claims, i.e., interpolate smoothly.

[TManikantan]

I think this issue can be closed "close as not planned" and label it invalid as we are going to solve the regFun3 .

[beutlich]

I still like to keep it open. It is possible to close/address/reference more than one issue with a single PR.

[casella]

The original problem was with a certain model. The root cause is regFun3 not behaving, but this ticket remains valid IMHO.