Issue with current implementation of Stream connectors

[Francesco Casella]

I am investigating the root causes of #3118. One interesting test case is ThermoPower.Test.WaterComponents.TestValves. There is one fairly large nonlinear system of equations, this is the dumpSimCode output:

108:  (NONLINEAR) index:1 jacobian: false
		crefs: $cse16 [Real] , $cse12 [Real] , $cse17 [Real] , $cse6 [Real] , V5 [model .ValveLiq$V5] .fluidState [record .ThermodynamicState$V5$fluidState] .p [Real]
	79: V5.dp=V5.fluidState.p - SinkP1.p0[Real ]
	80: V3.dp=SourceP2.p0 - V5.fluidState.p[Real ]
	81: V2.dp=SourceP1.p0 - V5.fluidState.p[Real ]
	82: V3.rho=Modelica.Media.Water.IF97_Utilities.rho_props_ph(SourceP2.p0, DIVISION($cse6 * SourceP2.h + $cse12 * SinkP3.h, $cse6 + $cse12), Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(SourceP2.p0, DIVISION($cse6 * SourceP2.h + $cse12 * SinkP3.h, $cse6 + $cse12), 0, 0))[Real ]
	83: V3.w=homotopy(ThermoPower.Water.ValveLiq$V3.FlowChar(CloseValves.y) * V3.Av * sqrt(V3.rho) * ThermoPower.Water.ValveLiq$V3.sqrtR(V3.dp, 3000.0), CloseValves.y * V3.wnom * DIVISION(SourceP2.p0 - V5.fluidState.p, V3.thetanom * V3.dpnom))[Real ]
	84: $cse2=max(V3.w, 1e-015)[Real (start = -V3.wnom, max = if V3.allowFlowReversal then 9.999999999999999e+059 else 0.0, quantity = "MassFlowRate.WaterIF97", min = -100000.0, unit = "kg/s")]
	85: $cse7=max(-V3.w, 1e-015)[Real (start = V3.wnom, min = if V3.allowFlowReversal then -9.999999999999999e+059 else 0.0, quantity = "MassFlowRate.WaterIF97", max = 100000.0, unit = "kg/s")]
	86: V3.outlet.h_outflow=DIVISION($cse6 * SourceP2.h + $cse12 * SinkP3.h, $cse6 + $cse12)[Real ]
	87: V2.outlet.h_outflow=DIVISION($cse16 * SinkP2.h + $cse17 * SourceP1.h, $cse16 + $cse17)[Real ]
	88: V2.rho=Modelica.Media.Water.IF97_Utilities.rho_props_ph(SourceP1.p0, DIVISION($cse16 * SinkP2.h + $cse17 * SourceP1.h, $cse16 + $cse17), Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(SourceP1.p0, DIVISION($cse16 * SinkP2.h + $cse17 * SourceP1.h, $cse16 + $cse17), 0, 0))[Real ]
	89: V2.w=homotopy(ThermoPower.Water.ValveLiq$V2.FlowChar(CloseValves.y) * V2.Av * sqrt(V2.rho) * ThermoPower.Water.ValveLiq$V2.sqrtR(V2.dp, 5000.0), CloseValves.y * V2.wnom * DIVISION(SourceP1.p0 - V5.fluidState.p, V2.thetanom * V2.dpnom))[Real ]
	90: V5.w=V2.w + V3.w[Real ]
	91: $cse11=max(-V5.w, 1e-015)[Real (start = V5.wnom, min = if V5.allowFlowReversal then -9.999999999999999e+059 else 0.0, quantity = "MassFlowRate.WaterIF97", max = 100000.0, unit = "kg/s")]
	92: V2.inlet.h_outflow=DIVISION($cse11 * SinkP1.h + $cse2 * V3.outlet.h_outflow, $cse11 + $cse2)[Real ]
	93: $cse1=max(V2.w, 1e-015)[Real (start = -V2.wnom, max = if V2.allowFlowReversal then 9.999999999999999e+059 else 0.0, quantity = "MassFlowRate.WaterIF97", min = -100000.0, unit = "kg/s")]
	94: V3.inlet.h_outflow=DIVISION($cse1 * V2.outlet.h_outflow + $cse11 * SinkP1.h, $cse1 + $cse11)[Real ]
	95: V4.rho=Modelica.Media.Water.IF97_Utilities.rho_props_ph(SourceP2.p0, DIVISION($cse6 * SourceP2.h + $cse7 * V3.inlet.h_outflow, $cse6 + $cse7), Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(SourceP2.p0, DIVISION($cse6 * SourceP2.h + $cse7 * V3.inlet.h_outflow, $cse6 + $cse7), 0, 0))[Real ]
	96: V4.w=homotopy(ThermoPower.Water.ValveLiq$V4.FlowChar(OpenRelief.y) * V4.Av * sqrt(V4.rho) * ThermoPower.Water.ValveLiq$V4.sqrtR(V4.dp, 8000.0), OpenRelief.y * V4.wnom * DIVISION(SourceP2.p0 - SinkP3.p0, V4.thetanom * V4.dpnom))[Real ]
	97: SourceP2.flange.m_flow=(-V4.w) - V3.w[Real ]
	98: $cse21=max(-V2.w, 1e-015)[Real (start = V2.wnom, min = if V2.allowFlowReversal then -9.999999999999999e+059 else 0.0, quantity = "MassFlowRate.WaterIF97", max = 100000.0, unit = "kg/s")]
	99: V1.rho=Modelica.Media.Water.IF97_Utilities.rho_props_ph(SourceP1.p0, DIVISION($cse21 * V2.inlet.h_outflow + $cse17 * SourceP1.h, $cse21 + $cse17), Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(SourceP1.p0, DIVISION($cse21 * V2.inlet.h_outflow + $cse17 * SourceP1.h, $cse21 + $cse17), 0, 0))[Real ]
	100: V1.w=homotopy(ThermoPower.Water.ValveLiq$V1.FlowChar(OpenRelief.y) * V1.Av * sqrt(V1.rho) * ThermoPower.Water.ValveLiq$V1.sqrtR(V1.dp, 9000.0), OpenRelief.y * V1.wnom * DIVISION(SourceP1.p0 - SinkP2.p0, V1.thetanom * V1.dpnom))[Real ]
	101: SourceP1.flange.m_flow=(-V2.w) - V1.w[Real ]
	102: V5.rho=Modelica.Media.Water.IF97_Utilities.rho_props_ph(V5.fluidState.p, DIVISION($cse1 * V2.outlet.h_outflow + $cse2 * V3.outlet.h_outflow, $cse1 + $cse2), Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(V5.fluidState.p, DIVISION($cse1 * V2.outlet.h_outflow + $cse2 * V3.outlet.h_outflow, $cse1 + $cse2), 0, 0))[Real ]
	107: max(-SourceP2.flange.m_flow, 1e-015) - $cse6 (RESIDUAL)
	106: homotopy(ThermoPower.Water.ValveLiq$V5.FlowChar(CloseLoad.y) * V5.Av * sqrt(V5.rho) * ThermoPower.Water.ValveLiq$V5.sqrtR(V5.dp, 4000.0), CloseLoad.y * V5.wnom * DIVISION(V5.fluidState.p - SinkP1.p0, V5.thetanom * V5.dpnom)) - V5.w (RESIDUAL)
	105: max(-V4.w, 1e-015) - $cse12 (RESIDUAL)
	104: max(-V1.w, 1e-015) - $cse16 (RESIDUAL)
	103: max(-SourceP1.flange.m_flow, 1e-015) - $cse17 (RESIDUAL)

