後方乱気流

シミュレータとOSSのメモ

SALOMEのShaperモジュール + Python3で攪拌翼をモデリングする

勝手にAdventカレンダー16日目の記事です。

動機

SALOMEは, ほぼすべての機能をSALOME内蔵のPythonインタープリタから実行可能です。
簡易3D-CAD機能を提供するShaperモジュールも, 例外ではありません。
サイズ違いのモデルを何通りか作成したい場合, PythonスクリプトからSALOME/Shaperを実行できれば省力化できると思われます。
今回はPythonからShaperを操作し, STLファイルを生成するまでの流れをまとめました。

やりたいこと

Shaperモジュールで攪拌翼, 形状が簡単なパドル翼を生成します。
パドル翼のスパン径や取付角度を自由に変更できるようPythonスクリプトに落とし込みます。

f:id:junkroom:20201216211832p:plain
Fig. 1 sketch
f:id:junkroom:20201216211841p:plain
Fig. 2 result

環境

DockerでUbuntu環境を構築しています。

  • ベースOS : Windows10 Home Edition 64 bit

  • 仮想環境 : Docker / WSL2

    • (Ubuntu 18.04上にSALOMEと依存パッケージ導入済みコンテナを使用)
  • SALOME : v9.5.0
  • Python : v3.6.5(SALOME同梱), v3.6.9(システムPython)

Pythonスクリプトを書いてみる

作業ディレクトリにwing.pyを作成します。

cd dir/to/working-dir
touch wing.py

Salome実行に必要なライブラリをインポートします。
wing.pyを開いて以下を書き込みます。

import salome
import salome_notebook
from SketchAPI import *
from salome.shaper import model
import GEOM
from salome.geom import geomBuilder

攪拌翼クラス Paddle とShaperモジュールの初期化メソッド群を用意します。

class Paddle:
    def __init__():
        pass

    def set_partset(self):
    """
    パーツセットを生成する。
    """
        self.partSet = model.moduleDocument()

    def create_new_part(self, part_name):
    """
    ShaperモジュールのNew Partと等価
    Args:
      part_name (str): パーツ名
    """
        self.parts[part_name] = model.addPart(self.partSet)
        self.parts[part_name].setName(part_name)
        self.part_docs[part_name] = self.parts[part_name].document()

    def create_new_sketch(self, part_name, i, plane='YOZ'):
    """
    ShaperモジュールのNew Sketchと等価
    Args:
      part_name (str): パーツ名
      i (int): パーツ名の末尾につける数字
      plane (str): スケッチ面
    """
        sketch_name = '{}_{}'.format(part_name, i)
        self.sketch[sketch_name] = model.addSketch(
            self.part_docs[part_name],
            model.defaultPlane(plane))
        self.sketch[sketch_name].setName(sketch_name)

ShaperモジュールはPartSetに各Partを登録し、Part内にSketchを保持する構造になっています。
新規PartSet生成は model.moduleDocument()
Partの登録は model.addPart()
Sketchの生成は model.addSketch() といったメソッドで行えます。
戻り値は、生成したオブジェクトクラスのインスタンスです。
インスタンスが持つメソッドやインスタンス自体を使ってオブジェクトを操作するので、戻り値をPaddleクラスのインスタンス変数に保持しておきましょう。

