noEvent(x>0) fails to guard against sqrt(x) for x<0

[mwetter]

Hi

The regression test

Buildings.Fluid.FixedResistances.Examples.FixedResistancesExplicit

fails with "Model error: Argument of sqrt(res22.dp) was -5 should be >= 0"

This is triggered from

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

According to the Modelica specification, this is the correct implementation. However, when reformulating as

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;

the example works. Hence, this looks like the implementation of noEvent is wrong in OpenModelica.

This fix makes the following tests work:

simulate(Buildings.Fluid.Examples.Performance.Example1v1);
simulate(Buildings.Fluid.Examples.Performance.Example1v2);
simulate(Buildings.Fluid.Examples.Performance.Example2);
simulate(Buildings.Fluid.HeatExchangers.Examples.DryCoilCounterFlowPControl);
simulate(Buildings.Fluid.FixedResistances.Examples.FixedResistancesExplicit);
simulate(Buildings.Fluid.Examples.ResistanceVolumeFlowReversal);

[sjoelund.se]

Well... dp can be <0, right? The algorithm says it checks against dp_turbulent, which may be non-zero I assume?

[mwetter]

dp is smaller than zero in the elseif branch, which is why there is a negative sign in sqrt(-dp).

It is true that we check against dp_turbulent, but because always dp_turbulent > 0, this guarantees

m_flow := if noEvent(dp>dp_turbulent /* > 0 */) then k*sqrt(dp)
             elseif noEvent(dp<-dp_turbulent /* < 0 */) then -k*sqrt(-dp)
             else (k^....

Hence, our implementation looks fine to me.

[rfranke]

The code generation just takes away noEvent, see template daeExpCall in OMCompiler/Compiler/Template/CodegenCFunctions.tpl:

  case CALL(path=IDENT(name="noEvent"), expLst={e1}) then
    daeExp(e1, context, &preExp, &varDecls, &auxFunction)

This means that either the backend has to generate specific code for expressions inside noEvent, or the code generation is wrong.

[sjoelund.se]

Replying to rfranke:

This means that either the backend has to generate specific code for expressions inside noEvent, or the code generation is wrong.

The backend produces special code for relations that do not use noEvent.

[sjoelund.se]

Replying to mwetter:

dp is smaller than zero in the elseif branch, which is why there is a negative sign in sqrt(-dp).

It is true that we check against dp_turbulent, but because always dp_turbulent > 0, this guarantees

m_flow := if noEvent(dp>dp_turbulent /* > 0 */) then k*sqrt(dp)
             elseif noEvent(dp<-dp_turbulent /* < 0 */) then -k*sqrt(-dp)
             else (k^....

Hence, our implementation looks fine to me.

It's not the elseif branch being taken: Model error: Argument of sqrt(res22.dp) was -5 should be >= 0 suggests it is the if-branch and res22.dp_turbulent < -5. Needs some additional looking into.

[mwetter]

Replying to sjoelund.se:

Replying to mwetter:

dp is smaller than zero in the elseif branch, which is why there is a negative sign in sqrt(-dp).

It is true that we check against dp_turbulent, but because always dp_turbulent > 0, this guarantees

m_flow := if noEvent(dp>dp_turbulent /* > 0 */) then k*sqrt(dp)
             elseif noEvent(dp<-dp_turbulent /* < 0 */) then -k*sqrt(-dp)
             else (k^....

Hence, our implementation looks fine to me.

It's not the elseif branch being taken: Model error: Argument of sqrt(res22.dp) was -5 should be >= 0 suggests it is the if-branch and res22.dp_turbulent < -5. Needs some additional looking into.

Yes, it is the if branch. Of course m_flow := if noEvent(dp>dp_turbulent and dp>0) then k*sqrt(dp) ... does not work either (and the additional test is not needed as dp>0 whenever dp>dp_turbulent).

[casella]

The reports in ​https://test.openmodelica.org/libraries/Buildings/BuildModelRecursive.html indicate that the failure no longer takes place, possibly because of the implemented workaround.

I set up a simple test case to reproduce the problem

model Test
  Real x(start = 1, fixed = true);
  Real y;
  Real k = 3;
equation 
  der(x) = -2;
  y = if noEvent(x > 0) then k*sqrt(x) else 0;
end Test;

which worked correctly without generating events, using OpenModelica-v1.13.0-dev-63-g453e1c7-64bit.

Feel free to reopen this ticket if the problem shows up again elsewhere.