Non-linear system causes stack protection issues

[Martin Sjölund]

​http://openmodelica.org/index.php?option=com_agora&task=topic&id=526&Itemid=87
{{{package Medium

function h_pt

input Real p;
input Real t;
output Real h;

algorithm

h := t;

end h_pt;

function d_ph

input Real p;
input Real h;
output Real d;

algorithm

d := h;

end d_ph;

end Medium;

connector Port

flow Real mf;
Real p;
stream Real h;

end Port;

model Pipe

Port port_a, port_b;
Real loopPressure;

protected

Real density;
Real outputPressure;

function func

input Real inputPressure;
input Real density;
output Real outputPressure;
output Real loopPressure;

algorithm

loopPressure := inputPressure;
outputPressure := inputPressure;

end func;

equation

density = Medium.d_ph(port_a.p, port_a.h);
(outputPressure, loopPressure) = func(port_a.p, density);

port_a.mf + port_b.mf = 0;
port_b.p = outputPressure;
port_a.h = inStream(port_b.h);
port_b.h = inStream(port_a.h);

end Pipe;

model Split

Port port_a, port_b, port_c;

equation

connect(port_a, port_b);
connect(port_a, port_c);

end Split;

model Join

Port port_a, port_b, port_c;
final parameter Real epsFlow = 1e-15;

equation

port_a.mf + port_b.mf + port_c.mf = 0;
port_b.p = min(port_a.p, port_c.p);
port_a.h = (max(port_b.mf, epsFlow) * inStream(port_b.h) + max(port_c.mf, epsFlow) * inStream(port_c.h)) / (max(port_b.mf, epsFlow) + max(port_c.mf, epsFlow));
port_b.h = (max(port_a.mf, epsFlow) * inStream(port_a.h) + max(port_c.mf, epsFlow) * inStream(port_c.h)) / (max(port_a.mf, epsFlow) + max(port_c.mf, epsFlow));
port_c.h = (max(port_a.mf, epsFlow) * inStream(port_a.h) + max(port_b.mf, epsFlow) * inStream(port_b.h)) / (max(port_a.mf, epsFlow) + max(port_b.mf, epsFlow));

end Join;

model Boundary

Port port;
Real t;
Real mf;
Real p;

equation

port.mf = mf;
port.p = p;
port.h = Medium.h_pt(p, t);

end Boundary;

model Test

Split[numberOfSplits - 1] split;
Join[numberOfSplits - 1] join;
Pipe[numberOfSplits] pipe;

Boundary source;
Boundary sink;

parameter Integer numberOfSplits = 2;
parameter Integer numberOfLoops = 10;

equation

source.mf = -10;
source.p = 100000;
source.t = 300;
sink.t = 300;

Grid:
split
--+--+---
| | | pipe
--+--+---
join

connect(source.port, split[1].port_a);

for i in 1:numberOfSplits - 1 loop

connect(split[i].port_c, pipe[i].port_a);
connect(pipe[i].port_b, join[i].port_c);

split[i].port_c.mf = -source.port.mf / numberOfSplits;

end for;

for i in 1:numberOfSplits - 2 loop

connect(split[i].port_b, split[i + 1].port_a);
connect(join[i + 1].port_b, join[i].port_a);

end for;

connect(split[numberOfSplits - 1].port_b, pipe[numberOfSplits].port_a);
connect(pipe[numberOfSplits].port_b, join[numberOfSplits - 1].port_a);

connect(join[1].port_b, sink.port);

end Test; }}}

{{{$ MODELICAUSERCFLAGS="-g -fstack-protector" omc +s +simCodeTarget=Dump a.mo | grep residual | wc -l
16}}}

{{{$ MODELICAUSERCFLAGS="-g -fstack-protector" omc +s +simCodeTarget=Dump a.mo | grep nonlinear | tr ,
n | cut -d: -f2 | tr -d " " | wc -l
12}}}

So we have 16 residual expressions written to 12 memory locations. WRONG WRONG WRONG.

