diff --git a/y3prediction/prediction.py b/y3prediction/prediction.py index 02fdb90..7dd9d08 100644 --- a/y3prediction/prediction.py +++ b/y3prediction/prediction.py @@ -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,)) @@ -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) @@ -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( @@ -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: ####") @@ -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 @@ -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, @@ -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( @@ -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) diff --git a/y3prediction/pyro_mechs/uiuc_const_gamma.py b/y3prediction/pyro_mechs/uiuc_const_gamma.py index ff29ad0..61aa604 100644 --- a/y3prediction/pyro_mechs/uiuc_const_gamma.py +++ b/y3prediction/pyro_mechs/uiuc_const_gamma.py @@ -3,6 +3,7 @@ """ +from warnings import warn import numpy as np @@ -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. @@ -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): @@ -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): @@ -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] ), ]) diff --git a/y3prediction/pyro_mechs/uiuc_sharp.py b/y3prediction/pyro_mechs/uiuc_sharp.py index d1d8728..c0bbd62 100644 --- a/y3prediction/pyro_mechs/uiuc_sharp.py +++ b/y3prediction/pyro_mechs/uiuc_sharp.py @@ -3,6 +3,7 @@ """ +from warnings import warn import numpy as np @@ -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. @@ -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): @@ -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): @@ -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] ), ]) diff --git a/y3prediction/unstart.py b/y3prediction/unstart.py index 4c2e497..aaa414d 100644 --- a/y3prediction/unstart.py +++ b/y3prediction/unstart.py @@ -111,12 +111,10 @@ def __call__(self, dcoll, x_vec, eos, *, time=0.0): mass = eos.get_density(self._pres_bulk, self._temp_bulk, self._y_bulk) + zeros velocity = self._vel_bulk + mom = mass*velocity energy = mass*(eos.get_internal_energy(self._temp_bulk, self._y_bulk) + 0.5*np.dot(velocity, velocity)) - mom = mass*velocity - energy = (energy + np.dot(mom, mom)/(2.0*mass)) - return make_conserved( dim=self._dim, mass=mass, @@ -128,20 +126,32 @@ def __call__(self, dcoll, x_vec, eos, *, time=0.0): def inlet_smoothing_func(self, x_vec, sigma): actx = x_vec[0].array_context - x0 = -0.013 - x1 = 0.013 - smth_bottom = smooth_step(actx, sigma*(x_vec[0] - x0)) - smth_top = smooth_step(actx, -sigma*(x_vec[0] - x1)) - return smth_bottom*smth_top + if self._dim == 2: + x0 = -0.013 + x1 = 0.013 + smth_bottom = smooth_step(actx, sigma*(x_vec[0] - x0)) + smth_top = smooth_step(actx, -sigma*(x_vec[0] - x1)) + return smth_bottom*smth_top + else: + r1 = 0.013 + radius = actx.np.sqrt((x_vec[1])**2 + (x_vec[2])**2) + smth_radial = smooth_step(actx, -sigma*(radius - r1)) + return smth_radial def outlet_smoothing_func(self, x_vec, sigma): actx = x_vec[0].array_context - x0 = -.2 - x1 = .2 - smth_bottom = smooth_step(actx, sigma*(x_vec[0] - x0)) - smth_top = smooth_step(actx, -sigma*(x_vec[0] - x1)) - return smth_bottom*smth_top + if self._dim == 2: + x0 = -.2 + x1 = .2 + smth_bottom = smooth_step(actx, sigma*(x_vec[0] - x0)) + smth_top = smooth_step(actx, -sigma*(x_vec[0] - x1)) + return smth_bottom*smth_top + else: + r1 = 0.2 + radius = actx.np.sqrt((x_vec[1])**2 + (x_vec[2])**2) + smth_radial = smooth_step(actx, -sigma*(radius - r1)) + return smth_radial def add_inlet(self, cv, pressure, temperature, x_vec, eos, *, time=0.0): """Create the solution state at locations *x_vec*. @@ -172,10 +182,14 @@ def add_inlet(self, cv, pressure, temperature, x_vec, eos, *, time=0.0): pres_inlet = self._pres_inlet # initial discontinuity location - y0 = -0.325 + if self._dim == 2: + y0 = -0.325 + dist = y0 - x_vec[1] + else: + x0 = -0.325 + dist = x0 - x_vec[0] # now solve for T, P, velocity - dist = y0 - x_vec[1] xtanh = self._disc_sigma*dist weight = 0.5*(1.0 - actx.np.tanh(xtanh)) pressure = pres_inlet + (pressure - pres_inlet)*weight @@ -185,7 +199,6 @@ def add_inlet(self, cv, pressure, temperature, x_vec, eos, *, time=0.0): # modify the temperature in the near wall region to match the # isothermal boundaries - sigma = self._temp_sigma if sigma > 0: wall_temperature = self._temp_wall @@ -244,10 +257,14 @@ def add_outlet(self, cv, pressure, temperature, x_vec, eos, *, time=0.0): pres_outlet = self._pres_outlet # initial discontinuity location - y0 = 0.825 + if self._dim == 2: + y0 = 0.825 + dist = x_vec[1] - y0 + else: + x0 = 1.1 + dist = x_vec[0] - x0 # now solve for T, P, velocity - dist = x_vec[1] - y0 xtanh = 50*dist weight = 0.5*(1.0 - actx.np.tanh(xtanh)) pressure = pres_outlet + (pressure - pres_outlet)*weight diff --git a/y3prediction/utils.py b/y3prediction/utils.py index b8a69b7..c8dc152 100644 --- a/y3prediction/utils.py +++ b/y3prediction/utils.py @@ -508,6 +508,7 @@ def __init__(self, *, dim, T0, P0, mass_frac, mach, gamma, self._x0 = smooth_x0 self._y0 = smooth_y0 self._z0 = smooth_z0 + self._r0 = smooth_r0 self._x1 = smooth_x1 self._y1 = smooth_y1 self._z1 = smooth_z1