いわて駐在研究日誌

OpenCAE、電子工作、R/C等、徒然なるままに

Dakotaの並列計算機能について(その2)

以下のような簡単な問題(Beale's function)を最適化で解いてみる。設計変数がx1, x2で、目的関数(単一)がf、制約条件(単一)がc>=0である。最小値の答えはf(3, 0.5)=0。

See: Test functions for optimization - Wikipedia, the free encyclopedia

#!/bin/python
from sys import argv
from math import *
import numpy as np

def f(x1, x2):
    return (1.5-x1+x1*x2)**2+(2.25-x1+x1*x2**2)**2+(2.625-x1+x1*x2**3)**2

def c(x1, x2):
    return x1**2-x2

if __name__ == '__main__':
    x1 = float(argv[1])
    x2 = float(argv[2])
    print f(x1,x2)
    print c(x1,x2)

これの最適解は、

 

①Gradientベースの最適化(conmin_mfd method)

勾配を用いる最適化の場合は、逐次的に最適値を探索する。勾配には、あらかじめ数式などで求まっている場合と数値的な微分を用いる場合があるが、後者が一般的。

Dakotaで実装されている一番一般的な勾配法として、conmin_mfdやconmin_frcg ( conjugate gradient optimization method)がある。ここでは前者を利用する。勾配はNumerical gradientである。初期値は(1,0)とする。

1)シリアル計算で実行した場合

<<<<< Function evaluation summary: 10 total (8 new, 2 duplicate)
<<<<< Best parameters          =
                      5.0524767374e-01 x1
                      2.5527521182e-01 x2
<<<<< Best objective function  =
                      8.9519566068e+00
<<<<< Best constraint values   =
                     -1.2678746941e-13
<<<<< Best data captured at function evaluation 6


<<<<< Iterator conmin_mfd completed.
<<<<< Environment execution completed.
DAKOTA execution time in seconds:
  Total CPU        =       0.02 [parent =   0.010999, child =   0.009001]
  Total wall clock =    1.26761

real    0m1.362s
user    0m1.528s
sys    0m0.221s

 2)evaluation_concurrency=4を指定した場合

<<<<< Function evaluation summary: 10 total (8 new, 2 duplicate)
<<<<< Best parameters          =
                      5.0524767374e-01 x1
                      2.5527521182e-01 x2
<<<<< Best objective function  =
                      8.9519566068e+00
<<<<< Best constraint values   =
                     -1.2678746941e-13
<<<<< Best data captured at function evaluation 6


<<<<< Iterator conmin_mfd completed.
<<<<< Environment execution completed.
DAKOTA execution time in seconds:
  Total CPU        =       0.02 [parent =   0.011998, child =   0.008002]
  Total wall clock =    1.26666

real    0m1.367s
user    0m1.534s
sys    0m0.234s

 3)mpirun -np 4で起動した場合(evaluation_concurrency指定なし)

<<<<< Function evaluation summary: 10 total (8 new, 2 duplicate)
<<<<< Best parameters          =
                      5.0524767374e-01 x1
                      2.5527521182e-01 x2
<<<<< Best objective function  =
                      8.9519566068e+00
<<<<< Best constraint values   =
                     -1.2678746941e-13
<<<<< Best data captured at function evaluation 6


<<<<< Iterator conmin_mfd completed.
<<<<< Environment execution completed.
DAKOTA master processor execution time in seconds:
  Total CPU        =       0.02 [parent   =   0.018997, child =   0.001003]
  Total wall clock =    1.42358 [MPI_Init = 1.69277e-05, run   =    1.42357]

 

という結果になった(3つとも当然ながら計算結果は一致し、勾配ベースはやはり計算が速い!)が、正解には到達していない。

 原理的に逐次計算を必要とする(数値的勾配を用いる場合)は、並列計算の意味があまりなさそうである(もしかしたらもっと良い方法があるのかもしれないが....)。

※ 勾配ベースの最適化は、実際には滑らかな応答局面でないと収束しなかったり、

局所最適解に落ちつくなどがあるので、流体計算などで適用できる場合は限られるかもしれない。

 

②Non-Gradientベースの最適化(moga method)

MOGA(Multi Objective GA)で確認してみる。mogaの設定は以下の通りで、最適解の可能性のあるものを1つ出力する。このあたりのパラメーターの設定をどうすればよいか、指針が欲しいところ(現状は、ネットで調べた設定をそのまま利用)。

