緩和係数の最適化

2021年6月20日

はじめに

OptunaBayesian Optimization を用いて OpenFOAM のケースの緩和係数の最適化を試みる。

使用バージョン

  • OpenFOAM v2012
  • Python 3

計算ケース

計算ケースには simpleFoam のチュートリアルケースである pitzDaily を用いる。

Optuna

インストール

$ pip install optuna

テスト

test-op.py

import optuna

def objective(trial):
    x = trial.suggest_uniform("x", 0, 10)
    return (x**2 - 2)**2

study = optuna.create_study()

#optuna.logging.set_verbosity(optuna.logging.INFO)
optuna.logging.set_verbosity(optuna.logging.WARNING)

study.optimize(objective, n_trials=100)
study.best_params

実行

$ python test-op.py

2 の平方根の値が出てくれば OK。

緩和係数の最適化

run-op.py

import sys
import subprocess as sp
import re
import optuna


def run_command_with_output(command):
    proc = sp.Popen(command,
            shell=True, stdout=sp.PIPE, stderr=sp.PIPE)

    while True:
        line = proc.stdout.readline()
        if not line:
            break
        line = line.decode().rstrip("\n")
        yield line

    if proc.wait() != 0:
        while True:
            line = proc.stderr.readline()
            if not line:
                break
            line = line.decode().rstrip("\n")
            yield line
        raise OSError("Command execution faied: %s" % command)


def write_fvSolution(ap, au, at):
    with open("system/fvSolution", "w") as f:
        f.write("""FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.1;
        smoother        GaussSeidel;
    }

    "(U|k|epsilon|omega|f|v2)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent      yes;

    residualControl
    {
        p               1e-2;
        U               1e-3;
        "(k|epsilon|omega|f|v2)" 1e-3;
    }
}

relaxationFactors
{
    fields
    {
""")

        f.write("        p %f;\n" % ap)

        f.write("""    }
    equations
    {
""")

        f.write("        U %f;\n" % au)
        f.write('        ".*" %f;\n' %at)

        f.write("""    }
}
""")


def run(ap, au, at):
    print(ap, au, at)
    write_fvSolution(ap, au, at)
    time = 0
    for line in run_command_with_output("simpleFoam"):
        r = re.search("Foam::error", line)
        if r is not None:
            return 1000
        r = re.search("ExecutionTime = (.+?) ", line)
        if r is not None:
            time = float(r.group(1))
    return time


def objective(trial):
    ap = trial.suggest_uniform("ap", 0.1, 1)
    au = trial.suggest_uniform("au", 0.1, 1)
    at = trial.suggest_uniform("at", 0.1, 1)
    return run(ap, au, at)


study = optuna.create_study()

optuna.logging.set_verbosity(optuna.logging.WARNING)

study.optimize(objective, n_trials=100)
print(study.best_params)

ここで、圧力の緩和係数を ap、速度の緩和係数を au、それ以外 (乱流統計量) の緩和係数を at としている。それぞれ探索範囲は 0.1〜1 としている。

実行

$ python run-op.py

結果

{'ap': 0.6937429666961807, 'au': 0.9501927541779454, 'at': 0.7819762799229192}

Bayesian Optimization

インストール

$ pip install bayesian-optimization

テスト

test-bo.py

from bayes_opt import BayesianOptimization

def objective(x):
    return -(x**2 - 2)**2

optimizer = BayesianOptimization(
        f=objective,
        pbounds={'x': (0, 2)},
        random_state=1
)

optimizer.maximize(init_points=3, n_iter=20)
print(optimizer.max)

緩和係数の最適化

run-bo.py

import sys
import subprocess as sp
import re
from bayes_opt import BayesianOptimization


def run_command_with_output(command):
    proc = sp.Popen(command,
            shell=True, stdout=sp.PIPE, stderr=sp.PIPE)

    while True:
        line = proc.stdout.readline()
        if not line:
            break
        line = line.decode().rstrip("\n")
        yield line

    if proc.wait() != 0:
        while True:
            line = proc.stderr.readline()
            if not line:
                break
            line = line.decode().rstrip("\n")
            yield line
        raise OSError("Command execution faied: %s" % command)


def write_fvSolution(ap, au, at):
    with open("system/fvSolution", "w") as f:
        f.write("""FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.1;
        smoother        GaussSeidel;
    }

    "(U|k|epsilon|omega|f|v2)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent      yes;

    residualControl
    {
        p               1e-2;
        U               1e-3;
        "(k|epsilon|omega|f|v2)" 1e-3;
    }
}

relaxationFactors
{
    fields
    {
""")

        f.write("        p %f;\n" % ap)

        f.write("""    }
    equations
    {
""")

        f.write("        U %f;\n" % au)
        f.write('        ".*" %f;\n' %at)

        f.write("""    }
}
""")


def run(ap, au, at):
    write_fvSolution(ap, au, at)
    time = 0
    for line in run_command_with_output("simpleFoam"):
        r = re.search("Foam::error", line)
        if r is not None:
            return 1000
        r = re.search("ExecutionTime = (.+?) ", line)
        if r is not None:
            time = float(r.group(1))
    return time


def objective(ap, au, at):
    return -run(ap, au, at)


optimizer = BayesianOptimization(
        f=objective,
        pbounds={"ap": (0.1, 1), "au": (0.1, 1), "at": (0.1, 1)},
        random_state=1
)

optimizer.maximize(init_points=10, n_iter=100)
print(optimizer.max)

結果

{'target': -5.57, 'params': {'ap': 0.9821591542366661, 'at': 0.8978768367274595, 'au': 0.9    428528721941154}}

比較

オリジナル (ap = 1、au = at = 0.9)

ExecutionTime = 5.97 s  ClockTime = 6 s
SIMPLE solution converged in 280 iterations

Optuna (ap = 0.694、au = 0.950、at = 0.782)

ExecutionTime = 5.2 s  ClockTime = 5 s
SIMPLE solution converged in 215 iterations

Bayesian Optimization (ap = 0.982、au = 0.943、at = 0.898)

ExecutionTime = 5.4 s  ClockTime = 5 s
SIMPLE solution converged in 217 iterations