From cfbc0293b477e37f527c73cd69d94d882591a81f Mon Sep 17 00:00:00 2001 From: Mike Anderson Date: Thu, 15 Aug 2024 07:05:01 -0700 Subject: [PATCH] updated unstart init for 3d --- y3prediction/prediction.py | 115 ++++++++++++++++++++++++++----------- y3prediction/unstart.py | 49 +++++++++++----- y3prediction/utils.py | 1 + 3 files changed, 116 insertions(+), 49 deletions(-) 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/unstart.py b/y3prediction/unstart.py index 7e5c765..aaa414d 100644 --- a/y3prediction/unstart.py +++ b/y3prediction/unstart.py @@ -126,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*. @@ -170,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 @@ -183,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 @@ -242,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