Skip to content

Commit

Permalink
Merge branch 'main' of github.com:illinois-ceesd/drivers_y3-prediction
Browse files Browse the repository at this point in the history
  • Loading branch information
@anderson2981
anderson2981 committed Aug 15, 2024
2 parents 732d9cb + 6d0eb3b commit 93f6419
Show file tree
Hide file tree
Showing 5 changed files with 223 additions and 128 deletions.
115 changes: 81 additions & 34 deletions y3prediction/prediction.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2449,7 +2449,10 @@ def get_data_from_hdf5(group):
# init params
disc_location = np.zeros(shape=(dim,))
fuel_location = np.zeros(shape=(dim,))
disc_location[1] = shock_loc_x
if dim == 2:
disc_location[1] = shock_loc_x
else:
disc_location[0] = shock_loc_x
fuel_location[1] = 10000.
plane_normal = np.zeros(shape=(dim,))

Expand Down Expand Up @@ -2518,9 +2521,14 @@ def get_data_from_hdf5(group):

sos = math.sqrt(inlet_gamma*pres_inflow/rho_inflow)

vel_inflow[1] = inlet_mach*sos
theta = 0.0
if dim == 2:
vel_inflow[1] = inlet_mach*sos
theta = np.pi/2.
else:
vel_inflow[0] = inlet_mach*sos

plane_normal = np.zeros(shape=(dim,))
theta = np.pi/2.
plane_normal[0] = np.cos(theta)
plane_normal[1] = np.sin(theta)
plane_normal = plane_normal/np.linalg.norm(plane_normal)
Expand All @@ -2532,7 +2540,8 @@ def get_data_from_hdf5(group):
print(f"\tinlet temperature {temp_inflow}")
print(f"\tinlet pressure {pres_inflow}")
print(f"\tinlet rho {rho_inflow}")
print(f"\tinlet velocity {vel_inflow[1]}")
print(f"\tinlet velocity x {vel_inflow[0]}")
print(f"\tinlet velocity y {vel_inflow[1]}")
#print(f"final inlet pressure {pres_inflow_final}")

