OpenModelica устраняет ошибку инициализации PDE

Question

Я пытаюсь использовать OpenModelica для численного решения очень простого PDE du/dx=du/dt с граничным условием u(0,t)=t^2 и u_x(0,t)=0. Я написал код ниже:

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;

Он компилируется, однако не завершает симуляцию с ошибкой ниже:

Blocstdout | OMEditInfo |

C: / Users /.../ AppData / Local / Temp / OpenModelica / OMEdit / pdetest_1.exe -port = 50450 -logFormat = xmltcp -override = startTime = 0, stopTime = 1, stepSize = 0.002, допуск =1e-6, решатель =dassl,outputFormat=mat,variableFilter=.* -r=pdetest_1_res.mat -jacobian= ColoredNumeric -w -lv = LOG_STATS

kquote LOG_INIT | ошибка |

Проблема инициализации противоречива из-за следующего уравнения: 0!= 0,000204061 = u[4]

стандартный вывод | предупреждение |

Ошибка в инициализации. Сохранение результатов и выход.
Используйте -lv=LOG_INIT -w для получения дополнительной информации.

стандартный вывод | ошибка |

Процесс моделирования не удался. Выход с кодом -1.

Буду признателен, если вы поможете мне узнать, в чем проблема и как я могу ее решить?

Answer 1

Хорошо, во-первых, очень грустно видеть, что сообщество Modelica настолько оцепенело по этому вопросу. Здесь, в SO или на форуме OpenModelica, есть дюжина вопросов, связанных с PDE, но не многие с правильным ответом. Я решил сделать это Github-репо, собирая все соответствующие материалы, которые я мог найти по всему Интернету, так что, по крайней мере, другим людям не придется обдумывать рабочий пример.

Но о коде выше. Код почти в порядке, и проблема заключается в физике проблемы. Я задал вопрос в вычислительной науке и получил очень хороший ответ.

Рабочий код:

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;

Я использовал код в этой презентации Яна Силара, чтобы решить эту проблему. и я также упомянул код в примере 4 упомянутого репозитория github.