1. はじめに

　今回は、最近注目? のJetson Nanoを用いて、非線形性をもつ偏微分方程式として知られているGray-Scott Modelを、2次元空間において時間発展させる数値計算（数値シミュレーション）を行います。　Jetson Nano自体は15,000円くらいもあればACアダプタ等一式を秋月電子などですぐに入手できます。購入しやすい価格帯なのにもかかわらず、GPUが搭載されているので、ちょこっと計算して、その数値計算結果を少ない資源で手軽にいい感じに描画したい場合などの開発勉強用途に、もしかしたら良いのかもしれません。

2. 描画環境の構築（openFrameworksの導入）

・Jetson NanoのUbuntu環境のセットアップなどは、調べれば色々見つかります。

・ここでは、C++ベースの豊富な描画環境が比較的カンタンに扱えるといった理由で、Ubuntuセットアップ済みのJetson NanoにopenFrameworks(oF)を導入してみます。oFのJetson Nanoへの導入は、こちらの記事を参照してください。英語でずらっと書かれているので、わかり難い部分もあるかもしれませんが、基本的にこの通りでoFの導入はばっちりできます。

・Jetson Nano上でoFの導入が完了後、動作確認のため、oFに元々入っているサンプル(例として、3DPrimitivesExamples)を実行するとこんな感じで描画されます。なお、下記コマンドで make & プログラム実行です。

~/of_v20xxxxx/examples/3d/3DPrimitivesExample$ make && make run

【補足】最初、openFrameworks内のインクルードファイル?の多くをコンパイル&ビルドするのでちょっと時間がかかります。なので、~/of_v20xxxxx/examples/3d/3DPrimitivesExample$ make && make run -j3のようにして実行すると処理が速くなります。

・なお、UbuntuにおけるC/C++の開発環境の作り方は、この記事が分かりやすいです。（※今回のoFを導入する前述の手順を踏んだ時点で必要な諸々がインストール済みとなり、もしかしたら、$sudo apt-get install build-essentialのコマンドは要らないかもしれませんが、とりあえずひと通り実行しました）

・あとは、vim, nanoエディタ等々でコードを直接編集したり書いたりすればOK、なはずです。

3. Turing不安定性（拡散不安定性）について

　「Turing不安定性」「拡散不安定性」「反応拡散方程式」などのキーワードで検索するといろいろ面白そうな文献が出てきます。一応ここでも、北海道大学の文献1を引用して、さわりを説明をしておこうと思います。　ところで、Turing不安定性（拡散不安定性）とは...

　生物は発生の過程で、1 個の卵（細胞）が分裂を繰り返して分化し、性質の異なる細胞集団群へと成長する。Turing はこのような空間一様状態から安定な空間非一様状態へと移行する現象に対して、1952 年に拡散不安定性という概念を提唱した。彼は物質の時間変化を反応 (局所的な物質の生成・消滅) と拡散(近隣領域との平均化) の和で記述し、拡散項がなければ定数解を持つ系が、拡散が加わることによって空間一様解が不安定化し、安定な空間非一様解が出現する場合があることを示した。この物質 (状態) の時間変化を反応と拡散の和で記述する枠組みは、反応拡散モデルと呼ばれる。

　自然界には多種多様なパターン（模様）が存在する。身近なものには、貝殻や動物の模様などがあるが、これらの特徴のひとつは自己組織的にパターンが形成されることである。後天的に動物の体表パターンの一部を欠損させた場合にも、周りの状況に応じてパターンが修復されることが実験的にも確認されている。近年、これらのパターン形成が上述の反応拡散モデルによってよく再現されることがわかってきた。中でも Turing モデルと呼ばれるモデルでは、パラメータによって斑点や縞柄などのパターンが安定な状態として出現することが知られている。

　次は、文献2からこの Gray-Scott Modelを具体的に以下に引用します。なお、ここで反応拡散現象がz方向に一様であると仮定すると、ラプラシアン△はx, yに関する２次元の２階微分となり、

以上が今回扱う Gray-Scott Model です。非線形項を含んだ方程式で記述されている事がわかります。写真の魚の模様のようなパターンが徐々に自然と生成されていく様子を、この後、コンピュータ(Jetson Nano)上でシミュレートします。その前に、コンピュータで計算を行えるようにするために、連続であるこのモデル式を「離散化」します。

4. Gray-Scott Modelの離散化

・「Gray Scott Model」でググると上位に出てきたこちらのサイトを参考にすると良いです。離散化した２階微分の項が漸化式として書いてあります。（各種パラメータを変化させたときの振る舞いをWeb上で確認できるサイトでもあり、とても助かります）

・差分式の誘導の仕方は書いてないですが、他にもグーグルで出てくるので、離散化手法の詳細などは調べてみて下さい。

