From 1a18be82391ad4c18ed5b5eb4c5f9632180e5de7 Mon Sep 17 00:00:00 2001 From: anderson2981 <41346365+anderson2981@users.noreply.github.com> Date: Sun, 12 May 2024 11:52:31 -0500 Subject: [PATCH] Update limiter (#61) * updated to test fixes to meshmode and gmsh-interop * replaced nodal to modal transformation when finding element averages --- smoke_test_quads/run.sh | 2 +- y3prediction/prediction.py | 68 +++++++++++++------------------------- 2 files changed, 24 insertions(+), 46 deletions(-) diff --git a/smoke_test_quads/run.sh b/smoke_test_quads/run.sh index 5193c15..e24422b 100755 --- a/smoke_test_quads/run.sh +++ b/smoke_test_quads/run.sh @@ -1,3 +1,3 @@ #!/bin/bash #mpirun -n 2 python -u -O -m mpi4py driver.py -i run_params.yaml --log --casename=prediction-eager -mpirun -n 2 python -u -O -m mpi4py driver.py -i run_params.yaml --casename=prediction-eager +mpirun -n 2 python -u -m mpi4py driver.py -i run_params.yaml --casename=prediction-eager diff --git a/y3prediction/prediction.py b/y3prediction/prediction.py index 8325f43..59a55c7 100644 --- a/y3prediction/prediction.py +++ b/y3prediction/prediction.py @@ -490,58 +490,33 @@ def limit_fluid_state(dcoll, cv, temperature_seed, gas_model, dd): return make_obj_array([cv_lim, pressure, temperature]) -def element_average(dcoll, field, dd, - positivity_preserving=False): +def element_average(dcoll, dd, field, volumes=None): # Compute cell averages of the state actx = field.array_context + cell_avgs = op.elementwise_integral(dcoll, dd, field) + if volumes is None: + volumes = abs(op.elementwise_integral( + dcoll, dd, actx.zeros_like(field) + 1.0)) - def cancel_polynomials(grp): - return actx.from_numpy( - np.asarray([1 if sum(mode_id) == 0 - else 0 for mode_id in grp.mode_ids()])) + return cell_avgs/volumes - dd_nodal = dd - dd_modal = dd_nodal.with_discr_tag(DISCR_TAG_MODAL) - modal_map = dcoll.connection_from_dds(dd_nodal, dd_modal) - nodal_map = dcoll.connection_from_dds(dd_modal, dd_nodal) - - modal_discr = dcoll.discr_from_dd(dd_modal) - modal_field = modal_map(field) - - # cancel the ``high-order"" polynomials p > 0 and keep the average - filtered_modal_field = DOFArray( - actx, - tuple(actx.einsum("ej,j->ej", - vec_i, - cancel_polynomials(grp), - arg_names=("vec", "filter"), - tagged=(FirstAxisIsElementsTag(),)) - for grp, vec_i in zip(modal_discr.groups, modal_field)) - ) - - # convert back to nodal to have the average at all points - cell_avgs = nodal_map(filtered_modal_field) - - return cell_avgs - - -def _element_average_cv(dcoll, cv, dd): +def _element_average_cv(dcoll, dd, cv, volumes=None): nspecies = cv.nspecies dim = cv.dim - density = element_average(dcoll, cv.mass, dd) + density = element_average(dcoll, dd, cv.mass, volumes) momentum = make_obj_array([ - element_average(dcoll, cv.momentum[i], dd) + element_average(dcoll, dd, cv.momentum[i], volumes) for i in range(dim)]) - energy = element_average(dcoll, cv.energy, dd) + energy = element_average(dcoll, dd, cv.energy, volumes) species_mass = None if nspecies > 0: species_mass = make_obj_array([ - element_average(dcoll, cv.species_mass[i], dd) + element_average(dcoll, dd, cv.species_mass[i], volumes) for i in range(0, nspecies)]) # make a new CV with the limited variables @@ -654,6 +629,8 @@ def limit_fluid_state_lv(dcoll, cv, temperature_seed, gas_model, dd, nspecies = cv.nspecies dim = cv.dim toler = 1.e-13 + element_vols = abs(op.elementwise_integral(dcoll, dd, + actx.zeros_like(cv.mass) + 1.0)) print_stuff = False index = 118 @@ -686,7 +663,7 @@ def limit_fluid_state_lv(dcoll, cv, temperature_seed, gas_model, dd, #rho_lim = elem_avg_cv.mass*0.1 rho_lim = 1.e-6 - cell_avgs = element_average(dcoll, cv.mass, dd) + cell_avgs = element_average(dcoll, dd, cv.mass, volumes=element_vols) #cell_avgs = elem_avg_cv.mass #print(f"rho_avg {cell_avgs}") @@ -777,8 +754,9 @@ def limit_fluid_state_lv(dcoll, cv, temperature_seed, gas_model, dd, mmin = 0. #cell_avgs = elem_avg_cv.species_mass_fractions[i] - cell_avgs = element_average(dcoll, - cv_update_rho.species_mass_fractions[i], dd) + cell_avgs = element_average( + dcoll, dd, cv_update_rho.species_mass_fractions[i], + volumes=element_vols) cell_avgs = actx.np.where(actx.np.greater(cell_avgs, mmin), cell_avgs, mmin) @@ -934,7 +912,7 @@ def limit_fluid_state_lv(dcoll, cv, temperature_seed, gas_model, dd, #entropy = actx.np.log(pressure_updated/cv_updated.mass**1.4) #print(f"{entropy=}") - elem_avg_cv = _element_average_cv(dcoll, cv_updated, dd) + elem_avg_cv = _element_average_cv(dcoll, dd, cv_updated, volumes=element_vols) elem_avg_temp = gas_model.eos.temperature( cv=elem_avg_cv, temperature_seed=temperature_seed) elem_avg_pres = gas_model.eos.pressure( @@ -3388,7 +3366,7 @@ def get_mesh_data(): def get_mesh_data(): from meshmode.mesh.io import read_gmsh mesh_construction_kwargs = { - "force_positive_orientation": False, + "force_positive_orientation": True, "skip_element_orientation_test": True} mesh, tag_to_elements = read_gmsh( mesh_filename, force_ambient_dim=dim, @@ -5058,18 +5036,16 @@ def _target_interface_state_func(**kwargs): uncoupled_fluid_boundaries, bndry_mapping) # check the boundary condition coverage - #from meshmode.mesh import check_bc_coverage + from meshmode.mesh import check_bc_coverage try: bound_list = [] for bound in list(uncoupled_fluid_boundaries.keys()): bound_list.append(bound.tag) #print(f"{uncoupled_fluid_boundaries=}") print(f"{bound_list=}") - """ check_bc_coverage(mesh=dcoll.discr_from_dd(dd_vol_fluid).mesh, boundary_tags=bound_list, incomplete_ok=False) - """ except (ValueError, RuntimeError): print(f"{uncoupled_fluid_boundaries=}") raise SimulationConfigurationError( @@ -6772,7 +6748,7 @@ def unfiltered_rhs(t, state): if use_injection_source is True: fluid_rhs = fluid_rhs + \ injection_source(x_vec=fluid_nodes, cv=cv, - eos=gas_model.eos, time=t)/current_dt + eos=gas_model.eos, time=t) if use_ignition > 0: fluid_rhs = fluid_rhs + \ @@ -6972,6 +6948,8 @@ def my_rhs(t, state): # Work around long compile issue by computing and filtering RHS in separate # compiled functions rhs_state = unfiltered_rhs_compiled(t, state) + #rhs_data_precompute = precompute_rhs_compiled(t, state) + #rhs_state = compute_rhs_compiled(t, rhs_data_precompute) # Use a spectral filter on the RHS if use_rhs_filter: