【OpenFOAM】移流方程式で色々な離散化スキーム(解の振る舞いと計算時間)

こんにちは（@t_kun_kamakiri）

前回の記事では移流方程式における移流項の離散化スキームを変えて結果がどう変わるかを検討しましたが、今回もその続きです。

前回の記事でも離散化スキームについて書きましたのでご参考ください。

モデルの入手

移流方程式のモデルは前回の記事で作成しましたのでそちらをご参考ください。
サンプルモデルも用意したのでご利用ください。

前回の記事

【OpenFOAM】移流方程式で離散化スキームの勉強をする
【OpenFOAM】移流方程式で色々な離散化スキームをPythonで自動実行して試してみた
モデル入手[追記]※最終モデルをgithubにアップしました。

モデルを入手して本記事の内容に沿って編集していくと再現できるでしょう。

フォルダ構成

フォルダ構成は以下のようになっています。

├── main.py
├── orgCase
│   ├── 0
│   │   ├── T
│   │   └── U
│   ├── Allrun
│   ├── constant
│   │   ├── polyMesh
│   │   │   ├── boundary
│   │   │   ├── faces
│   │   │   ├── neighbour
│   │   │   ├── owner
│   │   │   └── points
│   │   └── transportProperties
│   ├── post.foam
│   └── system
│       ├── blockMeshDict
│       ├── controlDict
│       ├── fvSchemes
│       ├── fvSolution
│       ├── sampleDict
│       └── setFieldsDict
├── resultDir
└── schemeTest

├── main.py

├── orgCase

│ ├── 0

│ │ ├── T

│ │ └── U

│ ├── Allrun

│ ├── constant

│ │ ├── polyMesh

│ │ │ ├── boundary

│ │ │ ├── faces

│ │ │ ├── neighbour

│ │ │ ├── owner

│ │ │ └── points

│ │ └── transportProperties

│ ├── post.foam

│ └── system

│ ├── blockMeshDict

│ ├── controlDict

│ ├── fvSchemes

│ ├── fvSolution

│ ├── sampleDict

│ └── setFieldsDict

├── resultDir

└── schemeTest

orgCase：元のケースファイル
schemeTest：スキームの一覧
Allrun：計算実行スクリプト
main.py：計算の自動実行スクリプト

スカラー輸送方程式（scalarTransportFoam）に計算時間を追加

scalarTransportFoamは以下の偏微分方程式を解く非定常のソルバーです。

$\begin{align*} \frac{\partial T}{\partial t}+\nabla\cdot(\boldsymbol{u}T)-\nabla\cdot (D_{T}\nabla T)=\boldsymbol{S}_{T}\tag{1} \end{align*}$

あわせて読みたい

scalarTransportFoam

本記事では(1)を少し簡単にするために、

生成項を0： $\boldsymbol{S}_{T}=0$
拡散項を0： $D_{T}=0$
流速は一定： $\boldsymbol{u}=(0.1, 0, 0)$

と設定するようにします。

$\Delta x = \frac{1}{100} = 0.01$ [m]
$\Delta t = 0.01$ [s]

なのでクーラン数 $C=\frac{u\Delta t}{\Delta x }=0.1$ としています。
クーラン数1だとSuperBeeでエラー出たので少し小さめにしておきました。

初期状態は以下のようになっています。

理論解は以下のように形を変えずに速度に乗って進行します。

しかし、微分方程式を離散化して代数的に解くためどうしても近似的な解になってしまいます。どのような離散化スキームを使うとどういう結果になるのかを知ることが本記事の目的であります。

もう一つの目的は、離散化スキームによって計算時間が異なってくるため計算時間を調査することです。しかし、困ったことにscalarTransportFoamはデフォルトで計算時間を出力する関数が使われていないため、計算時間がわかりませんでした。

そこで計算時間を知るために以下の方法を考えました。

