Design Optimization of a Thrust Nozzle

theta_cone = math.rad($theta_cone) -- nominal straight-cone angle
theta_init = math.rad($theta_init) -- starting angle for thrust nozzle
alpha = math.rad($alpha)  -- angle for setting b2 in Bezier curve
beta = math.rad($beta)  -- angle for setting b3 in Bezier curve

#! /usr/bin/env python3
# optimize.py
# Automate the running of the flow simulations while searching
# for the optimum parameters (angles) that define the nozzle shape.
#
# To monitor the progress of the optimizer, you can run the command:
# tail -f progress.txt
# to see the result of each objective evaluation.
#
# PJ, 2018-03-04, take bits from nenzfr
#     2021-12-22, update to accommodate loads-file changes, use nelmin
#                 and run jobs in their own directories

import sys, os
DGDINST = os.path.expandvars("$HOME/dgdinst")
sys.path.append(DGDINST)

import shutil, shlex, subprocess, queue
import string, math
import time
from gdtk.numeric.nelmin import minimize, NelderMeadMinimizer

start_time = time.time()
progress_file = open("progress.txt", 'w')
progress_file.write("# job wall_clock params[0] params[1] params[2] params[3] thrustx\n")
progress_file.flush()
# Each of the Eilmer simulations will be run in its own directory, identified by job number,
# so that it can be run independently of all other simulations.
# We will use a queue to regulate the specification of this job number.
obj_eval_queue = queue.Queue()
obj_eval_queue.put(0)

def run_command(cmdText, jobDir, logFile=None):
    """
    Run the command as a subprocess in directory jobDir.
    If a logFile name is provided, capture stdout+stderr and write to that file,
    else just let those streams go to the console.
    """
    # Flush before using subprocess to ensure output
    # text is in the right order.
    sys.stdout.flush()
    if (type(cmdText) is list):
        args = cmdText
    else:
        args = shlex.split(cmdText)
    # print("jobDir:", jobDir, "cmdText:", " ".join(args)) # debug
    captureFlag = True if logFile else False
    result = subprocess.run(args, cwd=jobDir, check=True, capture_output=captureFlag)
    if captureFlag:
        f = open(logFile, 'w')
        f.write(result.stdout.decode())
        f.write(result.stderr.decode())
        f.close()
    return result.returncode

def prepare_input_script(substituteDict, jobDir):
    """
    Prepare the actual input file for Eilmer4 from a template
    which has most of the Lua input script in place and just
    a few place-holders that need to be substituted for actual
    values.
    """
    fp = open("nozzle.template.lua", 'r')
    text = fp.read()
    fp.close()
    template = string.Template(text)
    text = template.substitute(substituteDict)
    fp = open(jobDir+"/nozzle.lua", 'w')
    fp.write(text)
    fp.close()
    return

def run_simulation(param_dict, jobDir):
    """
    Prepare and run a simulation in it own directory.
    We do this so that several simulations may run concurrently
    and so that we can easily find the data files later.
    """
    if not os.path.exists(jobDir): os.mkdir(jobDir)
    shutil.copy('ideal-air.inp', jobDir)
    logFile = jobDir+'.log'
    run_command('prep-gas ideal-air.inp ideal-air-gas-model.lua', jobDir, logFile)
    prepare_input_script(param_dict, jobDir)
    run_command('e4shared --prep --job=nozzle', jobDir, logFile)
    run_command('e4shared --run --job=nozzle --max-cpus=4 --verbosity=1', jobDir, logFile)
    return

def post_simulation_files(jobDir):
    """
    Postprocess a simulation is its own directory.
    Although not used by the optimizer, this function may be handy
    when exploring the behaviour of the optimization procedure.
    """
    run_command('e4shared --post --job=nozzle --tindx-plot=all'+
                ' --vtk-xml --add-vars="mach,pitot,total-p,total-h"', jobDir)
    run_command('e4shared --post --job=nozzle --slice-list="1,0,:,0"'+
                ' --output-file="nozzle-throat.data"', jobDir)
    run_command('e4shared --post --job=nozzle --slice-list="1,$,:,0"'+
                ' --output-file="nozzle-exit.data"', jobDir)
    return

