計算コードの自己レビュー

最近行っていた拡散方程式の数値計算のコードの見直しがひと段落したので、振り返りをしてみる。
コードは、つまびらかにするわけにいかないので、具体的な記述はあいまいに書いています。

良かったこと

• オブジェクト化がうまくできた。材料をクラスとして、材料のプロパティとしてそれぞれの物理定数 (質量、衝突断面積、温度係数など) を持たせることで、引数渡しが簡単にできた。
• 基本的な物理定数を入力するのをやめて、ライブラリーから取り込んだ。

import scipy.constants as cnst
c = cnst.physical_constants['speed of light in vacuum'][0] ## m/s
e = cnst.physical_constants['elementary charge'][0]
me = cnst.physical_constants['electron mass'][0] 
ips0 = cnst.physical_constants['electric constant'][0] 
kb = cnst.physical_constants['Boltzmann constant'][0] # [J K-1] = [m2 kg s-2 K-1]
na = cnst.physical_constants['Avogadro constant'][0] ## [mol^[-1]]

とする。打ち間違いや、有効桁数のぶれが生じない。

• ファイル名を自動的に拾って、結果ファイル名の共通部分とする。

from sys import argv
from os.path import basename

fn_base = basename(argv[0])[:-3] # ".py" stripped
## (中略)

images = []
if fig_save:    
    images.append("{fn_base:s}_fplot.png".format(**locals()))
    pyl.savefig(images[-1])

• 画像ファイルを含めて、ひとつの tiddly wiki の記事を生成する

addtiddler.py というライブラリを使っている。 (このファイルは、ネット上で見当たらなくなっている)

• 微分について、行列演算での表記を使った。

ベースは http://www.scientificpython.net/1/category/nonuniform/1.html
境界条件も含めて、コンパクトなコードになった。
バージョン 1 では全部場合分けしていたが、場合分け不要のコードにできた。
添え字 i,j,k を走査すると、各座標軸方向での終端にいるか否かはわかる。終端で処理も、行列表記の中で実現する。

• 対称境界条件を導入して、計算規模を小さくしたオブジェクトと、全モデルのオブジェクトがあるが、最初のクラス宣言だけが異なり、そのあとは、同じ振る舞いをするように作れた。
• スパース行列のライブラリを使った。pysparse1.2 scipy の sparse よりもかなり速い