Extending gas_out for a custom application

This script was written by Antony Cellier for a custom use case of gas_out in a cylindrical geometry.

First, two helper functions to convert an AVBP solution to a State are defined:

import copy as cp
import h5py as h5
import numpy as np
import hdfdict as hd
from ms_thermo.state import State
from ms_thermo.species import build_thermo_from_avbp


def load_mesh_and_solution(fname, mname):
    """Read HDF5 mesh and solution files"""

    print("Reading solution in ", fname)
    sol_h5 = h5.File(fname, 'r')

    sol = dict()
    for key1 in sol_h5.keys():
        sol_int=dict()
        for key2 in sol_h5[key1].keys():
            ARR_tot = np.array(sol_h5[key1][key2])
            sol_int[key2] = ARR_tot
        sol[key1] = sol_int

    print("Reading mesh in ", mname)
    mesh_h5 = h5.File(mname, 'r')

    mesh = dict()
    for key1 in mesh_h5.keys():
        mesh_int=dict()
        for key2 in mesh_h5[key1].keys():
            ARR_tot = np.array(mesh_h5[key1][key2])
            mesh_int[key2] = ARR_tot
        mesh[key1] = mesh_int

    return mesh, sol

def build_data_from_avbp(mesh, sol, species_db=None):
    """Build coordinates and state from AVBP"""

    state = State.from_cons(
        sol["GaseousPhase"]["rho"],
        sol["GaseousPhase"]["rhoE"],
        sol["RhoSpecies"], species_db=species_db)

    x_coor = mesh["Coordinates"]["x"]
    y_coor = mesh["Coordinates"]["y"]

    try:
        z_coor = mesh["Coordinates"]["z"]
        coor = np.stack((x_coor, y_coor, z_coor), axis=-1)
    except KeyError:
        coor = np.stack((x_coor, y_coor, np.zeros(x_coor.shape)), axis=-1)

    return coor, state

The ms-thermo State is created:

sname = 'path/to/species/db'
fname = 'path/to/solution'
mname = 'path/to/mesh'
species_db = build_thermo_from_avbp(sname)
MESH, INIT = ms_to.load_mesh_and_solution(fname, mname)
COOR, STATE = ms_to.build_data_from_avbp(MESH, INIT, species_db=species_db)

Then, a custom alteration mask adapted to the geometry is introduced:

def where_alter_cyl(COOR, axis, l_min, l_max, r_cyl):
    """Give the list of points to alter in the solution
    Avoid applying alteration to the entire domain"""

    ALTER = []

    for k in range(COOR.shape[0]):
        x = COOR[k, axis]
        d = [0, 1, 2]
        d.remove(axis)
        r = np.sqrt(COOR[k, d[0]]**2 + COOR[k, d[1]]**2)

        # Lower Cylinder
        if x >= l_min and x <= l_max and r < r_cyl :
            ALTER += [k]

    return ALTER

axis = ...
x_0_min = ...
x_0_max = ...
r_0 = ...
ALTER = where_alter_cyl(COOR, axis, x_0_min, x_0_max, r_0)

Temperature, pressure, species mass fractions are generated from Riemann solutions and altered with gas_out. For this case, velocity vectors are also generated and altered using the same mask:

T, P, Y, U = ...

def create_new_vectors(STATE, COOR, ALTER, T, P, Y, U):
    ''' Apply the local changes on the points to alter
        Build the array on the entire initial solution '''

    #Empty
    T_new = STATE.temperature
    P_new = STATE.pressure
    U_new = np.zeros((COOR.shape[0],))

    Y_new = STATE.mass_fracs

    #Fill
    for k,I in enumerate(ALTER):

        T_new[I] = T[k]
        P_new[I] = P[k]
        U_new[I] = U[k]

        for spec in Y:
            Y_new[spec][I] = Y[spec][k]

    return T_new, P_new, Y_new, U_new

T_new, P_new, Y_new, U_new = create_new_vectors(STATE, COOR, ALTER, T, P, Y, U)

The State is updated with the new values:

def alter_state_riemann(state, temp_new=None, press_new=None, y_new=None, verbose=False):
    ''' Fill the initial state with Riemann solutions'''

    log = str()

    if y_new is not None:

        state_species = list(state.mass_fracs.keys())

        for spec in y_new:
            if spec not in state_species:
                msgerr = "\n\n Species mismatch:"
                msgerr += "\n"+ str(spec) + " is not part of the mixture:\n"
                msgerr +=  "/".join(state_species)
                msgerr +=  "\nCheck your input file pretty please..."
                raise RuntimeWarning(msgerr)

        for spec in state.mass_fracs:
            if spec in y_new:
                state.mass_fracs[spec] = y_new[spec]

    if temp_new is not None:
        state.temperature = temp_new

    if press_new is not None:
        state.pressure = press_new

    if verbose:
        log += ("Filling of the Riemann solution DONE")

    return log, state

    _, STATE_FILL = alter_state(STATE, temp_new=T_new, press_new=P_new, y_new=Y_new)

Finally, the AVBP solution is updated with the new state. The kinetic energy from the velocity vector is added to the total energy.

def save_data_for_avbp(state, sol, fname, movement=False, U=None, V=None, W=None):
    """Update the full solution with the state parameters"""

    sol_new = cp.deepcopy(sol)
    sol_new["GaseousPhase"]["rho"] = state.rho

    if movement:
        if W is None:
            KIN = 1/2 * (np.power(U, 2) + np.power(V, 2))
        else:
            KIN = 1/2 * (np.power(U, 2) + np.power(V, 2) + np.power(W, 2))

        sol_new["GaseousPhase"]["rhoE"] = state.rho * (state.energy + KIN)
        sol_new["GaseousPhase"]["rhou"] = state.rho * U
        sol_new["GaseousPhase"]["rhov"] = state.rho * V

    else:
        sol_new["GaseousPhase"]["rhoE"] = state.rho * state.energy

    for spec in sol_new["RhoSpecies"]:
        sol_new["RhoSpecies"][spec] = state.rho * state.mass_fracs[spec]

    try:
        sol_new["Additionals"]["temperature"] = state.temperature
        sol_new["Additionals"]["pressure"] = state.pressure
    except KeyError:
        print("-Additionals- group is not part of the solution.")

    print("Saving solution in ", fname)
    with h5.File(fname, 'w') as fout:
        hd.dump(sol_new, fout)

V_new = np.zeros(U_new.shape)
save_data_for_avbp(STATE_FILL, INIT, fname, movement=True, U=U_new, V=V_new)