続いて攪拌翼生成メソッド create_wing() を作ります。

    def create_wing(self, h, t, w, m, angle, n_wing, p_name=None):
    """
    Shaperモジュール内にパドル翼を生成する.

    Args:
      h (float): 攪拌翼の高さ
      t (float): 攪拌翼の厚み
      w (float): 攪拌翼の幅(1枚分). 攪拌槽の半径以下に設定する.
      m (float): 攪拌翼の下端から攪拌槽底部までの距離
      angle (float): 攪拌翼の取付角度
      n_wing (int): パドルの枚数
      p_name (str): パーツ名
    """
        if p_name is None:
            p_name = 'wing'

        s_name = '{}_{}'.format(p_name, 1)
        self.create_new_part(p_name)
        self.create_new_sketch(p_name, i=1, plane='XOZ')
        lines = [0] * 5
        lines[0] = self.sketch[s_name].addLine(-t / 2, m, -t / 2, m + h)
        lines[1] = self.sketch[s_name].addLine(-t / 2, m + h, t / 2, m + h)
        lines[2] = self.sketch[s_name].addLine(t / 2, m + h, t / 2, m)
        lines[3] = self.sketch[s_name].addLine(t / 2, m, -t / 2, m)
        lines[4] = self.sketch[s_name].addLine(-t / 2,
                                               m + h / 2, t / 2, m + h / 2)
        lines[4].setAuxiliary(True)
        point = self.sketch[s_name].addPoint(0, (m + h) / 2)
        point.setAuxiliary(True)
        for i in range(4):
            if i != 3:
                self.sketch[s_name].setCoincident(lines[i].endPoint(),
                                                  lines[i + 1].startPoint())
            else:
                self.sketch[s_name].setCoincident(lines[i].endPoint(),
                                                  lines[0].startPoint())
        self.sketch[s_name].setCoincident(lines[4].startPoint(),
                                          lines[0].result())
        self.sketch[s_name].setCoincident(lines[4].endPoint(),
                                          lines[2].result())
        self.sketch[s_name].setMiddlePoint(lines[0].result(),
                                           lines[4].startPoint())
        self.sketch[s_name].setMiddlePoint(lines[2].result(),
                                           lines[4].endPoint())
        self.sketch[s_name].setCoincident(point.coordinates(),
                                          lines[4].result())
        self.sketch[s_name].setMiddlePoint(lines[4].result(),
                                           point.coordinates())
        self.sketch[s_name].setParallel(lines[1].result(), lines[3].result())
        self.sketch[s_name].setParallel(lines[0].result(), lines[2].result())
        self.sketch[s_name].setPerpendicular(lines[0].result(),
                                             lines[1].result())
        proj_z = self.sketch[s_name].addProjection(
            model.selection('EDGE', 'PartSet/OZ'),
            False)
        oz_axis = proj_z.createdFeature()
        self.sketch[s_name].setCoincident(point.coordinates(),
                                          oz_axis.result())
        self.sketch[s_name].setLength(lines[0], h)
        self.sketch[s_name].setLength(lines[1], t)
        proj_x = self.sketch[s_name].addProjection(
            model.selection('EDGE', 'PartSet/OX'),
            False)
        ox_axis = proj_x.createdFeature()
        self.sketch[s_name].setDistance(
            point.coordinates(), ox_axis.result(), (m + h) / 2)
        self.sketch[s_name].setAngle(
            lines[0].result(), ox_axis.result(), angle)
        model.do()
        extrusion = model.addExtrusion(self.part_docs[p_name],
                                       [model.selection('COMPOUND', s_name)],
                                       model.selection(), w, 0)
        extrusion.setName('wing')
        extrusion.result().setName('wing')
        angular_copy = model.addMultiRotation(self.part_docs[p_name],
                                              [model.selection(
                                                  'SOLID', 'wing')],
                                              model.selection(
                                                  'EDGE', 'PartSet/OZ'),
                                              n_wing,
                                              keepSubResults=True)
        angular_copy.setName('wings')
        angular_copy.result().setName('wings')
        lines[0] = self.sketch[s_name].addLine(-t / 2, m, -t / 2, m + h)

Sketchに線を追加するメソッドがsketch.addLine()です。 引数は、(始点のx座標, 始点のz座標, 終点のx座標, 終点のz座標)です(スケッチ面がXOZ面の場合)。 戻り値は変数やリスト、辞書の値として受けとれます。

        lines[4].setAuxiliary(True)

