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)