FreeFEM 4.6 でMUMPSで並列計算で求解する [FreeFEM]

FreeFEMのDocumentationのMUMPSに関する記述は (少なくともRelease4.6では) 古いモジュール "MUMPS_FreeFem" に関するものなので、最新のモジュール "MUMPS" の使い方はExample (examples/mpi/MUMPS.edp) を参照すること、という話。

varf lap(u,v) = ... として弱形式を記述して、matrix A = lap(Vh,Vh) 等として全体剛性行列を得る。このあと、load "MUMPS" を実行した上でset (A, solver=sparsesolver, master=xx) とすればソルバーとして直接法のMUMPSを使用するよう設定でき、x[] = A^-1 * b[] 等として求解できる。逐次版でもOpenMPを用いた並列化 (少なくともWindows版の"MUMPS_seq"ではそのようにビルドされているよう) でマルチコアで計算できるし、MPI版ならもちろん並列で求解できる。

ただMPI版で気になるのが、全体剛性行列はどのような値が用いられるのか。結論としては、set()の引数masterの値で動作を切り換えられる。

1) master=0の場合 (centralized)
係数行列はランク0のプロセスの行列Aの値を使う。他のプロセスの値は使用しない。このとき、MUMPSにはICNTL(18)=0 (the input matrix is centralized on the host) が与えられる。また、連立一次方程式の右辺bも同様に、ランク0のプロセスで与えた値を用いる。解xもランク0のプロセスのみが持つ。
つまり、ランク0が逐次版と全く同様に動作するように作成しておけば、求解部分のみを並列計算で行ってくれる。シンプル。

2) master=-1の場合 (distributed)
係数行列は全てのプロセスの (ローカルな) 行列Aの値の総和となる。このとき、MUMPSにはICNTL(18)=3 (user directly provides the distributed matrix, pattern and entries, ...) が与えられる。また、連立一次方程式の右辺bも同様に、全プロセスの値の合計を用いる (FreeFEMのプラグイン側でMPI_Reduceで求めてMUMPSに渡す)。解xも全てのプロセスにブロードキャストされる。

Example中のdistributed版では、change(Th, fregion=nuTriangle%mpisize) でregion番号をプロセス数を上限とする番号に付け直し、int2d(Th, mpirank) (...) としてregion番号が自分のランクと等しい要素のみをローカルな行列Aに計算している。FreeFEM側で各プロセスは自分の持ち分のみ行列要素を計算し、それを足し合わせる処理はMUMPSに任せていることになる。なるほど。

他にもMUMPSのset() のオプションにはstrategy, rinfoやinfoなどがあるようだけど、そちらは未調査です。

※以上、FreeFEM v4.6時点の話です。次にAPIが大変更されるまでは適用可能と思いますが。

参考:plugin/mpi/MUMPS.cpp
MUMPS Documentation

11月23日 17時10分
fartrip




コメント(0) 

11年 [日記]

ふと思い立って、ずいぶん前にニコ動に登録したまま放置していたAbyss in Heavenを観まして。観終わってマイリス登録の日付を見たら、2009年08月11日でした。

何と。

2020年3月27日25時58分
fartrip

コメント(0) 

gmshで有限要素法用のメッシュを作るときのTips [電子工作・プログラミング]

昔は直方体一つ作るにも随分な作業量が必要で必要でしたが、ver.3からOpenCASCADE kernelが使えるようになり、boxコマンド一発で直方体が作れるようになりました。
GUIはSalomeなどの方が親切なようですが、いったん出来上がったモデルの寸法をコマンドラインから変えたりできてそれはそれで便利なので、gmshを少しだけ使ってみました。
そこで、使う上で便利(というかまともなユーザビリティを得るために必要)と思ったことをメモします。

☆いくつかの要素が一体に形成されたメッシュが欲しい
 → Box, Cylinderなどを配置した後、最後にCoherence; もしくはBooleanFragmentsを行う。
   今回使ってみようと思った動機であり、各要素の重なった面・領域で分割をしつつ重複した面を結合し、一体の結合した領域が得られる。