lines[4]は、wing_angleで翼を傾けるときに回転中心を拘束するための補助線です。 setAuxiliary(True)で補助線にしておきます。

        for i in range(4):
            if i != 3:
                self.sketch[s_name].setCoincident(lines[i].endPoint(),
                                                  lines[i + 1].startPoint())

addLine()で線を追加した直後は、線同士がつながっていません。
setCoincident()メソッドで線の末端同士を拘束します。
拘束したい線分は一筆書きで書いておくとfor文で一気に処理できます。

ここでlinesがいくつかのattributeを持つことが分かります。
startPoint()が始点、endPoint()が終点、result()が線そのものを指します。
ここでは線分しか扱いませんが、ArcやCircleを操作する場合はresults()[1]で円弧の線を取り出せます。

後は線分の角度や長さを拘束し、自由度を下げていきます。
setMiddlePoint() : 線分の中点に別の点を拘束 setParallel() : 2つの線分を平行に拘束 setPerpendicular() : 2つの線分を垂直に拘束 setAngle() : 2つの線分間の角度を拘束 setLength() : 線分の長さを拘束 setDistance() : 点と線分までの長さを拘束

攪拌翼をSTLファイルに出力させる関数 export_paddle() を実装します。

    def export_paddle(
            self,
            wing_h, wing_t, wing_w, wing_m, wing_angle, n_wing,
            file_path=None, suffix=0):
        """
        パドル翼を生成する.

        Args:
            wing_h (float): 攪拌翼の高さ
            wing_t (float): 攪拌翼の厚み
            wing_w (float): 攪拌翼の幅(1枚分). 攪拌槽の半径以下に設定する.
            wing_m (float): 攪拌翼の下端から攪拌槽底部までの距離
            wing_angle (float): 攪拌翼の取付角度
            n_wing (int): パドルの枚数
            file_path
        """
        wing_name = f'ag_wing{suffix:04}.stl'
        if file_path is None:
            file_path = os.path.abspath('./')
        wing_path = os.path.join(file_path, wing_name)

        model.begin()
        self.set_partset()
        self.create_wing(wing_h, wing_t, wing_w, wing_m, wing_angle, n_wing)
        model.exportToXAO(
            self.part_docs['wing'], '/tmp/wings.xao',
            model.selection('COMPOUND', 'wings'), 'XAO')
        model.end()

        geompy = geomBuilder.New()
        (imported, wings, [], [], []) = geompy.ImportXAO(
            '/tmp/wings.xao')
        geompy.ExportSTL(wings, wing_path, True, 0.001, True)
model.begin()

は、おまじないです。
GUIでShaperモジュールを選択するのと同じ効果だと思います。
これまでに作成した関数群を実行し、パーツセットを作成、パーツを新規作成、スケッチング、押し出し、回転コピーまで進めます。

model.exportToXAO()

この関数はパーツをXAOファイルに出力します。
次段がXAOファイルが読み書きできるソフトウェアならば、ここで終了です。
私はSTLファイルとして出力したいので、GeometryモジュールにXAOファイルを読み込ませて、STLファイルに変換します。

geomBuilder.New()

でGeometryモジュールのハンドラを取得して、

geompy.ImportXAO()

メソッドを通してShaperで作成したXAOファイルをGeometryモジュール上に取り込みます。 GUI上でExport To Geomボタンを押すのと同じです。

geompy.ExportSTL()

を実行すれば読み込んだモデルをSTLファイルとしてエクスポートできます。 第3引数以降はよく分かりません。

お疲れさまでした。

(おまけ1)単独実行できるようにしておく。

pyファイルの末尾に

if __name__ == 'main':
    paddle = Paddle()
    paddle.export_paddle(hogehoge)

を追加しておきます。
SALOMEをヘッドレスモードで起動し、このファイルを渡せば翼が出力されると思います。

(おまけ2)SALOMEをヘッドレスモードで使う

シェルを開いてSALOMEインストールディレクトリに移動しておく。

./salome -t
./salome shell wing_build.py
./salome killall

