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

[Michael Wetter]

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);

[Martin Sjölund]

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

[Michael Wetter]

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.

[Rüdiger Franke]

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.

[Martin Sjölund]

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.

[Martin Sjölund]

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.

[Michael Wetter]

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).

[Francesco 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.