15 FEniCS で遊ぶ (もがく?)

(2017/7/15 開始)

整理したページは後で作るとして、やったことの記録。

学期末が近づいてめちゃ忙しいのだけれど、 東大数理の集中講義をのぞくことが出来て、 以前から興味は持っていた FEniCS とおつきあいすることに。

劉先生は、集中講義受講者に WWW で学習環境を提供してくれたのだけど、 後々のことを考えると、自分のマシンで直接使えると良いかなと考えて、 Mac で動かせる方法を探す。

色々やり方はあるみたい。 具体的なことは、FEniCS の WWW サイトの Download を読むと良い。

  1. Ubuntu (Linux) で使う。
  2. Docker という Linux コンテナー・ プラットーフォーム (というので合ってるのかな?) で使う。
  3. Anaconda で使う。
  4. ソースからインストールする。
  5. Homebrew で使う (使えると聞いた気がするけれど、本当かな?)。

研究室だったら、1 (Ubuntu) というのもありかな、と思ったのだけど、 自宅では、古めの MacBook Air で作業することが多いので、他の選択肢を検討する。

FEniCS のチュートリアル文書では、何となく 2 (Docker) を推している感じがして、 それを試してみた。でもそうしても、グラフィックスが簡単には出ないみたい。 追加でまた何か (きわめつけは ParaView, そのうち使ってみたいけれど、 今そこまで手を出せないなあ) 入れましょう、みたい。 ………パス!!

3 は、Anaconda という、Python 環境というのか、 Python にもとづくデータ・サイエンス・プラットホームというのか、 そういうのを利用して、FEniCS のもろもろをインストールして使いましょう、 ということらしい。

最初、Anaconda のサイトを見たら、一体何なのだろう?? よそのページを見ないと正体が分かりにくいのって何か変だよなあ。

気を取り直して。

macOS には、もともと Python (2.7系列) 入っているし、 MacPorts でも何故か Python 入ってしまっているし (Mac によっては、2.7 系列と 3.x 系列、両方が入っているぞ)、 この上また別枠で Python 入れるのかと、ちょっとバカバカしい感じがしたけれど、 FEniCS は 3.x で使うべきものらしいから、2.7 系列とダブルのは仕方がないし、 Ubuntu (OS) や Docker (OS よりは軽いけれどコンテナー) を使うことを考えれば、 そんなに大したことがないかと気を取り直す。

たくさん、気を取り直す必要があります。

Anaconda Download for Mac から 3.6 のパッケージを持ってきてインストールする。

ホームディレクトリィの下に anaconda3 というディレクトリィが出来て、 ~/.profile
export PATH=/Users/ユーザー名/anaconda3/bin:$PATH
のような PATH の設定をする。

それから FEniCS のインストール。
ターミナルでこうして FEniCS のインストール
conda create -n fenicsproject -c conda-forge fenics

ターミナルで、こうやって FEniCS を使える状態にする
source activate fenicsproject
(anaconda3/bin/activate というシェル・スクリプトを起動するみたい)

面倒なのでエイリアスでも定義するか。
~/.profile などに別名定義の設定をする
alias fenics='source activate fenicsproject'

(個人的には、普段 bash でなくて tcsh を使っているので、 tcsh 向きの activate 作らないといけないかなあ。)


これで準備OKかと思ったけれど、 matplotlib, mshr の追加インストールが必要のようだ。 注意しないといけないのは、この追加インストールは、 fenicsproject 環境の中で実行する必要がある(らしい)ことである。 「この Mac に matplotlib をインストールしたはずなのにおかしいな?」 と言いながら、インストールしたことが続いて気が付いた。

fenicsproject 環境で、matplotlib と mshr をインストール
conda install -c conda-forge matplotlib
conda install -c conda-forge mshr

以下は、FEniCS のチュートリアル の例題プログラムを実行した記録である。