・最終的にプログラムとして記述する際の数式が書いてあるので、これを参考にし、これを今回はoFのコードに落としこみます。

5. oFにおける数値計算＆描画コード【全文】

・oFでは、主に３つのファイルが必要となります。①main.cpp　②ofApp.h　③ofApp.cpp　の以上３つです。

・~/of_v20xxxxx/apps/myapps/emptyExample のディレクトリに空のプロジェクトフォルダ(つまりemptyExample)があるはずです。それをコピーしてどこか別の場所にバックアップとして保存しておきます。

・基本的に、oFをエラーなく動かすには、プロジェクトの階層構造の決まり？を守る必要があるらしく、その階層構造を崩さないようにするため、このemptyExampleの中のsrcフォルダ内にある①main.cpp　②ofApp.h　③ofApp.cpp　の３ファイルを丸ごと入れ替える、あるいは直接編集して開発するのが定番のようです。（プロジェクトファイル一式を自動生成してくれるProjectGeneratorというツールもあるようですがここでは割愛します）

・よくは分かってないのですが、oFのofVboというクラスメソッドを使うと、簡単にグラフィクスカード(≒GPU?)で描画できるようです。詳しくは公式リファレンスを参照ください。※またこちらのスライド資料も見ながら実装しました。

・以下に、３つのファイルのコード全文を記載します。空間はx, y方向ともに等間隔で刻むようにしました。今回の設定はWIDTH=200, HEIGHT=200とし、200x200の正方形領域を計算します。一応、これでひと通り動きました。

main.cpp

#include "ofMain.h"
#include "ofApp.h"
int main(){
    ofSetupOpenGL(300,300,OF_WINDOW);
    ofRunApp(new ofApp());
}

ofApp.h

#pragma once
#include "ofMain.h"
class ofApp : public ofBaseApp {
public:
    void setup();
    void update();
    void draw();

    void keyPressed(int key);
    void keyReleased(int key);
    void mouseMoved(int x, int y);
    void mouseDragged(int x, int y, int button);
    void mousePressed(int x, int y, int button);
    void mouseReleased(int x, int y, int button);
    void mouseEntered(int x, int y);
    void mouseExited(int x, int y);
    void windowResized(int w, int h);
    void dragEvent(ofDragInfo dragInfo);
    void gotMessage(ofMessage msg);
    void Init();

    int msec;

    // WIDTH x HEIGHT の正方形領域の指定
    static const int WIDTH = 200;
    static const int HEIGHT = 200;
    static const int NUM_PARTICLES = WIDTH * HEIGHT;

    // インスタンス生成
    ofEasyCam cam;
    ofVbo myBuff;

    // ３次元ベクトルで頂点情報を格納する
    ofVec3f myVerts[NUM_PARTICLES];

    // 頂点の色情報を格納する
    ofFloatColor myColor[NUM_PARTICLES];

    float de =  0.10/2; //空間ステップ 
    float dt = 0.050*2; //時間ステップ

    //各種パラメータ
    float Du = 0.00050; 
    float Dv = 0.0040;
    float F = 0.09;
    float k = 0.06;

    float Lapl_u;
    float Lapl_v;
    float Diff_u;
    float Diff_v;

    float t = 0;

    float u[WIDTH + 1][HEIGHT + 1];
    float v[WIDTH + 1][HEIGHT + 1];

};

ofApp.cpp

#include "ofApp.h"
//--------------------------------------------------------------
void ofApp::setup() {
    ofSetFrameRate(50);
    cam.setDistance(180);

    // 頂点情報を初期化
    for (int i = 0; i < WIDTH; i++) {
        for (int j = 0; j < HEIGHT; j++) {
            myVerts[i * WIDTH + j].set(i - WIDTH / 2, j - HEIGHT / 2, 0);
            myColor[i * WIDTH + j].set(0.0, ofRandom(0.0, 1.0), 0.0, 1.0);
        }
    }


    //初期値揺らぎの設定
    for (int j = 0; j < HEIGHT + 1; j++) {
        for (int i = 0; i < WIDTH + 1; i++) {
            u[i][j] = ofRandom(0.0, 1.0);
            v[i][j] = ofRandom(0.0, 1.0);
        }
    }

    //境界条件（ノイマン条件）の設定
    for (int j = 0; j < HEIGHT + 1; j++) {
        u[0][j] = u[1][j];
        u[WIDTH][j] = u[WIDTH - 1][j];
        v[0][j] = v[1][j];
        v[WIDTH][j] = v[WIDTH - 1][j];
    }

    //バッファに位置と色の情報を設定
    myBuff.setVertexData(myVerts, NUM_PARTICLES, GL_DYNAMIC_DRAW);
    myBuff.setColorData(myColor, NUM_PARTICLES, GL_DYNAMIC_DRAW);

}


