Portable, Extensible Toolkit for Scientific Computation/PETScの実践的な応用

はじめに

本章では、PETScを実際の科学技術計算問題に適用した具体的な事例をさらに詳しく紹介します。実践的な問題を通じて、PETScの活用方法を深く学びます。取り上げる事例は以下の通りです。

1. 流体力学シミュレーション
2. 構造解析
3. 最適化問題
4. データ同化

流体力学シミュレーション

問題設定

流体力学シミュレーションでは、ナビエ-ストークス方程式を解いて流体の挙動を解析します。以下の方程式を解きます。 $\frac{\partial 𝐮}{\partial t}+(𝐮\cdot \mathrm{\nabla })𝐮=-\mathrm{\nabla }p+\nu {\mathrm{\nabla }}^2𝐮$ ここで、$𝐮$ は速度場、$p$ は圧力、$\nu$ は動粘度です。

PETScによる解法

PETScの非線形ソルバー（SNES）と時間ステッピングソルバー（TS）を使用して、ナビエ-ストークス方程式を解きます。

例: 流体力学シミュレーションのコード

#include "petsc.h"
#include "petscsnes.h"
#include "petscts.h"

PetscErrorCode FormFunction(SNES snes, Vec U, Vec F, void *ctx) {
    // 非線形関数の実装
    Mat J;
    SNESGetJacobian(snes, &J, NULL, NULL, NULL);
    MatMult(J, U, F);
    return 0;
}

int main(int argc, char **argv) {
    PetscInitialize(&argc, &argv, NULL, NULL);

    SNES snes;
    TS ts;
    Vec U, F;
    Mat J;

    // 初期化と設定
    SNESCreate(PETSC_COMM_WORLD, &snes);
    TSCreate(PETSC_COMM_WORLD, &ts);
    VecCreate(PETSC_COMM_WORLD, &U);
    VecCreate(PETSC_COMM_WORLD, &F);
    MatCreate(PETSC_COMM_WORLD, &J);

    // 行列とベクトルの設定
    // ... (省略)

    // 非線形ソルバーの設定
    SNESSetFunction(snes, F, FormFunction, NULL);
    SNESSetFromOptions(snes);

    // 時間ステッピングの設定
    TSSetProblemType(ts, TS_NONLINEAR);
    TSSetRHSFunction(ts, NULL, FormRHSFunction, NULL);
    TSSetFromOptions(ts);

    // ソルバーの実行
    SNESSolve(snes, NULL, U);
    TSSolve(ts, U);

    // 終了処理
    SNESDestroy(&snes);
    TSDestroy(&ts);
    VecDestroy(&U);
    VecDestroy(&F);
    MatDestroy(&J);
    PetscFinalize();
    return 0;
}

構造解析

問題設定

構造解析では、弾性方程式を解いて物体の変形や応力を解析します。以下の方程式を解きます。 $\mathrm{\nabla }\cdot \sigma =𝐟$ ここで、$\sigma$ は応力テンソル、$𝐟$ は外力です。

PETScによる解法

PETScの線形ソルバー（KSP）を使用して、弾性方程式を解きます。

例: 構造解析のコード

#include "petsc.h"
#include "petscksp.h"

int main(int argc, char **argv) {
    PetscInitialize(&argc, &argv, NULL, NULL);

    KSP ksp;
    Vec U, F;
    Mat A;

    // 初期化と設定
    KSPCreate(PETSC_COMM_WORLD, &ksp);
    VecCreate(PETSC_COMM_WORLD, &U);
    VecCreate(PETSC_COMM_WORLD, &F);
    MatCreate(PETSC_COMM_WORLD, &A);

    // 行列とベクトルの設定
    // ... (省略)

    // 線形ソルバーの設定
    KSPSetOperators(ksp, A, A);
    KSPSetFromOptions(ksp);

    // ソルバーの実行
    KSPSolve(ksp, F, U);

    // 終了処理
    KSPDestroy(&ksp);
    VecDestroy(&U);
    VecDestroy(&F);
    MatDestroy(&A);
    PetscFinalize();
    return 0;
}

最適化問題

問題設定

最適化問題では、与えられた制約条件下で目的関数を最小化または最大化します。以下の最適化問題を解きます。 $\underset{𝐱}{\min }f(𝐱)$ ここで、$f(𝐱)$ は目的関数です。

PETScによる解法

PETScの最適化ツール（TAO）を使用して、最適化問題を解きます。

例: 最適化問題のコード

#include "petsc.h"
#include "petsctao.h"

PetscErrorCode FormObjective(Tao tao, Vec X, PetscReal *f, void *ctx) {
    // 目的関数の実装
    Vec g;
    VecDuplicate(X, &g);
    VecNorm(X, NORM_2, f);
    VecDestroy(&g);
    return 0;
}

int main(int argc, char **argv) {
    PetscInitialize(&argc, &argv, NULL, NULL);

    Tao tao;
    Vec X;

    // 初期化と設定
    TaoCreate(PETSC_COMM_WORLD, &tao);
    VecCreate(PETSC_COMM_WORLD, &X);

    // ベクトルの設定
    // ... (省略)

    // 最適化ソルバーの設定
    TaoSetObjective(tao, FormObjective, NULL);
    TaoSetFromOptions(tao);

    // ソルバーの実行
    TaoSolve(tao);

    // 終了処理
    TaoDestroy(&tao);
    VecDestroy(&X);
    PetscFinalize();
    return 0;
}

データ同化

問題設定

データ同化では、観測データと数値モデルを組み合わせて、システムの状態を推定します。以下の最適化問題を解きます。 $\underset{𝐱}{\min }\Vert 𝐲-H(𝐱){\Vert }^2$ ここで、$𝐲$ は観測データ、$H(𝐱)$ は観測演算子です。

PETScによる解法

PETScの最適化ツール（TAO）を使用して、データ同化問題を解きます。

例: データ同化のコード

#include "petsc.h"
#include "petsctao.h"

PetscErrorCode FormObjective(Tao tao, Vec X, PetscReal *f, void *ctx) {
    // 目的関数の実装
    Vec g;
    VecDuplicate(X, &g);
    VecNorm(X, NORM_2, f);
    VecDestroy(&g);
    return 0;
}

int main(int argc, char **argv) {
    PetscInitialize(&argc, &argv, NULL, NULL);

    Tao tao;
    Vec X;

    // 初期化と設定
    TaoCreate(PETSC_COMM_WORLD, &tao);
    VecCreate(PETSC_COMM_WORLD, &X);

    // ベクトルの設定
    // ... (省略)

    // 最適化ソルバーの設定
    TaoSetObjective(tao, FormObjective, NULL);
    TaoSetFromOptions(tao);

    // ソルバーの実行
    TaoSolve(tao);

    // 終了処理
    TaoDestroy(&tao);
    VecDestroy(&X);
    PetscFinalize();
    return 0;
}

練習問題

1. 流体力学シミュレーション: 簡単な流れ（例: ポアズイユ流れ）をシミュレートするプログラムを作成してください。
2. 構造解析: 梁の変形を解析するプログラムを作成し、変形状態を可視化してください。
3. 最適化問題: 二次関数の最小化問題を解くプログラムを作成し、最適解を表示してください。
4. データ同化: 簡単なデータ同化問題を解くプログラムを作成し、推定結果を評価してください。

まとめ

本章では、PETScを実際の科学技術計算問題に適用した具体的な事例を詳しく紹介しました。流体力学シミュレーション、構造解析、最適化問題、データ同化を通じて、PETScの実践的な活用方法を学びました。これらの事例を参考に、自身の問題にPETScを適用してみてください。