チュートリアルの例題プログラムを
(fenicsproject) bash-3.2$ python ft01_poisson.py
のように実行して、plot() するところまで来ても絵が出ないので、 https://fenicsproject.org/qa/13050/running-jupyterを参考に、プログラム先頭付近に次の3行を書き込む。
from dolfin import *
parameters['plotting_backend'] == 'matplotlib'
import matplotlib.pyplot as plt
(あれ?最初の2行は必要ないのかな???)

それで plot() でグラフィックスを表示させた直後に、 plt.show() を実行する。

例えば、ft01_poisson.py だと次のような修正を施す。
修正前と修正後を diff で比較
(fenicsproject) bash-3.2$ diff -c  fenics-samples/ft01_poisson.py ../fenics-samples/ft01_poisson.py
*** fenics-samples/ft01_poisson.py      2017-07-10 22:03:24.000000000 +0900
--- ../fenics-samples/ft01_poisson.py   2017-07-15 22:39:49.000000000 +0900
***************
*** 12,17 ****
--- 12,22 ----
  from __future__ import print_function
  from fenics import *

+ # https://fenicsproject.org/qa/13050/running-jupyter
+ from dolfin import *
+ parameters['plotting_backend'] == 'matplotlib'
+ import matplotlib.pyplot as plt
+
  # Create mesh and define function space
  mesh = UnitSquareMesh(8, 8)
  V = FunctionSpace(mesh, 'P', 1)
***************
*** 38,43 ****
--- 43,49 ----
  # Plot solution and mesh
  plot(u)
  plot(mesh)
+ plt.show()

  # Save solution to file in VTK format
  vtkfile = File('poisson/solution.pvd')
(fenicsproject) bash-3.2$

ターミナルでこうして実行
(fenicsproject) bash-3.2$ python ft01_poisson.py

とりあえず、この最初のサンプル・プログラム ft01_poisson.py は動いた。

ft02_poisson_membrane.py はもともと matplotlib が import されていて、無修正で動いた。