☆境界や領域の番号を見たい
 → Tools->Options->Geometry->Visibility->Surface labelsやCurve labelsのチェックをONにする。Label typeをElementary tagsにする。

☆境界番号の文字が小さくて見えない!
 → Tools->Options->General->Aspect->Default fontの値を大きくする

☆境界番号の文字の色が見えない!
 → Tools->Options->Geometry->Color の中の色を変更する。SurfaceやVolumesを変えたいことが多いと思う。

☆コンソールの文字が小さい!
 → 起動時、gmsh -fontsize 18などとする

☆立体をどの面からみているか分からなくなった
 → Tools->Options->General->Aspect->Projection mode をperspectiveにする

☆メッシュの向きを知りたい
 → メッシュ作成した後、Tools->Options->Mesh->Visibility->Normals and tangentsの値を大きくする

等々。気付いたら追記していきたいと思います。

11月15日21時37分
fartrip

コメント(0) 

続・電池の交換できないスマホ [電子工作・プログラミング]

以前に書いていた電池の交換できないスマホ、実は紆余曲折があり未だ使っています。しかし先日、また具合が悪くなりました。着信の白いLEDランプがずっと点滅しており、スリープしない状態に。LEDが眩しいのはひっくり返しておけば気にならないのですが、電池が1日持たないのは大変困ります。しかも、電源を切っていても電池の消耗は変わらない状態。普通に使っていた状態であれば、ついに故障もしくは電池の寿命かと思いあきらめるところです。

が、このスマホは既に一度分解の憂き目にあった猛者、半分ダメもと、半分確信を持って分解して電池のコネクタを外して再度装着したところ、LEDの点滅が収まり電池の持ちも復活しました。いくらgoogle先生に聞いても、上記の症状で分解せよとは一言も出てきませんでしたが、この方法で無事に直すことができました。やはり少し昔の電子機器、完全にバッテリーを外しての再起動が必須になるシーンもあるようです。いい加減古いのでもうすぐ新しいものに買い換えるつもりですが、バッテリーを交換できない機種だと、三たび同じ目に遭わないかどうか、少し心配もしているところです。

4月3日24時7分
fartrip
コメント(0) 

電池の交換できないスマホ [電子工作・プログラミング]

ある朝、気づいたらスマホの画面が真っ暗。以前からよくあったが、どうやら電源が勝手に落ちたらしい。いつもの通り電源を投入しようとしたら、なぜか起動せず。少しバタバタしていたので、帰宅後に電源をつなぎながらトライするも、やはり起動せず。よく見ると、充電を示す赤いランプが点滅している。どうも電池周りのエラーに見えて、google先生に聞いてみると、「電池を入れ直す」の対策が。このスマホ、電池が交換できないタイプなんです。。。

試してみてだめだったこと
・電源ボタン12秒長押し、その後一度離してさらに12秒長押し
・電源ボタンとマナーモードのボタンを同時に12秒長押し

試してみたら直ったこと
・分解し、電池と基盤をつなぐコネクタを外し、悪戦苦闘の末再度組み直す

原因は不明です。が、サルベージができたことは大変ほっとしました。
もうかなり経つし、次のにしないとなぁ。。。

11月25日26時26分
fartrip
コメント(0) 

Canon iP2700をLinux(CentOS 7 64bit)で使うのは [電子工作・プログラミング]

(8/12追記)
下手なことをせず、USBをつなぐだけの自動認識でできました。調べると、system-config-printerがGutenprintドライバというのを探して入れてくれるようです。最初、純正のドライバを予め用意してからつなごうと考えていたら、ドツボにはまり、結局できなかったという話が↓に書いてあります。いやもう、本当に徒労だったと言う他ありません。

----------(以下、古い内容です)---------------------------------------------------------
できないので、仮想マシンとかを使いましょう、という話。

数年前から我が家にあるCanonのインクジェットプリンターiP2700を、ふとメインマシンのCentOS 7 64bit版から使いたいと思ってしまったのが運の尽き。ずいぶん前に購入してすぐ、当時のfedoraから印刷できたので、難無くできるだろうと思っていたら、苦労の挙句Linuxからはできなかったので、トライしたことと愚痴でも綴ろうかと。

