ベイズ最適化のfluentへの適用

pip install ansys-fluent-core
pip install scikit-optimize

import ansys.fluent.core as pyfluent
import os
import glob
import re

from math import sqrt
import numpy as np
np.random.seed(237)
import matplotlib.pyplot as plt
from skopt.plots import plot_gaussian_process
from skopt.plots import plot_convergence
from skopt import gp_minimize

count = 0

################################################################
def run_fluent(x):
    global count
    count += 1
    solver = pyfluent.launch_fluent(
        precision="double",
        version="2d",
        processor_count=1,
        mode="solver",
    )
    tui = solver.tui
    trnfile = "log-{:02d}-{:f}.trn".format(count, x)
    solver.file.start_transcript(file_name=trnfile)
    solver.file.read_case(file_name="test.cas.h5")
    solver.file.read_data(file_name="test.dat.h5")

    solver.setup.materials.fluid['air'] = {
        "viscosity" : {
            "option" : "non-newtonian-power-law",
            "non_newtonian_power_law" : {
                "option" : "shear-rate-dependent",
                "maximum_viscosity" : 0.,
                "minimum_viscosity" : 0.,
                "power_law_index" : x,
                "consistency_index" : 100.
            }
        }
    }

    solver.solution.initialization.initialization_type = "standard"
    solver.solution.initialization.standard_initialize()
    solver.solution.run_calculation.iterate(iter_count=100)

    xyfile = "uy-{:02d}-{:f}.xy".format(count, x)
    tui.plot.plot("yes",xyfile,"yes","no","no","y-velocity","yes","1","0","outlet","()")
#   ident = "test-{:02d}-{:f}".format(count, x)
#   solver.file.write_case(file_name=ident+".cas.h5")
#   solver.file.write_data(file_name=ident+".dat.h5")
    solver.file.stop_transcript()


################################################################
def extract_data_from_file(filename):
    positions = []
    velocities = []

    with open(filename, 'r') as file:
        lines = file.readlines()

        # Regex to match data lines
        data_pattern = re.compile(r'^\s*([\d\.]+)\s+([\d\.Ee+-]+)\s*$')

        # Flag to indicate when we are inside the data section
        inside_data = False

        for line in lines:
            if '((xy/key/label' in line:
                inside_data = True
                continue
            if inside_data:
                match = data_pattern.match(line)
                if match:
                    position = float(match.group(1))
                    velocity = float(match.group(2))
                    positions.append(position)
                    velocities.append(velocity)
                elif line.strip() == ')':
                    break

    return velocities

################################################################
def cleanup_files():
    # 削除したいファイルのパターンリスト
    file_patterns = [
        'log-*.trn',
        'uy-*.xy',
        'test-*.h5'
    ]

    # 全てのパターンにマッチするファイルを格納するリスト
    all_files = []

    # 各パターンに対してファイルを取得し、all_filesリストに追加
    for pattern in file_patterns:
        all_files.extend(glob.glob(pattern))

    # 重複を除外（必要に応じて）
    all_files = list(set(all_files))

    # ファイルを一括削除
    for file_path in all_files:
        try:
            os.remove(file_path)
            print(f"{file_path} は正常に削除されました。")
        except FileNotFoundError:
            print(f"{file_path} が見つかりません。")
        except PermissionError:
            print(f"{file_path} の削除に必要な権限がありません。")
        except Exception as e:
            print(f"{file_path} の削除中にエラーが発生しました: {e}")


################################################################
def calculate_objective(vel):
    target = (  # <-- semi-circular profile (outlet)
    0.0,
    0.004178490302970204,
    0.006395205576424277,
    0.007981295678046266,
    0.00858,
    0.007981295678046266,
    0.006395205576424275,
    0.004178490302970204,
    0.0,
    )

    n = len(target)-2  # exclude two end points where uy=0
    err = 0
    for i in range(1,n+1):
        err += (target[i]-vel[i])**2
    return sqrt(err/n)