def neg_thrust(tindx, jobDir):
    """
    Read the loads file and return the x-component of thrust.

    Input:
    tindx : integer specifying which loads file to inspect
    jobDir : the dedicated directory for the simulation files.
    """
    fileName = jobDir+'/loads/t%04d/b0001.t%04d.loads.dat' % (tindx, tindx)
    print("Estimating thrust from loads file ", fileName) # debug
    f = open(fileName, 'r')
    thrustx = 0.0
    for line in f.readlines():
        items = line.strip().split()
        if items[0] == '#': continue
        if len(items) < 10: break
        dA = float(items[6]) # 7:area per radian for axisymmetric simulation
        nx = float(items[3]) # 4:n.x
        p = float(items[9]) # 10:p
        thrustx = thrustx + 2*math.pi*dA*nx*p
    print("thrustx=", thrustx) # debug
    return thrustx

def objective(params, *args):
    """
    Given a list of parameter values, run a simulation and compute the thrust.
    Since the thrust is in the negative x-direction, large negative values are good.
    The minimizer will drive toward good values.
    """
    global start_time # Immutable.
    global progress_file # Shared output stream, hopefully not too messed up.
    global obj_eval_queue # Thread-safe.
    job_number = obj_eval_queue.get()
    job_number += 1
    obj_eval_queue.put(job_number)
    print("Start job number:", job_number)
    pdict = {"theta_init":params[0], "alpha":params[1],
             "beta":params[2], "theta_cone":params[3]}
    jobDir = 'job-%04d' % job_number
    run_simulation(pdict, jobDir)
    # Note that, if we run an Eilmer simulation several times,
    # there may be several loads files, indexed in order of creation.
    # We want the most recent, if we are in the optimization process.
    f = open(jobDir+'/loads/nozzle-loads.times', 'r')
    tindx = 0
    for line in f.readlines():
        items = line.strip().split()
        if items[0] == '#': continue
        if len(items) < 2: break
        tindx = int(items[0])
    thrustx = neg_thrust(tindx, jobDir)
    progress_file.write("%d %.1f %.4f %.4f %.4f %.4f %.2f\n" %
                        (job_number, time.time()-start_time, params[0],
                         params[1], params[2], params[3], thrustx))
    progress_file.flush() # so that we can see the results as simulations run
    return thrustx

def main():
    """
    This script was built in stages.
    The if-statements are for testing the functions as the script
    was being developed. They might be still useful for exploring.
    """
    if 0:
        print("Let's run a simulation.")
        pdict = {"theta_init":30.0, "alpha":0.0, "beta":0.0, "theta_cone":30.0}
        run_simulation(pdict, "job-0000")
    if 0:
        print("Compute thrust from previously run simulation.")
        print("thrust=", neg_thrust(0, 'job-0000'))
    if 0:
        print("Evaluate objective function.")
        params = [30.0, 0.0, 0.0, 30.0] # [theta_init, alpha, beta, theta_cone]
        obj_eval_number = 1
        objv = objective(params)
        print("objective value=", objv)
    if 1:
        print("Let the optimizer take control and run the numerical experiment.")
        x0 = [30.0, 0.0, 0.0, 30.0] # [theta_init, alpha, beta, theta_cone] degrees
        result = minimize(objective, x0, [2.0, 2.0, 2.0, 2.0],
                          options={'tol':1.0e-4, 'P':2, 'maxfe':60, 'n_workers':2})
        print('optimized result:')
        print('  x=', result.x)
        print('  fx=', result.fun)
        print('  convergence-flag=', result.success)
        print('  number-of-fn-evaluations=', result.nfe)
        print('  number-of-restarts=', result.nrestarts)
        print('  vertices=', [str(v) for v in result.vertices])
    #
    print("overall calculation time:", time.time()-start_time)
    return

# Let's actually do some work...
main()
progress_file.close()

optimized result:
  x= [29.678627533594543, -6.232721263775233, 5.956211556232702, 18.49358284183389]
  fx= -2220.791439092127
  convergence-flag= False
  number-of-fn-evaluations= 71
  number-of-restarts= 0
  vertices= ['Vertex(x=[29.67862753 -6.23272126  5.95621156 18.49358284], f=-2220.791439092127)', 'Vertex(x=[29.27443431 -6.21183191  5.78719484 18.03676979], f=-2220.769113094416)', 'Vertex(x=[29.28534293 -6.17598784  6.12646338 17.88696286], f=-2220.4961864985025)', 'Vertex(x=[29.4983822  -6.31313774  6.90884094 18.14101969], f=-2220.414532010557)', 'Vertex(x=[29.48168364 -5.51494113  5.86079823 18.19524534], f=-2220.3968307587793)']
overall calculation time: 4692.7158970832825