-tオプションでSALOMEサーバーをヘッドレスモードで起動します。
2行目でSALOMEのPythonインタプリタへpyファイルを渡しています。
pyファイルに実行時引数を渡したいときは ''' bash ./salome shell wing_build.py args:hoge ''' のように args:をつけてpyファイルの後に置きます。
翼の形状パラメータを流し込むときに使えます。
スクリプトの実行が終わったらkillallコマンドでサーバーをシャットダウンします。

軽量コンテナsystemd-nspawnによるOpenFOAM-v1912導入

目的と趣旨

OpenFOAM-v1912をXubuntu20.04上に導入します。
今回はXubuntuに採用されているシステム管理デーモンsystemd配下の、systemd-nspawnを使用してUbuntu18.04コンテナを作成し、OpenFOAMのインストールを試みます。
ParaViewは ホスト側(Xubuntu20.04)のXサーバに接続して描画できないか挑戦してみました。

なぜsystemd-nspawn?

  • 動作が軽い。
  • ホスト/コンテナ両方でsystemd-networkdが起動していると、自動でホスト-コンテナ通信を確立してくれる。
  • Xubuntuにはsystemdが標準でインストールされているので、セットアップが楽。
  • OpenFOAMのコンパイルにしか使わないライブラリとか諸々をホスト側から隔離しておきたい。
  • ミスってもコンテナごと消し飛ばせばOK!
  • 将来的にはAnsibleなどのサーバー自動設定ツールでゴニョゴニョできるかも...?

なぜDockerを使わない?

  • 筆者のスキル不足が主要因です。
  • Linux上で使い捨てコンテナが作れれば何でも良かったのです。
  • ParaViewなどopenGL周りでトラブルが多発し、嫌気が差したのも理由の一つです。

今回の環境

導入したいもの

  • OpenFOAM-v1912
  • ParaView v5.6.3 (OpenFOAM-v1912 に付属)

注意

  • 筆者のための作業履歴 兼 備忘録です。手順は最適化されていません。
  • root権限で作業する場面が多々あります。コマンドの動作可否、結果について筆者は何ら保証しません。
  • 本記事通りに作業しても、OpenFOAM等が動作する保証もありません。
  • 特にParaview周りは鬼門だと思います。

作業

コンテナ作成

まずはsystemd-nspawnのためのコンテナを作成します。
UbuntuDebian系なのでdebootstrapで環境構築が可能です。

### systemd-nspawnに必要なソフトウェアをインストール
sudo apt install systemd-container
sudo apt install debootstrap

### コンテナのルートディレクトリを作成
### Ubuntu18.04を導入する。
mkdir container && cd container
mkdir Ubuntu1804
sudo debootstrap --arch=amd64 bionic Ubuntu1804 http://archive.ubuntu.com/ubuntu

### コンテナを起動してrootパスワードを設定し、ユーザー(user_nameは任意の文字列)を作成する。
### user_nameをsudoerに任命し終了する。
sudo systemd-nspawn -D ~/container/Ubuntu1804
passwd
adduser user_name
gpasswd -a user_name sudo
### ここまでできたら、ctrl + ] x3回でコンテナにkillシグナルを送り、終了する。

X11の整備

コンテナを起動して x window system環境を整えます。

### bオプションを付けてコンテナを起動する。
### bind-roオプションでコンテナ内のX11アプリケーションをホストの画面に出力させる。(ソケット通信)
### 先程作成したユーザー名とパスワードでログインする。
sudo systemd-nspawn -b -D ~/container/Ubuntu1804 --bind-ro /tmp/.X11-unix:/tmp/.X11-unix

### ダウンローダとテキストエディタをインストール
sudo apt install wget vim

### X11アプリケーションをホストの画面に出力させる。
### DISPLAY変数が設定されていないので、.profileに追記する。
vim ~/.profile
'''~/.profile
(前略)
export DISPLAY=:0.0
'''

### x window systemを整備する。
sudo apt install x11-utils x11-apps