system/controlDictにcodeを埋め込む
→上手くいかなかったので断念（自分がわかっていないだけなのでわかる方いたらコメント頂けると助かります_(._.)_）
追記：「time -f ‘,%U’ scalarTransportFoam > log.scalarTransportFoam」でコマンドを実行するのに使用したユーザーCPU時間が計測できるようです。
なので②のようにCPU計測のための関数を追加しなくても良かったです…
カスタマイズする
関数「runTime.printExecutionTime(Info); //Add」を追加するだけなのでこちらを採用します。

②の方法を採用しソルバにちょこっと関数を追加することで計算時間を出力するようにしました。

では、以下のコマンドでソルバをコピーします。

cp -r $FOAM_APP/solvers/basic/scalarTransportFoam MyscalarTransportFoam

1	cp -r $FOAM_APP/solvers/basic/scalarTransportFoam MyscalarTransportFoam

今回は使い捨てコードでローカルにコピーしても動作するコードなのでローカルにコピーして編集します。

フォルダ移動します。

cd MyscalarTransportFoam

1	cd MyscalarTransportFoam

フォルダ構成は以下のようになっています。

.
├── Make
│   ├── files
│   └── options
├── scalarTransportFoam.C
└── createFields.H

├── Make

│ ├── files

│ └── options

├── scalarTransportFoam.C

└── createFields.H

ファイルの名前を変えます。

mv scalarTransportFoam.C MyscalarTransportFoam.C

1	mv scalarTransportFoam.C MyscalarTransportFoam.C

MyscalarTransportFoam.Cの時間のループの中に「runTime.printExecutionTime(Info); //Add」を追加して計算時間を出力するようにします。

MyscalarTransportFoam.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    scalarTransportFoam

Group
    grpBasicSolvers

Description
    Passive scalar transport equation solver.

    \heading Solver details
    The equation is given by:

    \f[
        \ddt{T} + \div \left(\vec{U} T\right) - \div \left(D_T \grad T \right)
        = S_{T}
    \f]

    Where:
    \vartable
        T       | Passive scalar
        D_T     | Diffusion coefficient
        S_T     | Source
    \endvartable

    \heading Required fields
    \plaintable
        T       | Passive scalar
        U       | Velocity [m/s]
    \endplaintable

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "fvOptions.H"
#include "simpleControl.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Passive scalar transport equation solver."
    );

    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMesh.H"

    simpleControl simple(mesh);

    #include "createFields.H"

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nCalculating scalar transport\n" << endl;

    #include "CourantNo.H"

    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        while (simple.correctNonOrthogonal())
        {
            fvScalarMatrix TEqn
            (
                fvm::ddt(T)
              + fvm::div(phi, T)
              - fvm::laplacian(DT, T)
             ==
                fvOptions(T)
            );

            TEqn.relax();
            fvOptions.constrain(TEqn);
            TEqn.solve();
            fvOptions.correct(T);
        }

        runTime.write();
        runTime.printExecutionTime(Info); //Add
    }

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

