Skip to content

Commit

Permalink
updated unstart init for 3d
Browse files Browse the repository at this point in the history
  • Loading branch information
@anderson2981
anderson2981 committed Aug 15, 2024
1 parent 7fa2c66 commit cfbc029
Show file tree
Hide file tree
Showing 3 changed files with 116 additions and 49 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
49 changes: 34 additions & 15 deletions y3prediction/unstart.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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*.
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions y3prediction/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit cfbc029

Please sign in to comment.