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

[Jens Frenkel]

After a closer look at the equations the matching algorithm seems to produce a bad solution because some variables should be solved as function inputs. Reordering the system by hand shows a sulution without solving variables as function inputs:

split\[1].port_a.h = 2e-016 * join[1].port_a.h + 1.0 * join[1].port_c.h
join[1].port_b.h = 1.0 * split[1].port_b.h + 6.66666666666667e-017 * split[1].port_c.h
(pipe[1].outputPressure, pipe[1].loopPressure) := Pipe.func(100000.0,pipe[1].density)
(pipe[2].outputPressure, pipe[2].loopPressure) := Pipe.func(100000.0,pipe[2].density)
sink.port.h = Medium.h_pt(sink.p,300.0)
sink.p = min(pipe[2].outputPressure,pipe[1].outputPressure)
pipe[1].density = Medium.d_ph(100000.0,join[1].port_c.h)
pipe[2].density = Medium.d_ph(100000.0,join[1].port_a.h)
split[1].port_b.h = 0.666666666666667 * source.port.h + 0.333333333333333 * join[1].port_c.h
join[1].port_c.h = 1.0 * split[1].port_b.h + 6.66666666666667e-017 * sink.port.h
split[1].port_c.h = 1.0 * source.port.h + 1e-016 * join[1].port_a.h
join[1].port_a.h = 0.5 * (sink.port.h + split[1].port_c.h)

Looks like the function transforming the equations to residual form fails.

[Jens Frenkel]

Now with tearing the system is reduced to size 4 (from size 12), so I asume the problem is solved.

torn nonlinear Equationsystem:
1: pipe[2].outputPressure:VARIABLE(protected = true ) .Test, .Pipe, .Real type: Real [2]
2: pipe[2].loopPressure:VARIABLE() .Test, .Pipe, .Real type: Real [2]
3: split[1].port_c.h:VARIABLE() .Test, .Split, .Port, .Real type: Real [1]
4: pipe[1].outputPressure:VARIABLE(protected = true ) .Test, .Pipe, .Real type: Real [2]
5: pipe[1].loopPressure:VARIABLE() .Test, .Pipe, .Real type: Real [2]
6: sink.p:VARIABLE() .Test, .Boundary, .Real type: Real
7: split[1].port_b.h:VARIABLE() .Test, .Split, .Port, .Real type: Real [1]
8: sink.port.h:VARIABLE() .Test, .Boundary, .Port, .Real type: Real
1: join[1].port_c.h:VARIABLE() .Test, .Join, .Port, .Real type: Real [1]
2: join[1].port_a.h:VARIABLE() .Test, .Join, .Port, .Real type: Real [1]
3: pipe[1].density:VARIABLE(protected = true ) .Test, .Pipe, .Real type: Real [2]
4: pipe[2].density:VARIABLE(protected = true ) .Test, .Pipe, .Real type: Real [2]
1 : pipe[2].outputPressure := $_start(pipe[2].outputPressure);

pipe[2].loopPressure := $_start(pipe[2].loopPressure);
(pipe[2].outputPressure, pipe[2].loopPressure) := Pipe.func(100000.0, pipe[2].density);

2 : split[1].port_c.h = 1.0 * source.port.h + 1e-16 * join[1].port_a.h
3 : pipe[1].outputPressure := $_start(pipe[1].outputPressure);

pipe[1].loopPressure := $_start(pipe[1].loopPressure);
(pipe[1].outputPressure, pipe[1].loopPressure) := Pipe.func(100000.0, pipe[1].density);

4 : sink.p = min(pipe[2].outputPressure, pipe[1].outputPressure)
5 : split[1].port_b.h = 0.666666666666667 * source.port.h + 0.333333333333333 * join[1].port_c.h
6 : sink.port.h = Medium.h_pt(sink.p, 300.0)
1 : join[1].port_c.h = 1.0 * split[1].port_b.h + 6.66666666666667e-17 * sink.port.h
2 : pipe[1].density = Medium.d_ph(100000.0, join[1].port_c.h)
3 : join[1].port_a.h = 0.5 * (sink.port.h + split[1].port_c.h)
4 : pipe[2].density = Medium.d_ph(100000.0, join[1].port_a.h)

[Martin Sjölund]

Sounds good