#############################################################################
def polyfunc(x):
    """
    approximated objective function. this function is used to plot only.
    as a result, the calculated minimum value is not
    exactly match with actual objective
    """
    if x[0] <= 0.5:
        val = \
        -2.494e-3 * x[0]**6 \
        +5.564e-2 * x[0]**5 \
        -8.284e-2 * x[0]**4 \
        +4.235e-2 * x[0]**3 \
        -5.643e-3 * x[0]**2 \
        -4.136e-4 * x[0]    \
        +5.241e-4
    else:
        val = \
        -1.313e-3 * x[0]**6 \
        +6.371e-3 * x[0]**5 \
        -1.230e-2 * x[0]**4 \
        +1.196e-2 * x[0]**3 \
        -6.446e-3 * x[0]**2 \
        +3.064e-3 * x[0]    \
        -1.028e-4
    return val


#############################################################################
def objective_func(x):
    run_fluent(x[0])
    xyfile = "uy-{:02d}-{:f}.xy".format(count, x[0])
    vel = extract_data_from_file(xyfile)
    err = calculate_objective(vel)
    print("calculate_objective has finished. count={} x={} err={}".format(count,x[0],err))
    return err


#############################################################################
#  main routine start
#############################################################################

noise_level = 0.0

cleanup_files()

n_calls = 10

res = gp_minimize(objective_func,      # the function to minimize
                  [(0.05, 0.95)],      # the bounds on each dimension of x
                  acq_func="LCB",      # the acquisition function
                  kappa=1.96,
                  n_calls=n_calls,     # the number of evaluations of f
                  n_initial_points=0,  # the number of random initialization points
                  x0=[[0.3], [0.8]],
                  random_state=11)     # the random seed


#############################################################################
# Accordingly, the approximated minimum is found to be:

"x^*=%.4f, f(x^*)=%.4f" % (res.x[0], res.fun)

#############################################################################
# For further inspection of the results, attributes of the `res` named tuple
# provide the following information:
#
# - `x` [float]: location of the minimum.
# - `fun` [float]: function value at the minimum.
# - `models`: surrogate models used for each iteration.
# - `x_iters` [array]:
#   location of function evaluation for each iteration.
# - `func_vals` [array]: function value for each iteration.
# - `space` [Space]: the optimization space.
# - `specs` [dict]: parameters passed to the function.

print(res)

with open('result.dat', 'w') as f:
    for key, value in res.items():
        f.write(f"{key}: {value}\n")

#############################################################################
# Together these attributes can be used to visually inspect the results of the
# minimization, such as the convergence trace or the acquisition function at
# the last iteration:

ax = plot_convergence(res);
ax.plot()
#plt.show()
plotfile = "convergence.png"
plt.savefig(plotfile)
plt.close()


#############################################################################

plt.rcParams["figure.figsize"] = (8, 8)

for n_iter in range(n_calls-2):
    # Plot objective function.
    plt.subplot(2, 1, 1)
    ax = plot_gaussian_process(res, n_calls=n_iter,
#                              objective=objective_func,
                               objective=polyfunc,
                               noise_level=noise_level,
                               show_legend=True, show_title=True,
                               show_mu=True,
                               show_observations=True,
                               show_next_point=False, show_acq_func=False)
    ax.set_ylabel("Objective")
    ax.set_ylim(0.0004, 0.0014)
    ax.set_xlim(0.0, 1.0)

    # Plot acquisition function.
    plt.subplot(2, 1, 2)
    ax = plot_gaussian_process(res, n_calls=n_iter,
                               show_legend=True, show_title=False,
                               show_mu=False, show_acq_func=True,
                               show_observations=False,
                               show_next_point=True)
    ax.set_ylabel("LCB")
    ax.set_xlabel("n-value")
    ax.set_ylim(-0.0012, -0.0000)
    ax.set_xlim(0.0, 1.0)
    ax.plot()
#   plt.show()
    plotfile = "plot-{:02d}.png".format(n_iter)
    plt.savefig(plotfile)
    plt.close()
    


#############################################################################
#  final result

plt.rcParams["figure.figsize"] = (6, 4)

ax = plot_gaussian_process(res,
#                         objective=objective_func,
                          objective=polyfunc,
                          noise_level=noise_level)

ax.plot()
#plt.show()
plotfile = "final.png"
plt.savefig(plotfile)
plt.close()