diff --git a/y3prediction/actii_y3_cav5.py b/y3prediction/actii_y3_cav5.py index 3588fa2..16233ae 100644 --- a/y3prediction/actii_y3_cav5.py +++ b/y3prediction/actii_y3_cav5.py @@ -670,9 +670,10 @@ def add_injection(self, cv, pressure, temperature, x_vec, eos, *, time=0.0): xpos = x_vec[0] ypos = x_vec[1] + actx = xpos.array_context + zpos = actx.np.zeros_like(xpos) if self._dim == 3: zpos = x_vec[2] - actx = xpos.array_context zeros = actx.np.zeros_like(xpos) ones = zeros + 1.0 @@ -862,9 +863,10 @@ def add_injection_upstream(self, cv, pressure, temperature, xpos = x_vec[0] ypos = x_vec[1] + actx = xpos.array_context + zpos = actx.np.zeros_like(xpos) if self._dim == 3: zpos = x_vec[2] - actx = xpos.array_context zeros = actx.np.zeros_like(xpos) ones = zeros + 1.0 diff --git a/y3prediction/actii_y3_cav8.py b/y3prediction/actii_y3_cav8.py index 702ee20..0584469 100644 --- a/y3prediction/actii_y3_cav8.py +++ b/y3prediction/actii_y3_cav8.py @@ -131,11 +131,12 @@ def __call__(self, dcoll, x_vec, eos, *, time=0.0): xpos = x_vec[0] ypos = x_vec[1] + actx = xpos.array_context + zpos = actx.np.zeros_like(xpos) if self._dim == 3: zpos = x_vec[2] - ytop = 0*x_vec[0] - actx = xpos.array_context - zeros = 0*xpos + ytop = actx.zeros_like(xpos) + zeros = actx.zeros_like(xpos) ones = zeros + 1.0 mach = zeros @@ -670,9 +671,10 @@ def add_injection(self, cv, pressure, temperature, x_vec, eos, *, time=0.0): xpos = x_vec[0] ypos = x_vec[1] + actx = xpos.array_context + zpos = actx.np.zeros_like(xpos) if self._dim == 3: zpos = x_vec[2] - actx = xpos.array_context zeros = actx.np.zeros_like(xpos) ones = zeros + 1.0 @@ -862,9 +864,10 @@ def add_injection_upstream(self, cv, pressure, temperature, xpos = x_vec[0] ypos = x_vec[1] + actx = xpos.array_context + zpos = actx.np.zeros_like(xpos) if self._dim == 3: zpos = x_vec[2] - actx = xpos.array_context zeros = actx.np.zeros_like(xpos) ones = zeros + 1.0 @@ -1395,9 +1398,10 @@ def add_injection(self, cv, pressure, temperature, x_vec, eos, *, time=0.0): xpos = x_vec[0] ypos = x_vec[1] + actx = xpos.array_context + zpos = actx.np.zeros_like(xpos) if self._dim == 3: zpos = x_vec[2] - actx = xpos.array_context zeros = actx.np.zeros_like(xpos) ones = zeros + 1.0 @@ -1550,9 +1554,10 @@ def add_injection_upstream(self, cv, pressure, temperature, xpos = x_vec[0] ypos = x_vec[1] + actx = xpos.array_context + zpos = actx.np.zeros_like(xpos) if self._dim == 3: zpos = x_vec[2] - actx = xpos.array_context zeros = actx.np.zeros_like(xpos) ones = zeros + 1.0 diff --git a/y3prediction/prediction.py b/y3prediction/prediction.py index 694b1a9..b5fe839 100644 --- a/y3prediction/prediction.py +++ b/y3prediction/prediction.py @@ -1756,10 +1756,10 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): inv_flux_type = "Chandrashekar for single gas or passive species.\n" flux_msg = flux_msg + inv_flux_type else: - if inv_num_flux == "rusanov": - inviscid_numerical_flux_func = inviscid_facial_flux_rusanov - flux_msg = flux_msg + "Rusanov\n" - elif inv_num_flux == "hll": + inviscid_numerical_flux_func = inviscid_facial_flux_rusanov + flux_msg = flux_msg + "Rusanov\n" + + if inv_num_flux == "hll": inviscid_numerical_flux_func = inviscid_facial_flux_hll flux_msg = flux_msg + "HLL\n" @@ -1786,14 +1786,14 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): mw_ar = 39.948 univ_gas_const = 8314.59 - if gas_mat_prop == 0: - # working gas: O2/N2 # - # O2 mass fraction 0.273 - # gamma = 1.4 - # cp = 37.135 J/mol-K, - # rho= 1.977 kg/m^3 @298K - gamma = 1.4 - mw = mw_o2*mf_o2 + mw_n2*(1.0 - mf_o2) + # working gas: O2/N2 # + # O2 mass fraction 0.273 + # gamma = 1.4 + # cp = 37.135 J/mol-K, + # rho= 1.977 kg/m^3 @298K + gamma = 1.4 + mw = mw_o2*mf_o2 + mw_n2*(1.0 - mf_o2) + if gas_mat_prop == 1: # working gas: Ar # # O2 mass fraction 0.273 @@ -1818,11 +1818,11 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): Pr = 0.75 # viscosity @ 400C, Pa-s - if gas_mat_prop == 0: - # working gas: O2/N2 # - mu_o2 = 3.76e-5 - mu_n2 = 3.19e-5 - mu = mu_o2*mf_o2 + mu_n2*(1-mu_o2) # 3.3456e-5 + # working gas: O2/N2 # + mu_o2 = 3.76e-5 + mu_n2 = 3.19e-5 + mu = mu_o2*mf_o2 + mu_n2*(1-mu_o2) # 3.3456e-5 + if gas_mat_prop == 1: # working gas: Ar # mu_ar = 4.22e-5 @@ -1943,9 +1943,8 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): # set the species names if eos_type == 0: - if nspecies == 0: - species_names = ["inert"] - elif nspecies == 2: + species_names = ["inert"] + if nspecies == 2: species_names = ["air", "fuel"] elif nspecies == 3: species_names = ["air", "fuel", "inert"] @@ -2039,10 +2038,10 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): error_message = "Unknown transport_type {}".format(transport_type) raise RuntimeError(error_message) - if transport_type == 0: - physical_transport_model = SimpleTransport( - viscosity=mu, thermal_conductivity=kappa, - species_diffusivity=species_diffusivity) + physical_transport_model = SimpleTransport( + viscosity=mu, thermal_conductivity=kappa, + species_diffusivity=species_diffusivity) + if transport_type == 1: physical_transport_model = PowerLawTransport( alpha=transport_alpha, beta=transport_beta, @@ -2070,6 +2069,14 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): # with transport and eos sorted out, build the gas model gas_model = GasModel(eos=eos, transport=transport_model) + # quiescent initialization + bulk_init = Uniform( + dim=dim, + velocity=np.zeros(shape=(dim,)), + pressure=pres_bkrnd, + temperature=temp_bkrnd + ) + # select the initialization case if init_case == "shock1d": @@ -2899,9 +2906,8 @@ def injection_ramp_pressure(t): f"{vel_injection_upstream[1]}") print("#### Simluation initialization data: ####\n") - if actii_init_case == "cav8": - from y3prediction.actii_y3_cav8 import InitACTIIRamp - else: + from y3prediction.actii_y3_cav8 import InitACTIIRamp + if actii_init_case == "cav5": error_message = "Ramping init not fully implemented for cav5 config" bulk_init = InitACTIIRamp( @@ -3583,6 +3589,7 @@ def assign_fluid_boundaries(all_boundaries, bndry_mapping): = bndry_mapping[bndry_type] return all_boundaries + dd_vol_wall = None if use_wall: dd_vol_wall = DOFDesc(VolumeDomainTag("wall"), DISCR_TAG_BASE) wall_nodes = force_evaluation(actx, actx.thaw(dcoll.nodes(dd_vol_wall))) @@ -4060,6 +4067,7 @@ def filter_wall_rhs(rhs): logger.info("Positivity-preserving limiter enabled:") def my_limiter_func(cv, temperature_seed, gas_model, dd): + limiter_func = None if use_species_limiter == 1: limiter_func = limit_fluid_state( dcoll, cv, temperature_seed, gas_model, dd) @@ -4119,12 +4127,6 @@ def update_fluid_state(cv, dv, tv): from mirgecom.gas_model import ViscousFluidState return ViscousFluidState(cv, dv, tv) - def _create_wall_dependent_vars(wv): - return wall_model.dependent_vars(wv) - - create_wall_dependent_vars_compiled = actx.compile( - _create_wall_dependent_vars) - def get_temperature_update(cv, temperature): y = cv.species_mass_fractions e = gas_model.eos.internal_energy(cv)/cv.mass @@ -4133,6 +4135,98 @@ def get_temperature_update(cv, temperature): get_temperature_update_compiled = actx.compile(get_temperature_update) + # + # Setup the wall model + # + if use_wall: + def experimental_kappa(temperature): + return ( + 1.766e-10 * temperature**3 + - 4.828e-7 * temperature**2 + + 6.252e-4 * temperature + + 6.707e-3) + + def puma_kappa(mass_loss_frac): + return ( + 0.0988 * mass_loss_frac**2 + - 0.2751 * mass_loss_frac + + 0.201) + + def puma_effective_surface_area(mass_loss_frac): + # Original fit function: -1.1012e5*x**2 - 0.0646e5*x + 1.1794e5 + # Rescale by x==0 value and rearrange + return 1.1794e5 * ( + 1 + - 0.0547736137 * mass_loss_frac + - 0.9336950992 * mass_loss_frac**2) + + def _get_wall_kappa_fiber(mass, temperature): + mass_loss_frac = ( + (wall_insert_rho - mass)/wall_insert_rho + * wall_insert_mask) + scaled_insert_kappa = ( + experimental_kappa(temperature) + * puma_kappa(mass_loss_frac) + / puma_kappa(0)) + return ( + scaled_insert_kappa * wall_insert_mask + + wall_surround_kappa * wall_surround_mask) + + def _get_wall_kappa_inert(mass, temperature): + return ( + wall_insert_kappa * wall_insert_mask + + wall_surround_kappa * wall_surround_mask) + + def _get_wall_effective_surface_area_fiber(mass): + mass_loss_frac = ( + (wall_insert_rho - mass)/wall_insert_rho + * wall_insert_mask) + return ( + puma_effective_surface_area(mass_loss_frac) * wall_insert_mask) + + def _mass_loss_rate_fiber(mass, ox_mass, temperature, eff_surf_area): + actx = mass.array_context + alpha = ( + (0.00143+0.01*actx.np.exp(-1450.0/temperature)) + / (1.0+0.0002*actx.np.exp(13000.0/temperature))) + k = alpha*actx.np.sqrt( + (univ_gas_const*temperature)/(2.0*np.pi*mw_o2)) + return (mw_co/mw_o2 + mw_o/mw_o2 - 1)*ox_mass*k*eff_surf_area + + # inert + wall_model = WallModel( + heat_capacity=( + wall_insert_cp * wall_insert_mask + + wall_surround_cp * wall_surround_mask), + thermal_conductivity_func=_get_wall_kappa_inert) + + # non-porous + if wall_material == 1: + wall_model = WallModel( + heat_capacity=( + wall_insert_cp * wall_insert_mask + + wall_surround_cp * wall_surround_mask), + thermal_conductivity_func=_get_wall_kappa_fiber, + effective_surface_area_func=_get_wall_effective_surface_area_fiber, + mass_loss_func=_mass_loss_rate_fiber, + oxygen_diffusivity=wall_insert_ox_diff * wall_insert_mask) + # porous + elif wall_material == 2: + wall_model = WallModel( + heat_capacity=( + wall_insert_cp * wall_insert_mask + + wall_surround_cp * wall_surround_mask), + thermal_conductivity_func=_get_wall_kappa_fiber, + effective_surface_area_func=_get_wall_effective_surface_area_fiber, + mass_loss_func=_mass_loss_rate_fiber, + oxygen_diffusivity=wall_insert_ox_diff * wall_insert_mask) + + def _create_wall_dependent_vars(wv): + return wall_model.dependent_vars(wv) + + create_wall_dependent_vars_compiled = actx.compile( + _create_wall_dependent_vars) + if rank == 0: logger.info("Smoothness functions processing") @@ -4751,6 +4845,8 @@ def get_production_rates(cv, temperature): if rank == 0: logger.info("More gradient processing") + target_boundaries = {} + def grad_cv_operator_target(fluid_state, time): return grad_cv_operator(dcoll=dcoll, gas_model=gas_model, dd=dd_vol_fluid, @@ -4779,7 +4875,6 @@ def grad_t_operator_target(fluid_state, time): target_bndry_mapping["prescribed"] = DummyBoundary() target_bndry_mapping["isentropic_pressure_ramp"] = DummyBoundary() - target_boundaries = {} target_boundaries = assign_fluid_boundaries( target_boundaries, target_bndry_mapping) @@ -4811,92 +4906,6 @@ def grad_t_operator_target(fluid_state, time): smoothness_kappa=target_av_skappa, smoothness_d=target_av_sd) - # - # Setup the wall model - # - if use_wall: - def experimental_kappa(temperature): - return ( - 1.766e-10 * temperature**3 - - 4.828e-7 * temperature**2 - + 6.252e-4 * temperature - + 6.707e-3) - - def puma_kappa(mass_loss_frac): - return ( - 0.0988 * mass_loss_frac**2 - - 0.2751 * mass_loss_frac - + 0.201) - - def puma_effective_surface_area(mass_loss_frac): - # Original fit function: -1.1012e5*x**2 - 0.0646e5*x + 1.1794e5 - # Rescale by x==0 value and rearrange - return 1.1794e5 * ( - 1 - - 0.0547736137 * mass_loss_frac - - 0.9336950992 * mass_loss_frac**2) - - def _get_wall_kappa_fiber(mass, temperature): - mass_loss_frac = ( - (wall_insert_rho - mass)/wall_insert_rho - * wall_insert_mask) - scaled_insert_kappa = ( - experimental_kappa(temperature) - * puma_kappa(mass_loss_frac) - / puma_kappa(0)) - return ( - scaled_insert_kappa * wall_insert_mask - + wall_surround_kappa * wall_surround_mask) - - def _get_wall_kappa_inert(mass, temperature): - return ( - wall_insert_kappa * wall_insert_mask - + wall_surround_kappa * wall_surround_mask) - - def _get_wall_effective_surface_area_fiber(mass): - mass_loss_frac = ( - (wall_insert_rho - mass)/wall_insert_rho - * wall_insert_mask) - return ( - puma_effective_surface_area(mass_loss_frac) * wall_insert_mask) - - def _mass_loss_rate_fiber(mass, ox_mass, temperature, eff_surf_area): - actx = mass.array_context - alpha = ( - (0.00143+0.01*actx.np.exp(-1450.0/temperature)) - / (1.0+0.0002*actx.np.exp(13000.0/temperature))) - k = alpha*actx.np.sqrt( - (univ_gas_const*temperature)/(2.0*np.pi*mw_o2)) - return (mw_co/mw_o2 + mw_o/mw_o2 - 1)*ox_mass*k*eff_surf_area - - # inert - if wall_material == 0: - wall_model = WallModel( - heat_capacity=( - wall_insert_cp * wall_insert_mask - + wall_surround_cp * wall_surround_mask), - thermal_conductivity_func=_get_wall_kappa_inert) - # non-porous - elif wall_material == 1: - wall_model = WallModel( - heat_capacity=( - wall_insert_cp * wall_insert_mask - + wall_surround_cp * wall_surround_mask), - thermal_conductivity_func=_get_wall_kappa_fiber, - effective_surface_area_func=_get_wall_effective_surface_area_fiber, - mass_loss_func=_mass_loss_rate_fiber, - oxygen_diffusivity=wall_insert_ox_diff * wall_insert_mask) - # porous - elif wall_material == 2: - wall_model = WallModel( - heat_capacity=( - wall_insert_cp * wall_insert_mask - + wall_surround_cp * wall_surround_mask), - thermal_conductivity_func=_get_wall_kappa_fiber, - effective_surface_area_func=_get_wall_effective_surface_area_fiber, - mass_loss_func=_mass_loss_rate_fiber, - oxygen_diffusivity=wall_insert_ox_diff * wall_insert_mask) - ################################## # Set up the boundary conditions # ################################## @@ -5406,6 +5415,7 @@ def _sponge_source(sigma, cv, sponge_cv): logger.info("Viz & utilities processsing") # avoid making a second discretization if viz_order == order + wall_visualizer = None if viz_order == order: fluid_visualizer = make_visualizer(dcoll, volume_dd=dd_vol_fluid) if use_wall: @@ -6010,6 +6020,7 @@ def report_violators(ary, data_min, data_max): if np.any(mask): guilty_node_x = nodes_x[mask] guilty_node_y = nodes_y[mask] + guilty_node_z = None if dim == 3: guilty_node_z = nodes_z[mask] guilty_data = data[mask] @@ -6342,7 +6353,13 @@ def my_pre_step(step, t, dt, state): do_status = check_step(step=step, interval=nstatus) next_dump_number = step + dv = None + ts_field_fluid = None + ts_field_wall = None + cfl_fluid = 0. + cfl_wall = cfl_fluid if any([do_viz, do_restart, do_health, do_status]): + wv = None if not force_eval: fluid_state = force_evaluation(actx, fluid_state) #state = force_evaluation(actx, state) @@ -6358,15 +6375,12 @@ def my_pre_step(step, t, dt, state): t=t, dt=dt, cfl=current_cfl, t_final=t_final, constant_cfl=constant_cfl, fluid_dd=dd_vol_fluid) - ts_field_wall = None if use_wall: ts_field_wall, cfl_wall, dt_wall = my_get_timestep_wall( dcoll=dcoll, wv=wv, wall_kappa=wdv.thermal_conductivity, wall_temperature=wdv.temperature, t=t, dt=dt, cfl=current_cfl, t_final=t_final, constant_cfl=constant_cfl, wall_dd=dd_vol_wall) - else: - cfl_wall = cfl_fluid """ # adjust time for constant cfl, use the smallest timescale @@ -6564,6 +6578,7 @@ def boundary_flux(bdtag, bdry): off_axis_x = 1e-7 fluid_nodes_are_off_axis = actx.np.greater(fluid_nodes[0], off_axis_x) + wall_nodes_are_off_axis = None if use_wall: wall_nodes_are_off_axis = actx.np.greater(wall_nodes[0], off_axis_x) @@ -6766,6 +6781,10 @@ def unfiltered_rhs(t, state): dd=dd_vol_fluid, operator_states_quad=fluid_operator_states_quad, time=t, quadrature_tag=quadrature_tag) + smoothness_mu = actx.np.zeros_like(cv.mass) + smoothness_beta = actx.np.zeros_like(cv.mass) + smoothness_kappa = actx.np.zeros_like(cv.mass) + smoothness_d = actx.np.zeros_like(cv.mass) if use_av == 1: smoothness_mu = compute_smoothness( cv=cv, dv=fluid_state.dv, grad_cv=grad_fluid_cv) diff --git a/y3prediction/utils.py b/y3prediction/utils.py index 6b1cc91..b8a69b7 100644 --- a/y3prediction/utils.py +++ b/y3prediction/utils.py @@ -438,8 +438,7 @@ def __call__(self, sponge_field, x_vec, *, time=0.0): zeros = actx.np.zeros_like(xpos) x0 = zeros + self._x0 - if abs(self._direction) == 1: - coords = xpos + coords = xpos if abs(self._direction) == 2: coords = ypos