Also, many of the expressions are KNOWN ZERO!
{{{ residual: xloc[2];

residual: xloc[1];
residual: xloc[11] - Medium.d_ph(100000.0,xloc[10]);
residual: xloc[9] + (-0.9999999999999998 * source.port.h + -0.00000000000000009999999999999999 * xloc[10]);
residual: xloc[10] + -0.5 * (xloc[8] + xloc[9]);
residual: xloc[6] + (-0.9999999999999999 * xloc[7] + -0.00000000000000006666666666666667 * xloc[8]);
residual: xloc[7] + (-0.6666666666666666 * source.port.h + -0.3333333333333333 * xloc[6]);
residual: xloc[5] - Medium.d_ph(100000.0,xloc[6]);
residual: xloc[3];
residual: xloc[4];
residual: xloc[3];
residual: xloc[4];
residual: xloc[0] - min(xloc[2],xloc[3]);
residual: xloc[2];
residual: xloc[1];
residual: xloc[8] - Medium.h_pt(xloc[0],300.0);}}}

We can replace xloc[1]=xloc[2]=xloc[4]=xloc[3]=0 since this is a trivial solution (and the only correct one). We can then replace all these variables with the known constant (like removeSimpleExpressions), and re-run SCC algorithm since this may no longer be a SCC :)

In this case:
{{{ residual: xloc[11] - Medium.d_ph(100000.0,xloc[10]);

residual: xloc[9] + (-0.9999999999999998 * source.port.h + -0.00000000000000009999999999999999 * xloc[10]);
residual: xloc[10] + -0.5 * (xloc[8] + xloc[9]);
residual: xloc[6] + (-0.9999999999999999 * xloc[7] + -0.00000000000000006666666666666667 * xloc[8]);
residual: xloc[7] + (-0.6666666666666666 * source.port.h + -0.3333333333333333 * xloc[6]);
residual: xloc[5] - Medium.d_ph(100000.0,xloc[6]);
residual: xloc[0] - 0.0;
residual: xloc[8] - Medium.h_pt(xloc[0],300.0);}}}

In this case (re-run xloc[0]=0:
{{{ residual: xloc[11] - Medium.d_ph(100000.0,xloc[10]);

residual: xloc[9] + (-0.9999999999999998 * source.port.h + -0.00000000000000009999999999999999 * xloc[10]);
residual: xloc[10] + -0.5 * (xloc[8] + xloc[9]);
residual: xloc[6] + (-0.9999999999999999 * xloc[7] + -0.00000000000000006666666666666667 * xloc[8]);
residual: xloc[7] + (-0.6666666666666666 * source.port.h + -0.3333333333333333 * xloc[6]);
residual: xloc[5] - Medium.d_ph(100000.0,xloc[6]);
residual: xloc[8] - Medium.h_pt(0.0,300.0); /* Another known constant! */}}}

In this case (re-run xloc[8]=300.0:
{{{ residual: xloc[11] - Medium.d_ph(100000.0,xloc[10]);

residual: xloc[9] + (-0.9999999999999998 * source.port.h + -0.00000000000000009999999999999999 * xloc[10]);
residual: xloc[10] + -0.5 * (300.0 + xloc[9]);
residual: xloc[6] + (-0.9999999999999999 * xloc[7] + -0.00000000000000006666666666666667 * 300.0);
residual: xloc[7] + (-0.6666666666666666 * source.port.h + -0.3333333333333333 * xloc[6]);
residual: xloc[5] - Medium.d_ph(100000.0,xloc[6]);}}}

Much simpler to solve :) Feel free to split this into 2 different bugs if needed.

Re-sort:

Linear (?):
  residual: xloc[9] + (-0.9999999999999998 * source.port.h + -0.00000000000000009999999999999999 * xloc[10]);
  residual: xloc[10] + -0.5 * (300.0 + xloc[9]);
Scalar:
  residual: xloc[11] = Medium.d_ph(100000.0,xloc[10]);
Linear (?):
  residual: xloc[6] + (-0.9999999999999999 * xloc[7] + -0.00000000000000006666666666666667 * 300.0);
  residual: xloc[7] + (-0.6666666666666666 * source.port.h + -0.3333333333333333 * xloc[6]);
Scalar:
  residual: xloc[5] = Medium.d_ph(100000.0,xloc[6]);}}}

This equation system screams: Please make me faster, Willi and/or Jens