1. rpmのインストール
Canonのwebページから、Linux版rpm(但し32bit)が落とせます。何も考えず解凍し、install.shを叩くと、依存性のエラーがいくつも。libcほか諸々が無いと言ってくるのですが、一見はバージョン含めすでに入っているものばかり。よく考えると、32bit版が無いと言っていることに気付き、片っ端からyumで入れます。それで、最後に残ったのが「libpangox-1.0.so.0」。64bit版ならepelのpangox-compatで見つけたのですが、32bit版は見つけられず。あきらめました。

2. ソースのビルド
同じwebページから、GPLに従ってかソースコードが落とせます。解凍したディレクトリを見ると、何とREADMEもINSTALLも何も無し。Makefileはあるのでmakeしてみると、入った先のディレクトリにMakefileが無いと。実際それぞれの先のディレクトリを見てみると、INSTALLに./autogen.shしてmakeせよと指示が。いやはや。で、それぞれのディレクトリでmakeしていく訳ですが、当然の様にエラーが阻みます。今回やった事は、1. cngpijのbjcups.cのコンパイル時に685行目で「不完全型のポインタへの間接参照」が出たのでオプション「-D_IPP_PRIVATE_STRUCTURES」をつけることと、2. backendのcnij_backend_common.c に #include <cups/ppd.h> を追記する、の2つでした。

これでコンパイルはクリアしたのですが、backendnetのリンク時にlkbcnnet.soがinconpatibleである旨のエラーが出て、ビルド出来ず。で、調べてみると偉大な先人の方がいて、http://usagi.hatenablog.jp/entry/2012/03/06/030324 にて同様のことをされていました。そちらによれば、このビルド済みライブラリはx86用なので、とのこと。これで64bit環境で使うのはあきらめました。


そうなれば、無理してネイティブで動かす必要も無いわけで。かねてから使っていた、windowsの仮想環境でサクっと印刷して終了。ある意味、休日らしい大変贅沢な時間の使い方でした。無駄とも言います。

8月11日12時42分
fartrip

コメント(0) 

Microsoft MPI環境のテスト [電子工作・プログラミング]

WindowsでMPIを使う必要に迫られ、たまたまPCに昔入れていたMicrosoft HPC Pack2012のMPIをようやくテストしてみることに。結論としては(当たり前ですが)LAN内のPC2台でちゃんと動いてほっと一息。ポイントとしては、

1.使用するPCで、MS MPI(特にsmpd.exe)のバージョンを合わせる → 謎の通信エラーが出た
2.使用するPCで、msmpi.dllのバージョンを合わせる → stdout, stderrが開けないという謎のエラーが出た。思えば、msmpi.dllが内部で使用しているmsvcrtのバージョン違いのエラーなんじゃないかと。
3.OSのバージョン自体は違ってもOK → Windows 7 ProとWindows 8.1 Enterpriseで大丈夫だった。

以前に入れたHPC Packのインストーラを失くしていたので、適当にダウンロードしたら1.みたいなことに…。他の準備としては通常言われている「ディレクトリを同じ構造で用意しておくべし」「smpd.exe -p ポート番号」で起動しておくべし、など。mpiexec.exeが使うポート番号はどうやって設定しているのかわかりませんでしたが、どうやら8677を使うようです。

あと、実行時には

mpiexec.exe -host 2 localhost 1 nanigashi 1 a.exe

などとすると、ホストとプロセス数を設定できて便利です。プロセス数を1にして、-envオプションとかでOpenMPが効くかどうかはまだ試していないですが、実はこれが一番重要かもしれません。

(8月2日追記: できました。デフォルトではそのマシンの論理プロセッサ数でOpenMPスレッドが立ち上がり、-envオプションでOMP_NUM_THREADS環境変数を設定すれば、全プロセスでそれが反映されるようです。各プロセス毎にOpenMPスレッドが立ち上がってしまうので、むしろシーンに応じてOMP_NUM_THREADSを1にするように気をつける必要があるかもしれません。)

