概要

本記事は、OpenFOAMのsimpleFomaソルバー（定常乱流場解析）の残差をプロットするPython スクリプトを作成したメモである。

使用環境

使用環境は以下になる。
- OpenFOAM ver2312 - Python 3.10

OpenFOAMにおける設定

simpleFOAMの計算において、logファイルに残差が出力されるが、logファイルだけでは残差を抽出するには、少々フォーマットが煩雑で扱いにくい。幸い、OpemFOAMでは残差だけを抽出するような設定が準備されている。以下のように、/system/controlDictファイルに残差を出力するように命令を追記する。
continuityError1は連続式の残差（保存則）、residualsは速度、乱流諸量の残差の出力設定である。

functions
{
    continuityError1
    {
        type  continuityError;
        libs
        (
            fieldFunctionObjects
        ); // Mandatory enties (unmodifiable) // Optional entries (runtime modifiable)
        phi phi;  //Optional (inrerited) entries
        writePrecision 8;
        writeToFile yes;
        useUserTime yes;
        region region0;
        enabled yes;
        log yes;
        timeStart 0;
        timeEnd 500;
        executeControl timeStep;
        executeInterval 1;
        writeControl timeStep;
        writeInterval 50;
    }

    residuals
    {
        type    solverInfo;
        libs    ("libutilityFunctionObjects.so");
        fields  (U p k omega);
    }
}

残差の出力

残差は以下のとおり、/postProcessing/ディレクトリに出力される。
- 連続式の残差　/postProcessing/continuityError1/0/continuityError.dat - 流速、乱流諸量の残差　/postProcessing/residuals/0/solverInfo.dat

残差のプロット

上記で出力された残差データファイルを読み込み、プロットする。以下、プロットするのに作成したPython スクリプトを示す。

流速および乱流諸量の残差をプロットするPython スクリプト例。

import pandas as pd
import matplotlib.pyplot as plt
import os

# 残差ファイルの読み込み
def read_residual():
    df_residual = pd.read_table('./postProcessing/residuals/0/solverInfo.dat', skiprows=1)
    df_residual = pd.DataFrame(df_residual)
    return df_residual
    #print(df_residual.columns)

# 残差プロット
def residual_plot(df):
    init_res = [data for data in df.columns if "initial" in data]
    final_res = [data for data in df.columns if "final" in data]

    # initial residualプロット
    df.plot(x="# Time        ", y=init_res)
    plt.grid()
    plt.legend(loc='best', bbox_to_anchor=(1, 1))
    plt.yscale('log')
    plt.xlabel('time')
    plt.ylabel('initial residual')
    plt.show()

    # final residualプロット
    df.plot(x="# Time        ", y=final_res)
    plt.grid()
    plt.legend(loc='best', bbox_to_anchor=(1, 1))
    plt.yscale('log')
    plt.xlabel('time')
    plt.ylabel('final residual')
    plt.show()

df = read_residual()
residual_plot(df)

連続の式の残差（保存則）をプロットするPython スクリプト例。
連続の式の残差データファイルは、区切りスペースが統一されていないため、まずnumpyのloadtxt関数で読み込み、numpy配列とした後にpandasでデータフレーム化している。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os

# 連続の式の残差ファイルを読み込む
def read_continuity_error():
    df_continuity = np.loadtxt("./postProcessing/continuityError1/0/continuityError.dat")
    df_continuity = pd.DataFrame(df_continuity, columns=['Time','Local','Global','Cumulative'])
    return df_continuity

# 連続の式の残差プロット
def continuity_error_plot(df):
    x = ['Time']
    y = ['Local', 'Global', 'Cumulative']

    for i in range(len(y)):
        plt.plot(df[x[0]], df[y[i]])
        plt.grid()
        plt.legend(loc='best', bbox_to_anchor=(1, 1))
        plt.xlabel(x[0])
        plt.ylabel(y[i] + ' ' + 'continuity')
        plt.show()

df = read_continuity_error()
continuity_error_plot(df)