ft03_heat.py は、時間発展する系で、 plt.pause(0.01)plt.show() を併用した (参考: http://qiita.com/hausen6/items/b1b54f7325745ae43e47)。

(fenicsproject) bash-3.2$ diff -c ../tmp/fenics-samples/ft03_heat.py ft03_heat.py
*** ../tmp/fenics-samples/ft03_heat.py  2017-07-10 22:04:23.000000000 +0900
--- ft03_heat.py        2017-07-16 03:16:04.000000000 +0900
***************
*** 13,18 ****
--- 13,20 ----
  from __future__ import print_function
  from fenics import *
  import numpy as np
+ # add next line
+ import matplotlib.pyplot as plt

  T = 2.0            # final time
  num_steps = 10     # number of time steps
***************
*** 60,65 ****
--- 62,69 ----

      # Plot solution
      plot(u)
+     #
+     plt.pause(0.01)

      # Compute error at vertices
      u_e = interpolate(u_D, V)
***************
*** 70,73 ****
--- 74,78 ----
      u_n.assign(u)

  # Hold plot
+ plt.show()
  interactive()
(fenicsproject) bash-3.2$

ft04_heat_gaussian.py も同じ調子。 ft05_poisson_nonlinear.py, ft06_elasticity.py, ft07_navier_stokes_channel.py も問題ない。

ft08_navier_stokes_cylinder.py で、 前処理 hypre_amg がないと言われた。
*** Error:   Unable to solve linear system using Krylov iteration.
*** Reason:  Unknown preconditioner "hypre_amg". Use list_krylov_solver_preconditioners() to list available preconditioners().
*** Where:   This error was encountered inside KrylovSolver.cpp.

こうなった事情は、 The solvers and preconditioners change in DOLFIN 2017.1.0 に載っている。

実際どうなっているかは、今動いているものに尋ねるべきか。
(fenicsproject) bash-3.2$ python3
Python 3.6.1 | packaged by conda-forge | (default, May 23 2017, 14:31:56)
[GCC 4.2.1 Compatible Apple LLVM 6.1.0 (clang-602.0.53)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from fenics import *
>>> list_linear_solver_methods()
Solver method  |  Description
----------------------------------------------------------------------------
bicgstab       |  Biconjugate gradient stabilized method
cg             |  Conjugate gradient method
default        |  default linear solver
gmres          |  Generalized minimal residual method
minres         |  Minimal residual method
petsc          |  PETSc built in LU solver
richardson     |  Richardson method
tfqmr          |  Transpose-free quasi-minimal residual method
umfpack        |  UMFPACK (Unsymmetric MultiFrontal sparse LU factorization)
>>> list_krylov_solver_preconditioners()
Preconditioner  |  Description
----------------------------------------------------
default         |  default preconditioner
icc             |  Incomplete Cholesky factorization
ilu             |  Incomplete LU factorization
jacobi          |  Jacobi iteration
none            |  No preconditioner
petsc_amg       |  PETSc algebraic multigrid
sor             |  Successive over-relaxation
>>>

Solver と Preconditioner は、どういう組み合わせがいいんだろう?? とりあえず名前が近そうな petsc_amg で動いているようだけど……… ft08 は時間がかかるので、最後まで実行していない…

ft09_reaction_system.py は 前のプログラムが作ったファイルが必要だ。

ft10_poisson_extended.pyboxfield.py というのがなくなっていて、動かない。 ネットで boxfield.py を探し出して分かったのは、 boxfield.py は Python 2.x 用ということ。 print 文を関数に直しただけでは動かないみたい (まあ、動かない理由が boxfield.py にある、 と決まったわけではないけれど)。 これは困ったね。結局、チュートリアルのプログラムは保守されていないんだな。

(文化が違う感じがする。)

ft11_magnetostatics.py は動く。


一通り動いたので、推奨されている Jupyter Notebook を試してみる。

Jupyter Notebook     Jupyter notebook を実行するには、
(fenicsproject) bash-3.2$ jupyter notebook
あるいは後で説明する理由によって、次のようなオプションをつける。
(fenicsproject) bash-3.2$ jupyter notebook --matplotlib inline

jupyter notebook -matplotlib inline でも動くと書いてあるが、 出来ない。 (http://python.6.x6.nabble.com/iPython-3-0-dev-matplotlib-inline-td5085487.htmlに書いてある話なのかな?)

新しいNotebookを開き、
In[] = %load ft01_poisson.py \fbox{shift} + \fbox{return}
とすると、ft01_poisson.py が Notebook に読み込まれる。 そのセルで、 \fbox{shift} + \fbox{return} を入力すると実行される (「Mathematica みたいに」)。

matplotlib を使っている場合、 %matplotlib inline という (Python の) 注釈をプログラム先頭付近に書いておく必要がある。

これを書くのが面倒という場合は、設定をする。まず次のようにして設定ファイルの準備をする (空のファイルを作る)。
(fenicsproject) bash-3.2$ ipython profile create
これで ~/.ipython/profile_default/ipython_config.py というファイルが出来る。#c.InteractiveShellApp.exec_lines = [] という注釈行があるので、その下にでも、 '%matplotlib inline' というオプションを書き足す。
こんな感じ (3行目が追加したもの)
# lines of code to run at IPython startup.
# c.InteractiveShellApp.exec_lines = []
c.InteractiveShellApp.exec_lines = ['%matplotlib inline']

なるほど、Mathematica のノートブックのようなもの、なんですね。 確かに便利かもしれないので、しばらくは慣れるための時間を投入するのだろう。


(2017/7/22)
別名定義 (個人の趣味ですが)
alias fenics='source activate fenicsproject'
alias fenicsend='source deactivate fenicsproject'
alias note='jupyter notebook'

(2017/7/23) ちょっと気分転換に(現実逃避のため)触ってみる。

桂田 祐史
2018-11-04