<jupyter_start><jupyter_text>論文中の数式をSymPyで再現する 対象論文[Three-dimensional cascaded lattice Boltzmann method: Improved implementation and consistent forcing scheme](https://journals.aps.org/pre/pdf/10.1103/PhysRevE.97.053309) [参考] 論文要約```Phys. Rev. E 73, 066705 (2006)]で提案されたカスケード法または中心モーメントベースの格子ボルツマン法(CLBM)は、非常に優れた数値安定性を持っています。しかし、3次元CLBMのシミュレーションには2つの制約があります。第一に、従来の3次元CLBMの実装は、単一緩和時間（SRT）LBMに比べて煩雑な操作を必要とし、計算コストが非常に高くなります。第二に、一般的な力場を3次元CLBMに正確に組み込むことは困難である。本論文では、3DでCLBMを実装するための改良された手法を紹介します。その主な戦略は、簡略化された中心モーメントセットを採用し、一般的な多緩和時間（GMRT）フレームワークに基づいて中心モーメントベースの衝突演算子を実行することである。次に、最近提案されたCLBMのための一貫した強制スキーム［Fei and Luo, Phys. Rev. E 96, 053307 (2017)］を拡張して、一般的な力場を3D CLBMに組み込む。最近開発された非直交CLBM[Rosis, Phys. Rev. E 95, 013310 (2017)]と比較して、我々の実装は計算コストを大幅に削減することが証明された。非直交CLBMにおいて離散的な平衡分布関数を採用することの不整合を分析し、検証しています。ここで開発された3D CLBMと一貫した強制スキームは、いくつかのカノニカルな力駆動の流れの数値シミュレーションによって検証され、精度、収束性、ノンスリップルールとの整合性の点で非常に優れた特性が強調されている。最後に，ここで開発した3次元CLBMの技術は，3次元MRT-LBMの実装と実行をより効率的にするために応用できる。``` 再現内容論文中、**式(2)** で定義されている2つのベクトル、$$k_{nmp} = \langle f_i | e^m_x e^n_y e^p_y \rangle, \hspace{5mm} \tilde{k}_{nmp} = \langle f_i | (e_x - u_x)^m (e_y - u_y)^n (e_z - u_z)^p \rangle$$について、以下で定義される変換行列 $N$ の各要素を計算する。$$\tilde{k}_{nmp} = N k_{nmp}$$また、添字$\boldsymbol{nmp}$の順序は、次のようになる。$$k_{nmp} = \left[ k_{000},k_{100},k_{010},k_{001},k_{110},k_{101},k_{011},k_{200},k_{020},k_{002},k_{120},k_{102},k_{210},k_{201},k_{012},k_{021},k_{111},k_{220},k_{202},k_{022},k_{211},k_{121},k_{112},k_{122},k_{212},k_{221},k_{222} \right] \\\tilde{k}_{nmp} = \left[ \tilde{k}_{000},\tilde{k}_{100},\tilde{k}_{010},\tilde{k}_{001},\tilde{k}_{110},\tilde{k}_{101},\tilde{k}_{011},\tilde{k}_{200},\tilde{k}_{020},\tilde{k}_{002},\tilde{k}_{120},\tilde{k}_{102},\tilde{k}_{210},\tilde{k}_{201},\tilde{k}_{012},\tilde{k}_{021},\tilde{k}_{111},\tilde{k}_{220},\tilde{k}_{202},\tilde{k}_{022},\tilde{k}_{211},\tilde{k}_{121},\tilde{k}_{112},\tilde{k}_{122},\tilde{k}_{212},\tilde{k}_{221},\tilde{k}_{222} \right]$$  数式で使用するシンボルの定義- $u_x$, $u_y$, $u_z$- $e_x$, $e_y$, $e_z$<jupyter_code>ux, uy, uz = sp.symbols(r"u_x, u_y, u_z")
ux, uy, uz
ex, ey, ez = sp.symbols(r"e_x, e_y, e_z")
ex, ey, ez<jupyter_output><empty_output><jupyter_text>ベクトルの添字の順序を定義する<jupyter_code>k_list = "k000,k100,k010,k001,k110,k101,k011,k200,k020,k002,k120,\
k102,k210,k201,k012,k021,k111,k220,k202,k022,k211,k121,\
k112,k122,k212,k221,k222".split(',')
print(k_list)<jupyter_output>['k000', 'k100', 'k010', 'k001', 'k110', 'k101', 'k011', 'k200', 'k020', 'k002', 'k120', 'k102', 'k210', 'k201', 'k012', 'k021', 'k111', 'k220', 'k202', 'k022', 'k211', 'k121', 'k112', 'k122', 'k212', 'k221', 'k222']<jupyter_text>変換行列 $N$ の各要素を計算する<jupyter_code># 変換行列のサイズ
n_rows, n_column = (27, 27)

# 求めるべき変換行列 N
N = sp.zeros(n_rows, n_column)

for row, k in enumerate(k_list):
    m, n, p = int(k[1]), int(k[2]), int(k[3])
    k_expand = (ex - ux) ** m * (ey - uy) ** n * (ez - uz) ** p
    k_expand = sp.simplify(sp.expand(k_expand))
    for column, k2 in enumerate(k_list):
        m2, n2, p2 = int(k2[1]), int(k2[2]), int(k2[3])
        N[row, column] = k_expand.coeff(ex, n=m2).coeff(ey, n=n2).coeff(ez, n=p2)

N
# あまり意味がないが、逆行列も確認しておく
N.inv()<jupyter_output><empty_output><jupyter_text>変換行列 $N$ を使用したベクトルの変換をC言語として出力する<jupyter_code># まずシンボル定義
k = sp.MatrixSymbol(r'k_nmp', n_rows, 1)
k
k_tilde = N * k
k_tilde
for i in range(n_rows):
    code = sp.ccode(k_tilde[i], assign_to=f"k_tilde[{i}]")
    print(code)<jupyter_output>k_tilde[0] = k_nmp[0];
k_tilde[1] = -u_x*k_nmp[0] + k_nmp[1];
k_tilde[2] = -u_y*k_nmp[0] + k_nmp[2];
k_tilde[3] = -u_z*k_nmp[0] + k_nmp[3];
k_tilde[4] = u_x*u_y*k_nmp[0] - u_x*k_nmp[2] - u_y*k_nmp[1] + k_nmp[4];
k_tilde[5] = u_x*u_z*k_nmp[0] - u_x*k_nmp[3] - u_z*k_nmp[1] + k_nmp[5];
k_tilde[6] = u_y*u_z*k_nmp[0] - u_y*k_nmp[3] - u_z*k_nmp[2] + k_nmp[6];
k_tilde[7] = pow(u_x, 2)*k_nmp[0] - 2*u_x*k_nmp[1] + k_nmp[7];
k_tilde[8] = pow(u_y, 2)*k_nmp[0] - 2*u_y*k_nmp[2] + k_nmp[8];
k_tilde[9] = pow(u_z, 2)*k_nmp[0] - 2*u_z*k_nmp[3] + k_nmp[9];
k_tilde[10] = -u_x*pow(u_y, 2)*k_nmp[0] + 2*u_x*u_y*k_nmp[2] - u_x*k_nmp[8] + pow(u_y, 2)*k_nmp[1] - 2*u_y*k_nmp[4] + k_nmp[10];
k_tilde[11] = -u_x*pow(u_z, 2)*k_nmp[0] + 2*u_x*u_z*k_nmp[3] - u_x*k_nmp[9] + pow(u_z, 2)*k_nmp[1] - 2*u_z*k_nmp[5] + k_nmp[11];
k_tilde[12] = -pow(u_x, 2)*u_y*k_nmp[0] + pow(u_x, 2)*k_nmp[2] + 2*u_x*u_y*k_nmp[1] - 2*u_x*k_nmp[4] - u_y*k_nmp[7] + k_nmp[12];
k_tilde[13] = -pow(u_x, 2)*u_z*k_nmp[0] + pow(u_x, 2)*k_nmp[3] + 2*u[...]
{'fraction_non_alphanumeric': 0.23295070898041864, 'fraction_numerical': 0.11186135494035561, 'mean_word_length': 7.852589641434263, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 15, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '511752', 'license': 'MIT', 'n_tokens_mistral': 3129, 'n_tokens_neox': 2822, 'n_words': 267, 'path': 'notebooks/05-Metaprogramming.ipynb', 'repo_name': 'codemajin/Introduction-to-SymPy', 'size': '4443'}