method
  moga
  output silent
    seed = 10983
    final_solutions = 1
  max_function_evaluations = 5000
  initialization_type unique_random
  crossover_type shuffle_random
    num_offspring = 2 num_parents = 2
    crossover_rate = 0.8
  mutation_type replace_uniform
    mutation_rate = 0.1
  fitness_type domination_count
  replacement_type below_limit = 6
    shrinkage_percentage = 0.9
  convergence_type metric_tracker
    percent_change = 0.05 num_generations = 50

1)シリアル計算で実行した場合

<<<<< Function evaluation summary: 257 total (257 new, 0 duplicate)
<<<<< Best parameters          =
                     -9.7989947697e-01 x1
                      9.9330955813e-01 x2
<<<<< Best objective function  =
                      1.4384761666e+01
<<<<< Best constraint values   =
                     -3.3106573164e-02
<<<<< Best data captured at function evaluation 115


<<<<< Iterator moga completed.
<<<<< Environment execution completed.
DAKOTA execution time in seconds:
  Total CPU        =        0.2 [parent =    0.19597, child =    0.00403]
  Total wall clock =    47.6183

real    0m47.714s
user    0m41.059s
sys    0m5.870s

2)evaluation_concurrency=4を指定した場合

<<<<< Function evaluation summary: 257 total (257 new, 0 duplicate)
<<<<< Best parameters          =
                     -9.7989947697e-01 x1
                      9.9330955813e-01 x2
<<<<< Best objective function  =
                      1.4384761666e+01
<<<<< Best constraint values   =
                     -3.3106573164e-02
<<<<< Best data captured at function evaluation 115


<<<<< Iterator moga completed.
<<<<< Environment execution completed.
DAKOTA execution time in seconds:
  Total CPU        =       0.18 [parent =   0.187972, child =  -0.007972]
  Total wall clock =     22.097

real    0m22.193s
user    0m43.546s
sys    0m6.179s

3)mpirun -np 4で起動した場合(evaluation_concurrency指定なし)

<<<<< Function evaluation summary: 257 total (257 new, 0 duplicate)
<<<<< Best parameters          =
                     -9.7989947697e-01 x1
                      9.9330955813e-01 x2
<<<<< Best objective function  =
                      1.4384761666e+01
<<<<< Best constraint values   =
                     -3.3106573164e-02
<<<<< Best data captured at function evaluation 115


<<<<< Iterator moga completed.
<<<<< Environment execution completed.
DAKOTA master processor execution time in seconds:
  Total CPU        =       0.58 [parent   =   0.588911, child =  -0.008911]
  Total wall clock =    27.0884 [MPI_Init = 9.77516e-06, run   =    27.0884]

real    0m27.302s
user    1m20.837s
sys    0m7.109s

ということで、並列計算では計算時間は約半分になる....がやはり正解にたどり着かない!

 

(追記)

nonlinear constraintをなしにするとどちらも正解に近い値を出すことを確認できた。

Dakotaの並列計算機能について

ようやく復帰。Dakotaの並列計算機能について調べたのでまとめてみる。

概観と簡単なサマリー

Dakota自体から起動する"simulator"はシリアル計算でも並列計算でももちろん構わない(たとえばOFの並列計算など)ので、ここではDakota本体の並列計算について調べてみた。

 

ユーザーマニュアルのChapter 17 Parallel Computingを読んで理解した限りでは、Dakota本体としては、

  1. dakotaプロセス内からの複数のsimulation起動(応答関数の複数同時evaluation)
  2. MPIによる複数のdakotaプロセスの起動(multi node用?)

の機能があり、何をどう並列計算したいかにより、並列度を変えて使い分けることができる。とりあえず一番わかりやすいSingle-Level Parallelでは、

①asynchronous local

単一ノード内で、独立した設計変数による応答関数評価を同時に起動して計算効率化

②message passing with single simulation per node

MPIにより多ノードでdakotaを起動し、それぞれのdakotaプロセスで応答関数評価simlationを1つ起動する。masterノードがそれらの結果を取りまとめする

③Hybrid

①と②の組み合わせ

 となっている。

 

この中で、おそらく個々のsimlation負荷が小さいのであれば、①か③がよく、OF計算のような計算負荷が1つ1つ大きい場合には、②もしくは③がよいのかもしれない。

なお、並列プロセス実行にはforkを推奨とのこと。

 

【DACEでの例】

一番わかりやすいと思われる設計変数に対する応答関数の評価を独立に行うDACEでの並列計算の例を調べてみる。

具体例①:ラテン超方格を使ったDACEチュートリアル(ascynchronous evaluation_concurrency)

問題:Pythonを使った2設計変数2目的関数(objective.py)の最適評価

objective.py

#!/bin/python
from sys import argv
from math import *
import numpy as np

