How to extract values from algorithm section?

Question

I'm working with parabolic PDE's mixed with algebraic equation plus all these equation are coupled. I used Euler Method (Dassl is too slow) and big tolerance(for fast simulation) and recive error (for both types) in OM(1.9.1)(" residualFcn[some number]").The problem is that solver can't solve nonlinear system (mathematically, system is correct ). First question is what type of method use Euler method of integration in OM(explicit or implicit or Crank-Nicholson..or...)? So i tried to solve it numerically(explicit Euler Method (in code below "new[N]") (maybe the problem can be CFL condition.), but i have the problem (sample reconstruction for specific sample time). So, second question refers to reproducing values for specific sample time?! In code below there is array "a[3]". Idea is for every "ts(sample time)" reconstruct values in 3 nodes. How can this be done? How can i pass from algorithm section current values(i.e. for nodes) to equation section? I have all values in .txt file and i used Matlab to plot them, but i don't know how to pass it i.e. in "equation" section, or some another way. Same problem is for "new[N]", for specific sample time,plot function of nodes (N).

One more thing, if delta(t)/(delta(x))^2 >= 0.5 (delta(t) define the user, and refers to equation section, delta(x) is as in code below and spatial discretization is used in equation section (classical feedforward method)), is numerical stability satisfied? Same example but for algorithm section? Regards

Here is code:

model Euler1D
import Modelica.Utilities.*;
parameter Integer N=10; //50
parameter Real Lp=1e-6;
parameter Real deltax=1/(N-1)*Lp;
Real a[3];
Real old[N];
Real new[N];
Real b;
equation 
a[1]=if 
   (time>5) then 0 else time+5;
a[2]=time;
a[3]=2;
when 
(sample(0,1)) then
d=b;
end when;
algorithm 
// IN t=ts;
 when (sample(0,1)) then
 for i in 1:2 loop
  b:=a[i];
 Streams.print(String(time)+"   "+String(a[i])+ "     "+String(b), "C:/Some_Path/text.txt");
end for;
end when;

// Another problem
old[1]:=10;
old [N]:=0;
new[1]:=10;
new [N]:=0;
// Boundary
for i in 2:N-1 loop
old [i]:=10;
new[i]:=10;
end for;

for dx in deltax:deltax:Lp-deltax loop  // spatial discretization

for i in 2:N-1 loop
(new[i]):=(old[i]+0.5*(old[i + 1] +old[i-1]- 2*old[i]));
//def:=def+abs(new[i]-old[i]);
end for;
 for i in 2:N-1 loop
 old[i]:=new[i]; // switch the values
 end for;

 for i in 1:N loop
 Streams.print(String(time)+"   "+ String(new[2]), "C:/Some_Path/Anel_Nodes.txt");
 end for;

 annotation (uses(Modelica(version="3.2")));
 end Euler1D;

score 3 · Answer 1 · answered Feb 27 '13 at 15:16

Given your comments, this is how I would structure the model:

model Euler1D 
  parameter Integer N=10 "Spatial discretization";
  parameter Modelica.SIunits.Length L=1;
protected 
  parameter Modelica.SIunits.Length dx=L/(N-1) "Segment size";
  Real c[N] "Solution variable c";
  Real J[N] "Solution variable J";
  Real dc_dx "Spatial derivative at x==L";
  Real d2c_dx2[N] "Second spatial derivative of c";
initial equation 
  der(c[2:N-1]) = zeros(N-2);
equation 
  // Equations for spatial derivatives
  d2c_dx2[2:N-1] = { (c[i+1]-2*c[i]+c[i-1])/(dx^2) for i in 2:N-1};
  dc_dx = (c[N]-c[N-1])/dx;

  // Boundary conditions
  c[1] = 0 "Dirichlet B.C.";
  dc_dx = 1 "Neumann B.C.";

  // PDE
  J = c .* c "J = c^2";
  der(c) = d2c_dx2+J "PDE";
end Euler1D;

The initial equations avoid a big transient at the start of the simulation. So effectively, the model as I've posed it will give you a steady state solution. But you can make things time varying (e.g. the boundary conditions).

There could still be some errors in this, but it is hard for me to say because I don't have much intuition about this system and I don't have anyway to know if the values I've used are physically meaningful.

But hopefully it outlines how you could structure such a model. This model runs in Dymola. I did not try it in OpenModelica.

All that being said, you alway shave this issue with PDEs that there are implicit constraints that arise from the physics that the time integration scheme cannot explicitly see (e.g. the CFL condition). In general, variable time step algorithms have error estimation and they can indirectly see the instabilities introduced by violation of these implicit constraints. But they have to "fumble" through them. I am not sure to what extent this will be an issue with the model above.

If you want to add a little excitement to your solution, you can change the equation

  dc_dx = 1 "Neumann B.C.";

...to...

  dc_dx = 0.7+0.2*sin(time) "Neumann B.C.";

...and you will get a nice time varying solution. I noticed, when playing around with this, that if you make the gradient on the end too steep, the model goes unstable. Whether that is the true solution to the underlying mathematical equations or a numerical artifact due to the issue I mentioned above, I have no idea. So there are limits to what you can put into this model if you want a stable solution.

How to extract values from algorithm section?

1 Answers1