OpenModelicaで偏微分方程式を解く

# OpenModelicaで偏微分方程式を解く機能

いつ頃からか(多分OpenModelica1.12とか）OpenModelicaで偏微分方程式を解く機能が実装されました。
Modelica Conferenceなどでは発表がありましたが実はそんなに宣伝されていませんでした。
[勉強会](https://modelica-library.connpass.com/)でネタないかなーとOpenModelicaのドキュメントをあさっていると正式にUser Guideが整っていたのでとりあえず確認してみました。

偏微分方程式機能の名前はPDEModelicaという名称でModelica Specificationには無い機能をOpenModelicaで実装しているようです。

User Guide
https://openmodelica.org/doc/OpenModelicaUsersGuide/latest/pdemodelica.html

## 偏微分方程式機能(PDEModelica)で出来ること

偏微分方程式を解く、といっても解けるのは時間項ともう一つ任意独立変数の2変数偏微分方程式です。
三次元（x,y,z)を期待してはいけません。

また書いてはいませんが、元の論文からすると差分法を使っています。

# 設定方法

1.OpenModelicaを開いて「ツール」->「オプション」を開いてください。

2.「シミュレーション」をクリックし、「OMC Command Line Options」へ「--grammar=PDEModelica」と記入してください。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/70651d3d-2ecb-e9d8-b7dc-6ae1ddc2c954.png)


# 例題

## どんな方程式を解くか

以下のような方程式が例題として載っています

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/d26a54ab-8dd2-06e6-3021-03bc261eee89.png)
　
上記の方程式は図でいうと以下のようになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/c4e0b358-90f8-e374-70b6-547b9364db5b.png)

## Modelicaコード

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/34c4112b-8a8a-5b2d-ee88-f9f596de4675.png)

解説
3行目
parameter DomainLineSegment1D omega(L = 1, N = 100)  "domain";
DomainLineSegment1Dは領域を定義する型。
いわゆる独立変数
Lが長さ、Nが分割数

4行目
field Real u(domain = omega)                         "field";
fieldは解きたい場を定義する型
いわゆる従属変数

5-6行目
initial equation
  u = sin(2*pi*omega.x)                                "IC";
初期場の定義
DomainLineSegment1D内にxという変数があるのが分かる
ICはInitial Conditionの略

8行目
  der(u) + pder(u,x) = 0   indomain omega              "PDE";
解きたい偏微分方程式
indomain omega と方程式の後についているのが特徴的
決して記述ミスじゃないので消さないように
（最初消してしまって動きませんでした。。。）

9行目
  u = 0                    indomain omega.left         "BC";
境界条件

10行目
  u = extrapolateField(u)  indomain omega.right        "extrapolation";
境界条件を定義していない方の境界面を定義するためのおまじないと覚えておけばよいと思う
一応、User GUideの記述を引用する
> All fields that are spatially differentiated must have either BC or extrapolation at each boundary. This extrapolation should be done automatically by the compiler, but this has not been implemented yet. The current workaround is the usage of the extrapolateField() operator directly in the model.

要は本来コンパイラが自動的にやるべきことだが、まだ実装してないから自分で書いてね、ということだと思う

## 計算結果

以下のようにプロットウィンドウを操作するとアニメーションで結果が表示されます

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/61d6923b-b70f-c1c1-0657-2022e29575ad.png)

0.25sec
振動しているような
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/7cb7ccf3-5849-bdb9-bfe9-99000ff3847a.png)

0.50sec
やっぱり振動している
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/38c35da0-94ff-9207-844c-8f0d9c868421.png)

1.0sec
振動が消えず残っている
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/33e8a13a-72b7-460b-78a4-7db3bfbced98.png)

### 矩形波で計算

勉強会で矩形波で計算すると精度の評価を適切に出来るとご助言頂いた。
そこで以下のような初期値で計算を実施した。
綺麗な矩形にならなかったが。。。

初期状態
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/e26c3b5b-84ab-5413-e33e-f567d0121d04.png)

0.05sec　ああ・・・
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/4876fa00-7df3-0da2-fbe0-6cfcac15d3d9.png)

0.25sec

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/857f4d62-756d-fdca-a86a-904c833d089f.png)

0.50sec

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/b9bfb511-f2c4-96fb-5b76-72302fa77c82.png)

分割数N=500にして再計算

初期状態 今回は矩形波のようになっている
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/20907a98-7ec0-836e-1dd5-a489ce8646af.png)


0.05sec
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/e3575003-bb85-3bdb-6315-ff52005eaffd.png)

0.50sec
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244254/e9c779e4-e988-79e5-6b9e-704cdd8059ff.png)

精度が悪いようだ

# 所感
## 矩形波を試していない時

少し精度が気になるがまぁ解けたという感じ

## 矩形波を試した後

これは精度や計算的な不安定性からあまり使わない方がよさそう・・・

## 感想
元論文では中央差分、前進差分、後退差分程度が実装されていたと思うのでそこらへんが選べるようになるといいかもしれない
ただ、正直上記の機能程度なら自分で書いてしまえるのでなんとも。
(一次元非定常熱伝導方程式を実装した時の資料は以下
https://www.slideshare.net/ssusere33bfb/170624-room-andmbdandislideshare)

ただし、Modelica界隈が常微分だけでなく偏微分も取り組んでいることが分かる熱いメッセージだと思えば非常にわくわくする機能だと思う