def f1(x1, x2):
    return 2.0*sqrt(x1)

def f2(x1, x2):
    return x1-x1*x2+5.0

if __name__ == '__main__':
    x1 = float(argv[1])
    x2 = float(argv[2])
    print f1(x1,x2), f2(x1,x2)

dakota_lhs_par.in

# Dakota Input File: dakota.in

environment
#  graphics
  tabular_graphics_data
    tabular_graphics_file = 'objective.dat'

method
  dace oa_lhs
     seed = 5  
     samples = 121
model
  single

variables
  continuous_design = 2
    initial_point     2.0      1.5
    lower_bounds      1.0      1.0
    upper_bounds      4.0      2.0
    descriptors       'x1'     "x2"

interface
  fork
    asynchronous evaluation_concurrency = 4
  analysis_driver = 'simulator_script'
  parameters_file = 'params.in'
  results_file = 'results.out'
  work_directory directory_tag
    copy_files = 'templatedir/*'
  named 'workdir' file_save  directory_save
    aprepro

responses
  response_functions = 2
  no_gradients
  no_hessians

赤字の部分のようにevaluation_concurencyを指定すると同時に評価サンプルを複数実行してくれる。この場合、シングルで24.2s, evaluation_concurrency = 4 で7.5sとなった。

$ dakota dakota_lhs_par.in

[snip]

<<<<< Function evaluation summary: 121 total (121 new, 0 duplicate)

Simple Correlation Matrix among all inputs and outputs:
                       x1           x2 response_fn_1 response_fn_2
          x1  1.00000e+00
          x2  1.37441e-02  1.00000e+00
response_fn_1  9.96397e-01  1.52353e-02  1.00000e+00
response_fn_2 -4.95228e-01 -8.29560e-01 -4.94665e-01  1.00000e+00

Partial Correlation Matrix between input and output:
             response_fn_1 response_fn_2
          x1  9.96398e-01 -8.66506e-01
          x2  1.81689e-02 -9.47129e-01

Simple Rank Correlation Matrix among all inputs and outputs:
                       x1           x2 response_fn_1 response_fn_2
          x1  1.00000e+00
          x2  1.41173e-02  1.00000e+00
response_fn_1  1.00000e+00  1.41173e-02  1.00000e+00
response_fn_2 -4.56767e-01 -8.63711e-01 -4.56767e-01  1.00000e+00

Partial Rank Correlation Matrix between input and output:
             response_fn_1 response_fn_2
          x1  1.00000e+00 -8.82201e-01
          x2  1.16427e-10 -9.63760e-01


<<<<< Iterator dace completed.
<<<<< Environment execution completed.
DAKOTA execution time in seconds:
  Total CPU        =       0.09 [parent =   0.095986, child =  -0.005986]
  Total wall clock =    7.45185

※ analysis_concurrencyというキーワードもあるが、これを指定すると同時に評価する最大数を制限できる。

When asynchronous execution is enabled and each evaluation involves multiple analysis drivers, then the default behavior is to launch all drivers simultaneously. The analysis_concurrency keyword can be used to limit the number of concurrently run drivers.

 

具体例②:ラテン超方格を使ったDACEチュートリアル(mpirun on single node)

問題:Pythonを使った2設計変数2目的関数(objective.py)の最適評価(上記と同じ)

dakota_lhs.in (evaluation concurrencyの指定なしに注意)

# Dakota Input File: dakota.in
environment
#  graphics
  tabular_graphics_data
    tabular_graphics_file = 'objective.dat'

method
  dace oa_lhs
     seed = 5  
     samples = 100

model
  single

variables
  continuous_design = 2
    initial_point     2.0      1.5
    lower_bounds      1.0      1.0
    upper_bounds      4.0      2.0
    descriptors       'x1'     "x2"

interface
  fork

###    asynchronous evaluation_concurrency = 4
  analysis_driver = 'simulator_script'
  parameters_file = 'params.in'
  results_file = 'results.out'
  work_directory directory_tag
    copy_files = 'templatedir/*'
  named 'workdir' file_save  directory_save
    aprepro

responses
  response_functions = 2
  no_gradients
  no_hessians

$ mpirun -np 4 dakota dakota_lhs.in  > dakota_lhs_openmpi.out

→ dakotaプロセスが4つ起動する。前例では1つだけである。

[snip]

<<<<< Function evaluation summary: 121 total (121 new, 0 duplicate)

Simple Correlation Matrix among all inputs and outputs:
                       x1           x2 response_fn_1 response_fn_2
          x1  1.00000e+00
          x2  1.37441e-02  1.00000e+00
