剛体球の大規模数値計算をしたい人生だった

僕はあまり剛体球の計算の専門家では無いので、何か間違った事書いてたら教えてください。

剛体球シミュレーション

剛体球、英語ではhard sphereとか言われてますが、の計算とは、
「非弾性衝突する多数の粒子をどう計算するか」という事であり、
特に「2つあるいはそれ以上の数の粒子が衝突した際、その後の粒子の位置と速度はどのようになっているか」を計算機に教え込む事です。
こういったテクニックは、ゲームプログラミングや、粉体計算などの分野では頻出です。
ちなみに僕のフィールドでもある天文学とか惑星科学とかでもまぁまぁまぁまぁ使われていると思います。多分。

ゲームプログラムの世界では、複雑な形状の衝突判定をしなくてはならなかったりします。
そういう分野には営利的なアレコレが関わっているので、下手に実装すると特許に引っかかって怒られたりとか。
しかし天文屋さんは物は大体球なので、そんな頭のいい事やらずとも何とかなるわけです。

本稿は、(FDPS, https://github.com/FDPS/FDPS )という粒子計算加速ソフトウェア用いて、剛体球の数値計算を行うための個人用メモです。
個人用メモなんですけど、まぁ公開しとくと何かの訳にはたつかもね、というわけでQiitaに置いておく所存であります。

数式的な話

数学的な話をしましょ。
まず、粒子の座標を$\mathbf{x}$、速度を$\mathbf{v}$、質量を$m$、半径を$r$とします。
右下に小さい$i$をつける事で粒子のラベルにします。例えば、$m_i$なら$i$番目粒子の質量です。
これが他の粒子と接触しているかどうかを判定します。
今考えているのは球ですのでこれは簡単で、2つの粒子の相対距離が半径の和よりも小さかったら接触しているとみなします。
(面倒なんで絵は作って)ないです。
つまり、貫通深度$d_{ij}$と呼ばれる量を定義しまして、

d_{ij} = \max(r_i + r_j - |\mathbf{x}_j - \mathbf{x}_i|, 0)

とおきます。当然$d_{ij} > 0$なら衝突してるし、$0$ならしてません。
球って楽だなぁ。(球以外のものを扱った事があるかのような物言い)
さて、衝突後処理(拘束解消)ですが、色々あるような気もしますがとりあえず今回はRichardson (1994, MNRAS, 269, 493)に書いてある処方箋従います。

速度変化

やることはシンプルで、「衝突した方向に対して正反対の方向に跳ね返る」という挙動を実装してあげるだけです。
衝突がほんの一瞬の現象であり、瞬間的に速度が向きを変えるとするならば…それは撃力を加えていることに相当します。
なので「撃力法(impluse based method)」などと言われています。
まず、衝突面に対して垂直な法線ベクトルと平衡な接線ベクトルを考えましょ。
これは球の場合、当然2粒子の相対位置ベクトルを規格化した物であり、接線ベクトルなら垂直です。
なので、

\mathbf{n}_{ij} = \frac{\mathbf{x}_j - \mathbf{x}_i}{|\mathbf{x}_j - \mathbf{x}_i|}

となります。
次に、相対速度ベクトルの法線・接線方向成分ですが、これは単純に

\begin{eqnarray}
\mathbf{v}_{ij}^\mathrm{n} & = & \left[(\mathbf{v}_{j} - \mathbf{v}_i) \cdot \mathbf{n}_{ij} \right] \mathbf{n}_{ij}\\
\mathbf{v}_{ij}^\mathrm{t} & = & \mathbf{v}_{ij} - \mathbf{v}_{ij}^\mathrm{n}
\end{eqnarray}

となるだけです。
さて、法線方向の反発係数を$\epsilon_\mathrm{n}$、接線方向を$\epsilon_\mathrm{t}$と置き、衝突後の速度を$\tilde{\mathbf{v}}$と起きます。
衝突前後での速度の変化は、

\begin{eqnarray}
\tilde{\mathbf{v}}_{ij} & = & \tilde{\mathbf{v}}_{ij}^\mathrm{n} + \tilde{\mathbf{v}}_{ij}^\mathrm{t}\\
\tilde{\mathbf{v}}_{ij}^\mathrm{n} & = & - \epsilon_\mathrm{n} \mathbf{v}_{ij}^\mathrm{n}\\
\tilde{\mathbf{v}}_{ij}^\mathrm{t} & = & \epsilon_\mathrm{t} \mathbf{v}_{ij}^\mathrm{t}
\end{eqnarray}

となります。法線方向には反発するので符号が逆になるが、接線方向にはそうならない事に注意しましょう。

一方で、2粒子の持つ運動量は衝突前後で保存するとすれば、これは当然

m_i \tilde{\mathbf{v}}_i + m_j \tilde{\mathbf{v}}_j = m_i \mathbf{v}_i + m_j \mathbf{v}_j

となります。
反発後速度の式とこの式を合わせて、頑張って解きますと、

\begin{eqnarray}
\tilde{\mathbf{v}}_i & = & \mathbf{v}_i + (1 + \epsilon_\mathrm{n}) \frac{m_j}{m_i + m_j} \mathbf{v}_{ij}^\mathrm{n}\\
\tilde{\mathbf{v}}_j & = & \mathbf{v}_j - (1 + \epsilon_\mathrm{n}) \frac{m_i}{m_i + m_j} \mathbf{v}_{ij}^\mathrm{n}
\end{eqnarray}

ってなるんじゃないかな多分。
あと、しれっと$\epsilon_\mathrm{t} = 1$、つまり弾性衝突を仮定しました。
これは$\epsilon_\mathrm{t}$を考えるとなると、剛体球回転の効果を入れなければならないのでちょっと面倒だからです。
まぁあとで式だけ書きます。

位置変化

さて、上のセクションでは、衝突時に与えられる速度変化を見積もりましたが、
よくよく考えてみると粒子が衝突したその瞬間、粒子2つはめり込んでいる事になります。
これは物理的には不自然ですので、めり込みを解決する必要があります。
今回はめり込んだ場合、物理的な拘束条件の下で位置をめり込んでいないところまで更新する、という処理で解決します。
そこで、まずは粒子を、めり込みが$0$になる様な位置まで引き戻す手法を取ります。
Richardson (1994)では、

\begin{eqnarray}
\tilde{\mathbf{x}}_i & = & \mathbf{x}_i - \frac{d_{ij}}{2} \mathbf{n}_{ij}\\
\tilde{\mathbf{x}}_j & = & \mathbf{x}_j + \frac{d_{ij}}{2} \mathbf{n}_{ij}
\end{eqnarray}

という様にしていますが、これはmassが大きく違う2つの物がぶつかったにも関わらず、同じだけ引き離すというちょっと不自然な動作をしそうな気がしますので、
重心を保存させつつ

\begin{eqnarray}
\tilde{\mathbf{x}}_i & = & \mathbf{x}_i - \frac{m_j}{m_i + m_j} d_{ij} \mathbf{n}_{ij}\\
\tilde{\mathbf{x}}_j & = & \mathbf{x}_j + \frac{m_i}{m_i + m_j} d_{ij} \mathbf{n}_{ij}
\end{eqnarray}

としてもいいかもしれません。

ペナルティ法

引き離し法は、粒子同士が滅多に衝突しない系なら良いのですが、大量に粒子が存在して頻繁な衝突を引き起こすような稠密系だと、
2つの粒子に挟まれた粒子が高速で振動をするなど、あまりよろしくない挙動もあります。
そこで、例えば、ペナルティ法と呼ばれる手法では、貫通深度に比例した力を粒子に加えます。
この時、バネ定数として$k$と置くと、

\begin{eqnarray}
\frac{\mathrm{d}^2 \mathbf{x}_i}{\mathrm{d} t^2} = - k d_{ij}\frac{\mathbf{x}_j - \mathbf{x}_i}{|\mathbf{x}_j - \mathbf{x}_i|}
\end{eqnarray}

となります。
ただし、この手法は僕は使ったことがないので、どの$k$の値が正しいのか、よく知りません。
(でもこの方法ってもしかしてLagrangian書き下すと撃力法と言ってる事大して変わらない？)

多粒子系

上記のやり方は、2粒子間での衝突ですので、多粒子系では一つの$i$粒子に対して複数の$j$粒子が存在します。
なので、全粒子からの反発寄与を足し合わせまして、

\begin{eqnarray}
\tilde{\mathbf{x}}_i & = & \mathbf{x}_i - \sum_j \frac{m_j}{m_i + m_j} d_{ij} \mathbf{n}_{ij}\\
\tilde{\mathbf{v}}_i & = & \mathbf{v}_i + \sum_j (1 + \epsilon_\mathrm{n}) \frac{m_j}{m_i + m_j} \mathbf{v}_{ij}^\mathrm{n}
\end{eqnarray}

となります。

実装

各$i$粒子に関して、衝突によって誘起される運動量と位置の変化をバッファに足しこむという事を行い、それを全$i$粒子に関して行ったら、バッファをフラッシュするという方法を取ります。
今回、vec3<double>という3次元ベクタ用のクラスが用意されているものとします。
内積*と絶対値関数abs( )を用意しておきまして、あとは

struct{
    vec3<double> x, v, dx, dv;
    double m, r;
} ptcl[N];
for(int i = 0 ; i < N ; ++ i){
    ptcl[i].dx = ptcl[i].dv = 0;
    for(int j = 0 ; j < N ; ++ j){
        const vec3<double> xij = ptcl[j].x - ptcl[i].x;
        const vec3<double> vij = ptcl[j].v - ptcl[i].v;
        const vec3<double> nij = xij / abs(xij);
        if(abs(xij) < ptcl[i].r + ptcl[j].r){
            ptcl[i].dx += - ptcl[j].m / (ptcl[i].m + ptcl[j].m) * (ptcl[i].r + ptcl[j].r - abs(xij)) * nij;
            ptcl[i].dv += (1 + eps_n) * ptcl[j].m / (ptcl[i].m + ptcl[j].m) * (vij * nij) * nij;
        }
    }
}
for(int i = 0 ; i < N ; ++ i){
    ptcl[i].x += ptcl[i].dx;
    ptcl[i].v += ptcl[i].dv;
}

です。ね、簡単でしょう？

FDPS ver.

えー、はい、まぁその、結論から書くとこうなんですけど、長いな！解説は別記事にしよう！

#include <particle_simulator.hpp>
class Force{
    public:
    PS::F64vec pos;
    PS::F64vec vel;
    void clear(){
        pos = vel = 0;
    }
};

class FP{
    public:
    PS::S64    id;
    PS::F64    mass;
    PS::F64vec pos;
    PS::F64vec vel;
    PS::F64    rad;

    PS::F64vec getPos() const{
        return pos;
    }
    PS::F64 getCharge() const{
        return mass;
    }
    PS::F64 getRSearch() const{
        return 2.0 * rad;
    }

    void setPos(const PS::F64vec& pos){
        this->pos = pos;
    }

    void copyFromForce(const Force& diff){
        pos += diff.pos;
        vel += diff.vel;
    }
    void copyFromFP(const FP& fp){
        *this = fp;
    }
    void writeAscii(FILE* fp) const{
        fprintf(fp, "%lld\t%g\t%g\t%g\t%g\t%g\t%g\n", id, mass, rad, pos.x, pos.y, vel.x, vel.y);
    }
};

class CollisionReactance{
    static const PS::F64 eps_n = 0.1;
    public:
    void operator () (const FP* const ep_i, const PS::S32 Nip, const FP* const ep_j, const PS::S32 Njp, Force* const diff){
        for(PS::S32 i = 0 ; i < Nip ; ++ i){
            const FP& ith = ep_i[i];
            for(PS::S32 j = 0 ; j < Njp ; ++ j){
                const FP& jth = ep_j[j];
                const PS::F64vec dr = jth.pos - ith.pos;
                const PS::F64vec dv = jth.vel - ith.vel;
                const PS::F64    pd = ith.rad + jth.rad - sqrt(dr * dr);
                if(pd < 0 || dr * dr <= 0.0) continue;
                const PS::F64vec n = dr / sqrt(dr * dr);
                diff[i].pos += - jth.mass / (ith.mass + jth.mass) * pd * n;
                diff[i].vel += (1.0 + eps_n) * jth.mass / (ith.mass + jth.mass) * (dv * n) * n;
            }
        }
    }
};

template <class ThisPtcl> double GetGlobalTimestep(const PS::ParticleSystem<ThisPtcl>& ptcl){
    return 1./128.;
}

struct system_t{
    PS::F64 time, dt, end_time;
    unsigned int step;
    system_t() : time(0.0), step(0), end_time(0.0), dt(1.0e+30){
    }
};

template <class ThisPtcl> void SetIC(PS::ParticleSystem<FP>& ptcl, PS::DomainInfo& dinfo, system_t& sysinfo){
    ptcl.setNumberOfParticleLocal(0);
    //dinfo.setBoundaryCondition(PS::BOUNDARY_CONDITION_PERIODIC_XYZ);
    dinfo.setBoundaryCondition(PS::BOUNDARY_CONDITION_PERIODIC_XY);
    //dinfo.setPosRootDomain(PS::F64vec(0.0, 0.0, 0.0), PS::F64vec(1.0, 1.0, 1.0));
    dinfo.setPosRootDomain(PS::F64vec(0.0, 0.0), PS::F64vec(1.0, 1.0));
    sysinfo.end_time = 1.0;
    if(PS::Comm::getRank() != 0) return;
    ptcl.setNumberOfParticleLocal(256);
    for(int i = 0 ; i < ptcl.getNumberOfParticleLocal() ; ++ i){
        ptcl[i].pos.x = rand()/(double)(RAND_MAX);
        ptcl[i].pos.y = rand()/(double)(RAND_MAX);
        //ptcl[i].pos.z = rand()/(double)(RAND_MAX);
        ptcl[i].vel.x = 1.0 - 2.0 * rand()/(double)(RAND_MAX);
        ptcl[i].vel.y = 1.0 - 2.0 * rand()/(double)(RAND_MAX);
        //ptcl[i].vel.z = 1.0 - 2.0 * rand()/(double)(RAND_MAX);
        ptcl[i].vel *= 0.1;
        ptcl[i].mass = 1.0 / ptcl.getNumberOfParticleLocal();
        ptcl[i].rad = 0.01;
    }
}

template <class ThisPtcl> void OutputFileWithTimeInterval(PS::ParticleSystem<ThisPtcl>& ptcl, const system_t& sysinfo){
    const int NumberOfSnapshot = 10;

    static PS::F64 time = sysinfo.time;
    static PS::S64 step = sysinfo.step;
    if(sysinfo.time >= time){
        char filename[256];
        sprintf(filename, "result/%05d", step);
        ptcl.writeParticleAscii(filename, "%s_%05d_%05d.dat");
        if(PS::Comm::getRank() == 0){
            std::cout << "output " << filename << "." << std::endl;
        }
        time += sysinfo.end_time / NumberOfSnapshot;
        ++ step;
    }
}

int main(int argc, char *argv[]){
    PS::Initialize(argc, argv);
    PS::ParticleSystem<FP> ptcl;
    ptcl.initialize();
    PS::TreeForForceShort<Force, FP, FP>::Symmetry tree_coll;
    PS::DomainInfo dinfo;
    dinfo.initialize();

    system_t sysinfo;
    SetIC<FP>(ptcl, dinfo, sysinfo);
    tree_coll.initialize(ptcl.getNumberOfParticleLocal());

    for(sysinfo.time = 0, sysinfo.step = 0 ; sysinfo.time < sysinfo.end_time ; sysinfo.time += sysinfo.dt, ++ sysinfo.step){
        dinfo.decomposeDomainAll(ptcl);
        ptcl.exchangeParticle(dinfo);
        sysinfo.dt = GetGlobalTimestep<FP>(ptcl);
        tree_coll.calcForceAllAndWriteBack(CollisionReactance(), ptcl, dinfo);
        for(int i = 0 ; i < ptcl.getNumberOfParticleLocal() ; ++ i) ptcl[i].pos += ptcl[i].vel * sysinfo.dt;
        ptcl.adjustPositionIntoRootDomain(dinfo);
        if(PS::Comm::getRank() == 0) std::cout << "time = " << sysinfo.time << " (dt = " << sysinfo.dt << ")" << std::endl;
        OutputFileWithTimeInterval<FP>(ptcl, sysinfo);
    }

    PS::Finalize();
    return 0;
}