あと、HPC Pack2012にも関わらず、Visual Studio 2017のx64ネイティブ版のclでビルドしたのが動きました。このあたりの仕組みはよくわかっていませんが、簡単なテストプログラムだからかもしれません。

7月29日26時6分
fartrip

コメント(0) 

CentOS 7でParaView 5.1.2のビルド [電子工作・プログラミング]

できた。CentOS 7のyumで入るcmakeのバージョンが古かったので、それだけは別途ビルドして。Qtのdevel用ライブラリをyumで用意するのが必要なはずだったのだけど、僕の環境ではなぜかすでにQt4が使える環境だったようで、特にその作業は必要なし。

メモとしては、その後色々使いたいと思っていたので、ParaViewのビルド用のccmakeを実行するときにOpenGL2を使えるように「VTK_RENDERING_BACKEND」を「OpenGL2」にすることと(デフォルトは「OpenGL」になっていることに気づかなかった)、プラグインを作る練習もしたいと思っていたので、「PARAVIEW_INSTALL_DEVELOPMENT_FILES」をONにしておくことでしょうか。ファイルを唸るほどダウンロードするので、make時に-jオプションで使用コアを増やしておくと良い感じ。

でも、すんなり行った方だと思います。

6月17日23時5分
fartrip
コメント(0) 

SONY TA-F333ESXの電源ランプ [電子工作・プログラミング]

往年のオーディオアンプ SONY TA-F333ESXをひょんなことから譲り受けました。と言っても、僕自身はアンプ回りのことは全然わからないので、その有難みも今一つ分からないのですが。
で、電源投入してみると、カチッと音がして音はなるようになるのですが、電源ランプは赤いまま。google先生に聞いてみると、下記のページがヒットして、どうやら同様の症状に思います。
http://blogs.yahoo.co.jp/ta3113ko2006/7801685.html
調べてみると、4mm 12Vのムギ球だそうで、そのものは見つけられなかったのですが、とりあえずパイロ電子のT-3WT12Vをゲット。今度、換装できるかどうかやってみます。

でも、今になって
http://blogs.yahoo.co.jp/mrsa0204/246314.html
の記事を見つけて、330Ωの抵抗とLEDでも大丈夫なようですね。まぁ、ムギ球が切れるまで使えるかどうかわかりませんが。

3月12日23時30分
fartrip

コメント(0) 

Intel MKLでHelmholtz方程式とデータ補完 [電子工作・プログラミング]

今や無料で使えてしまうIntelの数値計算ライブラリMKL (Math Kernel Library) のPoisson SolverとData Fitting Funtionsを使ってみました。

前者は簡単な形のLaplace方程式・Poisson方程式・Helmholtz方程式を解けるサービスルーチンで、内部的にはフーリエ変換を使っているよう。確かに波数空間で展開をしてPoisson方程式を解く、というのを大学のころやった気がしますが、理解が追いついていません(が、いまググってみたところ、グリーン関数を使って出した解を計算しているように思えます)。
ソルバーとしては比較的簡単に使えて、結構速いように思えます。ただ、なにぶん問題の制限が強くて。右辺のソース項こそ空間的な分布を持たせられますが、左辺の係数は全領域で一定である必要があります。しかも変数は実数のみ、境界条件はDirichlet条件かNeumann条件のみです。あまりに適用し辛いのですが、メジャーな用途は何になるのでしょうか。座標系として、直交座標の他に球殻上の座標が選べるので、例えば地球上の大気拡散のシミュレーションなどでしょうか。マニアックです。

後者は1次元関数の補完をトライしてみました。Python(のSciPy)だと大変使いやすそうな補完ルーチンがあるのですが、C言語で使えるものは意外と見当たらず、GSLくらいでしょうか。ただ、GSLもSciPyの元ネタも使うのは一見面倒で。ここで、MKLでも用意されていることに気づいたので、使ってみた次第です。やりたいことは確かにできましたが、パラメータの設定が微妙に面倒くさいような…?
でもまぁ、Windows環境下でライブラリのビルドをする手間がないことを思うと、使いやすいのかもしれません。

2月19日25時4分
fartrip
コメント(0)