The code shows several instances of max(someFlowRate, 1e-015), which is apparently the simplest implementation of the positiveMax() operator defined in Sect. 15.2 of the Modelica Specification. The spec is a bit vague as it refers to "a small flow" for the threshold eps1.

The underlying idea is that when flow rates are so small that their contributions to mass and energy balances are negligible, they can be approximated in the mixing equations to avoid singularities in the computation of stream variables. On the other hand, picking a flow threshold which is actually too small could lead to unnecessary rounding errors or convergence problems, if the stream variables are involved in bigger nonlinear systems.

The engineering models I am aware of have flow rates in SI units with order of magnitude between 1e-4 and 1e4, so that 1e-15 seems way to small even in the case of the smallest possible flows (e.g. in pharmaceutical plant applications).

I propose the following actions:

• change the default value of eps1 to 1e-7, which is still low enough for 99.9% of applications, while safely distant from machine precision
• allow to change the default of eps1 by a compiler flag, to allow advanced users find out the best value in case of critical problems
• use eps1*flowVariable.nominal as the threshold, in order to allow proper problem-dependent scaling

In fact, the last action will only be fully effective if proper nominal attribute propagation is also implemented for alias equations such as a = +/-b and a+/-b = 0, but I would implement it anyway even if this latter feature is still not available, because it might be provided later in the future.

[Adrian Pop]

I'm assigning this to Per. Hopefully he's not in holidays, but if he is I can probably do the change myself.

[Per Östlund]

Ok, I'll fix it tomorrow. The reason I chose 1e-15 was because the spec didn't say what eps1 was, so I simply used Modelica.Constants.eps. But I can easily change it to be 1e-7 by default and add a flag to set it with (name suggestions are welcome). I'll have to investigate the nominal stuff, but I think we propagate nominal values to the connection handling, so that shouldn't be a problem either.

[Francesco Casella]

Great, thanks!

I'm impressed by the reaction times of the development team. You guys rock!

[Per Östlund]

Fixed in ca563fc/OMCompiler. 1e-7 is now used by default, or 1e-7*flow.nominal if the flow variable has a nominal value. There's also a new flag --flowThreshold that can be used to change the threshold.