response_fn_1  9.96397e-01  1.52353e-02  1.00000e+00
response_fn_2 -4.95228e-01 -8.29560e-01 -4.94665e-01  1.00000e+00

Partial Correlation Matrix between input and output:
             response_fn_1 response_fn_2
          x1  9.96398e-01 -8.66506e-01
          x2  1.81689e-02 -9.47129e-01

Simple Rank Correlation Matrix among all inputs and outputs:
                       x1           x2 response_fn_1 response_fn_2
          x1  1.00000e+00
          x2  1.41173e-02  1.00000e+00
response_fn_1  1.00000e+00  1.41173e-02  1.00000e+00
response_fn_2 -4.56767e-01 -8.63711e-01 -4.56767e-01  1.00000e+00

Partial Rank Correlation Matrix between input and output:
             response_fn_1 response_fn_2
          x1  1.00000e+00 -8.82201e-01
          x2  1.16427e-10 -9.63760e-01

 

<<<<< Iterator dace completed.
<<<<< Environment execution completed.
DAKOTA master processor execution time in seconds:
  Total CPU        =       0.21 [parent   =   0.206969, child =   0.003031]
  Total wall clock =    7.38266 [MPI_Init = 9.05991e-06, run   =    7.38266]

※ 計算時間はevaluation_concurrency指定の場合とほぼ変わらず7.4s。

 ということで、マルチノードならmpi、シングルノードならevaluation concurrency指定でよさそうである。

 

→ 次は、gradien base, Non-gradient baseの最適化の例について調査すべし。

最近の状況

子育て中心で、なかなかブログを更新できておりません。

定時を過ぎたらばたばたと迎えにいって、おむつ替えて、遊びの相手して、風呂入れて、ミルク飲ませて、寝かしつける、で8時過ぎ。

飯を食って風呂入ると、疲れて眠くなり、机で気絶してたりする。少なくとも家ではまとまった研究とか読書の時間がとれないのは、要反省ですな。

しかも、夜中に起きる子供の相手すると2時間はすぐに経つし、朝5時には起きるし、慢性的に寝不足な感じです。どうもいかんなー。

パスフレーズなしのSSHログイン(OpenMPI)と、NFSマウントの設定

パスフレーズなしのSSHログイン設定メモ

 

(1)クライアント(ログインユーザはwaku)でssh2のRSA暗号鍵を生成

$ ssh-keygen -t rsa

パスフレーズを聞いてくるが、なしでよければそのままリターンでOK

 

(2)公開鍵(*.pub)のほうを接続したいサーバにコピーする。

楽をしたければ、以下のコマンド1つでOK。

$ ssh-copy-id -i ~/.ssh/id_rsa.pub waku@SEREVER

もしくは以下のように、クライアントからscpなどをつかって公開鍵をコピーし、その後、サーバー側で設定する。

※ サーバ自身も計算ノードに使う場合は、サーバ自身にたいしてもやっておくこと。

※ 同様に、各クライアント→サーバ、各クライアント→他のクライアントもパスフレーズなしで認証できるようにしておく必要がある(重要)。

@クライアント

$ scp ~/.ssh/id_dsa.pub waku@SERVER

@サーバ(コピーした公開鍵をauthoized_keysに追記)

$ cat ~/id_dsa.pub >> authorized_keys server

$ chmod 600 authorized_keys

$ rm ~/id_dsa.pub

 

(3)パスワードなしでsshログインできるか確認

$ ssh waku@ensis13
Last login: Mon Jul 13 08:55:12 2015 from ensis10
[waku@ensis13 ~]$

 

NFS設定メモ

 ※ mpirunには--preload_binaryというオプションがあるので、sshでホスト名を使ったログインができていれば、いちいちNFSマウントしていなくても実行バイナリを配ってくれるらしい。

※ 今回は真面目に、ジョブサーバの特定ディレクトリをNFSで公開し、計算ノードからNFSクライアントとしてマウントすることにする。

 

(1)使用条件

NFSサーバ(ジョブサーバ、マウントされる側)

ensis10: /home/waku/exports

NFSクライアント(マウントする側)

ensis{12/13/14}: /home/waku/exports

 

(2)NFSサーバの設定

/etc/exportsに以下を記述

/home/waku/exports  192.168.0.12(rw,no_root_squash)    192.168.0.13(rw,no_root_squash)        192.168.0.14(rw,no_root_squash)

exportsを確認

exportfs -v

/home/waku/exports
        192.168.0.12(rw,wdelay,no_root_squash,no_subtree_check)
/home/waku/exports
        192.168.0.13(rw,wdelay,no_root_squash,no_subtree_check)