//--------------------------------------------------------------
void ofApp::update() {
    
for (int i = 1; i < WIDTH; i++) {
    for (int j = 1; j < HEIGHT; j++) {
        //拡散項
        Lapl_u = ( (u[i - 1][j] - 2 * u[i][j] + u[i + 1][j]) + (u[i][j - 1] - 2 * u[i][j] + u[i][j + 1]) ) / (de*de);
        Lapl_v = ( (v[i - 1][j] - 2 * v[i][j] + v[i + 1][j]) + (v[i][j - 1] - 2 * v[i][j] + v[i][j + 1]) ) / (de*de);
        Diff_u = Du * Lapl_u;
        Diff_v = Dv * Lapl_v;
        float fu = std::pow(u[i][j], 2)*v[i][j] - (F + k) * u[i][j];
        float fv = -std::pow(u[i][j], 2)*v[i][j] + F * (1.0 - v[i][j]);
    
        u[i][j] = u[i][j] + dt * ((fu)+Diff_u);
        v[i][j] = v[i][j] + dt * ((fv)+Diff_v);

        float val = ofMap(u[i][j], 0.0, 1.0, 100 / 255, 1.0);
        myColor[i * WIDTH + j].set( val, 0.5647, 0.5647, 1.0);
    }
}
//全点における色情報の更新
myBuff.updateColorData(myColor,NUM_PARTICLES);
t = t + dt;

}

//--------------------------------------------------------------
void ofApp::draw() {
    // カメラ開始
    cam.begin();
    // 頂点の位置をドットで表示
    glPointSize(3.0);
    myBuff.draw(GL_POINTS, 0, NUM_PARTICLES);
    // カメラ停止
    cam.end();
    // FPS,Time の表示
    string info;

    info += " FPS = " + ofToString(ofGetFrameRate(), 2) + "\n" + "Time = " + ofToString(t, 2) + "\n" ;
    ofDrawBitmapString(info, 30, 30);
    ofDrawBitmapString(info, 30, 30);

}

void ofApp::Init() {
    //初期値揺らぎの設定
    for (int j = 0; j < HEIGHT + 1; j++) {
        for (int i = 0; i < WIDTH + 1; i++) {
            u[i][j] = ofRandom(0.0, 1.0);
            v[i][j] = ofRandom(0.0, 1.0);
        }
    }

    //境界条件（ノイマン条件）の設定
    for (int j = 0; j < HEIGHT + 1; j++) {
        u[0][j] = u[1][j];
        u[WIDTH][j] = u[WIDTH - 1][j];
        v[0][j] = v[1][j];
        v[WIDTH][j] = v[WIDTH - 1][j];
    }
    t = 0;
    
}

//--------------------------------------------------------------
void ofApp::keyPressed(int key) {
    if (key == 'r') {  
          Init();
    }
}

//--------------------------------------------------------------
void ofApp::keyReleased(int key) {

}

//--------------------------------------------------------------
void ofApp::mouseMoved(int x, int y) {

}

//--------------------------------------------------------------
void ofApp::mouseDragged(int x, int y, int button) {

}

//--------------------------------------------------------------
void ofApp::mousePressed(int x, int y, int button) {

}

//--------------------------------------------------------------
void ofApp::mouseReleased(int x, int y, int button) {

}

//--------------------------------------------------------------
void ofApp::mouseEntered(int x, int y) {

}

//--------------------------------------------------------------
void ofApp::mouseExited(int x, int y) {

}

//--------------------------------------------------------------
void ofApp::windowResized(int w, int h) {

}

//--------------------------------------------------------------
void ofApp::gotMessage(ofMessage msg) {

}

//--------------------------------------------------------------
void ofApp::dragEvent(ofDragInfo dragInfo) {

}

6. 数値実験結果

・ランダムな微小摂動を初期条件（初期値）として与えて、その初期状態から計算を開始すると、系が時間発展するにつれて、徐々に縞模様や斑点状のパターンが形成されてゆくことがわかります。

・今回は例として、WIDTH x HEIGHT = 200 x 200 の小さめの空間領域において計算しました。このくらいの計算領域であれば、Jetson Nanoの描画能力でも、50近くのFPSが常に出ていることがわかります。

・各種パラメータ(Du, Dv, F, k)やタイムステップ dt を変えてみたり、もっと空間領域を大きくするなどして、数値シミュレーションしながら遊んでみて下さい。

・プログラム実行中、キーボード入力 "r" で初期リセットがかかります。やってみてください。

7. まとめ

　Jetson Nanoに数値計算&描画環境としてC++ベースのopenFrameworksを導入し、その中でGray-Scott Modelを例にとって、Turing不安定性（拡散不安定性）の2次元数値シミュレーションを試行しました。比較的お手軽に、非線形性をもつ偏微分方程式の数値的解析のイントロダクションが行えました。