From bb960f6ac4502a6dd900742947944a13a641a284 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Fri, 10 Dec 2021 10:31:23 -0600 Subject: [PATCH 01/20] Remove .join() calls in Euler operator --- mirgecom/euler.py | 45 +++++++++++++++++++++++++++++++-------------- 1 file changed, 31 insertions(+), 14 deletions(-) diff --git a/mirgecom/euler.py b/mirgecom/euler.py index 26886c259..5fabccefa 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -57,11 +57,11 @@ inviscid_flux, inviscid_facial_flux ) -from grudge.eager import ( - interior_trace_pair, + +from grudge.trace_pair import ( + local_interior_trace_pair, cross_rank_trace_pairs ) -from grudge.trace_pair import TracePair from mirgecom.fluid import make_conserved from mirgecom.operators import div_operator @@ -100,20 +100,37 @@ def euler_operator(discr, eos, boundaries, cv, time=0.0): Agglomerated object array of DOF arrays representing the RHS of the Euler flow equations. """ + # Compute volume contributions inviscid_flux_vol = inviscid_flux(discr, eos, cv) + + # Compute interface contributions inviscid_flux_bnd = ( - inviscid_facial_flux(discr, eos=eos, cv_tpair=interior_trace_pair(discr, cv)) - + sum(inviscid_facial_flux( - discr, eos=eos, cv_tpair=TracePair( - part_tpair.dd, interior=make_conserved(discr.dim, q=part_tpair.int), - exterior=make_conserved(discr.dim, q=part_tpair.ext))) - for part_tpair in cross_rank_trace_pairs(discr, cv.join())) - + sum(boundaries[btag].inviscid_boundary_flux(discr, btag=btag, cv=cv, - eos=eos, time=time) - for btag in boundaries) + # Rank-local contributions + inviscid_facial_flux( + discr, + eos=eos, + cv_tpair=local_interior_trace_pair(discr, cv) + ) + # Cross-rank contributions + + sum( + inviscid_facial_flux( + discr, + eos=eos, + cv_tpair=part_tpair + ) for part_tpair in cross_rank_trace_pairs(discr, cv) + ) + # Boundary condition contributions + + sum( + boundaries[btag].inviscid_boundary_flux( + discr, + btag=btag, + cv=cv, + eos=eos, + time=time + ) for btag in boundaries + ) ) - q = -div_operator(discr, inviscid_flux_vol.join(), inviscid_flux_bnd.join()) - return make_conserved(discr.dim, q=q) + return -div_operator(discr, inviscid_flux_vol, inviscid_flux_bnd) def inviscid_operator(discr, eos, boundaries, q, t=0.0): From 0450f5c7c3422e30cfdbab91fb971f9fa93f8bf9 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Fri, 10 Dec 2021 11:54:55 -0600 Subject: [PATCH 02/20] Compute componentwise norms for array containers --- mirgecom/simutil.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index 7ac50ccdd..e0cbb22c3 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -48,6 +48,13 @@ import numpy as np import grudge.op as op +from arraycontext import map_array_container, flatten + +from functools import partial + +from meshmode.dof_array import DOFArray + + logger = logging.getLogger(__name__) @@ -290,7 +297,18 @@ def compare_fluid_solutions(discr, red_state, blue_state): """ actx = red_state.array_context resid = red_state - blue_state - return [actx.to_numpy(discr.norm(v, np.inf)) for v in resid.join()] + resid_errs = actx.to_numpy( + flatten(componentwise_norms(discr, resid, order=np.inf), actx)) + + return resid_errs.tolist() + + +def componentwise_norms(discr, fields, order=np.inf): + if not isinstance(fields, DOFArray): + return map_array_container( + partial(componentwise_norms, discr, order=order), fields) + + return discr.norm(fields, order) def generate_and_distribute_mesh(comm, generate_mesh): From a29209d471f862d892064c802404e664954530f4 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Fri, 10 Dec 2021 12:02:48 -0600 Subject: [PATCH 03/20] Document componentwise norms --- mirgecom/simutil.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index e0cbb22c3..7696278d3 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -12,6 +12,7 @@ -------------------- .. autofunction:: compare_fluid_solutions +.. autofunction:: componentwise_norms .. autofunction:: check_naninf_local .. autofunction:: check_range_local @@ -304,6 +305,11 @@ def compare_fluid_solutions(discr, red_state, blue_state): def componentwise_norms(discr, fields, order=np.inf): + """Return the *order*-norm for each component of *fields*. + + .. note:: + This is a collective routine and must be called by all MPI ranks. + """ if not isinstance(fields, DOFArray): return map_array_container( partial(componentwise_norms, discr, order=order), fields) From 72948f28dea51a52bbce028fa219535397afc353 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Fri, 10 Dec 2021 12:51:35 -0600 Subject: [PATCH 04/20] Start updating unit tests --- test/test_euler.py | 12 ++++++------ test/test_simutil.py | 7 +++++-- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/test/test_euler.py b/test/test_euler.py index 821e9b43e..c0a22d978 100644 --- a/test/test_euler.py +++ b/test/test_euler.py @@ -234,7 +234,7 @@ def test_vortex_rhs(actx_factory, order): discr, eos=IdealSingleGas(), boundaries=boundaries, cv=vortex_soln, time=0.0) - err_max = actx.to_numpy(discr.norm(inviscid_rhs.join(), np.inf)) + err_max = actx.to_numpy(discr.norm(inviscid_rhs, np.inf)) eoc_rec.add_data_point(1.0 / nel_1d, err_max) logger.info( @@ -294,7 +294,7 @@ def test_lump_rhs(actx_factory, dim, order): expected_rhs = lump.exact_rhs(discr, cv=lump_soln, time=0) err_max = actx.to_numpy( - discr.norm((inviscid_rhs-expected_rhs).join(), np.inf)) + discr.norm((inviscid_rhs-expected_rhs), np.inf)) if err_max > maxxerr: maxxerr = err_max @@ -371,7 +371,7 @@ def test_multilump_rhs(actx_factory, dim, order, v0): print(f"expected_rhs = {expected_rhs}") err_max = actx.to_numpy( - discr.norm((inviscid_rhs-expected_rhs).join(), np.inf)) + discr.norm((inviscid_rhs-expected_rhs), np.inf)) if err_max > maxxerr: maxxerr = err_max @@ -437,7 +437,7 @@ def _euler_flow_stepper(actx, parameters): def write_soln(state, write_status=True): dv = eos.dependent_vars(cv=state) expected_result = initializer(nodes, t=t) - result_resid = (state - expected_result).join() + result_resid = state - expected_result maxerr = [np.max(np.abs(result_resid[i].get())) for i in range(dim + 2)] mindv = [np.min(dvfld.get()) for dvfld in dv] maxdv = [np.max(dvfld.get()) for dvfld in dv] @@ -493,7 +493,7 @@ def rhs(t, q): cv = rk4_step(cv, t, dt, rhs) cv = make_conserved( - dim, q=filter_modally(discr, "vol", cutoff, frfunc, cv.join()) + dim, q=filter_modally(discr, "vol", cutoff, frfunc, cv) ) t += dt @@ -506,7 +506,7 @@ def rhs(t, q): maxerr = max(write_soln(cv, False)) else: expected_result = initializer(nodes, time=t) - maxerr = actx.to_numpy(discr.norm((cv - expected_result).join(), np.inf)) + maxerr = actx.to_numpy(discr.norm(cv - expected_result, np.inf)) logger.info(f"Max Error: {maxerr}") if maxerr > exittol: diff --git a/test/test_simutil.py b/test/test_simutil.py index a590abcd9..a80cb7fdb 100644 --- a/test/test_simutil.py +++ b/test/test_simutil.py @@ -29,6 +29,7 @@ from arraycontext import ( # noqa thaw, + flatten, pytest_generate_tests_for_pyopencl_array_context as pytest_generate_tests ) @@ -108,7 +109,7 @@ def test_basic_cfd_healthcheck(actx_factory): def test_analytic_comparison(actx_factory): """Quick test of state comparison routine.""" from mirgecom.initializers import Vortex2D - from mirgecom.simutil import compare_fluid_solutions + from mirgecom.simutil import compare_fluid_solutions, componentwise_norms actx = actx_factory() nel_1d = 4 @@ -134,7 +135,9 @@ def test_analytic_comparison(actx_factory): cv = make_conserved(dim, mass=mass, energy=energy, momentum=mom) resid = vortex_soln - cv - expected_errors = [actx.to_numpy(discr.norm(v, np.inf)) for v in resid.join()] + + expected_errors = actx.to_numpy( + flatten(componentwise_norms(discr, resid, order=np.inf), actx)).tolist() errors = compare_fluid_solutions(discr, cv, cv) assert max(errors) == 0 From b9e714a01e2a0edae0b3753e83ea392282eee0b9 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Sun, 12 Dec 2021 13:45:49 -0600 Subject: [PATCH 05/20] Purge join from tests --- mirgecom/filter.py | 8 +++++++- test/test_filter.py | 5 +++-- test/test_fluid.py | 13 +++++-------- test/test_lazy.py | 18 ++++++------------ test/test_operators.py | 29 +++++++++++------------------ test/test_restart.py | 2 +- test/test_viscous.py | 17 +++++++++-------- 7 files changed, 42 insertions(+), 50 deletions(-) diff --git a/mirgecom/filter.py b/mirgecom/filter.py index 2737970b9..1b31a199a 100644 --- a/mirgecom/filter.py +++ b/mirgecom/filter.py @@ -201,7 +201,13 @@ def filter_modally(dcoll, dd, cutoff, mode_resp_func, field): dd_modal = dof_desc.DD_VOLUME_MODAL discr = dcoll.discr_from_dd(dd) - assert isinstance(field, DOFArray) + from arraycontext import map_array_container + from functools import partial + if not isinstance(field, DOFArray): + return map_array_container( + partial(filter_modally, dcoll, dd, cutoff, mode_resp_func), field + ) + actx = field.array_context modal_map = dcoll.connection_from_dds(dd, dd_modal) diff --git a/test/test_filter.py b/test/test_filter.py index 85814fb9c..2bbee78f5 100644 --- a/test/test_filter.py +++ b/test/test_filter.py @@ -186,13 +186,14 @@ def test_filter_function(actx_factory, dim, order, do_viz=False): # the filter unharmed. from mirgecom.initializers import Uniform initr = Uniform(dim=dim) - uniform_soln = initr(t=0, x_vec=nodes).join() + uniform_soln = initr(t=0, x_vec=nodes) from mirgecom.filter import filter_modally filtered_soln = filter_modally(discr, "vol", cutoff, frfunc, uniform_soln) soln_resid = uniform_soln - filtered_soln - max_errors = [actx.to_numpy(discr.norm(v, np.inf)) for v in soln_resid] + from mirgecom.simutil import componentwise_norms + max_errors = componentwise_norms(discr, soln_resid, np.inf) tol = 1e-14 diff --git a/test/test_fluid.py b/test/test_fluid.py index 9cac5d7aa..4def2ff25 100644 --- a/test/test_fluid.py +++ b/test/test_fluid.py @@ -77,8 +77,7 @@ def test_velocity_gradient_sanity(actx_factory, dim, mass_exp, vel_fac): cv = make_conserved(dim, mass=mass, energy=energy, momentum=mom) from grudge.op import local_grad - grad_cv = make_conserved(dim, - q=local_grad(discr, cv.join())) + grad_cv = local_grad(discr, cv) grad_v = velocity_gradient(discr, cv, grad_cv) @@ -123,8 +122,7 @@ def test_velocity_gradient_eoc(actx_factory, dim): cv = make_conserved(dim, mass=mass, energy=energy, momentum=mom) from grudge.op import local_grad - grad_cv = make_conserved(dim, - q=local_grad(discr, cv.join())) + grad_cv = local_grad(discr, cv) grad_v = velocity_gradient(discr, cv, grad_cv) def exact_grad_row(xdata, gdim, dim): @@ -178,8 +176,7 @@ def test_velocity_gradient_structure(actx_factory): cv = make_conserved(dim, mass=mass, energy=energy, momentum=mom) from grudge.op import local_grad - grad_cv = make_conserved(dim, - q=local_grad(discr, cv.join())) + grad_cv = local_grad(discr, cv) grad_v = velocity_gradient(discr, cv, grad_cv) tol = 1e-11 @@ -233,8 +230,8 @@ def test_species_mass_gradient(actx_factory, dim): cv = make_conserved(dim, mass=mass, energy=energy, momentum=mom, species_mass=species_mass) from grudge.op import local_grad - grad_cv = make_conserved(dim, - q=local_grad(discr, cv.join())) + grad_cv = local_grad(discr, cv) + from mirgecom.fluid import species_mass_fraction_gradient grad_y = species_mass_fraction_gradient(discr, cv, grad_cv) diff --git a/test/test_lazy.py b/test/test_lazy.py index 75e645e0f..11ba007ff 100644 --- a/test/test_lazy.py +++ b/test/test_lazy.py @@ -22,7 +22,7 @@ import numpy as np from functools import partial -from pytools.obj_array import make_obj_array, obj_array_vectorize +from pytools.obj_array import make_obj_array, obj_array_vectorize # noqa import pyopencl as cl import pyopencl.tools as cl_tools import pyopencl.array as cla # noqa @@ -75,19 +75,13 @@ def get_discr(order): # Mimics math.isclose for state arrays def _isclose(discr, x, y, rel_tol=1e-9, abs_tol=0, return_operands=False): - def componentwise_norm(a): - from mirgecom.fluid import ConservedVars - if isinstance(a, ConservedVars): - return componentwise_norm(a.join()) - from arraycontext import get_container_context_recursively - actx = get_container_context_recursively(a) - return obj_array_vectorize(lambda b: actx.to_numpy(discr.norm(b, np.inf)), a) - - lhs = componentwise_norm(x - y) + + from mirgecom.simutil import componentwise_norms + lhs = componentwise_norms(discr, x - y, np.inf) rhs = np.maximum( rel_tol * np.maximum( - componentwise_norm(x), - componentwise_norm(y)), + componentwise_norms(discr, x, np.inf), + componentwise_norms(discr, y, np.inf)), abs_tol) is_close = np.all(lhs <= rhs) diff --git a/test/test_operators.py b/test/test_operators.py index d60a8ec80..761d32f97 100644 --- a/test/test_operators.py +++ b/test/test_operators.py @@ -217,13 +217,12 @@ def sym_eval(expr, x_vec): test_data = test_func(nodes) exact_grad = grad_test_func(nodes) - def inf_norm(x): - return actx.to_numpy(discr.norm(x, np.inf)) + from mirgecom.simutil import componentwise_norms + from arraycontext import flatten + + err_scale = max(flatten(componentwise_norms(discr, exact_grad, np.inf), + actx)) - if isinstance(test_data, ConservedVars): - err_scale = inf_norm(exact_grad.join()) - else: - err_scale = inf_norm(exact_grad) if err_scale <= 1e-16: err_scale = 1 @@ -236,18 +235,12 @@ def inf_norm(x): test_data_int_tpair, boundaries) from mirgecom.operators import grad_operator - if isinstance(test_data, ConservedVars): - test_grad = make_conserved( - dim=dim, q=grad_operator(discr, test_data.join(), - test_data_flux_bnd.join()) - ) - grad_err = inf_norm((test_grad - exact_grad).join())/err_scale - else: - test_grad = grad_operator(discr, test_data, test_data_flux_bnd) - grad_err = inf_norm(test_grad - exact_grad)/err_scale - - print(f"{test_grad=}") - eoc.add_data_point(actx.to_numpy(h_max), grad_err) + test_grad = grad_operator(discr, test_data, test_data_flux_bnd) + grad_err = \ + max(flatten(componentwise_norms(discr, test_grad - exact_grad, np.inf), + actx)) / err_scale + + eoc.add_data_point(actx.to_numpy(h_max), actx.to_numpy(grad_err)) assert ( eoc.order_estimate() >= order - 0.5 diff --git a/test/test_restart.py b/test/test_restart.py index 2b6b6dcdb..71ec62a8d 100644 --- a/test/test_restart.py +++ b/test/test_restart.py @@ -76,4 +76,4 @@ def test_restart_cv(actx_factory, nspecies): restart_data = read_restart_data(actx, rst_filename) resid = test_state - restart_data["state"] - assert actx.to_numpy(discr.norm(resid.join(), np.inf)) == 0 + assert actx.to_numpy(discr.norm(resid, np.inf)) == 0 diff --git a/test/test_viscous.py b/test/test_viscous.py index d30e37cdf..926c57553 100644 --- a/test/test_viscous.py +++ b/test/test_viscous.py @@ -83,7 +83,7 @@ def test_viscous_stress_tensor(actx_factory, transport_model): mom = mass * velocity cv = make_conserved(dim, mass=mass, energy=energy, momentum=mom) - grad_cv = make_conserved(dim, q=op.local_grad(discr, cv.join())) + grad_cv = op.local_grad(discr, cv) if transport_model: tv_model = SimpleTransport(bulk_viscosity=1.0, viscosity=0.5) @@ -203,8 +203,7 @@ def inf_norm(x): cv_flux_bnd = _elbnd_flux(discr, cv_flux_interior, cv_flux_boundary, cv_int_tpair, boundaries) from mirgecom.operators import grad_operator - grad_cv = make_conserved(dim, q=grad_operator(discr, cv.join(), - cv_flux_bnd.join())) + grad_cv = grad_operator(discr, cv, cv_flux_bnd) xp_grad_cv = initializer.exact_grad(x_vec=nodes, eos=eos, cv_exact=cv) xp_grad_v = 1/cv.mass * xp_grad_cv.momentum @@ -244,9 +243,11 @@ def inf_norm(x): assert inf_norm(heat_flux - xp_heat_flux) < 2e-8 # verify diffusive mass flux is zilch (no scalar components) - from mirgecom.viscous import diffusive_flux - j = diffusive_flux(discr, eos, cv, grad_cv) - assert len(j) == 0 + # Don't call for non-multi CV + if cv.has_multispecies: + from mirgecom.viscous import diffusive_flux + j = diffusive_flux(discr, eos, cv, grad_cv) + assert len(j) == 0 xp_e_flux = np.dot(xp_tau, cv.velocity) - xp_heat_flux xp_mom_flux = xp_tau @@ -319,7 +320,7 @@ def test_species_diffusive_flux(actx_factory): cv = make_conserved(dim, mass=mass, energy=energy, momentum=mom, species_mass=species_mass) - grad_cv = make_conserved(dim, q=op.local_grad(discr, cv.join())) + grad_cv = op.local_grad(discr, cv) mu_b = 1.0 mu = 0.5 @@ -392,7 +393,7 @@ def test_diffusive_heat_flux(actx_factory): cv = make_conserved(dim, mass=mass, energy=energy, momentum=mom, species_mass=species_mass) - grad_cv = make_conserved(dim, q=op.local_grad(discr, cv.join())) + grad_cv = op.local_grad(discr, cv) mu_b = 1.0 mu = 0.5 From a5476c2b565b3b12f550ddcc6143f2fc34af29b4 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Sun, 12 Dec 2021 15:01:01 -0600 Subject: [PATCH 06/20] Update tests to subvert discr.norm of emtpy array --- mirgecom/simutil.py | 16 +++++++++++++++- test/test_euler.py | 13 ++++++------- test/test_restart.py | 3 ++- test/test_viscous.py | 2 +- 4 files changed, 24 insertions(+), 10 deletions(-) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index 7696278d3..3eae215a2 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -13,6 +13,7 @@ .. autofunction:: compare_fluid_solutions .. autofunction:: componentwise_norms +.. autofunction:: componentwise_max_error .. autofunction:: check_naninf_local .. autofunction:: check_range_local @@ -313,8 +314,21 @@ def componentwise_norms(discr, fields, order=np.inf): if not isinstance(fields, DOFArray): return map_array_container( partial(componentwise_norms, discr, order=order), fields) + if len(fields) > 0: + return discr.norm(fields, order) + else: + return 0 + + +def max_component_norm(discr, fields, order=np.inf): + """Return the *order*-norm for each component of *fields*. - return discr.norm(fields, order) + .. note:: + This is a collective routine and must be called by all MPI ranks. + """ + actx = fields.array_context + return max(actx.to_numpy(flatten( + componentwise_norms(discr, fields, order), actx))) def generate_and_distribute_mesh(comm, generate_mesh): diff --git a/test/test_euler.py b/test/test_euler.py index c0a22d978..f93cbe70c 100644 --- a/test/test_euler.py +++ b/test/test_euler.py @@ -53,6 +53,7 @@ pytest_generate_tests_for_pyopencl_array_context as pytest_generate_tests) +from mirgecom.simutil import max_component_norm from grudge.shortcuts import make_visualizer from mirgecom.inviscid import get_inviscid_timestep @@ -234,7 +235,8 @@ def test_vortex_rhs(actx_factory, order): discr, eos=IdealSingleGas(), boundaries=boundaries, cv=vortex_soln, time=0.0) - err_max = actx.to_numpy(discr.norm(inviscid_rhs, np.inf)) + err_max = max_component_norm(discr, inviscid_rhs, np.inf) + eoc_rec.add_data_point(1.0 / nel_1d, err_max) logger.info( @@ -293,8 +295,7 @@ def test_lump_rhs(actx_factory, dim, order): ) expected_rhs = lump.exact_rhs(discr, cv=lump_soln, time=0) - err_max = actx.to_numpy( - discr.norm((inviscid_rhs-expected_rhs), np.inf)) + err_max = max_component_norm(discr, inviscid_rhs-expected_rhs, np.inf) if err_max > maxxerr: maxxerr = err_max @@ -492,9 +493,7 @@ def rhs(t, q): write_soln(state=cv) cv = rk4_step(cv, t, dt, rhs) - cv = make_conserved( - dim, q=filter_modally(discr, "vol", cutoff, frfunc, cv) - ) + cv = filter_modally(discr, "vol", cutoff, frfunc, cv) t += dt istep += 1 @@ -506,7 +505,7 @@ def rhs(t, q): maxerr = max(write_soln(cv, False)) else: expected_result = initializer(nodes, time=t) - maxerr = actx.to_numpy(discr.norm(cv - expected_result, np.inf)) + maxerr = max_component_norm(discr, cv-expected_result, np.inf) logger.info(f"Max Error: {maxerr}") if maxerr > exittol: diff --git a/test/test_restart.py b/test/test_restart.py index 71ec62a8d..55894c538 100644 --- a/test/test_restart.py +++ b/test/test_restart.py @@ -76,4 +76,5 @@ def test_restart_cv(actx_factory, nspecies): restart_data = read_restart_data(actx, rst_filename) resid = test_state - restart_data["state"] - assert actx.to_numpy(discr.norm(resid, np.inf)) == 0 + from mirgecom.simutil import max_component_norm + assert max_component_norm(discr, resid, np.inf) == 0 diff --git a/test/test_viscous.py b/test/test_viscous.py index 926c57553..99d776a0f 100644 --- a/test/test_viscous.py +++ b/test/test_viscous.py @@ -244,7 +244,7 @@ def inf_norm(x): # verify diffusive mass flux is zilch (no scalar components) # Don't call for non-multi CV - if cv.has_multispecies: + if cv.nspecies > 0: from mirgecom.viscous import diffusive_flux j = diffusive_flux(discr, eos, cv, grad_cv) assert len(j) == 0 From 35ff60fdbcb3f49a62574eeec68cf22e89c6aa1a Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Sun, 12 Dec 2021 15:12:11 -0600 Subject: [PATCH 07/20] Correct/add docstring for max_component_norm --- mirgecom/simutil.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index 3eae215a2..89c593230 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -13,7 +13,7 @@ .. autofunction:: compare_fluid_solutions .. autofunction:: componentwise_norms -.. autofunction:: componentwise_max_error +.. autofunction:: max_component_norm .. autofunction:: check_naninf_local .. autofunction:: check_range_local @@ -321,7 +321,7 @@ def componentwise_norms(discr, fields, order=np.inf): def max_component_norm(discr, fields, order=np.inf): - """Return the *order*-norm for each component of *fields*. + """Return the max *order*-norm over the components of *fields*. .. note:: This is a collective routine and must be called by all MPI ranks. From 82c21f7df55a553caaea77157dc3c4e9a9f233db Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Sun, 12 Dec 2021 16:22:57 -0600 Subject: [PATCH 08/20] Udpate lazy test to handle join purge. --- test/test_lazy.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/test/test_lazy.py b/test/test_lazy.py index 11ba007ff..5b3472bcc 100644 --- a/test/test_lazy.py +++ b/test/test_lazy.py @@ -77,11 +77,14 @@ def get_discr(order): def _isclose(discr, x, y, rel_tol=1e-9, abs_tol=0, return_operands=False): from mirgecom.simutil import componentwise_norms - lhs = componentwise_norms(discr, x - y, np.inf) + from arraycontext import flatten + actx = x.array_context + lhs = actx.to_numpy(flatten(componentwise_norms(discr, x - y, np.inf), actx)) + rhs = np.maximum( rel_tol * np.maximum( - componentwise_norms(discr, x, np.inf), - componentwise_norms(discr, y, np.inf)), + actx.to_numpy(flatten(componentwise_norms(discr, x, np.inf), actx)), + actx.to_numpy(flatten(componentwise_norms(discr, y, np.inf), actx))), abs_tol) is_close = np.all(lhs <= rhs) From 8a592ab13d256bb34aec75c0a591c366a5a53fb0 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Sun, 12 Dec 2021 19:16:22 -0600 Subject: [PATCH 09/20] Use custom production branch temporarily --- .ci-support/production-testing-env.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.ci-support/production-testing-env.sh b/.ci-support/production-testing-env.sh index 130a3c53e..956014435 100755 --- a/.ci-support/production-testing-env.sh +++ b/.ci-support/production-testing-env.sh @@ -8,7 +8,7 @@ set -x # The production capability may be in a CEESD-local mirgecom branch or in a # fork, and is specified through the following settings: # -# export PRODUCTION_BRANCH="" # The base production branch to be installed by emirge +export PRODUCTION_BRANCH="production-purge-join" # The base production branch to be installed by emirge # export PRODUCTION_FORK="" # The fork/home of production changes (if any) # # Multiple production drivers are supported. The user should provide a ':'-delimited From 3c4c24506a83a38032b2d9e204436bfd32588342 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Mon, 13 Dec 2021 07:56:58 -0600 Subject: [PATCH 10/20] Undo production env customization --- .ci-support/production-testing-env.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.ci-support/production-testing-env.sh b/.ci-support/production-testing-env.sh index 956014435..130a3c53e 100755 --- a/.ci-support/production-testing-env.sh +++ b/.ci-support/production-testing-env.sh @@ -8,7 +8,7 @@ set -x # The production capability may be in a CEESD-local mirgecom branch or in a # fork, and is specified through the following settings: # -export PRODUCTION_BRANCH="production-purge-join" # The base production branch to be installed by emirge +# export PRODUCTION_BRANCH="" # The base production branch to be installed by emirge # export PRODUCTION_FORK="" # The fork/home of production changes (if any) # # Multiple production drivers are supported. The user should provide a ':'-delimited From 9e8824ac1157fda9bbd4212a68c932b4626c3414 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Mon, 13 Dec 2021 08:01:48 -0600 Subject: [PATCH 11/20] Use productions version of Euler operator --- mirgecom/euler.py | 53 ++++++++++++----------------------------------- 1 file changed, 13 insertions(+), 40 deletions(-) diff --git a/mirgecom/euler.py b/mirgecom/euler.py index 5fabccefa..324ae8338 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -58,15 +58,12 @@ inviscid_facial_flux ) -from grudge.trace_pair import ( - local_interior_trace_pair, - cross_rank_trace_pairs -) -from mirgecom.fluid import make_conserved +from grudge.trace_pair import interior_trace_pairs from mirgecom.operators import div_operator -def euler_operator(discr, eos, boundaries, cv, time=0.0): +def euler_operator(discr, eos, boundaries, cv, time=0.0, + dv=None): r"""Compute RHS of the Euler flow equations. Returns @@ -101,44 +98,20 @@ def euler_operator(discr, eos, boundaries, cv, time=0.0): flow equations. """ # Compute volume contributions - inviscid_flux_vol = inviscid_flux(discr, eos, cv) - + inviscid_flux_vol = inviscid_flux(discr, eos.pressure(cv), cv) + interior_cv = interior_trace_pairs(discr, cv) # Compute interface contributions inviscid_flux_bnd = ( - # Rank-local contributions - inviscid_facial_flux( - discr, - eos=eos, - cv_tpair=local_interior_trace_pair(discr, cv) - ) - # Cross-rank contributions - + sum( - inviscid_facial_flux( - discr, - eos=eos, - cv_tpair=part_tpair - ) for part_tpair in cross_rank_trace_pairs(discr, cv) - ) - # Boundary condition contributions - + sum( - boundaries[btag].inviscid_boundary_flux( - discr, - btag=btag, - cv=cv, - eos=eos, - time=time - ) for btag in boundaries - ) + # Interior faces + + sum(inviscid_facial_flux(discr, eos=eos, cv_tpair=part_tpair) + for part_tpair in interior_cv) + # Domain boundaries + + sum(boundaries[btag].inviscid_divergence_flux(discr, btag=btag, cv=cv, + eos=eos, time=time) + for btag in boundaries) ) - return -div_operator(discr, inviscid_flux_vol, inviscid_flux_bnd) - -def inviscid_operator(discr, eos, boundaries, q, t=0.0): - """Interface :function:`euler_operator` with backwards-compatible API.""" - from warnings import warn - warn("Do not call inviscid_operator; it is now called euler_operator. This" - "function will disappear August 1, 2021", DeprecationWarning, stacklevel=2) - return euler_operator(discr, eos, boundaries, make_conserved(discr.dim, q=q), t) + return -div_operator(discr, inviscid_flux_vol, inviscid_flux_bnd) # By default, run unitless From b3e589144e1d418ce7692f99f29bcca173fa7c90 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Mon, 13 Dec 2021 08:12:27 -0600 Subject: [PATCH 12/20] Revert to main-compatible Euler API --- mirgecom/euler.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mirgecom/euler.py b/mirgecom/euler.py index 324ae8338..78d993bd5 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -98,7 +98,7 @@ def euler_operator(discr, eos, boundaries, cv, time=0.0, flow equations. """ # Compute volume contributions - inviscid_flux_vol = inviscid_flux(discr, eos.pressure(cv), cv) + inviscid_flux_vol = inviscid_flux(discr, eos, cv) interior_cv = interior_trace_pairs(discr, cv) # Compute interface contributions inviscid_flux_bnd = ( @@ -106,7 +106,7 @@ def euler_operator(discr, eos, boundaries, cv, time=0.0, + sum(inviscid_facial_flux(discr, eos=eos, cv_tpair=part_tpair) for part_tpair in interior_cv) # Domain boundaries - + sum(boundaries[btag].inviscid_divergence_flux(discr, btag=btag, cv=cv, + + sum(boundaries[btag].inviscid_boundary_flux(discr, btag=btag, cv=cv, eos=eos, time=time) for btag in boundaries) ) From 6e2e4d3e1178de7190ffb17ec21cf4d7732c5137 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Mon, 13 Dec 2021 10:34:41 -0600 Subject: [PATCH 13/20] Remove extraneous argument to Euler operator --- mirgecom/euler.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/mirgecom/euler.py b/mirgecom/euler.py index 78d993bd5..c21da5ea6 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -62,8 +62,7 @@ from mirgecom.operators import div_operator -def euler_operator(discr, eos, boundaries, cv, time=0.0, - dv=None): +def euler_operator(discr, eos, boundaries, cv, time=0.0): r"""Compute RHS of the Euler flow equations. Returns From 2c8c3e816ba882b44a3bfe729e8c838f82ffdc1c Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Mon, 13 Dec 2021 15:07:52 -0600 Subject: [PATCH 14/20] Clean up max expression for util Co-authored-by: Matt Smith --- mirgecom/simutil.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index 89c593230..0fa76f8df 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -327,8 +327,8 @@ def max_component_norm(discr, fields, order=np.inf): This is a collective routine and must be called by all MPI ranks. """ actx = fields.array_context - return max(actx.to_numpy(flatten( - componentwise_norms(discr, fields, order), actx))) + return actx.to_numpy(actx.np.max( + componentwise_norms(discr, fields, order), actx)) def generate_and_distribute_mesh(comm, generate_mesh): From d3275d73a1b87ab7ddedf8d52d4d817e2248020b Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Mon, 13 Dec 2021 17:09:07 -0600 Subject: [PATCH 15/20] Revert to previous expression --- mirgecom/simutil.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index 0fa76f8df..89c593230 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -327,8 +327,8 @@ def max_component_norm(discr, fields, order=np.inf): This is a collective routine and must be called by all MPI ranks. """ actx = fields.array_context - return actx.to_numpy(actx.np.max( - componentwise_norms(discr, fields, order), actx)) + return max(actx.to_numpy(flatten( + componentwise_norms(discr, fields, order), actx))) def generate_and_distribute_mesh(comm, generate_mesh): From 0824668e17ab09910d13504beec76993f87b6369 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Mon, 13 Dec 2021 17:13:36 -0600 Subject: [PATCH 16/20] Remove crufty check on empty species flux --- test/test_viscous.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/test/test_viscous.py b/test/test_viscous.py index 99d776a0f..2d6e7a70c 100644 --- a/test/test_viscous.py +++ b/test/test_viscous.py @@ -242,13 +242,6 @@ def inf_norm(x): xp_heat_flux = -kappa*xp_grad_t assert inf_norm(heat_flux - xp_heat_flux) < 2e-8 - # verify diffusive mass flux is zilch (no scalar components) - # Don't call for non-multi CV - if cv.nspecies > 0: - from mirgecom.viscous import diffusive_flux - j = diffusive_flux(discr, eos, cv, grad_cv) - assert len(j) == 0 - xp_e_flux = np.dot(xp_tau, cv.velocity) - xp_heat_flux xp_mom_flux = xp_tau from mirgecom.viscous import viscous_flux From a2e7449d4d1535fc1717211458a1a453fccf77c3 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Tue, 14 Dec 2021 18:08:38 -0600 Subject: [PATCH 17/20] Make comment about dumb work-around --- mirgecom/simutil.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index 89c593230..48dee31b9 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -316,7 +316,8 @@ def componentwise_norms(discr, fields, order=np.inf): partial(componentwise_norms, discr, order=order), fields) if len(fields) > 0: return discr.norm(fields, order) - else: + else: + # FIXME: This work-around for #575 can go away after #569 return 0 From b4df6b8f93a89e3b73ef1946164e41388e26cac1 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Tue, 14 Dec 2021 18:13:29 -0600 Subject: [PATCH 18/20] deflakeate --- mirgecom/simutil.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index 48dee31b9..f310b7a7f 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -316,7 +316,7 @@ def componentwise_norms(discr, fields, order=np.inf): partial(componentwise_norms, discr, order=order), fields) if len(fields) > 0: return discr.norm(fields, order) - else: + else: # FIXME: This work-around for #575 can go away after #569 return 0 From 74558fa31d5f7a4807117a4e4c527b7e6fceef5a Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Tue, 14 Dec 2021 22:34:38 -0600 Subject: [PATCH 19/20] Update filter docs and remove decorator --- mirgecom/filter.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/mirgecom/filter.py b/mirgecom/filter.py index 1b31a199a..da09d9526 100644 --- a/mirgecom/filter.py +++ b/mirgecom/filter.py @@ -48,9 +48,13 @@ import numpy as np import grudge.dof_desc as dof_desc +from arraycontext import map_array_container + +from functools import partial + from meshmode.dof_array import DOFArray + from pytools import keyed_memoize_in -from pytools.obj_array import obj_array_vectorized_n_args def exponential_mode_response_function(mode, alpha, cutoff, nfilt, filter_order): @@ -163,7 +167,6 @@ def apply_spectral_filter(actx, modal_field, discr, cutoff, ) -@obj_array_vectorized_n_args def filter_modally(dcoll, dd, cutoff, mode_resp_func, field): """Stand-alone procedural interface to spectral filtering. @@ -189,26 +192,23 @@ def filter_modally(dcoll, dd, cutoff, mode_resp_func, field): Mode below which *field* will not be filtered mode_resp_func: Modal response function returns a filter coefficient for input mode id - field: :class:`numpy.ndarray` - DOFArray or object array of DOFArrays + field: :class:`mirgecom.fluid.ConservedVars` + An array container containing the relevant field(s) to filter. Returns ------- - result: numpy.ndarray - Filtered version of *field*. + result: :class:`mirgecom.fluid.ConservedVars` + An array container containing the filtered field(s). """ - dd = dof_desc.as_dofdesc(dd) - dd_modal = dof_desc.DD_VOLUME_MODAL - discr = dcoll.discr_from_dd(dd) - - from arraycontext import map_array_container - from functools import partial if not isinstance(field, DOFArray): return map_array_container( partial(filter_modally, dcoll, dd, cutoff, mode_resp_func), field ) actx = field.array_context + dd = dof_desc.as_dofdesc(dd) + dd_modal = dof_desc.DD_VOLUME_MODAL + discr = dcoll.discr_from_dd(dd) modal_map = dcoll.connection_from_dds(dd, dd_modal) nodal_map = dcoll.connection_from_dds(dd_modal, dd) From 4951228875776a84d8a26e2b7eceee896b8e7852 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Wed, 15 Dec 2021 05:58:23 -0600 Subject: [PATCH 20/20] Fix up inviscid boundary flux interface per upstream changes. --- mirgecom/euler.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mirgecom/euler.py b/mirgecom/euler.py index c21da5ea6..057224a67 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -105,8 +105,8 @@ def euler_operator(discr, eos, boundaries, cv, time=0.0): + sum(inviscid_facial_flux(discr, eos=eos, cv_tpair=part_tpair) for part_tpair in interior_cv) # Domain boundaries - + sum(boundaries[btag].inviscid_boundary_flux(discr, btag=btag, cv=cv, - eos=eos, time=time) + + sum(boundaries[btag].inviscid_divergence_flux(discr, btag=btag, cv=cv, + eos=eos, time=time) for btag in boundaries) )