OpenModelica solving a PDE Initialization error

Question

I'm trying to use OpenModelica to numerically solve a very simple PDE du/dx=du/dt with boundary condition u(0,t)=t^2 and u_x(0,t)=0. I have written the code below:

model pdetest_1

    parameter Real L=1;
    parameter Integer N=100;
    parameter Real dx=L/(N-1);
    parameter Real[N] x=array(i*dx for i in 0:N-1);

    Real u[N],ux[N];

initial equation

    for i in 1:N loop
      u[i]=0;
    end for;

equation
    u[1]=(time)^2;
    ux[1]=0;

    for i in 2:N loop
      u[i]=u[i-1]+dx*ux[i-1];
      der(u[i])=ux[i];
    end for;

end pdetest_1;

It does compile however it does not finish the simulation quitting with the error below:

Blocstdout | OMEditInfo |

C:/Users/.../AppData/Local/Temp/OpenModelica/OMEdit/pdetest_1.exe -port=50450 -logFormat=xmltcp -override=startTime=0,stopTime=1,stepSize=0.002,tolerance=1e-6,solver=dassl,outputFormat=mat,variableFilter=.* -r=pdetest_1_res.mat -jacobian=coloredNumerical -w -lv=LOG_STATS

kquote LOG_INIT | error |

The initialization problem is inconsistent due to the following equation: 0 != 0.000204061 = u[4]

stdout | warning |

Error in initialization. Storing results and exiting.
Use -lv=LOG_INIT -w for more information.

stdout | error |

Simulation process failed. Exited with code -1.

I would appreciate if you could help me know what is the problem and how I can solve it?

Answer 1

Ok, first of all it is very sad to see that the Modelica community has been so numb about this subject. There are a dozen of PDE related questions here in SO or in OpenModelica forum, not many with a proper answer. I decided to make this Github repo collecting all the relevant materials I could find all over the internet, so at least other people wont have to wonder around for a working example.

But about the code above. The code is almost fine, and the issue lies in the physics of the problem. I asked the question in computational science and got a very good answer.

The working code is:

model pdetest_1
  parameter Real L = 1;
  parameter Integer N = 100;
  parameter Real dx = L / (N - 1);
  parameter Real c = 1;
  Real u[N], ux[N];
initial equation
  for i in 1:N loop
    u[i] = 0;
  end for;
equation
  if c>0 then
    u[N] = time ^ 2;
    ux[N] = 0;
    for i in 1:N-1 loop
      u[i] = u[i + 1] - dx * ux[i];
      der(u[i]) = c*ux[i];
    end for;
  else
    u[1] = time ^ 2;
    ux[1] = 0;
    for i in 2:N loop
      u[i] = u[i - 1] + dx * ux[i];
      der(u[i]) = c*ux[i];
    end for;
  end if;
end pdetest_1;

I used the code in this presentation by Jan Silar to solve the issue. and I also mentioned the code in the example 4 of the said github repo.