### 設定がうまく行っているかxeyesでテストしてみる。
### 表示されないときはホスト側に戻り、xauth, xhostコマンドでlocalhostからの接続を許可する。
xeyes

目玉が表示されればグラフィック周りはひとまず大丈夫です。

OpenFOAM-v1912の導入の前準備

### リポジトリが"http://archive.ubuntu.com/ubuntu bionic main"しか登録されていない。
### /etc/apt/sources.list を開いて適宜追加する(要root権限)。今回は以下の様に追加した。
sudoedit /etc/apt/sources.list
''' sources.list
deb http://archive.ubuntu.com/ubuntu bionic main universe
deb http://archive.ubuntu.com/ubuntu bionic-security main universe
deb http://archive.ubuntu.com/ubuntu bionic-updates main universe
deb-src http://archive.ubuntu.com/ubuntu bionic main universe
deb-src http://archive.ubuntu.com/ubuntu bionic-security main universe
deb-src http://archive.ubuntu.com/ubuntu bionic-updates main universe
'''
### OpenFOAMソースファイルとThird Partyライブラリ群のダウンロード&解凍する。
mkdir ~/OpenFOAM && cd ~/OpenFOAM
wget https://sourceforge.net/projects/openfoam/files/v1912/OpenFOAM-v1912.tgz
wget https://sourceforge.net/projects/openfoam/files/v1912/ThirdParty-v1912.tgz
tar -xzf OpenFOAM-v1912.tgz && tar -xzf ThirdParty-v1912.tgz

### コンパイルに必要なライブラリやツールをインストールする。
### (参考) https://openfoamwiki.net/index.php/Installation/Linux/OpenFOAM-v1806/Ubuntu
sudo apt install build-essential flex bison cmake zlib1g-dev libboost-system-dev libboost-thread-dev
sudo apt install libopenmpi-dev openmpi-bin gnuplot libreadline-dev libncurses-dev libxt-dev
sudo apt install qt5-default libqt5x11extras5-dev libqt5help5 qtdeclarative5-dev qttools5-dev
sudo apt install libqtwebkit-dev freeglut3-dev libqt5opengl5-dev texinfo
sudo apt install libscotch-dev libcgal-dev python python-dev

apt install でxxxxx is not foundエラーが出るときは /etc/apt/sources.listを再確認してください。
"http://jp.archive.ubuntu.com/ubuntu/"の方がDL速いかもしれません。

OpenFOAMをコンパイル

ようやくコンパイルです。コア数の指定を忘れずに...(1敗)

### OpenFOAMをコンパイル
source ~/OpenFOAM/OpenFOAM-v1912/etc/bashrc
foam
./Allwmake > log.Allwmake 2>&1 &

### マシンのCPUコア数とRAMに余力がある方は、jオプションによる並列処理も可能です。
### ./Allwmake -j 8 > log.Allwmake 2>&1 &
tail +F log.Allwmake

で進行状況を確認できます。
1スレッドコンパイルで最大メモリ使用量は約2.5 GBでした。

ParaViewのコンパイル準備

気持ちが先走ってOpenFOAMを既にコンパイルしてしまいました。
本来はParaViewのコンパイルが先のようです。

cd ~/OpenFOAM/ThirdParty-v1912

### ubuntu18.04の/bin/shはdashのシンボリックリンクになっており、
### そのままではmakeParaViewがうまく動作しない。
### 修正方法は以下の通り(shebangの変更か/bin/shのリンク修正のどちらかで対応する。)
### makeParaView1行目を#/bin/bashに変更する。
vim makeParaView
'''$WM_THIRD_PARTY_DIR/makeParaView
#!/bin/bash
(後略)
'''
### /bin/shのシンボリックリンクをdashからbashに変更
### 選択肢はNoを選ぶ。
sudo dpkg-reconfigure dash

ようやくParaViewコンパイル