bulk_init = PlanarDiscontinuityMulti(
Expand Down Expand Up @@ -2650,7 +2659,10 @@ def inlet_ramp_pressure(t):

sos = math.sqrt(inlet_gamma*pres_inflow/rho_inflow)

vel_inflow[1] = inlet_mach*sos
if dim == 2:
vel_inflow[1] = inlet_mach*sos
else:
vel_inflow[0] = inlet_mach*sos

if rank == 0:
print("#### Simluation initialization data: ####")
Expand All @@ -2661,7 +2673,8 @@ def inlet_ramp_pressure(t):
print(f"\tinlet pressure begin {inlet_ramp_beginP}")
print(f"\tinlet pressure end {inlet_ramp_endP}")
print(f"\tinlet rho {rho_inflow}")
print(f"\tinlet velocity {vel_inflow[1]}")
print(f"\tinlet velocity x {vel_inflow[0]}")
print(f"\tinlet velocity y {vel_inflow[1]}")
#print(f"final inlet pressure {pres_inflow_final}")

from y3prediction.unstart import InitUnstartRamp
Expand Down Expand Up @@ -5002,21 +5015,41 @@ def inflow_ramp_pressure(t):
p_fun=inflow_ramp_pressure)
else:
normal_dir[0] = 1
inflow_state = IsentropicInflow(
dim=dim,
temp_wall=temp_wall,
temp_sigma=temp_sigma,
vel_sigma=vel_sigma,
smooth_y0=-0.013,
smooth_y1=0.013,
normal_dir=normal_dir,
gamma=gamma,
nspecies=nspecies,
mass_frac=y,
T0=total_temp_inflow,
P0=ramp_beginP,
mach=inlet_mach,
p_fun=inflow_ramp_pressure)
if dim == 2:
inflow_state = IsentropicInflow(
dim=dim,
temp_wall=temp_wall,
temp_sigma=temp_sigma,
vel_sigma=vel_sigma,
smooth_y0=-0.013,
smooth_y1=0.013,
normal_dir=normal_dir,
gamma=gamma,
nspecies=nspecies,
mass_frac=y,
T0=total_temp_inflow,
P0=ramp_beginP,
mach=inlet_mach,
p_fun=inflow_ramp_pressure)
else:
smooth_r0 = fluid_nodes
smooth_r0[1] = actx.np.zeros_like(fluid_nodes[0])
smooth_r0[2] = actx.np.zeros_like(fluid_nodes[0])
inflow_state = IsentropicInflow(
dim=dim,
temp_wall=temp_wall,
temp_sigma=temp_sigma,
vel_sigma=vel_sigma,
smooth_r0=smooth_r0,
smooth_r1=0.013,
normal_dir=normal_dir,
gamma=gamma,
nspecies=nspecies,
mass_frac=y,
T0=total_temp_inflow,
P0=ramp_beginP,
mach=inlet_mach,
p_fun=inflow_ramp_pressure)
else:
inflow_state = IsentropicInflow(
dim=dim,
Expand Down Expand Up @@ -5162,7 +5195,7 @@ def _target_interface_state_func(**kwargs):
print(f"{bound_list=}")
check_bc_coverage(mesh=dcoll.discr_from_dd(dd_vol_fluid).mesh,
boundary_tags=bound_list,
incomplete_ok=False)
incomplete_ok=True)
except (ValueError, RuntimeError):
print(f"{uncoupled_fluid_boundaries=}")
raise SimulationConfigurationError(
Expand Down Expand Up @@ -5360,18 +5393,32 @@ def _sponge_sigma(sponge_field, x_vec):
return sponge_field

elif init_case == "unstart" or init_case == "unstart_ramp":
sponge_init_inlet = InitSponge(x0=inlet_sponge_x0,
thickness=inlet_sponge_thickness,
amplitude=sponge_amp,
direction=-2)
sponge_init_outlet = InitSponge(x0=outlet_sponge_x0,
thickness=outlet_sponge_thickness,
amplitude=sponge_amp,
direction=2)
sponge_init_top = InitSponge(x0=top_sponge_x0,
thickness=top_sponge_thickness,
amplitude=sponge_amp,
direction=1)
if dim == 2:
sponge_init_inlet = InitSponge(x0=inlet_sponge_x0,
thickness=inlet_sponge_thickness,
amplitude=sponge_amp,
direction=-2)
sponge_init_outlet = InitSponge(x0=outlet_sponge_x0,
thickness=outlet_sponge_thickness,
amplitude=sponge_amp,
direction=2)
sponge_init_top = InitSponge(x0=top_sponge_x0,
thickness=top_sponge_thickness,
amplitude=sponge_amp,
direction=1)
else:
sponge_init_inlet = InitSponge(x0=inlet_sponge_x0,
thickness=inlet_sponge_thickness,
amplitude=sponge_amp,
direction=-1)
sponge_init_outlet = InitSponge(x0=outlet_sponge_x0,
thickness=outlet_sponge_thickness,
amplitude=sponge_amp,
direction=1)
#sponge_init_top = InitSponge(x0=top_sponge_x0,
#thickness=top_sponge_thickness,
#amplitude=sponge_amp,
#direction=1)

def _sponge_sigma(sponge_field, x_vec):
sponge_field = sponge_init_outlet(sponge_field=sponge_field, x_vec=x_vec)
Expand Down
91 changes: 53 additions & 38 deletions y3prediction/pyro_mechs/uiuc_const_gamma.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
"""


from warnings import warn
import numpy as np


Expand Down Expand Up @@ -83,8 +84,22 @@ def __init__(self, usr_np=np):
self.species_names = ['C2H4', 'O2', 'CO2', 'CO', 'H2O', 'H2', 'N2']
self.species_indices = {'C2H4': 0, 'O2': 1, 'CO2': 2, 'CO': 3, 'H2O': 4, 'H2': 5, 'N2': 6}

self.wts = np.array([28.054, 31.998, 44.009, 28.009999999999998, 18.015, 2.016, 28.014])
self.iwts = 1/self.wts
self.molecular_weights = np.array([28.054, 31.998, 44.009, 28.009999999999998, 18.015, 2.016, 28.014])
self.inv_molecular_weights = 1/self.molecular_weights

@property
def wts(self):
warn("Thermochemistry.wts is deprecated and will go away in 2024. "
"Use molecular_weights instead.", DeprecationWarning, stacklevel=2)

return self.molecular_weights

@property
def iwts(self):
warn("Thermochemistry.iwts is deprecated and will go away in 2024. "
"Use inv_molecular_weights instead.", DeprecationWarning, stacklevel=2)

return self.inv_molecular_weights

def _pyro_zeros_like(self, argument):
# FIXME: This is imperfect, as a NaN will stay a NaN.
Expand Down Expand Up @@ -135,13 +150,13 @@ def species_index(self, species_name):

def get_specific_gas_constant(self, mass_fractions):
return self.gas_constant * (
+ self.iwts[0]*mass_fractions[0]
+ self.iwts[1]*mass_fractions[1]
+ self.iwts[2]*mass_fractions[2]
+ self.iwts[3]*mass_fractions[3]
+ self.iwts[4]*mass_fractions[4]
+ self.iwts[5]*mass_fractions[5]
+ self.iwts[6]*mass_fractions[6]
+ self.inv_molecular_weights[0]*mass_fractions[0]
+ self.inv_molecular_weights[1]*mass_fractions[1]
+ self.inv_molecular_weights[2]*mass_fractions[2]
+ self.inv_molecular_weights[3]*mass_fractions[3]
+ self.inv_molecular_weights[4]*mass_fractions[4]
+ self.inv_molecular_weights[5]*mass_fractions[5]
+ self.inv_molecular_weights[6]*mass_fractions[6]
)

def get_density(self, p, temperature, mass_fractions):
Expand All @@ -156,39 +171,39 @@ def get_pressure(self, rho, temperature, mass_fractions):

def get_mix_molecular_weight(self, mass_fractions):
return 1/(
+ self.iwts[0]*mass_fractions[0]
+ self.iwts[1]*mass_fractions[1]
+ self.iwts[2]*mass_fractions[2]
+ self.iwts[3]*mass_fractions[3]
+ self.iwts[4]*mass_fractions[4]
+ self.iwts[5]*mass_fractions[5]
+ self.iwts[6]*mass_fractions[6]
+ self.inv_molecular_weights[0]*mass_fractions[0]
+ self.inv_molecular_weights[1]*mass_fractions[1]
+ self.inv_molecular_weights[2]*mass_fractions[2]
+ self.inv_molecular_weights[3]*mass_fractions[3]
+ self.inv_molecular_weights[4]*mass_fractions[4]
+ self.inv_molecular_weights[5]*mass_fractions[5]
+ self.inv_molecular_weights[6]*mass_fractions[6]
)

def get_concentrations(self, rho, mass_fractions):
return self._pyro_make_array([
self.iwts[0] * rho * mass_fractions[0],
self.iwts[1] * rho * mass_fractions[1],
self.iwts[2] * rho * mass_fractions[2],
self.iwts[3] * rho * mass_fractions[3],
self.iwts[4] * rho * mass_fractions[4],
self.iwts[5] * rho * mass_fractions[5],
self.iwts[6] * rho * mass_fractions[6],
self.inv_molecular_weights[0] * rho * mass_fractions[0],
self.inv_molecular_weights[1] * rho * mass_fractions[1],
self.inv_molecular_weights[2] * rho * mass_fractions[2],
self.inv_molecular_weights[3] * rho * mass_fractions[3],
self.inv_molecular_weights[4] * rho * mass_fractions[4],
self.inv_molecular_weights[5] * rho * mass_fractions[5],
self.inv_molecular_weights[6] * rho * mass_fractions[6],
])

def get_mole_fractions(self, mix_mol_weight, mass_fractions):
return self._pyro_make_array([
self.iwts[0] * mass_fractions[0] * mix_mol_weight,
self.iwts[1] * mass_fractions[1] * mix_mol_weight,
self.iwts[2] * mass_fractions[2] * mix_mol_weight,
self.iwts[3] * mass_fractions[3] * mix_mol_weight,
self.iwts[4] * mass_fractions[4] * mix_mol_weight,
self.iwts[5] * mass_fractions[5] * mix_mol_weight,
self.iwts[6] * mass_fractions[6] * mix_mol_weight,
self.inv_molecular_weights[0] * mass_fractions[0] * mix_mol_weight,
self.inv_molecular_weights[1] * mass_fractions[1] * mix_mol_weight,
self.inv_molecular_weights[2] * mass_fractions[2] * mix_mol_weight,
self.inv_molecular_weights[3] * mass_fractions[3] * mix_mol_weight,
self.inv_molecular_weights[4] * mass_fractions[4] * mix_mol_weight,
self.inv_molecular_weights[5] * mass_fractions[5] * mix_mol_weight,
self.inv_molecular_weights[6] * mass_fractions[6] * mix_mol_weight,
])

def get_mass_average_property(self, mass_fractions, spec_property):
return sum([mass_fractions[i] * spec_property[i] * self.iwts[i]
return sum([mass_fractions[i] * spec_property[i] * self.inv_molecular_weights[i]
for i in range(self.num_species)])

def get_mixture_specific_heat_cp_mass(self, temperature, mass_fractions):
Expand Down Expand Up @@ -336,31 +351,31 @@ def get_species_mass_diffusivities_mixavg(self, pressure, temperature,
])
return self._pyro_make_array([
temp_pres*self.usr_np.where(self.usr_np.greater(denom[0], zeros),
(mmw - mole_fracs[0] * self.wts[0])/(mmw * denom[0]),
(mmw - mole_fracs[0] * self.molecular_weights[0])/(mmw * denom[0]),
bdiff_ij[0, 0]
),
temp_pres*self.usr_np.where(self.usr_np.greater(denom[1], zeros),
(mmw - mole_fracs[1] * self.wts[1])/(mmw * denom[1]),
(mmw - mole_fracs[1] * self.molecular_weights[1])/(mmw * denom[1]),
bdiff_ij[1, 1]
),
temp_pres*self.usr_np.where(self.usr_np.greater(denom[2], zeros),
(mmw - mole_fracs[2] * self.wts[2])/(mmw * denom[2]),
(mmw - mole_fracs[2] * self.molecular_weights[2])/(mmw * denom[2]),
bdiff_ij[2, 2]
),
temp_pres*self.usr_np.where(self.usr_np.greater(denom[3], zeros),
(mmw - mole_fracs[3] * self.wts[3])/(mmw * denom[3]),
(mmw - mole_fracs[3] * self.molecular_weights[3])/(mmw * denom[3]),
bdiff_ij[3, 3]
),
temp_pres*self.usr_np.where(self.usr_np.greater(denom[4], zeros),
(mmw - mole_fracs[4] * self.wts[4])/(mmw * denom[4]),
(mmw - mole_fracs[4] * self.molecular_weights[4])/(mmw * denom[4]),
bdiff_ij[4, 4]
),
temp_pres*self.usr_np.where(self.usr_np.greater(denom[5], zeros),
(mmw - mole_fracs[5] * self.wts[5])/(mmw * denom[5]),
(mmw - mole_fracs[5] * self.molecular_weights[5])/(mmw * denom[5]),
bdiff_ij[5, 5]
),
temp_pres*self.usr_np.where(self.usr_np.greater(denom[6], zeros),
(mmw - mole_fracs[6] * self.wts[6])/(mmw * denom[6]),
(mmw - mole_fracs[6] * self.molecular_weights[6])/(mmw * denom[6]),
bdiff_ij[6, 6]
),
])
Expand Down
Loading

0 comments on commit 93f6419

Please sign in to comment.