/home/waku/exports
        192.168.0.14(rw,wdelay,no_root_squash,no_subtree_check)

nfs, rpcbindの設定・起動

# yum install nfs* rpcbind

# chkconfig rpcbind on

# chkconfig nfslock on

# chkconfig nfs on

# /etc/init.d/nfs start

 

(3)NFSクライアントの設定

# /etc/rc.d/init.d/rpcbind start

# chkconfig rpcbind on

# /etc/rc.d/init.d/netfs start

# chkconfig netfs on

$ mkdir /home/waku/exports

手動マウントしてみる

# mount -t nfs 192.168.0.10:/home/waku/exports /home/waku/exports

うまくいけば、

# vi /etc/fstab

192.168.0.10:/home/waku/exports /home/waku/exports nfs _netdev,rw,async,hard,intr 0 0

(注)ここでサーバ、クライアントは同じGID/UIDを持つものとしているので(事前に同じGID/UIDに変更している)、同じユーザwakuで読み書きできるはずである。ただし、NFSv4(SL6/CentOS6)のバグ?のようなものがあり、nobody.nobodyでマッピングされてしまうことがある。対処方法は、以下のようにすれば良い。

# vi /etc/idmapd.conf

  Domain = localdomain (記述を追記)

# /etc/rc.d/init.d/rpcimapd restart

#  nfsidmap -c (nfsidのキャッシュクリア)

# /etc/rcd./init.d/netfs restart

nfsv4 mounts files as nobody

 

 

OF+外部ライブラリ(物性値ライブラリ)その3

 水・水蒸気のライブラリであるfreesteamとOFの併用テスト

1)まず、freesteamをビルドして(sconsが必要)、.soライブラリ”libfreesteam.so.1”を作成しておく。OFのビルド環境で実行する。

tar xf freesteam-2.1.tar.bz2

cd freesteam-2.1

scons

scons test && ./test. (テスト)

 2)続いて、IAPWS-IF97-OFのライブラリ”libIAPWSRangeThermo.so”のビルドを行う。OFのビルド環境で実行する。

git clone https://github.com/romansCode/IAPWS-IF97-OF.git

cd IAPWS-IF97-OF

wclean

wmake libso

 

 

VAWT翼型最適化(2)

メッシャーの問題がひと段落したので、いよいよ計算という感じであるが、ここから泥臭い作業が続くのだろうなと覚悟しておく。

 

その前に、Dakotaで計算する場合、特定のCFD計算が発散したりしてエラーになる場合、それを飛ばして次の計算に移る仕組みを調べておく必要がありそう(リスタートとは別物・・・リスタート機能もちゃんと調べないとな)。

 

で、ちょっと調べてみた。

おそらくこのあたりだろうと思うが、以下のキーワードを設定すると良さそう。

Dakota Reference Manual: failure_capture

interface>analysis_drivers>failure_capture

  • abort (the default)
  • retry
  • recover
  • continuation

デフォルトではabortになっているので、各設計変数での計算が失敗した場合は、そこでdakotaも終了ということになる。retryはもう一度analysisを実行。recoverは、一緒に指定した”reals for specifying the dummy function values”をダミーで出力して次のanalysisへ進む。cotinuationは良く分からないが、それまでに成功した設計変数から、失敗した設計変数のanalysisに近い設計変数を生成し再トライするような感じである。

ということで、recoverでダミー目的関数値をしておくことで後からそれらを弾けば良さそうである。

 

続いてリスタートであるが、dakotaは黙っていても、dakota.rstというリスタートファイルを作成してくれる。ファイル名を指定したければ、

dakota -i dakota.in -write_restart my_restart_file

このファイルを読み込む場合は、

dakota -i dakota.in -read_restart my_restart_file

で良いとのこと。-write_restartと -read_restartで同時にファイルを指定すると、-write_restartの指定にappendする。

さらに、

dakota -i dakota.in -r dakota.rst -s 50 -w dakota_new.rst

とすると、dakota.rstの最初の50個のevaluationを読み込み、計算を続行する。

 

cfMeshで機能追加を希望するところ

cfMeshをいろいろ使っていて、やっぱり非常に使いやすい。

ちょっと不具合があってtwitterでぼやいていたら、わざわざ中の人が声を掛けてくれて直してくれたし(正確には、自分のSTLデータの造り方が想定外だったようだけど)、非常に好感触である。

 

さらに使い勝手向上のためには、

  • MultiRegionへの対応
  • cellZoneを指定可(patchのrenameはできるので)

してくれると嬉しいかなーと思いました。まぁ、既存のやり方で十分対応できるんですが...。