### 勝手にマルチスレッド処理してくれるようです。
./makeParaView > log.makePV 2>&1 &
less +F log.makePV

ParaViewプラグインコンパイル

### ParaviewをOpenFOAMより後にコンパイルしたせいかPVReader系がうまく動かない。
### エラー文に従って対処する。
cd $WM_PROJECT_DIR/applications/utilities/postProcessing/graphics/PVReaders
./Allwclean
./Allwmake
wmRefresh

On-Screen Mesa化

openGL周りのエラーに嫌気が差していたので、この際mesaに頼ることにしました。
結果的にParaViewが動くようになりましたが、他のホスト環境(AMDNVIDIA GPU)でどうなるかは不明です。
(あとはレンダリングの負荷がどうなるか...)

### (参考) https://openfoamwiki.net/index.php/Installation/On-Screen_Mesa/Ubuntu#Ubuntu_18.04
sudo apt install mesa-utils libglu1-mesa-dev scons llvm-dev python-pip libyaml-dev
pip install prettytable Mako pyaml dateutils --upgrade
sudo apt build-dep mesa

cd $WM_THIRD_PARTY_DIR
sudo apt source mesa
### mesaソースディレクトリのオーナーがrootになっていたので、
### ログイン中のユーザーに所有権を移した。(user_name, group_nameは自分のコンテナ環境に合わせて変更する。)
cd mesa-19.2.8 && sudo chown -hR user_name:group_name .
scons build=release texture_float=yes libgl-xlib > log.makeMesa 2>&1
cp -vr build/linux-x86_64/gallium/targets/libgl-xlib/* $ParaView_DIR/lib/

ここはノートラブルでした。
ダメな場合はlog.makeMesaの末尾を読んで対処してください。

動作チェック

適当にチュートリアルを実行して、完走するか確かめてみてください。

mkdir ~/Documents/openFoam/ && cd ~/Documents/openFoam/
tut
cp -r ./incompressible/simpleFoam/pitzDaily ~/Documents/openFoam/
cd - 
cd pitzDaily
blockMesh > log.blockMesh 2>&1
simpleFoam > log.simpleFoam 2>&1
paraFoam

筆者の環境では、OpenFOAM/ParaViewどちらも問題なく動きました。
追々、decomposeParやsnappyHexMeshの動作チェックもしていきたいです。

参考

systemd-nspawn - ArchWiki

OpenCFD Release OpenFOAM® v1912

Installation/On-Screen Mesa/Ubuntu - OpenFOAMWiki

Psi4の導入

Psi4とは

psicode.org

Psi4は第一原理計算を実行するためのオープンソースソフトウェアです。 フリー使用可能な計算ソフトとして、GAMESS-USやFireflyが有名ですが、INPUTファイルの複雑さやOUTPUTファイルの互換性の前に心を折られることもしばしばあります。 一方、Psi4はPython向けAPIが充実しており、Jupyter Labなどを通してインタラクティブな操作が可能です。 機械学習やケモインフォマティクスに興味を持っている方々の飛び道具として、お手頃なソフトだと思います。

環境

OS : Ubuntu 18.04 LTS
Python環境 : Ubuntuシステム標準のもの
プリプロセスソフト : Avogadro 1.2.0
ポストプロセスソフト : wxMacMolPlt 7.7

AvogadroやwxMacMolPltは計算対象分子のモデリングや計算結果の可視化に使用します。 Avogadroはaptコマンドなどでインストール可能です。 wxMacMolPltは公式サイトの手順で同様にインストールできました。
https://brettbode.github.io/wxmacmolplt/downloads.html

Psi4のインストール方法

Psi4公式の方法に従ってインストールするのが良いでしょう。

admiring-tesla-08529a.netlify.com

ご自身の環境に合わせてボタンをクリックしていくと、ページ下部に、端末で実行すべきコマンドが表示されます。

シェルを再起動した後、psi4 --testが完走すれば、Psi4導入完了です。お疲れ様でした。