id	summary	reporter	owner	description	type	status	priority	milestone	component	version	resolution	keywords	cc
1692	Non-linear system causes stack protection issues	Martin Sjölund	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"	defect	accepted	high		Backend				Martin Sjölund Frenkel TUD Willi Braun