/*---------------------------------------------------------------------------*\

========= |

\\ / F ield | OpenFOAM: The Open Source CFD Toolbox

\\ / O peration |

\\ / A nd | www.openfoam.com

\\/ M anipulation |

-------------------------------------------------------------------------------

License

This file is part of OpenFOAM.

OpenFOAM is free software: you can redistribute it and/or modify it

under the terms of the GNU General Public License as published by

the Free Software Foundation, either version 3 of the License, or

(at your option) any later version.

OpenFOAM is distributed in the hope that it will be useful, but WITHOUT

ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License

for more details.

You should have received a copy of the GNU General Public License

along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.

Application

scalarTransportFoam

Group

grpBasicSolvers

Description

Passive scalar transport equation solver.

\heading Solver details

The equation is given by:

\f[

\ddt{T} + \div \left(\vec{U} T\right) - \div \left(D_T \grad T \right)

= S_{T}

\f]

Where:

\vartable

T | Passive scalar

D_T | Diffusion coefficient

S_T | Source

\endvartable

\heading Required fields

\plaintable

T | Passive scalar

U | Velocity [m/s]

\endplaintable

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"

#include "fvOptions.H"

#include "simpleControl.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])

{

argList::addNote

(

"Passive scalar transport equation solver."

);

#include "addCheckCaseOptions.H"

#include "setRootCaseLists.H"

#include "createTime.H"

#include "createMesh.H"

simpleControl simple(mesh);

#include "createFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nCalculating scalar transport\n" << endl;

#include "CourantNo.H"

while (simple.loop())

{

Info<< "Time = " << runTime.timeName() << nl << endl;

while (simple.correctNonOrthogonal())

{

fvScalarMatrix TEqn

(

fvm::ddt(T)

+ fvm::div(phi, T)

- fvm::laplacian(DT, T)

fvOptions(T)

);

TEqn.relax();

fvOptions.constrain(TEqn);

TEqn.solve();

fvOptions.correct(T);

}

runTime.write();

runTime.printExecutionTime(Info); //Add

}

Info<< "End\n" << endl;

return 0;

}

// ************************************************************************* //

そして実行ファイルを出力する先を指定してコマンドに名前を付けます。

Make/files

MyscalarTransportFoam.C

EXE = $(FOAM_USER_APPBIN)/MyscalarTransportFoam

MyscalarTransportFoam.C

EXE = $(FOAM_USER_APPBIN)/MyscalarTransportFoam

以下のコマンドコンパイルします。

wmake

wmake

これで「$FOAM_USER_APPBIN」にMyscalarTransportFoamという実行ファイルが生成されています。ちなみに環境変数は「$FOAM_USER_APPBIN=/home/ユーザ名/OpenFOAM/kamakiri-v2012/platforms/linux64Gcc63DPInt32Opt/bin」となっています。

では、計算時間を出力してくれるかを確認しましょう。
フォルダを移動して、

cd ../orgCase

1	cd ../orgCase

以下のコマンドで計算を実行します。

MyscalarTransportFoam

1	MyscalarTransportFoam

計算のログに「ExecutionTime = 0.05 s 」などが出ていればオッケーです。

Time = 0.01

DILUPBiCGStab:  Solving for T, Initial residual = 0.00552486, Final residual = 6.37066e-19, No Iterations 1
ExecutionTime = 0.05 s  ClockTime = 0 s

Time = 0.02

DILUPBiCGStab:  Solving for T, Initial residual = 0.00558659, Final residual = 7.33487e-19, No Iterations 1
ExecutionTime = 0.05 s  ClockTime = 0 s

Time = 0.03

DILUPBiCGStab:  Solving for T, Initial residual = 0.00562634, Final residual = 1.54832e-18, No Iterations 1
ExecutionTime = 0.05 s  ClockTime = 0 s
(以下省略)

Time = 0.01

DILUPBiCGStab: Solving for T, Initial residual = 0.00552486, Final residual = 6.37066e-19, No Iterations 1

ExecutionTime = 0.05 s ClockTime = 0 s

Time = 0.02

DILUPBiCGStab: Solving for T, Initial residual = 0.00558659, Final residual = 7.33487e-19, No Iterations 1

ExecutionTime = 0.05 s ClockTime = 0 s

Time = 0.03

DILUPBiCGStab: Solving for T, Initial residual = 0.00562634, Final residual = 1.54832e-18, No Iterations 1

ExecutionTime = 0.05 s ClockTime = 0 s

(以下省略)

計算時間の桁数はもう少し多くほしいですが・・・・

scalarTransportFoamの処理時間を知りたかったけど記述がなかったから、
runTime.printExecutionTime(Info);
を追加してみたけどコードを追っていくと結局<ctime>を使っていて10[ms]程度の分解能しかないというところまで理解した_(._.)_ https://t.co/Hb7mWab7WF

— カマキリ技術士勉強中 (@t_kun_kamakiri) February 19, 2023

複数回計算させて平均を取るかもう少し大きなモデルで解析をするかなど考えないといけませんね。

ひとまず正常に動作しているようであれば、初期状態に戻しておきましょう。
以下のコマンドでメッシュ状態はそのままで計算した結果だけを削除することができます。

foamListTimes -rm

1	foamListTimes -rm

では、Pythonコードのあるフォルダに戻ります。

cd ..

cd ..

前回の記事よりも多く色々なスキームを試した

schemeTestというファイルにスキームの一覧を記載します。
前回の記事よりスキームは多くなっています。

schemeTest

upwind 
linear
cubic
QUICK
Minmod
midPoint
SFCD
SuperBee
vanLeer
vanAlbada
UMIST
MUSCL
LUST grad(T)
linearUpwind grad(T)
linearUpwind faceLimited grad(T)
linearUpwind cellLimited grad(T)
linearUpwind cellLimited2 grad(T)
linearUpwind cellLimited31 grad(T)
linearUpwind cellLimited32 grad(T)
linearUpwind cellLimited33 grad(T)
quadraticUpwindFit 1
cubicUpwindFit 1
limitedLinear 1
limitedCubic 1
limitedLinear 0.2
limitedLinear 0.5
limitedLinear 1.0
limitedCubic 0.2
limitedCubic 0.5
limitedCubic 1.0
Gamma 0.2
Gamma 0.5
Gamma 1.0
limitWith QUICK Minmod
limitWith QUICK SuperBee
limitWith QUICK vanAlbada

upwind

linear

cubic

QUICK

Minmod

midPoint

SFCD

SuperBee

vanLeer

vanAlbada

UMIST

MUSCL

LUST grad(T)

linearUpwind grad(T)

linearUpwind faceLimited grad(T)

linearUpwind cellLimited grad(T)

linearUpwind cellLimited2 grad(T)

linearUpwind cellLimited31 grad(T)

linearUpwind cellLimited32 grad(T)

linearUpwind cellLimited33 grad(T)

quadraticUpwindFit 1

cubicUpwindFit 1

limitedLinear 1

limitedCubic 1

limitedLinear 0.2

limitedLinear 0.5

limitedLinear 1.0

limitedCubic 0.2

limitedCubic 0.5

limitedCubic 1.0

Gamma 0.2

Gamma 0.5

Gamma 1.0

limitWith QUICK Minmod

limitWith QUICK SuperBee

limitWith QUICK vanAlbada

グラフに計算時間を載せるようにした

今回は計算時間も測定してグラフに載せるようにしました。
計算時間はあくまで1回の計算をさせた場合での結果です。他のアプリの使用状態によってはCPUの使用状態も変わってくるので、本来は複数回計算させた際の平均値を取るのが良いでしょう。

from pathlib import Path
import shutil
import subprocess
import numpy as np
from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile
import matplotlib.pyplot as plt
import os

# スキームファイルを読み込む
def open_SchemeList():
    with open('./schemeTest') as f:
        scheme_list = [s.strip() for s in f.readlines()]
    return scheme_list
    
# baseをコピーする
def clone_file(orgCase, resultDir):
    "orgCaseをresultDirに重複しないようにコピー"
    baseName = Path(orgCase).name
 
    addName = None
    if addName!=None:
        baseName += f'_{addName}'
 
    Path(resultDir).mkdir(exist_ok=True, parents=True)
 
    n=0
    newCase = Path(resultDir) / f'{baseName}_{n}'
 
    while newCase.is_dir():
        n += 1
        newCase = Path(resultDir) / f'{baseName}_{n}'
 
    shutil.copytree(orgCase, newCase)
    
    return newCase

# スキームを入れ替える
def new_parameter(newScheme, newCase):
    fileScheme = ParsedParameterFile(newCase / 'system/fvSchemes')
    fileScheme.content['divSchemes']
    fileScheme.content['divSchemes']['div(phi,T)'] = f'Gauss {newScheme}'
    fileScheme.writeFile()
    
# Allrunスクリプトを実行する(postProcess)
def Allrun(newCase):
    "Allrunを実行"
    out_run = subprocess.check_output(['./Allrun'], cwd=newCase)
    return out_run.decode()

# 計算時間の出力
def ExecutionTime_func(newCase, resultDir):
    print(f'Create figure {Path(newCase)/"log.MyscalarTransportFoam"}')
    with open( Path(newCase)/'log.MyscalarTransportFoam') as f:
        log_list = [s.strip() for s in f.readlines() if ('ExecutionTime' in s) or ('End' in s)]

    ExecutionTime = log_list[-2:][0].split('s')[0] + 's'
    print(ExecutionTime)
    if log_list[-1] != 'End':
        print("計算途中ストップ")
    
    return ExecutionTime

# グラフにする
def graph(scheme, ExecutionTime, newCase, resultDir):
    time_list = os.listdir(Path(newCase) / f'postProcessing/samples')
    plt.figure()
    for time in time_list[::7]:
        temp_data = np.loadtxt(Path(newCase) / f'postProcessing/samples/{time}/x_T_T.xy').T
        plt.plot(temp_data[0], temp_data[1],label=f'time = {time}')    

    plt.grid()
    plt.yticks(np.arange(-0.5,2,0.25))
    plt.legend()
    plt.xlabel("X(m)", fontsize=12)
    plt.ylabel("Temperature(K)", fontsize=12)
    plt.text(0, 2,ExecutionTime)
    plt.title(f"Temperature Advection Equation({scheme})", fontsize=14)
    
    # 画像賦存用データのフォルダ作成
    result_graph_dir = Path(resultDir)/'savefig'
    if(os.path.isdir(result_graph_dir) == False):
        os.mkdir(result_graph_dir)
    
    # グラフを保存
    plt.savefig(f'{result_graph_dir}/{scheme}.png')
    plt.close()
    
# 全部のグラフを2列ずつにしてまとめる
def Allgraph_png(scheme_list, ExecutionTime_list):
    fig = plt.figure(figsize = (14, 50))
    for i, scheme in enumerate(scheme_list):
         #グラフを表示する領域を，figオブジェクトとして作成．
        col = len(scheme_list)//2 + 1 
        ax1 = fig.add_subplot(col,2,i+1)
        time_list = os.listdir(f'./resultDir/orgCase_{i}/postProcessing/samples')
        
        for time in time_list[::7]:
            temp_data = np.loadtxt(f'./resultDir/orgCase_{i}/postProcessing/samples/{time}/x_T_T.xy').T
            ax1.plot(temp_data[0], temp_data[1],label=f'time = {time}')
        
        # 計算時間
        ExecutionTime = ExecutionTime_list[i]
        plt.subplots_adjust(left=0.1, right=0.9, bottom=0.01, top=0.95, wspace=0.6, hspace=0.6)
        ax1.text(0, 1.5,ExecutionTime)
        ax1.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0, fontsize=8)  
        ax1.grid()
        ax1.set_yticks(np.arange(-0.5,2,0.25))
        ax1.set_xlabel("X(m)", fontsize=12)
        ax1.set_ylabel("Temperature(K)", fontsize=12)
        ax1.set_title(f"Temperature Advection Equation({scheme})", fontsize=12)


        # グラフを保存
        fig.savefig('./graph_Allschemes.png')
    plt.close()

#================== main ========================
orgCase = 'orgCase'
resultDir = 'resultDir'

# スキームのリスト
scheme_list = open_SchemeList()
print(f'スキームの数 : {len(scheme_list)}')

ExecutionTime_list = []
# 計算実行
for i, scheme in enumerate(scheme_list):
    print(f'{i+1} ======== {scheme} =========')
    newCase = clone_file(orgCase, resultDir)                # 関数を実行する
    new_parameter(scheme, newCase)                          # スキームを入れ替える
    Allrun(newCase)                                         # 計算を実行
    ExecutionTime = ExecutionTime_func(newCase, resultDir)  # 計算時間を取得
    ExecutionTime_list.append(ExecutionTime)
    graph(scheme, ExecutionTime, newCase, resultDir)         # グラフ化
    
# 全てのグラフをまとめる
Allgraph_png(scheme_list, ExecutionTime_list)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

from pathlib import Path

import shutil

import subprocess

import numpy as np

from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile

import matplotlib.pyplot as plt

import os

# スキームファイルを読み込む

def open_SchemeList():

with open('./schemeTest') as f:

scheme_list = [s.strip() for s in f.readlines()]

return scheme_list

# baseをコピーする

def clone_file(orgCase, resultDir):

"orgCaseをresultDirに重複しないようにコピー"

baseName = Path(orgCase).name

addName = None

if addName!=None:

baseName += f'_{addName}'

Path(resultDir).mkdir(exist_ok=True, parents=True)

n=0

newCase = Path(resultDir) / f'{baseName}_{n}'

while newCase.is_dir():

n += 1

newCase = Path(resultDir) / f'{baseName}_{n}'

shutil.copytree(orgCase, newCase)

return newCase

# スキームを入れ替える

def new_parameter(newScheme, newCase):

fileScheme = ParsedParameterFile(newCase / 'system/fvSchemes')

fileScheme.content['divSchemes']

fileScheme.content['divSchemes']['div(phi,T)'] = f'Gauss {newScheme}'

fileScheme.writeFile()

# Allrunスクリプトを実行する(postProcess)

def Allrun(newCase):

"Allrunを実行"

out_run = subprocess.check_output(['./Allrun'], cwd=newCase)

return out_run.decode()

# 計算時間の出力

def ExecutionTime_func(newCase, resultDir):

print(f'Create figure {Path(newCase)/"log.MyscalarTransportFoam"}')

with open( Path(newCase)/'log.MyscalarTransportFoam') as f:

log_list = [s.strip() for s in f.readlines() if ('ExecutionTime' in s) or ('End' in s)]

ExecutionTime = log_list[-2:][0].split('s')[0] + 's'

print(ExecutionTime)

if log_list[-1] != 'End':

print("計算途中ストップ")

return ExecutionTime

# グラフにする

def graph(scheme, ExecutionTime, newCase, resultDir):

time_list = os.listdir(Path(newCase) / f'postProcessing/samples')

plt.figure()

for time in time_list[::7]:

temp_data = np.loadtxt(Path(newCase) / f'postProcessing/samples/{time}/x_T_T.xy').T

plt.plot(temp_data[0], temp_data[1],label=f'time = {time}')

plt.grid()

plt.yticks(np.arange(-0.5,2,0.25))

plt.legend()

plt.xlabel("X(m)", fontsize=12)

plt.ylabel("Temperature(K)", fontsize=12)

plt.text(0, 2,ExecutionTime)

plt.title(f"Temperature Advection Equation({scheme})", fontsize=14)

# 画像賦存用データのフォルダ作成

result_graph_dir = Path(resultDir)/'savefig'

if(os.path.isdir(result_graph_dir) == False):

os.mkdir(result_graph_dir)

# グラフを保存

plt.savefig(f'{result_graph_dir}/{scheme}.png')

plt.close()

# 全部のグラフを2列ずつにしてまとめる

def Allgraph_png(scheme_list, ExecutionTime_list):

fig = plt.figure(figsize = (14, 50))

for i, scheme in enumerate(scheme_list):

#グラフを表示する領域を，figオブジェクトとして作成．

col = len(scheme_list)//2 + 1

ax1 = fig.add_subplot(col,2,i+1)

time_list = os.listdir(f'./resultDir/orgCase_{i}/postProcessing/samples')

for time in time_list[::7]:

temp_data = np.loadtxt(f'./resultDir/orgCase_{i}/postProcessing/samples/{time}/x_T_T.xy').T

ax1.plot(temp_data[0], temp_data[1],label=f'time = {time}')

# 計算時間

ExecutionTime = ExecutionTime_list[i]

plt.subplots_adjust(left=0.1, right=0.9, bottom=0.01, top=0.95, wspace=0.6, hspace=0.6)

ax1.text(0, 1.5,ExecutionTime)

ax1.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0, fontsize=8)

ax1.grid()

ax1.set_yticks(np.arange(-0.5,2,0.25))

ax1.set_xlabel("X(m)", fontsize=12)

ax1.set_ylabel("Temperature(K)", fontsize=12)

ax1.set_title(f"Temperature Advection Equation({scheme})", fontsize=12)

# グラフを保存

fig.savefig('./graph_Allschemes.png')

plt.close()

#================== main ========================

orgCase = 'orgCase'

resultDir = 'resultDir'

# スキームのリスト

scheme_list = open_SchemeList()

print(f'スキームの数 : {len(scheme_list)}')

ExecutionTime_list = []

# 計算実行

for i, scheme in enumerate(scheme_list):

print(f'{i+1} ======== {scheme} =========')

newCase = clone_file(orgCase, resultDir) # 関数を実行する

new_parameter(scheme, newCase) # スキームを入れ替える

Allrun(newCase) # 計算を実行

ExecutionTime = ExecutionTime_func(newCase, resultDir) # 計算時間を取得

ExecutionTime_list.append(ExecutionTime)

graph(scheme, ExecutionTime, newCase, resultDir) # グラフ化

# 全てのグラフをまとめる

Allgraph_png(scheme_list, ExecutionTime_list)

Pythonを実行します。

Pythonで「system/fvSchemes」の「divSchemes」スキームを変えています。

divSchemes
{
    default         none;
    div(phi,T) Gauss upwind grad(T);
}

divSchemes

{

default none;

div(phi,T) Gauss upwind grad(T);

}

結果はスキームごとに分けて「resultDir/savefig」に画像を作成し保存、全てのスキームでの結果を「graph_AllShemes.png」に画像を作成し保存しています。

pythonを実行します。

python3 main.py

1	python3 main.py

以下のように計算が進みます。

スキームの数 : 36
1 ======== linear =========
Create figure resultDir/orgCase_36/log.MyscalarTransportFoam
ExecutionTime = 3.43s
2 ======== cubic =========
Create figure resultDir/orgCase_37/log.MyscalarTransportFoam
ExecutionTime = 1.95 s
(以下省略)

スキームの数 : 36

1 ======== linear =========

Create figure resultDir/orgCase_36/log.MyscalarTransportFoam

ExecutionTime = 3.43s

2 ======== cubic =========

Create figure resultDir/orgCase_37/log.MyscalarTransportFoam

ExecutionTime = 1.95 s

(以下省略)

結果を確認

結果は「resultDir/savefig」に画像として保存するようにしていますのでじっくり眺めることができます。
全体のグラフの比較は「graph_Allschemes.png」で行うことができます。

graph_Allschemes.png

やはり複数回計算させた結果を計算時間とするべきでしょうね。なので、今回の計算時間は参考程度ということです。

アニメーション作成

アニメーションにした方がわかりやすいかと思ったのでアニメーションにしてみました。

移流方程式というナビエストークス方程式よりは単純な形をした一般解のある偏微分方程式ですが、移流項 $u\frac{\partial T}{\partial x}$ の部分のスキーム違いで結果が大きく変わりますね。
特に初期状態の矩形波の変化が激しい部分の微分を離散化が難しく、離散化スキームをどれにするかによって結果が変わるということがわかりました。これがもしsin波のような滑らかな初期状態だったり、拡散項があるとあまり結果に違いがでないかもしれません….(試す必要がありますね）

アニメーション作成用のコードは即興で作ったのでお粗末ですがご参考に…
先ほどのPythonコードに以下の関数を追加して実行すればOK。

from matplotlib.animation import ArtistAnimation
# 全部のグラフを2列ずつにしてまとめる
def AllAnimation(scheme_list, ExecutionTime_list):
    ax_list = []
    fig = plt.figure(figsize = (14, 50))
    plt.subplots_adjust(left=0.1, right=0.9, bottom=0.01, top=0.95, wspace=0.6, hspace=0.6)
    time_list = os.listdir(f'./resultDir/orgCase_{0}/postProcessing/samples')
    for i, scheme in enumerate(scheme_list):
        #グラフを表示する領域を，figオブジェクトとして作成．
        col = len(scheme_list)//2 + 1 
        ax_list.append(fig.add_subplot(col,2,i+1))
    
    frames = []
    for time in time_list[::2]:
        frame = []
        print(time)
        title = ax.text(0, 1.05, 'Time: {time}', fontsize=15) # 追加
        for i, scheme in enumerate(scheme_list):
            temp_data = np.loadtxt(f'./resultDir/orgCase_{i}/postProcessing/samples/{time}/x_T_T.xy').T
            
            frame.append(ax_list[i].plot(temp_data[0], temp_data[1], color="blue")[0])
            
            # 計算時間
            ExecutionTime = ExecutionTime_list[i]
            ax_list[i].text(0, 1.5,ExecutionTime)
            ax_list[i].grid()
            ax_list[i].set_yticks(np.arange(-0.5,2,0.25))
            ax_list[i].set_xlabel("X(m)", fontsize=12)
            ax_list[i].set_ylabel("Temperature(K)", fontsize=12)
            ax_list[i].set_title(f"Temperature Advection Equation({scheme})", fontsize=12)
       
        frames.append(frame)
            
    ani = ArtistAnimation(fig, frames, interval=1)
    ani.save("anim1-18.gif", writer="pillow")

# ========== main  ========== 
AllAnimation(scheme_list, ExecutionTime_list)

from matplotlib.animation import ArtistAnimation

# 全部のグラフを2列ずつにしてまとめる

def AllAnimation(scheme_list, ExecutionTime_list):

ax_list = []

fig = plt.figure(figsize = (14, 50))

plt.subplots_adjust(left=0.1, right=0.9, bottom=0.01, top=0.95, wspace=0.6, hspace=0.6)

time_list = os.listdir(f'./resultDir/orgCase_{0}/postProcessing/samples')

for i, scheme in enumerate(scheme_list):

#グラフを表示する領域を，figオブジェクトとして作成．

col = len(scheme_list)//2 + 1

ax_list.append(fig.add_subplot(col,2,i+1))

frames = []

for time in time_list[::2]:

frame = []

print(time)

title = ax.text(0, 1.05, 'Time: {time}', fontsize=15) # 追加

for i, scheme in enumerate(scheme_list):

temp_data = np.loadtxt(f'./resultDir/orgCase_{i}/postProcessing/samples/{time}/x_T_T.xy').T

frame.append(ax_list[i].plot(temp_data[0], temp_data[1], color="blue")[0])

# 計算時間

ExecutionTime = ExecutionTime_list[i]

ax_list[i].text(0, 1.5,ExecutionTime)

ax_list[i].grid()

ax_list[i].set_yticks(np.arange(-0.5,2,0.25))

ax_list[i].set_xlabel("X(m)", fontsize=12)

ax_list[i].set_ylabel("Temperature(K)", fontsize=12)

ax_list[i].set_title(f"Temperature Advection Equation({scheme})", fontsize=12)

frames.append(frame)

ani = ArtistAnimation(fig, frames, interval=1)

ani.save("anim1-18.gif", writer="pillow")

# ========== main ==========

AllAnimation(scheme_list, ExecutionTime_list)

まとめ

スキームは試して結果を眺めてみないとわからないので、本記事のような簡単な移流方程式で色々と試してみました。
今回は結果と計算時間だけを一覧で眺めれるようにしました。スキームによっては解が安定しなかったり、安定していても拡散する傾向にあったりします。離散化スキームの理論的な話は割愛しましたが、勉強するとなぜそのようなことが起こるのかがわかってきて理解が深まります。