From b9bb3f47375f442cce1e2446012f0b7f22aebb06 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Mon, 6 Jun 2022 09:24:02 -0500 Subject: [PATCH 01/21] Add @w-hagen artificial viscosity operator. --- .gitignore | 2 +- conda-env.yml | 1 + doc/operators/artificial_viscosity.rst | 4 + doc/operators/operators.rst | 1 + examples/README.md | 1 + examples/doublemach-mpi.py | 491 +++++++++++++++++++++++++ mirgecom/artificial_viscosity.py | 340 +++++++++++++++++ mirgecom/boundary.py | 77 +++- mirgecom/initializers.py | 118 ++++++ test/test_av.py | 265 +++++++++++++ 10 files changed, 1297 insertions(+), 3 deletions(-) create mode 100644 doc/operators/artificial_viscosity.rst create mode 100644 examples/doublemach-mpi.py create mode 100644 mirgecom/artificial_viscosity.py create mode 100644 test/test_av.py diff --git a/.gitignore b/.gitignore index 5c0c82d53..96005093a 100644 --- a/.gitignore +++ b/.gitignore @@ -46,5 +46,5 @@ test/nodal-dg .run-pylint.py # Emacs backups -.#* +.\#* \#* diff --git a/conda-env.yml b/conda-env.yml index 3d34f71b9..6b96edb87 100644 --- a/conda-env.yml +++ b/conda-env.yml @@ -25,3 +25,4 @@ dependencies: - pydocstyle - cantera - h5py * nompi_* # Make sure cantera does not pull in MPI through h5py +- gmsh diff --git a/doc/operators/artificial_viscosity.rst b/doc/operators/artificial_viscosity.rst new file mode 100644 index 000000000..9b583612a --- /dev/null +++ b/doc/operators/artificial_viscosity.rst @@ -0,0 +1,4 @@ +Artificial Viscosity +==================== + +.. automodule:: mirgecom.artificial_viscosity diff --git a/doc/operators/operators.rst b/doc/operators/operators.rst index c0b40c7fe..0907c286d 100644 --- a/doc/operators/operators.rst +++ b/doc/operators/operators.rst @@ -6,4 +6,5 @@ Operators wave-eq diffusion + artificial_viscosity gas-dynamics diff --git a/examples/README.md b/examples/README.md index fa3728ea1..9f0cd2861 100644 --- a/examples/README.md +++ b/examples/README.md @@ -15,5 +15,6 @@ unique features they exercise are as follows: - `sod-mpi.py`: Sod's shock case: Fluid test case with strong shock - `vortex-mpi.py`: Isentropic vortex advection: outflow boundaries, verification - `hotplate-mpi.py`: Isothermal BC verification (prescribed exact soln) +- `doublemach-mpi.py`: AV test case - `nsmix-mpi.py`: Viscous mixture w/Pyrometheus-based EOS - `poiseuille-mpi.py`: Poiseuille flow verification case \ No newline at end of file diff --git a/examples/doublemach-mpi.py b/examples/doublemach-mpi.py new file mode 100644 index 000000000..18931ee39 --- /dev/null +++ b/examples/doublemach-mpi.py @@ -0,0 +1,491 @@ +"""Demonstrate double mach reflection.""" + +__copyright__ = """ +Copyright (C) 2020 University of Illinois Board of Trustees +""" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import logging +import numpy as np +import pyopencl as cl +import pyopencl.tools as cl_tools +from functools import partial + +from arraycontext import thaw +from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from grudge.dof_desc import DTAG_BOUNDARY +from grudge.eager import EagerDGDiscretization +from grudge.shortcuts import make_visualizer + + +from mirgecom.euler import euler_operator +from mirgecom.artificial_viscosity import ( + av_laplacian_operator, + smoothness_indicator +) +from mirgecom.io import make_init_message +from mirgecom.mpi import mpi_entry_point +from mirgecom.integrators import rk4_step +from mirgecom.steppers import advance_state +from mirgecom.boundary import ( + AdiabaticNoslipMovingBoundary, + PrescribedFluidBoundary +) +from mirgecom.initializers import DoubleMachReflection +from mirgecom.eos import IdealSingleGas +from mirgecom.transport import SimpleTransport +from mirgecom.simutil import get_sim_timestep + +from logpyle import set_dt +from mirgecom.euler import extract_vars_for_logging, units_for_logging +from mirgecom.logging_quantities import ( + initialize_logmgr, + logmgr_add_many_discretization_quantities, + logmgr_add_device_name, + logmgr_add_device_memory_usage, + set_sim_state +) + +logger = logging.getLogger(__name__) + + +class MyRuntimeError(RuntimeError): + """Simple exception to kill the simulation.""" + + pass + + +def get_doublemach_mesh(): + """Generate or import a grid using `gmsh`. + + Input required: + doubleMach.msh (read existing mesh) + + This routine will generate a new grid if it does + not find the grid file (doubleMach.msh). + """ + from meshmode.mesh.io import ( + read_gmsh, + generate_gmsh, + ScriptSource, + ) + import os + meshfile = "doubleMach.msh" + if not os.path.exists(meshfile): + mesh = generate_gmsh( + ScriptSource(""" + x0=1.0/6.0; + setsize=0.025; + Point(1) = {0, 0, 0, setsize}; + Point(2) = {x0,0, 0, setsize}; + Point(3) = {4, 0, 0, setsize}; + Point(4) = {4, 1, 0, setsize}; + Point(5) = {0, 1, 0, setsize}; + Line(1) = {1, 2}; + Line(2) = {2, 3}; + Line(5) = {3, 4}; + Line(6) = {4, 5}; + Line(7) = {5, 1}; + Line Loop(8) = {-5, -6, -7, -1, -2}; + Plane Surface(8) = {8}; + Physical Surface('domain') = {8}; + Physical Curve('ic1') = {6}; + Physical Curve('ic2') = {7}; + Physical Curve('ic3') = {1}; + Physical Curve('wall') = {2}; + Physical Curve('out') = {5}; + """, "geo"), force_ambient_dim=2, dimensions=2, target_unit="M", + output_file_name=meshfile) + else: + mesh = read_gmsh(meshfile, force_ambient_dim=2) + + return mesh + + +@mpi_entry_point +def main(ctx_factory=cl.create_some_context, use_logmgr=True, + use_leap=False, use_profiling=False, use_overintegration=False, + casename=None, rst_filename=None, actx_class=None, lazy=False): + """Drive the example.""" + if actx_class is None: + raise RuntimeError("Array context class missing.") + + cl_ctx = ctx_factory() + + if casename is None: + casename = "mirgecom" + + from mpi4py import MPI + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + nparts = comm.Get_size() + + logmgr = initialize_logmgr(use_logmgr, + filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) + + if use_profiling: + queue = cl.CommandQueue( + cl_ctx, properties=cl.command_queue_properties.PROFILING_ENABLE) + else: + queue = cl.CommandQueue(cl_ctx) + + if lazy: + actx = actx_class(comm, queue, mpi_base_tag=12000) + else: + actx = actx_class(comm, queue, + allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)), + force_device_scalars=True) + + # Timestepping control + current_step = 0 + timestepper = rk4_step + t_final = 1.0e-3 + current_cfl = 0.1 + current_dt = 1.0e-4 + current_t = 0 + constant_cfl = False + + # Some i/o frequencies + nstatus = 10 + nviz = 100 + nrestart = 100 + nhealth = 1 + + rst_path = "restart_data/" + rst_pattern = ( + rst_path + "{cname}-{step:04d}-{rank:04d}.pkl" + ) + if rst_filename: # read the grid from restart data + rst_filename = f"{rst_filename}-{rank:04d}.pkl" + + from mirgecom.restart import read_restart_data + restart_data = read_restart_data(actx, rst_filename) + local_mesh = restart_data["local_mesh"] + local_nelements = local_mesh.nelements + global_nelements = restart_data["global_nelements"] + assert restart_data["nparts"] == nparts + else: # generate the grid from scratch + gen_grid = partial(get_doublemach_mesh) + from mirgecom.simutil import generate_and_distribute_mesh + local_mesh, global_nelements = generate_and_distribute_mesh(comm, gen_grid) + local_nelements = local_mesh.nelements + + from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD + from meshmode.discretization.poly_element import \ + default_simplex_group_factory, QuadratureSimplexGroupFactory + + order = 3 + discr = EagerDGDiscretization( + actx, local_mesh, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: default_simplex_group_factory( + base_dim=local_mesh.dim, order=order), + DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order + 1) + }, + mpi_communicator=comm + ) + nodes = thaw(discr.nodes(), actx) + + if use_overintegration: + quadrature_tag = DISCR_TAG_QUAD + else: + quadrature_tag = None # noqa + + dim = 2 + if logmgr: + logmgr_add_device_name(logmgr, queue) + logmgr_add_device_memory_usage(logmgr, queue) + logmgr_add_many_discretization_quantities(logmgr, discr, dim, + extract_vars_for_logging, units_for_logging) + + logmgr.add_watches([ + ("step.max", "step = {value}, "), + ("t_sim.max", "sim time: {value:1.6e} s\n"), + ("min_pressure", "------- P (min, max) (Pa) = ({value:1.9e}, "), + ("max_pressure", "{value:1.9e})\n"), + ("min_temperature", "------- T (min, max) (K) = ({value:1.9e}, "), + ("max_temperature", "{value:1.9e})\n"), + ("t_step.max", "------- step walltime: {value:6g} s, "), + ("t_log.max", "log walltime: {value:6g} s") + ]) + + # Solution setup and initialization + s0 = -6.0 + kappa = 1.0 + alpha = 2.0e-2 + # {{{ Initialize simple transport model + kappa_t = 1e-5 + sigma_v = 1e-5 + # }}} + + initializer = DoubleMachReflection() + + from mirgecom.gas_model import GasModel, make_fluid_state + transport_model = SimpleTransport(viscosity=sigma_v, + thermal_conductivity=kappa_t) + eos = IdealSingleGas() + gas_model = GasModel(eos=eos, transport=transport_model) + + def _boundary_state(discr, btag, gas_model, state_minus, **kwargs): + actx = state_minus.array_context + bnd_discr = discr.discr_from_dd(btag) + nodes = thaw(bnd_discr.nodes(), actx) + return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, + **kwargs), gas_model) + + boundaries = { + DTAG_BOUNDARY("ic1"): + PrescribedFluidBoundary(boundary_state_func=_boundary_state), + DTAG_BOUNDARY("ic2"): + PrescribedFluidBoundary(boundary_state_func=_boundary_state), + DTAG_BOUNDARY("ic3"): + PrescribedFluidBoundary(boundary_state_func=_boundary_state), + DTAG_BOUNDARY("wall"): AdiabaticNoslipMovingBoundary(), + DTAG_BOUNDARY("out"): AdiabaticNoslipMovingBoundary(), + } + + if rst_filename: + current_t = restart_data["t"] + current_step = restart_data["step"] + current_cv = restart_data["cv"] + if logmgr: + from mirgecom.logging_quantities import logmgr_set_time + logmgr_set_time(logmgr, current_step, current_t) + else: + # Set the current state from time 0 + current_cv = initializer(nodes) + current_state = make_fluid_state(cv=current_cv, gas_model=gas_model) + + visualizer = make_visualizer(discr, + discr.order if discr.dim == 2 else discr.order) + + initname = initializer.__class__.__name__ + eosname = eos.__class__.__name__ + init_message = make_init_message( + dim=dim, + order=order, + nelements=local_nelements, + global_nelements=global_nelements, + dt=current_dt, + t_final=t_final, + nstatus=nstatus, + nviz=nviz, + cfl=current_cfl, + constant_cfl=constant_cfl, + initname=initname, + eosname=eosname, + casename=casename, + ) + if rank == 0: + logger.info(init_message) + + def my_write_viz(step, t, state, dv, tagged_cells): + viz_fields = [("cv", state), + ("dv", dv), + ("tagged_cells", tagged_cells)] + from mirgecom.simutil import write_visfile + write_visfile(discr, viz_fields, visualizer, vizname=casename, + step=step, t=t, overwrite=True) + + def my_write_restart(step, t, state): + rst_fname = rst_pattern.format(cname=casename, step=step, rank=rank) + if rst_fname != rst_filename: + rst_data = { + "local_mesh": local_mesh, + "cv": state, + "t": t, + "step": step, + "order": order, + "global_nelements": global_nelements, + "num_parts": nparts + } + from mirgecom.restart import write_restart_file + write_restart_file(actx, rst_data, rst_fname, comm) + + def my_health_check(state, dv): + # Note: This health check is tuned s.t. it is a test that + # the case gets the expected solution. If dt,t_final or + # other run parameters are changed, this check should + # be changed accordingly. + health_error = False + from mirgecom.simutil import check_naninf_local, check_range_local + if check_naninf_local(discr, "vol", dv.pressure): + health_error = True + logger.info(f"{rank=}: NANs/Infs in pressure data.") + + from mirgecom.simutil import allsync + if allsync(check_range_local(discr, "vol", dv.pressure, .9, 18.6), + comm, op=MPI.LOR): + health_error = True + from grudge.op import nodal_max, nodal_min + p_min = actx.to_numpy(nodal_min(discr, "vol", dv.pressure)) + p_max = actx.to_numpy(nodal_max(discr, "vol", dv.pressure)) + logger.info(f"Pressure range violation ({p_min=}, {p_max=})") + + if check_naninf_local(discr, "vol", dv.temperature): + health_error = True + logger.info(f"{rank=}: NANs/INFs in temperature data.") + + if allsync( + check_range_local(discr, "vol", dv.temperature, 2.48e-3, 1.071e-2), + comm, op=MPI.LOR): + health_error = True + from grudge.op import nodal_max, nodal_min + t_min = actx.to_numpy(nodal_min(discr, "vol", dv.temperature)) + t_max = actx.to_numpy(nodal_max(discr, "vol", dv.temperature)) + logger.info(f"Temperature range violation ({t_min=}, {t_max=})") + + return health_error + + def my_pre_step(step, t, dt, state): + fluid_state = make_fluid_state(state, gas_model) + dv = fluid_state.dv + try: + + if logmgr: + logmgr.tick_before() + + from mirgecom.simutil import check_step + do_viz = check_step(step=step, interval=nviz) + do_restart = check_step(step=step, interval=nrestart) + do_health = check_step(step=step, interval=nhealth) + + if do_health: + from mirgecom.simutil import allsync + health_errors = allsync(my_health_check(state, dv), comm, + op=MPI.LOR) + if health_errors: + if rank == 0: + logger.info("Fluid solution failed health check.") + raise MyRuntimeError("Failed simulation health check.") + + if do_restart: + my_write_restart(step=step, t=t, state=state) + + if do_viz: + tagged_cells = smoothness_indicator(discr, state.mass, s0=s0, + kappa=kappa) + my_write_viz(step=step, t=t, state=state, dv=dv, + tagged_cells=tagged_cells) + + except MyRuntimeError: + if rank == 0: + logger.info("Errors detected; attempting graceful exit.") + tagged_cells = smoothness_indicator(discr, state.mass, s0=s0, + kappa=kappa) + my_write_viz(step=step, t=t, state=state, + tagged_cells=tagged_cells, dv=dv) + my_write_restart(step=step, t=t, state=state) + raise + + dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_cfl, t_final, constant_cfl) + + return state, dt + + def my_post_step(step, t, dt, state): + # Logmgr needs to know about EOS, dt, dim? + # imo this is a design/scope flaw + if logmgr: + set_dt(logmgr, dt) + set_sim_state(logmgr, dim, state, eos) + logmgr.tick_after() + return state, dt + + def my_rhs(t, state): + fluid_state = make_fluid_state(state, gas_model) + return ( + euler_operator(discr, state=fluid_state, time=t, + boundaries=boundaries, + gas_model=gas_model, quadrature_tag=quadrature_tag) + + av_laplacian_operator(discr, fluid_state=fluid_state, + boundaries=boundaries, + time=t, gas_model=gas_model, + alpha=alpha, s0=s0, kappa=kappa, + quadrature_tag=quadrature_tag) + ) + + current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_cfl, t_final, constant_cfl) + + current_step, current_t, current_cv = \ + advance_state(rhs=my_rhs, timestepper=timestepper, + pre_step_callback=my_pre_step, + post_step_callback=my_post_step, dt=current_dt, + state=current_state.cv, t=current_t, t_final=t_final) + + # Dump the final data + if rank == 0: + logger.info("Checkpointing final state ...") + current_state = make_fluid_state(current_cv, gas_model) + final_dv = current_state.dv + tagged_cells = smoothness_indicator(discr, current_cv.mass, s0=s0, kappa=kappa) + my_write_viz(step=current_step, t=current_t, state=current_cv, dv=final_dv, + tagged_cells=tagged_cells) + my_write_restart(step=current_step, t=current_t, state=current_cv) + + if logmgr: + logmgr.close() + elif use_profiling: + print(actx.tabulate_profiling_data()) + + finish_tol = 1e-16 + assert np.abs(current_t - t_final) < finish_tol + + +if __name__ == "__main__": + import argparse + casename = "doublemach" + parser = argparse.ArgumentParser(description=f"MIRGE-Com Example: {casename}") + parser.add_argument("--overintegration", action="store_true", + help="use overintegration in the RHS computations") + parser.add_argument("--lazy", action="store_true", + help="switch to a lazy computation mode") + parser.add_argument("--profiling", action="store_true", + help="turn on detailed performance profiling") + parser.add_argument("--log", action="store_true", default=True, + help="turn on logging") + parser.add_argument("--leap", action="store_true", + help="use leap timestepper") + parser.add_argument("--restart_file", help="root name of restart file") + parser.add_argument("--casename", help="casename to use for i/o") + args = parser.parse_args() + lazy = args.lazy + if args.profiling: + if lazy: + raise ValueError("Can't use lazy and profiling together.") + + from grudge.array_context import get_reasonable_array_context_class + actx_class = get_reasonable_array_context_class(lazy=lazy, distributed=True) + + logging.basicConfig(format="%(message)s", level=logging.INFO) + if args.casename: + casename = args.casename + rst_filename = None + if args.restart_file: + rst_filename = args.restart_file + + main(use_logmgr=args.log, use_leap=args.leap, use_profiling=args.profiling, + use_overintegration=args.overintegration, lazy=lazy, + casename=casename, rst_filename=rst_filename, actx_class=actx_class) + +# vim: foldmethod=marker diff --git a/mirgecom/artificial_viscosity.py b/mirgecom/artificial_viscosity.py new file mode 100644 index 000000000..3a429ce6a --- /dev/null +++ b/mirgecom/artificial_viscosity.py @@ -0,0 +1,340 @@ +r""":mod:`mirgecom.artificial_viscosity` Artificial viscosity for hyperbolic systems. + +Consider the following system of conservation laws of the form: + +.. math:: + + \partial_t \mathbf{Q} + \nabla\cdot\mathbf{F} = 0, + +where $\mathbf{Q}$ is the state vector and $\mathbf{F}$ is the vector of +fluxes. This module applies an artificial viscosity term by augmenting +the governing equations in the following way: + +.. math:: + + \partial_t \mathbf{Q} + \nabla\cdot\mathbf{F} = + \nabla\cdot{\varepsilon\nabla\mathbf{Q}}, + +where $\varepsilon$ is the artificial viscosity coefficient. +To evalutate the second order derivative numerically, the problem +is recast as a set of first order problems: + +.. math:: + + \partial_t{\mathbf{Q}} + \nabla\cdot\mathbf{F} &= \nabla\cdot\mathbf{R} \\ + \mathbf{R} &= \varepsilon\nabla\mathbf{Q} + +where $\mathbf{R}$ is an auxiliary variable, and the artificial viscosity +coefficient, $\varepsilon$, is spatially dependent and calculated using the +smoothness indicator of [Persson_2012]_: + +.. math:: + + S_e = \frac{\langle u_{N_p} - u_{N_{p-1}}, u_{N_p} - + u_{N_{p-1}}\rangle_e}{\langle u_{N_p}, u_{N_p} \rangle_e} + +where: + +- $S_e$ is the smoothness indicator +- $u_{N_p}$ is the solution in modal basis at the current polynomial order +- $u_{N_{p-1}}$ is the truncated modal represention to the polynomial order $p-1$ +- The $L_2$ inner product on an element is denoted $\langle \cdot,\cdot \rangle_e$ + +The elementwise viscosity is then calculated: + +.. math:: + + \varepsilon_e = + \begin{cases} + 0, & s_e < s_0 - \kappa \\ + \frac{1}{2}\left( 1 + \sin \frac{\pi(s_e - s_0)}{2 \kappa} \right ), + & s_0-\kappa \le s_e \le s_0+\kappa \\ + 1, & s_e > s_0+\kappa + \end{cases} + +where: + +- $\varepsilon_e$ is the element viscosity +- $s_e = \log_{10}{S_e} \sim 1/p^4$ is the smoothness indicator +- $s_0$ is a reference smoothness value +- $\kappa$ controls the width of the transition between 0 to 1 + +Smoothness Indicator Evaluation +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. autofunction:: smoothness_indicator + +AV RHS Evaluation +^^^^^^^^^^^^^^^^^ + +.. autofunction:: av_laplacian_operator +""" + +__copyright__ = """ +Copyright (C) 2020 University of Illinois Board of Trustees +""" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np + +from pytools import memoize_in, keyed_memoize_in +from functools import partial +from meshmode.dof_array import thaw, DOFArray + +from mirgecom.flux import num_flux_central +from mirgecom.operators import div_operator + +from grudge.trace_pair import ( + interior_trace_pairs, + tracepair_with_discr_tag +) + +from grudge.dof_desc import ( + DOFDesc, + as_dofdesc, + DD_VOLUME_MODAL, + DD_VOLUME +) + +import grudge.op as op + + +class _AVCVTag: + pass + + +class _AVRTag: + pass + + +def av_laplacian_operator(discr, boundaries, fluid_state, alpha, gas_model=None, + kappa=1., s0=-6., time=0, operator_states_quad=None, + grad_cv=None, quadrature_tag=None, boundary_kwargs=None, + indicator=None, divergence_numerical_flux=num_flux_central, + **kwargs): + r"""Compute the artificial viscosity right-hand-side. + + Computes the the right-hand-side term for artificial viscosity. + + .. math:: + + \mbox{RHS}_{\mbox{av}} = \nabla\cdot{\varepsilon\nabla\mathbf{Q}} + + Parameters + ---------- + fluid_state: :class:`mirgecom.gas_model.FluidState` + Fluid state object with the conserved and thermal state. + + boundaries: dict + Dictionary of boundary functions, one for each valid boundary tag + + alpha: float + The maximum artificial viscosity coefficient to be applied + + indicator: :class:`~meshmode.dof_array.DOFArray` + The indicator field used for locating where AV should be applied. If not + supplied by the user, then + :func:`~mirgecom.artificial_viscosity.smoothness_indicator` will be used + with fluid mass density as the indicator field. + + kappa + An optional argument that controls the width of the transition from 0 to 1, + $\kappa$. This parameter defaults to $\kappa=1$. + + s0 + An optional argument that sets the smoothness level to limit + on, $s_0$. Values in the range $(-\infty,0]$ are allowed, where $-\infty$ + results in all cells being tagged and 0 results in none. This parameter + defaults to $s_0=-6$. + + quadrature_tag + An optional identifier denoting a particular quadrature + discretization to use during operator evaluations. + The default value is *None*. + + boundary_kwargs: :class:`dict` + dictionary of extra arguments to pass through to the boundary conditions + + Returns + ------- + :class:`mirgecom.fluid.ConservedVars` + The artificial viscosity operator applied to *q*. + """ + cv = fluid_state.cv + actx = cv.array_context + dd_vol = DOFDesc("vol", quadrature_tag) + dd_faces = DOFDesc("all_faces", quadrature_tag) + + from warnings import warn + + if boundary_kwargs is not None: + warn("The AV boundary_kwargs interface is deprecated, please pass gas_model" + " and time directly.") + if gas_model is None: + gas_model = boundary_kwargs["gas_model"] + if "time" in boundary_kwargs: + time = boundary_kwargs["time"] + + interp_to_surf_quad = partial(tracepair_with_discr_tag, discr, quadrature_tag) + + def interp_to_vol_quad(u): + return op.project(discr, "vol", dd_vol, u) + + if operator_states_quad is None: + from mirgecom.gas_model import make_operator_fluid_states + operator_states_quad = make_operator_fluid_states( + discr, fluid_state, gas_model, boundaries, quadrature_tag) + + vol_state_quad, inter_elem_bnd_states_quad, domain_bnd_states_quad = \ + operator_states_quad + + # Get smoothness indicator based on mass component + if indicator is None: + indicator = smoothness_indicator(discr, fluid_state.mass_density, + kappa=kappa, s0=s0) + + if grad_cv is None: + from mirgecom.navierstokes import grad_cv_operator + grad_cv = grad_cv_operator(discr, gas_model, boundaries, fluid_state, + time=time, quadrature_tag=quadrature_tag, + operator_states_quad=operator_states_quad) + + # Compute R = alpha*grad(Q) + r = -alpha * indicator * grad_cv + + def central_flux_div(utpair): + dd = utpair.dd + normal = thaw(actx, discr.normal(dd)) + return op.project(discr, dd, dd.with_dtag("all_faces"), + # This uses a central vector flux along nhat: + # flux = 1/2 * (grad(Q)- + grad(Q)+) .dot. nhat + divergence_numerical_flux(utpair.int, utpair.ext)@normal) + + # Total flux of grad(Q) across element boundaries + r_bnd = ( + # Rank-local and cross-rank (across parallel partitions) contributions + + sum(central_flux_div(interp_to_surf_quad(tpair=tpair)) + for tpair in interior_trace_pairs(discr, r, tag=_AVRTag)) + + # Contributions from boundary fluxes + + sum(boundaries[btag].av_flux( + discr, + btag=as_dofdesc(btag).with_discr_tag(quadrature_tag), + diffusion=r) for btag in boundaries) + ) + + # Return the AV RHS term + return -div_operator(discr, dd_vol, dd_faces, interp_to_vol_quad(r), r_bnd) + + +def smoothness_indicator(discr, u, kappa=1.0, s0=-6.0): + r"""Calculate the smoothness indicator. + + Parameters + ---------- + u: meshmode.dof_array.DOFArray + The field that is used to calculate the smoothness indicator. + + kappa + An optional argument that controls the width of the transition from 0 to 1. + + s0 + An optional argument that sets the smoothness level to limit + on. Values in the range $(-\infty,0]$ are allowed, where $-\infty$ results in + all cells being tagged and 0 results in none. + + Returns + ------- + meshmode.dof_array.DOFArray + The elementwise constant values between 0 and 1 which indicate the smoothness + of a given element. + """ + if not isinstance(u, DOFArray): + raise ValueError("u argument must be a DOFArray.") + + actx = u.array_context + + @memoize_in(actx, (smoothness_indicator, "smooth_comp_knl")) + def indicator_prg(): + """Compute the smoothness indicator for all elements.""" + from arraycontext import make_loopy_program + from meshmode.transform_metadata import (ConcurrentElementInameTag, + ConcurrentDOFInameTag) + import loopy as lp + t_unit = make_loopy_program([ + "{[iel]: 0 <= iel < nelements}", + "{[idof]: 0 <= idof < ndiscr_nodes_in}", + "{[jdof]: 0 <= jdof < ndiscr_nodes_in}", + "{[kdof]: 0 <= kdof < ndiscr_nodes_in}" + ], + """ + result[iel,idof] = sum(kdof, vec[iel, kdof] \ + * vec[iel, kdof] \ + * modes_active_flag[kdof]) / \ + sum(jdof, vec[iel, jdof] \ + * vec[iel, jdof] \ + + 1.0e-12 / ndiscr_nodes_in) + """, + name="smooth_comp", + ) + return lp.tag_inames(t_unit, {"iel": ConcurrentElementInameTag(), + "idof": ConcurrentDOFInameTag()}) + + @keyed_memoize_in(actx, (smoothness_indicator, + "highest_mode"), + lambda grp: grp.discretization_key()) + def highest_mode(grp): + return actx.from_numpy( + np.asarray([1 if sum(mode_id) == grp.order + else 0 + for mode_id in grp.mode_ids()]) + ) + + # Convert to modal solution representation + modal_map = discr.connection_from_dds(DD_VOLUME, DD_VOLUME_MODAL) + uhat = modal_map(u) + + # Compute smoothness indicator value + indicator = DOFArray( + actx, + data=tuple( + actx.call_loopy( + indicator_prg(), + vec=uhat[grp.index], + modes_active_flag=highest_mode(grp) + )["result"] + for grp in discr.discr_from_dd("vol").groups + ) + ) + indicator = actx.np.log10(indicator + 1.0e-12) + + # Compute artificial viscosity percentage based on indicator and set parameters + yesnol = actx.np.greater(indicator, (s0 - kappa)) + yesnou = actx.np.greater(indicator, (s0 + kappa)) + sin_indicator = actx.np.where( + yesnol, + 0.5 * (1.0 + actx.np.sin(np.pi * (indicator - s0) / (2.0 * kappa))), + 0.0 * indicator, + ) + indicator = actx.np.where(yesnou, 1.0 + 0.0 * indicator, sin_indicator) + + return indicator diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py index 986218ea4..e4a8da5d6 100644 --- a/mirgecom/boundary.py +++ b/mirgecom/boundary.py @@ -50,7 +50,11 @@ from grudge.trace_pair import TracePair from mirgecom.viscous import viscous_facial_flux_central from mirgecom.flux import num_flux_central -from mirgecom.gas_model import make_fluid_state +from mirgecom.gas_model import ( + make_fluid_state, + project_fluid_state +) + from mirgecom.inviscid import inviscid_facial_flux_rusanov from abc import ABCMeta, abstractmethod @@ -237,9 +241,11 @@ class PrescribedFluidBoundary(FluidBoundary): .. automethod:: __init__ .. automethod:: inviscid_divergence_flux - .. automethod:: cv_gradient_flux .. automethod:: viscous_divergence_flux + .. automethod:: cv_gradient_flux .. automethod:: temperature_gradient_flux + .. automethod:: av_flux + .. automethod:: soln_gradient_flux """ def __init__(self, @@ -261,6 +267,8 @@ def __init__(self, boundary_gradient_cv_func=None, # Returns the boundary value for grad(temperature) boundary_gradient_temperature_func=None, + # For artificial viscosity - grad fluid soln on boundary + boundary_grad_av_func=None, ): """Initialize the PrescribedFluidBoundary and methods.""" self._bnd_state_func = boundary_state_func @@ -272,6 +280,11 @@ def __init__(self, self._viscous_flux_func = viscous_flux_func self._bnd_grad_cv_func = boundary_gradient_cv_func self._bnd_grad_temperature_func = boundary_gradient_temperature_func + self._av_num_flux_func = num_flux_central + self._bnd_grad_av_func = boundary_grad_av_func + + if not self._bnd_grad_av_func: + self._bnd_grad_av_func = self._identical_grad_av if not self._inviscid_flux_func and not self._bnd_state_func: from warnings import warn @@ -451,6 +464,39 @@ def viscous_divergence_flux(self, discr, btag, gas_model, state_minus, numerical_flux_func=numerical_flux_func, **kwargs) + # {{{ Boundary interface for artificial viscosity + + def _identical_grad_av(self, grad_av_minus, **kwargs): + return grad_av_minus + + def soln_gradient_flux(self, discr, btag, fluid_state, gas_model, **kwargs): + """Get the flux for solution gradient with AV API.""" + # project the conserved and thermal state to the boundary + fluid_state_minus = project_fluid_state(discr=discr, + src="vol", + tgt=btag, + gas_model=gas_model, + state=fluid_state) + # get the boundary flux for the grad(CV) + return self.cv_gradient_flux(discr=discr, btag=btag, + gas_model=gas_model, + state_minus=fluid_state_minus, + **kwargs) + + def av_flux(self, discr, btag, diffusion, **kwargs): + """Get the diffusive fluxes for the AV operator API.""" + grad_av_minus = discr.project("vol", btag, diffusion) + actx = grad_av_minus.mass[0].array_context + nhat = thaw(discr.normal(btag), actx) + grad_av_plus = self._bnd_grad_av_func( + discr=discr, btag=btag, grad_av_minus=grad_av_minus, **kwargs) + bnd_grad_pair = TracePair(btag, interior=grad_av_minus, + exterior=grad_av_plus) + num_flux = self._av_num_flux_func(bnd_grad_pair.int, bnd_grad_pair.ext)@nhat + return self._boundary_quantity(discr, btag, num_flux, **kwargs) + + # }}} + class DummyBoundary(PrescribedFluidBoundary): """Boundary type that assigns boundary-adjacent soln as the boundary solution.""" @@ -476,12 +522,14 @@ class AdiabaticSlipBoundary(PrescribedFluidBoundary): boundary conditions described in detail in [Poinsot_1992]_. .. automethod:: adiabatic_slip_state + .. automethod:: adiabatic_slip_grad_av """ def __init__(self): """Initialize AdiabaticSlipBoundary.""" PrescribedFluidBoundary.__init__( self, boundary_state_func=self.adiabatic_slip_state, + boundary_grad_av_func=self.adiabatic_slip_grad_av ) def adiabatic_slip_state(self, discr, btag, gas_model, state_minus, **kwargs): @@ -517,17 +565,38 @@ def adiabatic_slip_state(self, discr, btag, gas_model, state_minus, **kwargs): return make_fluid_state(cv=ext_cv, gas_model=gas_model, temperature_seed=t_seed) + def adiabatic_slip_grad_av(self, discr, btag, grad_av_minus, **kwargs): + """Get the exterior grad(Q) on the boundary.""" + # Grab some boundary-relevant data + dim, = grad_av_minus.mass.shape + actx = grad_av_minus.mass[0].array_context + nhat = thaw(discr.norm(btag), actx) + + # Subtract 2*wall-normal component of q + # to enforce q=0 on the wall + s_mom_normcomp = np.outer(nhat, + np.dot(grad_av_minus.momentum, nhat)) + s_mom_flux = grad_av_minus.momentum - 2*s_mom_normcomp + + # flip components to set a neumann condition + return make_conserved(dim, mass=-grad_av_minus.mass, + energy=-grad_av_minus.energy, + momentum=-s_mom_flux, + species_mass=-grad_av_minus.species_mass) + class AdiabaticNoslipMovingBoundary(PrescribedFluidBoundary): r"""Boundary condition implementing a noslip moving boundary. .. automethod:: adiabatic_noslip_state + .. automethod:: adiabatic_noslip_grad_av """ def __init__(self, wall_velocity=None, dim=2): """Initialize boundary device.""" PrescribedFluidBoundary.__init__( self, boundary_state_func=self.adiabatic_noslip_state, + boundary_grad_av_func=self.adiabatic_noslip_grad_av, ) # Check wall_velocity (assumes dim is correct) if wall_velocity is None: @@ -553,6 +622,10 @@ def adiabatic_noslip_state(self, discr, btag, gas_model, state_minus, **kwargs): tseed = state_minus.temperature if state_minus.is_mixture else None return make_fluid_state(cv=cv, gas_model=gas_model, temperature_seed=tseed) + def adiabatic_noslip_grad_av(self, grad_av_minus, **kwargs): + """Get the exterior solution on the boundary.""" + return(-grad_av_minus) + class IsothermalNoSlipBoundary(PrescribedFluidBoundary): r"""Isothermal no-slip viscous wall boundary. diff --git a/mirgecom/initializers.py b/mirgecom/initializers.py index f522248c0..d28b7edab 100644 --- a/mirgecom/initializers.py +++ b/mirgecom/initializers.py @@ -5,6 +5,7 @@ ^^^^^^^^^^^^^^^^^^^^^ .. autoclass:: Vortex2D .. autoclass:: SodShock1D +.. autoclass:: DoubleMachReflection .. autoclass:: Lump .. autoclass:: MulticomponentLump .. autoclass:: Uniform @@ -267,6 +268,123 @@ def __call__(self, x_vec, *, eos=None, **kwargs): momentum=mom) +class DoubleMachReflection: + r"""Implement the double shock reflection problem. + + The double shock reflection solution is crafted after [Woodward_1984]_ + and is defined by: + + .. math:: + + {\rho}(x < x_s(y,t)) &= \gamma \rho_j\\ + {\rho}(x > x_s(y,t)) &= \gamma\\ + {\rho}{V_x}(x < x_s(y,t)) &= u_p \cos(\pi/6)\\ + {\rho}{V_x}(x > x_s(y,t)) &= 0\\ + {\rho}{V_y}(x > x_s(y,t)) &= u_p \sin(\pi/6)\\ + {\rho}{V_y}(x > x_s(y,t)) &= 0\\ + {\rho}E(x < x_s(y,t)) &= (\gamma-1)p_j\\ + {\rho}E(x > x_s(y,t)) &= (\gamma-1), + + where the shock position is given, + + .. math:: + + x_s = x_0 + y/\sqrt{3} + 2 u_s t/\sqrt{3} + + and the normal shock jump relations are + + .. math:: + + \rho_j &= \frac{(\gamma + 1) u_s^2}{(\gamma-1) u_s^2 + 2} \\ + p_j &= \frac{2 \gamma u_s^2 - (\gamma - 1)}{\gamma+1} \\ + u_p &= 2 \frac{u_s^2-1}{(\gamma+1) u_s} + + The initial shock location is given by $x_0$ and $u_s$ is the shock speed. + + .. automethod:: __init__ + .. automethod:: __call__ + """ + + def __init__( + self, shock_location=1.0/6.0, shock_speed=4.0 + ): + """Initialize double shock reflection parameters. + + Parameters + ---------- + shock_location: float + initial location of shock + shock_speed: float + shock speed, Mach number + """ + self._shock_location = shock_location + self._shock_speed = shock_speed + + def __call__(self, x_vec, *, eos=None, time=0, **kwargs): + r""" + Create double mach reflection solution at locations *x_vec*. + + At times $t > 0$, calls to this routine create an advanced solution + under the assumption of constant normal shock speed *shock_speed*. + The advanced solution *is not* the exact solution, but is appropriate + for use as an exact boundary solution on the top and upstream (left) + side of the domain. + + Parameters + ---------- + time: float + Time at which to compute the solution + x_vec: numpy.ndarray + Nodal coordinates + eos: :class:`mirgecom.eos.GasEOS` + Equation of state class to be used in construction of soln (if needed) + """ + t = time + # Fail if numdim is other than 2 + if(len(x_vec)) != 2: + raise ValueError("Case only defined for 2 dimensions") + if eos is None: + eos = IdealSingleGas() + + gm1 = eos.gamma() - 1.0 + gp1 = eos.gamma() + 1.0 + gmn1 = 1.0 / gm1 + x_rel = x_vec[0] + y_rel = x_vec[1] + actx = x_rel.array_context + + # Normal Shock Relations + shock_speed_2 = self._shock_speed * self._shock_speed + rho_jump = gp1 * shock_speed_2 / (gm1 * shock_speed_2 + 2.) + p_jump = (2. * eos.gamma() * shock_speed_2 - gm1) / gp1 + up = 2. * (shock_speed_2 - 1.) / (gp1 * self._shock_speed) + + rhol = eos.gamma() * rho_jump + rhor = eos.gamma() + ul = up * np.cos(np.pi/6.0) + ur = 0.0 + vl = - up * np.sin(np.pi/6.0) + vr = 0.0 + rhoel = gmn1 * p_jump + rhoer = gmn1 * 1.0 + + xinter = (self._shock_location + y_rel/np.sqrt(3.0) + + 2.0*self._shock_speed*t/np.sqrt(3.0)) + sigma = 0.05 + xtanh = 1.0/sigma*(x_rel-xinter) + mass = rhol/2.0*(actx.np.tanh(-xtanh)+1.0)+rhor/2.0*(actx.np.tanh(xtanh)+1.0) + rhoe = (rhoel/2.0*(actx.np.tanh(-xtanh)+1.0) + + rhoer/2.0*(actx.np.tanh(xtanh)+1.0)) + u = ul/2.0*(actx.np.tanh(-xtanh)+1.0)+ur/2.0*(actx.np.tanh(xtanh)+1.0) + v = vl/2.0*(actx.np.tanh(-xtanh)+1.0)+vr/2.0*(actx.np.tanh(xtanh)+1.0) + + vel = make_obj_array([u, v]) + mom = mass * vel + energy = rhoe + .5*mass*np.dot(vel, vel) + + return make_conserved(dim=2, mass=mass, energy=energy, momentum=mom) + + class Lump: r"""Solution initializer for N-dimensional Gaussian lump of mass. diff --git a/test/test_av.py b/test/test_av.py new file mode 100644 index 000000000..4d73d4587 --- /dev/null +++ b/test/test_av.py @@ -0,0 +1,265 @@ +"""Test the artificial viscosity functions.""" + +__copyright__ = """ +Copyright (C) 2020 University of Illinois Board of Trustees +""" + + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import logging +import numpy as np +import pyopencl as cl +import pytest + +from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw +from meshmode.mesh import BTAG_ALL + +from mirgecom.artificial_viscosity import ( + av_laplacian_operator, + smoothness_indicator +) +from mirgecom.fluid import make_conserved +from mirgecom.gas_model import ( + GasModel, + make_fluid_state, + project_fluid_state +) +from mirgecom.eos import IdealSingleGas + +from grudge.eager import EagerDGDiscretization + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests, +) + +from pytools.obj_array import make_obj_array + +logger = logging.getLogger(__name__) + + +@pytest.mark.parametrize("dim", [1, 2, 3]) +@pytest.mark.parametrize("order", [1, 5]) +def test_tag_cells(ctx_factory, dim, order): + """Test tag_cells. + + Tests that the cells tagging properly tags cells + given prescribed solutions. + """ + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + nel_1d = 2 + tolerance = 1.e-16 + + def norm_indicator(expected, discr, soln, **kwargs): + return(discr.norm(expected-smoothness_indicator(discr, soln, **kwargs), + np.inf)) + + from meshmode.mesh.generation import generate_regular_rect_mesh + + mesh = generate_regular_rect_mesh( + a=(-1.0, )*dim, b=(1.0, )*dim, n=(nel_1d, ) * dim + ) + + discr = EagerDGDiscretization(actx, mesh, order=order) + nodes = thaw(actx, discr.nodes()) + nele = mesh.nelements + zeros = 0.0*nodes[0] + + # test jump discontinuity + soln = actx.np.where(nodes[0] > 0.0+zeros, 1.0+zeros, zeros) + err = norm_indicator(1.0, discr, soln) + + assert err < tolerance, "Jump discontinuity should trigger indicator (1.0)" + + # get meshmode polynomials + group = discr.discr_from_dd("vol").groups[0] + basis = group.basis() # only one group + unit_nodes = group.unit_nodes + modes = group.mode_ids() + order = group.order + + # loop over modes and check smoothness + for i, mode in enumerate(modes): + ele_soln = basis[i](unit_nodes) + soln[0].set(np.tile(ele_soln, (nele, 1))) + if sum(mode) == order: + expected = 1.0 + else: + expected = 0.0 + err = norm_indicator(expected, discr, soln) + assert err < tolerance, "Only highest modes should trigger indicator (1.0)" + + # Test s0 + s0 = -1. + eps = 1.0e-6 + + phi_n_p = np.sqrt(np.power(10, s0)) + phi_n_pm1 = np.sqrt(1. - np.power(10, s0)) + + # pick a polynomial of order n_p, n_p-1 + n_p = np.array(np.nonzero((np.sum(modes, axis=1) == order))).flat[0] + n_pm1 = np.array(np.nonzero((np.sum(modes, axis=1) == order-1))).flat[0] + + # create test soln perturbed around + # Solution above s0 + ele_soln = ((phi_n_p+eps)*basis[n_p](unit_nodes) + + phi_n_pm1*basis[n_pm1](unit_nodes)) + soln[0].set(np.tile(ele_soln, (nele, 1))) + err = norm_indicator(1.0, discr, soln, s0=s0, kappa=0.0) + assert err < tolerance, ( + "A function with an indicator >s0 should trigger indicator") + + # Solution below s0 + ele_soln = ((phi_n_p-eps)*basis[n_p](unit_nodes) + + phi_n_pm1*basis[n_pm1](unit_nodes)) + soln[0].set(np.tile(ele_soln, (nele, 1))) + err = norm_indicator(0.0, discr, soln, s0=s0, kappa=0.0) + assert err < tolerance, ( + "A function with an indicator tolerance, "s_e>s_0-kappa should trigger indicator" + + # lower bound + err = norm_indicator(1.0, discr, soln, s0=s0-(kappa+shift), kappa=kappa) + assert err < tolerance, "s_e>s_0+kappa should fully trigger indicator (1.0)" + err = norm_indicator(1.0, discr, soln, s0=s0-(kappa-shift), kappa=kappa) + assert err > tolerance, "s_e Date: Mon, 6 Jun 2022 16:48:35 -0500 Subject: [PATCH 02/21] eliminate some deprecated constructs --- examples/doublemach-mpi.py | 40 ++++++++++++++------------------------ mirgecom/discretization.py | 5 +++-- 2 files changed, 18 insertions(+), 27 deletions(-) diff --git a/examples/doublemach-mpi.py b/examples/doublemach-mpi.py index 18931ee39..172fc17a1 100644 --- a/examples/doublemach-mpi.py +++ b/examples/doublemach-mpi.py @@ -33,7 +33,6 @@ from arraycontext import thaw from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.dof_desc import DTAG_BOUNDARY -from grudge.eager import EagerDGDiscretization from grudge.shortcuts import make_visualizer @@ -60,7 +59,7 @@ from mirgecom.logging_quantities import ( initialize_logmgr, logmgr_add_many_discretization_quantities, - logmgr_add_device_name, + logmgr_add_cl_device_info, logmgr_add_device_memory_usage, set_sim_state ) @@ -189,22 +188,15 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, local_mesh, global_nelements = generate_and_distribute_mesh(comm, gen_grid) local_nelements = local_mesh.nelements - from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD - from meshmode.discretization.poly_element import \ - default_simplex_group_factory, QuadratureSimplexGroupFactory + from mirgecom.simutil import global_reduce as _global_reduce + global_reduce = partial(_global_reduce, comm=comm) + from mirgecom.discretization import create_discretization_collection order = 3 - discr = EagerDGDiscretization( - actx, local_mesh, - discr_tag_to_group_factory={ - DISCR_TAG_BASE: default_simplex_group_factory( - base_dim=local_mesh.dim, order=order), - DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order + 1) - }, - mpi_communicator=comm - ) + discr = create_discretization_collection(actx, local_mesh, order, comm) nodes = thaw(discr.nodes(), actx) + from grudge.dof_desc import DISCR_TAG_QUAD if use_overintegration: quadrature_tag = DISCR_TAG_QUAD else: @@ -212,7 +204,7 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, dim = 2 if logmgr: - logmgr_add_device_name(logmgr, queue) + logmgr_add_cl_device_info(logmgr, queue) logmgr_add_device_memory_usage(logmgr, queue) logmgr_add_many_discretization_quantities(logmgr, discr, dim, extract_vars_for_logging, units_for_logging) @@ -275,8 +267,7 @@ def _boundary_state(discr, btag, gas_model, state_minus, **kwargs): current_cv = initializer(nodes) current_state = make_fluid_state(cv=current_cv, gas_model=gas_model) - visualizer = make_visualizer(discr, - discr.order if discr.dim == 2 else discr.order) + visualizer = make_visualizer(discr, order) initname = initializer.__class__.__name__ eosname = eos.__class__.__name__ @@ -332,9 +323,8 @@ def my_health_check(state, dv): health_error = True logger.info(f"{rank=}: NANs/Infs in pressure data.") - from mirgecom.simutil import allsync - if allsync(check_range_local(discr, "vol", dv.pressure, .9, 18.6), - comm, op=MPI.LOR): + if global_reduce(check_range_local(discr, "vol", dv.pressure, .9, 18.6), + op="lor"): health_error = True from grudge.op import nodal_max, nodal_min p_min = actx.to_numpy(nodal_min(discr, "vol", dv.pressure)) @@ -345,9 +335,9 @@ def my_health_check(state, dv): health_error = True logger.info(f"{rank=}: NANs/INFs in temperature data.") - if allsync( + if global_reduce( check_range_local(discr, "vol", dv.temperature, 2.48e-3, 1.071e-2), - comm, op=MPI.LOR): + op="lor"): health_error = True from grudge.op import nodal_max, nodal_min t_min = actx.to_numpy(nodal_min(discr, "vol", dv.temperature)) @@ -370,9 +360,9 @@ def my_pre_step(step, t, dt, state): do_health = check_step(step=step, interval=nhealth) if do_health: - from mirgecom.simutil import allsync - health_errors = allsync(my_health_check(state, dv), comm, - op=MPI.LOR) + health_errors = \ + global_reduce(my_health_check(state, dv), op="lor") + if health_errors: if rank == 0: logger.info("Fluid solution failed health check.") diff --git a/mirgecom/discretization.py b/mirgecom/discretization.py index 9c356b6df..63a855aed 100644 --- a/mirgecom/discretization.py +++ b/mirgecom/discretization.py @@ -40,7 +40,7 @@ # when we want to change it. # TODO: Make this return an actual grudge `DiscretizationCollection` # when we are ready to change mirgecom to support that change. -def create_discretization_collection(actx, mesh, order): +def create_discretization_collection(actx, mesh, order, comm=None): """Create and return a grudge DG discretization object.""" from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD from grudge.eager import EagerDGDiscretization @@ -52,6 +52,7 @@ def create_discretization_collection(actx, mesh, order): discr_tag_to_group_factory={ DISCR_TAG_BASE: PolynomialWarpAndBlendGroupFactory(order), DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(3*order), - } + }, + mpi_communicator=comm ) return discr From a9743ee864eb88e9f8e7954f2ef4366c1615a736 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Tue, 7 Jun 2022 08:22:23 -0500 Subject: [PATCH 03/21] Sharpen the test comments and simplify just a bit. --- test/test_av.py | 86 +++++++++++++++++++++++-------------------------- 1 file changed, 41 insertions(+), 45 deletions(-) diff --git a/test/test_av.py b/test/test_av.py index 4d73d4587..bb62cbfea 100644 --- a/test/test_av.py +++ b/test/test_av.py @@ -62,8 +62,13 @@ def test_tag_cells(ctx_factory, dim, order): """Test tag_cells. - Tests that the cells tagging properly tags cells - given prescribed solutions. + This test checks that tag_cells properly tags cells near + discontinuous or nearly-discontinuous features in 1d test solutions. + The following tests/functions are used: + - Discontinuity detection/Heaviside step discontinuity + - Detection smoothness/Element basis functions for each mode + - Detection thresholding/(p, p-1) polynomials greater/less/equal to s0 + - Detection bounds checking for s_e = s_0 +/- delta """ cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) @@ -87,10 +92,9 @@ def norm_indicator(expected, discr, soln, **kwargs): nele = mesh.nelements zeros = 0.0*nodes[0] - # test jump discontinuity + # Test jump discontinuity soln = actx.np.where(nodes[0] > 0.0+zeros, 1.0+zeros, zeros) err = norm_indicator(1.0, discr, soln) - assert err < tolerance, "Jump discontinuity should trigger indicator (1.0)" # get meshmode polynomials @@ -157,11 +161,13 @@ def norm_indicator(expected, discr, soln, **kwargs): err = norm_indicator(0.0, discr, soln, s0=s0+kappa-shift, kappa=kappa) assert err > tolerance, "s_e>s_0-kappa should trigger indicator" - # lower bound + # upper bound err = norm_indicator(1.0, discr, soln, s0=s0-(kappa+shift), kappa=kappa) - assert err < tolerance, "s_e>s_0+kappa should fully trigger indicator (1.0)" + # s_e>s_0+kappa should fully trigger indicator (1.0) + assert err < tolerance err = norm_indicator(1.0, discr, soln, s0=s0-(kappa-shift), kappa=kappa) - assert err > tolerance, "s_e tolerance @pytest.mark.parametrize("dim", [1, 2, 3]) @@ -169,24 +175,25 @@ def norm_indicator(expected, discr, soln, **kwargs): def test_artificial_viscosity(ctx_factory, dim, order): """Test artificial_viscosity. - Tests the application on a few simple functions - to confirm artificial viscosity returns the analytical result. + Test AV operator for some test functions to verify artificial viscosity + returns the analytical result. + - 1d x^n (expected rhs = n*(n-1)*x^(n-2) + - x^2 + y^2 + z^2 (expected rhs = 2*dim) """ cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) nel_1d = 10 - tolerance = 1.e-8 + tolerance = 1.e-6 from meshmode.mesh.generation import generate_regular_rect_mesh mesh = generate_regular_rect_mesh( - a=(-1.0, )*dim, b=(1.0, )*dim, nelements_per_axis=(nel_1d, )*dim + a=(1.0, )*dim, b=(2.0, )*dim, nelements_per_axis=(nel_1d, )*dim ) discr = EagerDGDiscretization(actx, mesh, order=order) nodes = thaw(actx, discr.nodes()) - zeros = discr.zeros(actx) class TestBoundary: def cv_gradient_flux(self, disc, btag, state_minus, gas_model, **kwargs): @@ -215,40 +222,29 @@ def av_flux(self, disc, btag, diffusion, **kwargs): return disc.project(btag, "all_faces", flux_weak) boundaries = {BTAG_ALL: TestBoundary()} - # Uniform field return 0 rhs - soln = zeros + 1.0 - cv = make_conserved( - dim, - mass=soln, - energy=soln, - momentum=make_obj_array([soln for _ in range(dim)]), - species_mass=make_obj_array([soln for _ in range(dim)]) - ) - gas_model = GasModel(eos=IdealSingleGas()) - fluid_state = make_fluid_state(cv=cv, gas_model=gas_model) - rhs = av_laplacian_operator(discr, boundaries=boundaries, - gas_model=gas_model, - fluid_state=fluid_state, alpha=1.0, s0=-np.inf) - err = discr.norm(rhs, np.inf) - assert err < tolerance - - # Linear field return 0 rhs - soln = nodes[0] - cv = make_conserved( - dim, - mass=soln, - energy=soln, - momentum=make_obj_array([soln for _ in range(dim)]), - species_mass=make_obj_array([soln for _ in range(dim)]) - ) - fluid_state = make_fluid_state(cv=cv, gas_model=gas_model) - rhs = av_laplacian_operator(discr, boundaries=boundaries, - gas_model=gas_model, - fluid_state=fluid_state, alpha=1.0, s0=-np.inf) - err = discr.norm(rhs, np.inf) - assert err < tolerance - # Quadratic field return constant 2 + # Up to quadratic 1d + for nsol in range(3): + soln = nodes[0]**nsol + exp_rhs_1d = nsol * (nsol-1) * nodes[0]**(nsol - 2) + print(f"{nsol=},{soln=},{exp_rhs_1d=}") + cv = make_conserved( + dim, + mass=soln, + energy=soln, + momentum=make_obj_array([soln for _ in range(dim)]), + species_mass=make_obj_array([soln for _ in range(dim)]) + ) + gas_model = GasModel(eos=IdealSingleGas()) + fluid_state = make_fluid_state(cv=cv, gas_model=gas_model) + rhs = av_laplacian_operator(discr, boundaries=boundaries, + gas_model=gas_model, + fluid_state=fluid_state, alpha=1.0, s0=-np.inf) + print(f"{rhs=}") + err = discr.norm(rhs-exp_rhs_1d, np.inf) + assert err < tolerance + + # Quadratic field return constant 2*dim soln = np.dot(nodes, nodes) cv = make_conserved( dim, From 29115636100b5af2c5523a0eff6aba3e0ce47bca Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Tue, 7 Jun 2022 10:43:06 -0500 Subject: [PATCH 04/21] Replace 'an exact' --> 'a' MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andreas Klöckner --- mirgecom/initializers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mirgecom/initializers.py b/mirgecom/initializers.py index d28b7edab..4216560e1 100644 --- a/mirgecom/initializers.py +++ b/mirgecom/initializers.py @@ -327,7 +327,7 @@ def __call__(self, x_vec, *, eos=None, time=0, **kwargs): At times $t > 0$, calls to this routine create an advanced solution under the assumption of constant normal shock speed *shock_speed*. The advanced solution *is not* the exact solution, but is appropriate - for use as an exact boundary solution on the top and upstream (left) + for use as a boundary solution on the top and upstream (left) side of the domain. Parameters From e741a24cbe77705fb0d8cd0da87b595b60e6d245 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Fri, 10 Jun 2022 14:13:34 -0500 Subject: [PATCH 05/21] Move toward EagerDiscretization --> Discretization Collection --- mirgecom/boundary.py | 3 ++- mirgecom/diffusion.py | 43 +++++++++++++++++----------------- mirgecom/gas_model.py | 5 ++-- mirgecom/inviscid.py | 5 ++-- mirgecom/io.py | 8 +++---- mirgecom/logging_quantities.py | 16 ++++++++----- mirgecom/operators.py | 8 +++---- mirgecom/simutil.py | 2 +- mirgecom/viscous.py | 11 +++++---- mirgecom/wave.py | 15 ++++++------ 10 files changed, 63 insertions(+), 53 deletions(-) diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py index 986218ea4..872ee9dbf 100644 --- a/mirgecom/boundary.py +++ b/mirgecom/boundary.py @@ -48,6 +48,7 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from mirgecom.fluid import make_conserved from grudge.trace_pair import TracePair +import grudge.op as op from mirgecom.viscous import viscous_facial_flux_central from mirgecom.flux import num_flux_central from mirgecom.gas_model import make_fluid_state @@ -305,7 +306,7 @@ def _boundary_quantity(self, discr, btag, quantity, local=False, **kwargs): """Get a boundary quantity on local boundary, or projected to "all_faces".""" from grudge.dof_desc import as_dofdesc btag = as_dofdesc(btag) - return quantity if local else discr.project( + return quantity if local else op.project(discr, btag, btag.with_dtag("all_faces"), quantity) def _boundary_state_pair(self, discr, btag, gas_model, state_minus, **kwargs): diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py index f41743a53..1463193cd 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -41,6 +41,7 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.dof_desc import DOFDesc, as_dofdesc, DISCR_TAG_BASE from grudge.trace_pair import TracePair, interior_trace_pairs +import grudge.op as op def grad_flux(discr, u_tpair, *, quadrature_tag=DISCR_TAG_BASE): @@ -54,12 +55,12 @@ def grad_flux(discr, u_tpair, *, quadrature_tag=DISCR_TAG_BASE): normal_quad = thaw(discr.normal(dd_quad), actx) def to_quad(a): - return discr.project(dd, dd_quad, a) + return op.project(discr, dd, dd_quad, a) def flux(u, normal): return -u * normal - return discr.project(dd_quad, dd_allfaces_quad, flux( + return op.project(discr, dd_quad, dd_allfaces_quad, flux( to_quad(u_tpair.avg), normal_quad)) @@ -75,7 +76,7 @@ def diffusion_flux( normal_quad = thaw(discr.normal(dd_quad), actx) def to_quad(a): - return discr.project(dd, dd_quad, a) + return op.project(discr, dd, dd_quad, a) def flux(kappa, grad_u, normal): return -kappa * np.dot(grad_u, normal) @@ -87,7 +88,7 @@ def flux(kappa, grad_u, normal): to_quad(kappa_tpair.ext), to_quad(grad_u_tpair.ext), normal_quad) ) - return discr.project(dd_quad, dd_allfaces_quad, flux_tpair.avg) + return op.project(discr, dd_quad, dd_allfaces_quad, flux_tpair.avg) class DiffusionBoundary(metaclass=abc.ABCMeta): @@ -141,16 +142,16 @@ def __init__(self, value): def get_grad_flux( self, discr, dd, u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 - u_int = discr.project("vol", dd, u) + u_int = op.project(discr, "vol", dd, u) u_tpair = TracePair(dd, interior=u_int, exterior=2*self.value-u_int) return grad_flux(discr, u_tpair, quadrature_tag=quadrature_tag) def get_diffusion_flux( self, discr, dd, kappa, grad_u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 - kappa_int = discr.project("vol", dd, kappa) + kappa_int = op.project(discr, "vol", dd, kappa) kappa_tpair = TracePair(dd, interior=kappa_int, exterior=kappa_int) - grad_u_int = discr.project("vol", dd, grad_u) + grad_u_int = op.project(discr, "vol", dd, grad_u) grad_u_tpair = TracePair(dd, interior=grad_u_int, exterior=grad_u_int) return diffusion_flux( discr, kappa_tpair, grad_u_tpair, quadrature_tag=quadrature_tag) @@ -194,7 +195,7 @@ def __init__(self, value): def get_grad_flux( self, discr, dd, u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 - u_int = discr.project("vol", dd, u) + u_int = op.project(discr, "vol", dd, u) u_tpair = TracePair(dd, interior=u_int, exterior=u_int) return grad_flux(discr, u_tpair, quadrature_tag=quadrature_tag) @@ -208,10 +209,10 @@ def get_diffusion_flux( # spatially-varying kappa case (the other approach would result in a # grad_u_tpair that lives in the quadrature discretization; diffusion_flux # would need to be modified to accept such values). - kappa_int_quad = discr.project("vol", dd_quad, kappa) - value_quad = discr.project(dd, dd_quad, self.value) + kappa_int_quad = op.project(discr, "vol", dd_quad, kappa) + value_quad = op.project(discr, dd, dd_quad, self.value) flux_quad = -kappa_int_quad*value_quad - return discr.project(dd_quad, dd_allfaces_quad, flux_quad) + return op.project(discr, dd_quad, dd_allfaces_quad, flux_quad) class _DiffusionStateTag: @@ -268,10 +269,10 @@ def grad_operator(discr, boundaries, u, quadrature_tag=DISCR_TAG_BASE): dd_allfaces_quad = DOFDesc("all_faces", quadrature_tag) - return discr.inverse_mass( - discr.weak_grad(-u) - - # noqa: W504 - discr.face_mass( + return op.inverse_mass(discr, + op.weak_grad(discr, -u) + - 1.0 # noqa: W504 + * op.face_mass(discr, dd_allfaces_quad, sum( grad_flux(discr, u_tpair, quadrature_tag=quadrature_tag) @@ -383,13 +384,13 @@ def diffusion_operator(discr, *args, return_grad_u=False, **kwargs): grad_u = grad_operator(discr, boundaries, u, quadrature_tag=quadrature_tag) - kappa_quad = discr.project("vol", dd_quad, kappa) - grad_u_quad = discr.project("vol", dd_quad, grad_u) + kappa_quad = op.project(discr, "vol", dd_quad, kappa) + grad_u_quad = op.project(discr, "vol", dd_quad, grad_u) - diff_u = discr.inverse_mass( - discr.weak_div(dd_quad, -kappa_quad*grad_u_quad) - - # noqa: W504 - discr.face_mass( + diff_u = op.inverse_mass(discr, + op.weak_div(discr, dd_quad, -kappa_quad*grad_u_quad) + - 1. # noqa: W504 + * op.face_mass(discr, dd_allfaces_quad, sum( diffusion_flux( diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py index e6067b53e..a91b3b2e3 100644 --- a/mirgecom/gas_model.py +++ b/mirgecom/gas_model.py @@ -60,6 +60,7 @@ GasTransportVars ) from grudge.dof_desc import DOFDesc, as_dofdesc +import grudge.op as op from grudge.trace_pair import ( interior_trace_pairs, tracepair_with_discr_tag @@ -318,10 +319,10 @@ def project_fluid_state(discr, src, tgt, state, gas_model): Thermally consistent fluid state """ - cv_sd = discr.project(src, tgt, state.cv) + cv_sd = op.project(discr, src, tgt, state.cv) temperature_seed = None if state.is_mixture: - temperature_seed = discr.project(src, tgt, state.dv.temperature) + temperature_seed = op.project(discr, src, tgt, state.dv.temperature) return make_fluid_state(cv=cv_sd, gas_model=gas_model, temperature_seed=temperature_seed) diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py index 31ccd30c6..515aad2b2 100644 --- a/mirgecom/inviscid.py +++ b/mirgecom/inviscid.py @@ -41,6 +41,7 @@ import numpy as np from arraycontext import thaw +import grudge.op as op from mirgecom.fluid import make_conserved @@ -265,14 +266,14 @@ def inviscid_flux_on_element_boundary( from grudge.dof_desc import as_dofdesc def _interior_flux(state_pair): - return discr.project( + return op.project(discr, state_pair.dd, state_pair.dd.with_dtag("all_faces"), numerical_flux_func( state_pair, gas_model, thaw(discr.normal(state_pair.dd), state_pair.int.array_context))) def _boundary_flux(dd_bdry, boundary, state_minus): - return discr.project( + return op.project(discr, dd_bdry, dd_bdry.with_dtag("all_faces"), boundary.inviscid_divergence_flux( discr, dd_bdry, gas_model, state_minus=state_minus, diff --git a/mirgecom/io.py b/mirgecom/io.py index e6c7eb817..aac74db51 100644 --- a/mirgecom/io.py +++ b/mirgecom/io.py @@ -28,7 +28,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ - +from functools import partial +import grudge.op as op from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa @@ -54,9 +55,8 @@ def make_init_message(*, dim, order, dt, t_final, def make_status_message(*, discr, t, step, dt, cfl, dependent_vars): r"""Make simulation status and health message.""" dv = dependent_vars - from functools import partial - _min = partial(discr.nodal_min, "vol") - _max = partial(discr.nodal_max, "vol") + _min = partial(op.nodal_min, discr, "vol") + _max = partial(op.nodal_max, discr, "vol") statusmsg = ( f"Status: {step=} {t=}\n" f"------- P({_min(dv.pressure):.3g}, {_max(dv.pressure):.3g})\n" diff --git a/mirgecom/logging_quantities.py b/mirgecom/logging_quantities.py index 40dd2161f..70f401b9e 100644 --- a/mirgecom/logging_quantities.py +++ b/mirgecom/logging_quantities.py @@ -50,6 +50,8 @@ from typing import Optional, Callable import numpy as np +import grudge.op as oper + def initialize_logmgr(enable_logmgr: bool, filename: str = None, mode: str = "wu", @@ -102,10 +104,11 @@ def logmgr_add_device_memory_usage(logmgr: LogManager, queue: cl.CommandQueue): def logmgr_add_many_discretization_quantities(logmgr: LogManager, discr, dim, extract_vars_for_logging, units_for_logging): """Add default discretization quantities to the logmgr.""" - for op in ["min", "max", "L2_norm"]: + for reduction_op in ["min", "max", "L2_norm"]: for quantity in ["pressure", "temperature"]: logmgr.add_quantity(DiscretizationBasedQuantity( - discr, quantity, op, extract_vars_for_logging, units_for_logging)) + discr, quantity, reduction_op, extract_vars_for_logging, + units_for_logging)) for quantity in ["mass", "energy"]: logmgr.add_quantity(DiscretizationBasedQuantity( @@ -113,7 +116,8 @@ def logmgr_add_many_discretization_quantities(logmgr: LogManager, discr, dim, for d in range(dim): logmgr.add_quantity(DiscretizationBasedQuantity( - discr, "momentum", op, extract_vars_for_logging, units_for_logging, + discr, "momentum", reduction_op, extract_vars_for_logging, + units_for_logging, axis=d)) @@ -257,13 +261,13 @@ def __init__(self, discr: Discretization, quantity: str, op: str, from functools import partial if op == "min": - self._discr_reduction = partial(self.discr.nodal_min, "vol") + self._discr_reduction = partial(op.nodal_min, self.discr, "vol") self.rank_aggr = min elif op == "max": - self._discr_reduction = partial(self.discr.nodal_max, "vol") + self._discr_reduction = partial(op.nodal_max, self.discr, "vol") self.rank_aggr = max elif op == "L2_norm": - self._discr_reduction = partial(self.discr.norm, p=2) + self._discr_reduction = partial(op.norm, self.discr, p=2) self.rank_aggr = max else: raise ValueError(f"unknown operation {op}") diff --git a/mirgecom/operators.py b/mirgecom/operators.py index 240ee542e..fa398b5ce 100644 --- a/mirgecom/operators.py +++ b/mirgecom/operators.py @@ -56,9 +56,9 @@ def grad_operator(discr, dd_vol, dd_faces, u, flux): meshmode.dof_array.DOFArray or numpy.ndarray the dg gradient operator applied to *u* """ - return -discr.inverse_mass( + return -1.0 * (op.inverse_mass(discr, op.weak_local_grad(discr, dd_vol, u) - - op.face_mass(discr, dd_faces, flux) + - op.face_mass(discr, dd_faces, flux)) ) @@ -87,7 +87,7 @@ def div_operator(discr, dd_vol, dd_faces, v, flux): meshmode.dof_array.DOFArray or numpy.ndarray the dg divergence operator applied to vector-valued function(s) *v*. """ - return -discr.inverse_mass( + return -1.0 * (op.inverse_mass(discr, op.weak_local_div(discr, dd_vol, v) - - op.face_mass(discr, dd_faces, flux) + - op.face_mass(discr, dd_faces, flux)) ) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index 1cdd2d4cc..b9d7e498a 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -323,7 +323,7 @@ def componentwise_norms(discr, fields, order=np.inf): return map_array_container( partial(componentwise_norms, discr, order=order), fields) if len(fields) > 0: - return discr.norm(fields, order) + return op.norm(discr, fields, order) else: # FIXME: This work-around for #575 can go away after #569 return 0 diff --git a/mirgecom/viscous.py b/mirgecom/viscous.py index 8f5982ad7..1bf5d73d6 100644 --- a/mirgecom/viscous.py +++ b/mirgecom/viscous.py @@ -47,6 +47,8 @@ from meshmode.dof_array import DOFArray from arraycontext import thaw +import grudge.op as op + from mirgecom.fluid import ( velocity_gradient, species_mass_fraction_gradient, @@ -384,7 +386,6 @@ def viscous_flux_on_element_boundary( Time """ from grudge.dof_desc import as_dofdesc - from grudge.op import project dd_base = as_dofdesc("vol") @@ -392,7 +393,7 @@ def viscous_flux_on_element_boundary( # viscous fluxes across interior faces (including partition and periodic bnd) def _fvisc_divergence_flux_interior(state_pair, grad_cv_pair, grad_t_pair): - return discr.project( + return op.project(discr, state_pair.dd, state_pair.dd.with_dtag("all_faces"), numerical_flux_func( discr=discr, gas_model=gas_model, state_pair=state_pair, @@ -402,13 +403,13 @@ def _fvisc_divergence_flux_interior(state_pair, grad_cv_pair, grad_t_pair): def _fvisc_divergence_flux_boundary(dd_btag, boundary, state_minus): # Make sure we fields on the quadrature grid # restricted to the tag *btag* - return project( + return op.project( discr, dd_btag, dd_btag.with_dtag("all_faces"), boundary.viscous_divergence_flux( discr=discr, btag=dd_btag, gas_model=gas_model, state_minus=state_minus, - grad_cv_minus=project(discr, dd_base, dd_btag, grad_cv), - grad_t_minus=project(discr, dd_base, dd_btag, grad_t), + grad_cv_minus=op.project(discr, dd_base, dd_btag, grad_cv), + grad_t_minus=op.project(discr, dd_base, dd_btag, grad_t), time=time, numerical_flux_func=numerical_flux_func)) # }}} viscous flux helpers diff --git a/mirgecom/wave.py b/mirgecom/wave.py index d62b98f3f..8dc1803c9 100644 --- a/mirgecom/wave.py +++ b/mirgecom/wave.py @@ -34,6 +34,7 @@ from meshmode.dof_array import thaw from grudge.trace_pair import TracePair from grudge.eager import interior_trace_pair, cross_rank_trace_pairs +import grudge.op as op def _flux(discr, c, w_tpair): @@ -55,7 +56,7 @@ def _flux(discr, c, w_tpair): 0.5*normal*np.dot(normal, v.ext-v.int), ) - return discr.project(w_tpair.dd, "all_faces", c*flux_weak) + return op.project(discr, w_tpair.dd, "all_faces", c*flux_weak) class _WaveTag: @@ -82,19 +83,19 @@ def wave_operator(discr, c, w): u = w[0] v = w[1:] - dir_u = discr.project("vol", BTAG_ALL, u) - dir_v = discr.project("vol", BTAG_ALL, v) + dir_u = op.project(discr, "vol", BTAG_ALL, u) + dir_v = op.project(discr, "vol", BTAG_ALL, v) dir_bval = flat_obj_array(dir_u, dir_v) dir_bc = flat_obj_array(-dir_u, dir_v) return ( - discr.inverse_mass( + op.inverse_mass(discr, flat_obj_array( - -c*discr.weak_div(v), - -c*discr.weak_grad(u) + -c*op.weak_div(discr, v), + -c*op.weak_grad(discr, u) ) + # noqa: W504 - discr.face_mass( + op.face_mass(discr, _flux(discr, c=c, w_tpair=interior_trace_pair(discr, w)) + _flux(discr, c=c, w_tpair=TracePair(BTAG_ALL, interior=dir_bval, exterior=dir_bc)) From ce8e9e098df178f0295d7ae1e2623fe1b2989015 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Fri, 10 Jun 2022 14:28:31 -0500 Subject: [PATCH 06/21] Fix op --> oper instances. --- mirgecom/logging_quantities.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/mirgecom/logging_quantities.py b/mirgecom/logging_quantities.py index 70f401b9e..21c9b9101 100644 --- a/mirgecom/logging_quantities.py +++ b/mirgecom/logging_quantities.py @@ -112,7 +112,8 @@ def logmgr_add_many_discretization_quantities(logmgr: LogManager, discr, dim, for quantity in ["mass", "energy"]: logmgr.add_quantity(DiscretizationBasedQuantity( - discr, quantity, op, extract_vars_for_logging, units_for_logging)) + discr, quantity, reduction_op, extract_vars_for_logging, + units_for_logging)) for d in range(dim): logmgr.add_quantity(DiscretizationBasedQuantity( @@ -261,13 +262,13 @@ def __init__(self, discr: Discretization, quantity: str, op: str, from functools import partial if op == "min": - self._discr_reduction = partial(op.nodal_min, self.discr, "vol") + self._discr_reduction = partial(oper.nodal_min, self.discr, "vol") self.rank_aggr = min elif op == "max": - self._discr_reduction = partial(op.nodal_max, self.discr, "vol") + self._discr_reduction = partial(oper.nodal_max, self.discr, "vol") self.rank_aggr = max elif op == "L2_norm": - self._discr_reduction = partial(op.norm, self.discr, p=2) + self._discr_reduction = partial(oper.norm, self.discr, p=2) self.rank_aggr = max else: raise ValueError(f"unknown operation {op}") From 658183fc8f6c5ab46a303a9809b960f0d23cd31b Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Fri, 10 Jun 2022 14:45:26 -0500 Subject: [PATCH 07/21] Fix up weak_{grad, div} --> weak_local_{grad, div} --- mirgecom/diffusion.py | 5 +++-- mirgecom/wave.py | 4 ++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py index 1463193cd..7ec5615da 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -268,9 +268,10 @@ def grad_operator(discr, boundaries, u, quadrature_tag=DISCR_TAG_BASE): "Must be an instance of DiffusionBoundary.") dd_allfaces_quad = DOFDesc("all_faces", quadrature_tag) + dd_vol_quad = DOFDesc("vol", quadrature_tag) return op.inverse_mass(discr, - op.weak_grad(discr, -u) + op.weak_local_grad(discr, dd_vol_quad, -u) - 1.0 # noqa: W504 * op.face_mass(discr, dd_allfaces_quad, @@ -388,7 +389,7 @@ def diffusion_operator(discr, *args, return_grad_u=False, **kwargs): grad_u_quad = op.project(discr, "vol", dd_quad, grad_u) diff_u = op.inverse_mass(discr, - op.weak_div(discr, dd_quad, -kappa_quad*grad_u_quad) + op.weak_local_div(discr, dd_quad, -kappa_quad*grad_u_quad) - 1. # noqa: W504 * op.face_mass(discr, dd_allfaces_quad, diff --git a/mirgecom/wave.py b/mirgecom/wave.py index 8dc1803c9..6c89dc01c 100644 --- a/mirgecom/wave.py +++ b/mirgecom/wave.py @@ -91,8 +91,8 @@ def wave_operator(discr, c, w): return ( op.inverse_mass(discr, flat_obj_array( - -c*op.weak_div(discr, v), - -c*op.weak_grad(discr, u) + -c*op.weak_local_div(discr, "vol", v), + -c*op.weak_local_grad(discr, "vol", u) ) + # noqa: W504 op.face_mass(discr, From 884ebcdac56bc02a9c67b4e5f51e5a6b170dc556 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Fri, 10 Jun 2022 15:38:51 -0500 Subject: [PATCH 08/21] Use proper dofdescs --- examples/wave.py | 7 ++++--- mirgecom/wave.py | 5 +++-- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/examples/wave.py b/examples/wave.py index f575418d0..607af2951 100644 --- a/examples/wave.py +++ b/examples/wave.py @@ -32,7 +32,8 @@ from grudge.eager import EagerDGDiscretization from grudge.shortcuts import make_visualizer - +from grudge.op import nodal_min +from grudge.dof_desc import as_dofdesc from mirgecom.wave import wave_operator from mirgecom.integrators import rk4_step @@ -108,8 +109,8 @@ def main(use_profiling=False, use_logmgr=False, lazy_eval: bool = False): wave_speed = 1.0 from grudge.dt_utils import characteristic_lengthscales nodal_dt = characteristic_lengthscales(actx, discr) / wave_speed - from grudge.op import nodal_min - dt = actx.to_numpy(current_cfl * nodal_min(discr, "vol", nodal_dt))[()] + dt = actx.to_numpy(current_cfl * nodal_min(discr, as_dofdesc("vol"), + nodal_dt))[()] print("%d elements" % mesh.nelements) diff --git a/mirgecom/wave.py b/mirgecom/wave.py index 6c89dc01c..7d501b73b 100644 --- a/mirgecom/wave.py +++ b/mirgecom/wave.py @@ -33,6 +33,7 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from meshmode.dof_array import thaw from grudge.trace_pair import TracePair +from grudge.dof_desc import as_dofdesc from grudge.eager import interior_trace_pair, cross_rank_trace_pairs import grudge.op as op @@ -91,8 +92,8 @@ def wave_operator(discr, c, w): return ( op.inverse_mass(discr, flat_obj_array( - -c*op.weak_local_div(discr, "vol", v), - -c*op.weak_local_grad(discr, "vol", u) + -c*op.weak_local_div(discr, as_dofdesc("vol"), v), + -c*op.weak_local_grad(discr, as_dofdesc("vol"), u) ) + # noqa: W504 op.face_mass(discr, From c7a5677a15649094d7c5f075b6abab48157a857d Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Fri, 10 Jun 2022 19:32:55 -0500 Subject: [PATCH 09/21] Use base domain, not quad domain --- mirgecom/diffusion.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py index 7ec5615da..18eb31570 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -268,10 +268,9 @@ def grad_operator(discr, boundaries, u, quadrature_tag=DISCR_TAG_BASE): "Must be an instance of DiffusionBoundary.") dd_allfaces_quad = DOFDesc("all_faces", quadrature_tag) - dd_vol_quad = DOFDesc("vol", quadrature_tag) return op.inverse_mass(discr, - op.weak_local_grad(discr, dd_vol_quad, -u) + op.weak_local_grad(discr, as_dofdesc("vol"), -u) - 1.0 # noqa: W504 * op.face_mass(discr, dd_allfaces_quad, From 732f005b632bd303fb8bbc61669a95193ce1379e Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Sun, 12 Jun 2022 09:24:04 -0500 Subject: [PATCH 10/21] Excise unneeded array value arg to actx.np.where --- mirgecom/artificial_viscosity.py | 45 ++++++++++++++++++++++++-------- 1 file changed, 34 insertions(+), 11 deletions(-) diff --git a/mirgecom/artificial_viscosity.py b/mirgecom/artificial_viscosity.py index 3a429ce6a..52f72cdcc 100644 --- a/mirgecom/artificial_viscosity.py +++ b/mirgecom/artificial_viscosity.py @@ -314,27 +314,50 @@ def highest_mode(grp): uhat = modal_map(u) # Compute smoothness indicator value - indicator = DOFArray( - actx, - data=tuple( - actx.call_loopy( - indicator_prg(), - vec=uhat[grp.index], - modes_active_flag=highest_mode(grp) - )["result"] - for grp in discr.discr_from_dd("vol").groups + if actx.supports_nonscalar_broadcasting: + from meshmode.transform_metadata import DiscretizationDOFAxisTag + indicator = DOFArray( + actx, + data=tuple( + actx.tag_axis( + 1, + DiscretizationDOFAxisTag(), + actx.np.broadcast_to( + ((actx.einsum("ek,k->e", + uhat[grp.index]**2, + highest_mode(grp)) + / (actx.einsum("ej->e", + (uhat[grp.index]**2+(1e-12/grp.nunit_dofs)) + ))) + .reshape(-1, 1)), + uhat[grp.index].shape)) + for grp in discr.discr_from_dd("vol").groups + ) ) - ) + else: + indicator = DOFArray( + actx, + data=tuple( + actx.call_loopy( + indicator_prg(), + vec=uhat[grp.index], + modes_active_flag=highest_mode(grp) + )["result"] + for grp in discr.discr_from_dd("vol").groups + ) + ) + indicator = actx.np.log10(indicator + 1.0e-12) # Compute artificial viscosity percentage based on indicator and set parameters yesnol = actx.np.greater(indicator, (s0 - kappa)) yesnou = actx.np.greater(indicator, (s0 + kappa)) + saintly_value = 1.0 sin_indicator = actx.np.where( yesnol, 0.5 * (1.0 + actx.np.sin(np.pi * (indicator - s0) / (2.0 * kappa))), 0.0 * indicator, ) - indicator = actx.np.where(yesnou, 1.0 + 0.0 * indicator, sin_indicator) + indicator = actx.np.where(yesnou, saintly_value, sin_indicator) return indicator From 2f9c0d3134e6ef954ce50bec1395c1a0a035df31 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Mon, 13 Jun 2022 08:29:32 -0500 Subject: [PATCH 11/21] Customize production env. --- .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..6ddf8f69e 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="move-toward-dcoll-production" # 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 3721471a6294d509ca5af27cddc8ee61c625cb27 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 13 Jun 2022 14:47:47 -0500 Subject: [PATCH 12/21] Notation fix in div/grad --- examples/wave.py | 2 +- mirgecom/operators.py | 6 ++---- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/examples/wave.py b/examples/wave.py index 607af2951..26e9708aa 100644 --- a/examples/wave.py +++ b/examples/wave.py @@ -109,7 +109,7 @@ def main(use_profiling=False, use_logmgr=False, lazy_eval: bool = False): wave_speed = 1.0 from grudge.dt_utils import characteristic_lengthscales nodal_dt = characteristic_lengthscales(actx, discr) / wave_speed - dt = actx.to_numpy(current_cfl * nodal_min(discr, as_dofdesc("vol"), + dt = actx.to_numpy(current_cfl * nodal_min(discr, "vol", nodal_dt))[()] print("%d elements" % mesh.nelements) diff --git a/mirgecom/operators.py b/mirgecom/operators.py index fa398b5ce..e677c11b8 100644 --- a/mirgecom/operators.py +++ b/mirgecom/operators.py @@ -56,10 +56,9 @@ def grad_operator(discr, dd_vol, dd_faces, u, flux): meshmode.dof_array.DOFArray or numpy.ndarray the dg gradient operator applied to *u* """ - return -1.0 * (op.inverse_mass(discr, + return - op.inverse_mass(discr, op.weak_local_grad(discr, dd_vol, u) - op.face_mass(discr, dd_faces, flux)) - ) def div_operator(discr, dd_vol, dd_faces, v, flux): @@ -87,7 +86,6 @@ def div_operator(discr, dd_vol, dd_faces, v, flux): meshmode.dof_array.DOFArray or numpy.ndarray the dg divergence operator applied to vector-valued function(s) *v*. """ - return -1.0 * (op.inverse_mass(discr, + return - op.inverse_mass(discr, op.weak_local_div(discr, dd_vol, v) - op.face_mass(discr, dd_faces, flux)) - ) From 1bf3fe3a87151577be72a758e5cc1f111dc943be Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Mon, 13 Jun 2022 15:20:39 -0500 Subject: [PATCH 13/21] Sharpen and clean up some superfluous patterns. --- examples/wave.py | 1 - mirgecom/diffusion.py | 2 +- mirgecom/operators.py | 2 ++ mirgecom/wave.py | 5 ++--- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/wave.py b/examples/wave.py index 26e9708aa..1d15d494f 100644 --- a/examples/wave.py +++ b/examples/wave.py @@ -33,7 +33,6 @@ from grudge.eager import EagerDGDiscretization from grudge.shortcuts import make_visualizer from grudge.op import nodal_min -from grudge.dof_desc import as_dofdesc from mirgecom.wave import wave_operator from mirgecom.integrators import rk4_step diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py index 18eb31570..6b7286e38 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -270,7 +270,7 @@ def grad_operator(discr, boundaries, u, quadrature_tag=DISCR_TAG_BASE): dd_allfaces_quad = DOFDesc("all_faces", quadrature_tag) return op.inverse_mass(discr, - op.weak_local_grad(discr, as_dofdesc("vol"), -u) + op.weak_local_grad(discr, "vol", -u) - 1.0 # noqa: W504 * op.face_mass(discr, dd_allfaces_quad, diff --git a/mirgecom/operators.py b/mirgecom/operators.py index e677c11b8..b584bb5d8 100644 --- a/mirgecom/operators.py +++ b/mirgecom/operators.py @@ -56,6 +56,7 @@ def grad_operator(discr, dd_vol, dd_faces, u, flux): meshmode.dof_array.DOFArray or numpy.ndarray the dg gradient operator applied to *u* """ + # pylint: disable=invalid-unary-operand-type return - op.inverse_mass(discr, op.weak_local_grad(discr, dd_vol, u) - op.face_mass(discr, dd_faces, flux)) @@ -86,6 +87,7 @@ def div_operator(discr, dd_vol, dd_faces, v, flux): meshmode.dof_array.DOFArray or numpy.ndarray the dg divergence operator applied to vector-valued function(s) *v*. """ + # pylint: disable=invalid-unary-operand-type return - op.inverse_mass(discr, op.weak_local_div(discr, dd_vol, v) - op.face_mass(discr, dd_faces, flux)) diff --git a/mirgecom/wave.py b/mirgecom/wave.py index 7d501b73b..6c89dc01c 100644 --- a/mirgecom/wave.py +++ b/mirgecom/wave.py @@ -33,7 +33,6 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from meshmode.dof_array import thaw from grudge.trace_pair import TracePair -from grudge.dof_desc import as_dofdesc from grudge.eager import interior_trace_pair, cross_rank_trace_pairs import grudge.op as op @@ -92,8 +91,8 @@ def wave_operator(discr, c, w): return ( op.inverse_mass(discr, flat_obj_array( - -c*op.weak_local_div(discr, as_dofdesc("vol"), v), - -c*op.weak_local_grad(discr, as_dofdesc("vol"), u) + -c*op.weak_local_div(discr, "vol", v), + -c*op.weak_local_grad(discr, "vol", u) ) + # noqa: W504 op.face_mass(discr, From 82f617cb297f8206d9364b2fba0326901a93769d Mon Sep 17 00:00:00 2001 From: Luke Olson Date: Mon, 13 Jun 2022 15:36:57 -0500 Subject: [PATCH 14/21] add sin test for convergence --- test/test_av.py | 81 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/test/test_av.py b/test/test_av.py index bb62cbfea..7ef59436f 100644 --- a/test/test_av.py +++ b/test/test_av.py @@ -259,3 +259,84 @@ def av_flux(self, disc, btag, diffusion, **kwargs): fluid_state=fluid_state, alpha=1.0, s0=-np.inf) err = discr.norm(2.*dim-rhs, np.inf) assert err < tolerance + +@pytest.mark.parametrize("order", [1, 2, 3]) +@pytest.mark.parametrize("dim", [1, 2]) +def test_mms(ctx_factory, dim, order): + """Test artificial_viscosity. + + Test AV operator for some test functions to verify artificial viscosity + returns the analytical result. + """ + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + #nel_1d = 10 + olderr = None + print('\n') + for nel_1d in [8, 16, 32, 64]: + tolerance = 1.e-6 + + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(0.0, )*dim, b=(1.0, )*dim, nelements_per_axis=(nel_1d, )*dim + ) + + discr = EagerDGDiscretization(actx, mesh, order=order) + nodes = thaw(actx, discr.nodes()) + + class TestBoundary: + def cv_gradient_flux(self, disc, btag, state_minus, gas_model, **kwargs): + fluid_state_int = project_fluid_state(disc, "vol", btag, fluid_state, + gas_model) + cv_int = fluid_state_int.cv + from grudge.trace_pair import TracePair + bnd_pair = TracePair(btag, + interior=cv_int, + exterior=cv_int) + nhat = thaw(actx, disc.normal(btag)) + from mirgecom.flux import num_flux_central + from arraycontext import outer + # Do not project to "all_faces" as now we use built-in grad_cv_operator + return outer(num_flux_central(bnd_pair.int, bnd_pair.ext), nhat) + + def av_flux(self, disc, btag, diffusion, **kwargs): + nhat = thaw(actx, disc.normal(btag)) + diffusion_minus = discr.project("vol", btag, diffusion) + diffusion_plus = diffusion_minus + from grudge.trace_pair import TracePair + bnd_grad_pair = TracePair(btag, interior=diffusion_minus, + exterior=diffusion_plus) + from mirgecom.flux import num_flux_central + flux_weak = num_flux_central(bnd_grad_pair.int, bnd_grad_pair.ext)@nhat + return disc.project(btag, "all_faces", flux_weak) + + boundaries = {BTAG_ALL: TestBoundary()} + + # u = sin(pi x)sin(pi y)sin(pi z) + # rhs = -dim pi^2 u + soln = 1.0 + exp_rhs = -dim * np.pi**2 + for w in nodes: + soln *= actx.np.sin(w * np.pi) + exp_rhs *= actx.np.sin(w * np.pi) + + gas_model = GasModel(eos=IdealSingleGas()) + cv = make_conserved( + dim, + mass=soln, + energy=soln, + momentum=make_obj_array([soln for _ in range(dim)]), + species_mass=make_obj_array([soln for _ in range(dim)]) + ) + fluid_state = make_fluid_state(cv=cv, gas_model=gas_model) + rhs = av_laplacian_operator(discr, boundaries=boundaries, + gas_model=gas_model, + fluid_state=fluid_state, alpha=1.0, s0=-np.inf) + + err = discr.norm(rhs-exp_rhs, np.inf) + if olderr is not None: + rate = olderr / err + print(f'{dim=}, {nel_1d=}, {order=}, {rate=}, {err=}') + olderr = err From 68b70d259044531ef5ea00ef1e9f1378c0b08c49 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Tue, 14 Jun 2022 10:09:45 -0500 Subject: [PATCH 15/21] Add trig test from @lukeo slightly twisted. --- mirgecom/artificial_viscosity.py | 31 +++++++++++++ test/test_av.py | 75 +++++++++++++------------------- 2 files changed, 61 insertions(+), 45 deletions(-) diff --git a/mirgecom/artificial_viscosity.py b/mirgecom/artificial_viscosity.py index 52f72cdcc..d8ce2560d 100644 --- a/mirgecom/artificial_viscosity.py +++ b/mirgecom/artificial_viscosity.py @@ -59,6 +59,37 @@ - $s_0$ is a reference smoothness value - $\kappa$ controls the width of the transition between 0 to 1 + +Boundary Conditions +^^^^^^^^^^^^^^^^^^^ + +The artificial viscosity operator as currently implemented re-uses the fluid +solution gradient $\nabla{\mathbf{Q}}$ for the auxiliary equation: + +.. math:: + +\mathbf{R} = \varepsilon\nabla\mathbf{Q}_\text{fluid} + +As such, the fluid-system imposes the appropriate boundary solution $\mathbf{Q}^+$ +for the comptuation of $\nabla{\mathbf{Q}}$. This approach leaves the boundary +condition on $\mathbf{R}$ to be imposed by boundary treatment for the operator when +computing the divergence for the RHS, $\nabla \cdot \mathbf{R}$. + +Similar to the fluid boundary treatments; when no boundary conditions are imposed +on $\mathbf{R}$, the interior solution is simply extrapolated to the boundary, +(i.e., $\mathbf{R}^+ = \mathbf{R}^-). If such a boundary condition is imposed, +usually for selected components of $\mathbf{R}$, then such boundary conditions +are used directly: $\mathbf{R}^+ = \mathbf{R}_\text{bc}$. + +A central numerical flux is then employed to transmit the boundary condition to +the domain for the divergence operator: + +.. math:: + + \mathbf{R} \cdot \hat{mathbf{n}} = \frac{1}{2}\left(\frac{\mathbf{R}^- + + \mathbf{R}^+}\right) \cdot \hat{\mathbf{n}} + + Smoothness Indicator Evaluation ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/test/test_av.py b/test/test_av.py index 7ef59436f..8498ac154 100644 --- a/test/test_av.py +++ b/test/test_av.py @@ -260,9 +260,10 @@ def av_flux(self, disc, btag, diffusion, **kwargs): err = discr.norm(2.*dim-rhs, np.inf) assert err < tolerance -@pytest.mark.parametrize("order", [1, 2, 3]) + +@pytest.mark.parametrize("order", [2, 3]) @pytest.mark.parametrize("dim", [1, 2]) -def test_mms(ctx_factory, dim, order): +def test_trig(ctx_factory, dim, order): """Test artificial_viscosity. Test AV operator for some test functions to verify artificial viscosity @@ -272,55 +273,31 @@ def test_mms(ctx_factory, dim, order): queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) - #nel_1d = 10 - olderr = None - print('\n') - for nel_1d in [8, 16, 32, 64]: - tolerance = 1.e-6 + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + nel_1d_base = 4 + + for nfac in [1, 2, 4]: + nel_1d = nfac * nel_1d_base from meshmode.mesh.generation import generate_regular_rect_mesh mesh = generate_regular_rect_mesh( - a=(0.0, )*dim, b=(1.0, )*dim, nelements_per_axis=(nel_1d, )*dim + a=(0.0, )*dim, b=(1.0, )*dim, nelements_per_axis=(nel_1d, )*dim, + periodic=(True,)*dim ) discr = EagerDGDiscretization(actx, mesh, order=order) nodes = thaw(actx, discr.nodes()) - class TestBoundary: - def cv_gradient_flux(self, disc, btag, state_minus, gas_model, **kwargs): - fluid_state_int = project_fluid_state(disc, "vol", btag, fluid_state, - gas_model) - cv_int = fluid_state_int.cv - from grudge.trace_pair import TracePair - bnd_pair = TracePair(btag, - interior=cv_int, - exterior=cv_int) - nhat = thaw(actx, disc.normal(btag)) - from mirgecom.flux import num_flux_central - from arraycontext import outer - # Do not project to "all_faces" as now we use built-in grad_cv_operator - return outer(num_flux_central(bnd_pair.int, bnd_pair.ext), nhat) - - def av_flux(self, disc, btag, diffusion, **kwargs): - nhat = thaw(actx, disc.normal(btag)) - diffusion_minus = discr.project("vol", btag, diffusion) - diffusion_plus = diffusion_minus - from grudge.trace_pair import TracePair - bnd_grad_pair = TracePair(btag, interior=diffusion_minus, - exterior=diffusion_plus) - from mirgecom.flux import num_flux_central - flux_weak = num_flux_central(bnd_grad_pair.int, bnd_grad_pair.ext)@nhat - return disc.project(btag, "all_faces", flux_weak) - - boundaries = {BTAG_ALL: TestBoundary()} - - # u = sin(pi x)sin(pi y)sin(pi z) + boundaries = {} + + # u = sin(2 pi x)sin(2 pi y)sin(2 pi z) # rhs = -dim pi^2 u soln = 1.0 - exp_rhs = -dim * np.pi**2 + exp_rhs = -dim * 4. * np.pi**2 for w in nodes: - soln *= actx.np.sin(w * np.pi) - exp_rhs *= actx.np.sin(w * np.pi) + soln = soln * actx.np.sin(w * 2 * np.pi) + exp_rhs = exp_rhs * actx.np.sin(w * 2 * np.pi) gas_model = GasModel(eos=IdealSingleGas()) cv = make_conserved( @@ -335,8 +312,16 @@ def av_flux(self, disc, btag, diffusion, **kwargs): gas_model=gas_model, fluid_state=fluid_state, alpha=1.0, s0=-np.inf) - err = discr.norm(rhs-exp_rhs, np.inf) - if olderr is not None: - rate = olderr / err - print(f'{dim=}, {nel_1d=}, {order=}, {rate=}, {err=}') - olderr = err + err_rhs = actx.to_numpy(discr.norm(rhs-exp_rhs, np.inf)) + eoc_rec.add_data_point(1.0/nel_1d, err_rhs) + + logger.info( + f"{dim=}, {order=}" + f"Errors:\n{eoc_rec}" + ) + # For RHS-only checks we expect order - 1. + expected_order = order - 1 + assert ( + eoc_rec.order_estimate() >= expected_order - .5 + or eoc_rec.max_error() < 1e-9 + ) From cb23f1f28e1dbbabfb568a6483b92a27f2f6b686 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Tue, 14 Jun 2022 10:36:31 -0500 Subject: [PATCH 16/21] Sharpen docs, note limitations and AV replacement. --- mirgecom/artificial_viscosity.py | 4 ++-- test/test_av.py | 9 +++++++++ 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/mirgecom/artificial_viscosity.py b/mirgecom/artificial_viscosity.py index d8ce2560d..1aa9daef9 100644 --- a/mirgecom/artificial_viscosity.py +++ b/mirgecom/artificial_viscosity.py @@ -68,7 +68,7 @@ .. math:: -\mathbf{R} = \varepsilon\nabla\mathbf{Q}_\text{fluid} + \mathbf{R} = \varepsilon\nabla\mathbf{Q}_\text{fluid} As such, the fluid-system imposes the appropriate boundary solution $\mathbf{Q}^+$ for the comptuation of $\nabla{\mathbf{Q}}$. This approach leaves the boundary @@ -77,7 +77,7 @@ Similar to the fluid boundary treatments; when no boundary conditions are imposed on $\mathbf{R}$, the interior solution is simply extrapolated to the boundary, -(i.e., $\mathbf{R}^+ = \mathbf{R}^-). If such a boundary condition is imposed, +(i.e., $\mathbf{R}^+ = \mathbf{R}^-$). If such a boundary condition is imposed, usually for selected components of $\mathbf{R}$, then such boundary conditions are used directly: $\mathbf{R}^+ = \mathbf{R}_\text{bc}$. diff --git a/test/test_av.py b/test/test_av.py index 8498ac154..383a6cbbb 100644 --- a/test/test_av.py +++ b/test/test_av.py @@ -57,6 +57,15 @@ logger = logging.getLogger(__name__) +# NOTE: Testing of this av_laplacian_operator is currently +# pretty limited. This fact is somewhat indicative of the +# limitations and shortcomings of this operator. We intend +# to soon replace our shock-handling approach with one that is +# more robust in the presence of discontinuous coeffcients +# and in which we understand the required boundary conditions. +# Tracking the replacement endeavor: +# https://github.com/illinois-ceesd/mirgecom/issues/684 + @pytest.mark.parametrize("dim", [1, 2, 3]) @pytest.mark.parametrize("order", [1, 5]) def test_tag_cells(ctx_factory, dim, order): From bfa11437e0240c77e4b06bac28ad4eac978aaa1e Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Tue, 14 Jun 2022 10:50:47 -0500 Subject: [PATCH 17/21] debug latex --- mirgecom/artificial_viscosity.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mirgecom/artificial_viscosity.py b/mirgecom/artificial_viscosity.py index 1aa9daef9..4537be041 100644 --- a/mirgecom/artificial_viscosity.py +++ b/mirgecom/artificial_viscosity.py @@ -86,8 +86,8 @@ .. math:: - \mathbf{R} \cdot \hat{mathbf{n}} = \frac{1}{2}\left(\frac{\mathbf{R}^- - + \mathbf{R}^+}\right) \cdot \hat{\mathbf{n}} + \mathbf{R} \cdot \hat{mathbf{n}} = \frac{1}{2}\left(\mathbf{R}^- + + \mathbf{R}^+\right) \cdot \hat{\mathbf{n}} Smoothness Indicator Evaluation From c101096af1461770d60f98e2d8e590cdd2acffd8 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Tue, 14 Jun 2022 14:38:08 -0500 Subject: [PATCH 18/21] Add boundary handling test --- mirgecom/boundary.py | 20 +----- test/test_av.py | 164 +++++++++++++++++++++++++++++++++++++++---- 2 files changed, 152 insertions(+), 32 deletions(-) diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py index e4a8da5d6..42748ad01 100644 --- a/mirgecom/boundary.py +++ b/mirgecom/boundary.py @@ -50,10 +50,7 @@ from grudge.trace_pair import TracePair from mirgecom.viscous import viscous_facial_flux_central from mirgecom.flux import num_flux_central -from mirgecom.gas_model import ( - make_fluid_state, - project_fluid_state -) +from mirgecom.gas_model import make_fluid_state from mirgecom.inviscid import inviscid_facial_flux_rusanov @@ -245,7 +242,6 @@ class PrescribedFluidBoundary(FluidBoundary): .. automethod:: cv_gradient_flux .. automethod:: temperature_gradient_flux .. automethod:: av_flux - .. automethod:: soln_gradient_flux """ def __init__(self, @@ -469,20 +465,6 @@ def viscous_divergence_flux(self, discr, btag, gas_model, state_minus, def _identical_grad_av(self, grad_av_minus, **kwargs): return grad_av_minus - def soln_gradient_flux(self, discr, btag, fluid_state, gas_model, **kwargs): - """Get the flux for solution gradient with AV API.""" - # project the conserved and thermal state to the boundary - fluid_state_minus = project_fluid_state(discr=discr, - src="vol", - tgt=btag, - gas_model=gas_model, - state=fluid_state) - # get the boundary flux for the grad(CV) - return self.cv_gradient_flux(discr=discr, btag=btag, - gas_model=gas_model, - state_minus=fluid_state_minus, - **kwargs) - def av_flux(self, discr, btag, diffusion, **kwargs): """Get the diffusive fluxes for the AV operator API.""" grad_av_minus = discr.project("vol", btag, diffusion) diff --git a/test/test_av.py b/test/test_av.py index 383a6cbbb..41f54eec4 100644 --- a/test/test_av.py +++ b/test/test_av.py @@ -29,9 +29,12 @@ import numpy as np import pyopencl as cl import pytest - -from meshmode.array_context import PyOpenCLArrayContext -from meshmode.dof_array import thaw +from meshmode.array_context import ( # noqa + PyOpenCLArrayContext, + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests +) +from arraycontext import thaw from meshmode.mesh import BTAG_ALL from mirgecom.artificial_viscosity import ( @@ -41,8 +44,7 @@ from mirgecom.fluid import make_conserved from mirgecom.gas_model import ( GasModel, - make_fluid_state, - project_fluid_state + make_fluid_state ) from mirgecom.eos import IdealSingleGas @@ -97,7 +99,7 @@ def norm_indicator(expected, discr, soln, **kwargs): ) discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(discr.nodes(), actx) nele = mesh.nelements zeros = 0.0*nodes[0] @@ -202,25 +204,24 @@ def test_artificial_viscosity(ctx_factory, dim, order): ) discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(discr.nodes(), actx) class TestBoundary: + def cv_gradient_flux(self, disc, btag, state_minus, gas_model, **kwargs): - fluid_state_int = project_fluid_state(disc, "vol", btag, fluid_state, - gas_model) - cv_int = fluid_state_int.cv + cv_int = state_minus.cv from grudge.trace_pair import TracePair bnd_pair = TracePair(btag, interior=cv_int, exterior=cv_int) - nhat = thaw(actx, disc.normal(btag)) + nhat = thaw(disc.normal(btag), actx) from mirgecom.flux import num_flux_central from arraycontext import outer # Do not project to "all_faces" as now we use built-in grad_cv_operator return outer(num_flux_central(bnd_pair.int, bnd_pair.ext), nhat) def av_flux(self, disc, btag, diffusion, **kwargs): - nhat = thaw(actx, disc.normal(btag)) + nhat = thaw(disc.normal(btag), actx) diffusion_minus = discr.project("vol", btag, diffusion) diffusion_plus = diffusion_minus from grudge.trace_pair import TracePair @@ -296,7 +297,7 @@ def test_trig(ctx_factory, dim, order): ) discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(discr.nodes(), actx) boundaries = {} @@ -334,3 +335,140 @@ def test_trig(ctx_factory, dim, order): eoc_rec.order_estimate() >= expected_order - .5 or eoc_rec.max_error() < 1e-9 ) + + +# Box grid generator widget lifted from @majosm's diffusion tester +def _get_box_mesh(dim, a, b, n): + dim_names = ["x", "y", "z"] + boundary_tag_to_face = {} + for i in range(dim): + boundary_tag_to_face["-"+str(i+1)] = ["-"+dim_names[i]] + boundary_tag_to_face["+"+str(i+1)] = ["+"+dim_names[i]] + from meshmode.mesh.generation import generate_regular_rect_mesh + return generate_regular_rect_mesh(a=(a,)*dim, b=(b,)*dim, n=(n,)*dim, + boundary_tag_to_face=boundary_tag_to_face) + + +class _VortexSoln: + + def __init__(self, **kwargs): + self._dim = 2 + from mirgecom.initializers import Vortex2D + origin = np.zeros(2) + velocity = origin + 1 + self._initializer = Vortex2D(center=origin, velocity=velocity) + + def dim(self): + return self._dim + + def __call__(self, r, eos, **kwargs): + return self._initializer(x_vec=r, eos=eos) + + +class _TestSoln: + + def __init__(self, dim=2): + self._dim = dim + from mirgecom.initializers import Uniform + origin = np.zeros(self._dim) + velocity = origin + 1 + self._initializer = Uniform(dim=self._dim, velocity=velocity) + + def dim(self): + return self._dim + + def __call__(self, r, eos, **kwargs): + return self._initializer(x_vec=r, eos=eos) + + +# This test makes sure that the fluid boundaries work mechanically for AV, +# and give an expected result when using the AV boundary interface. +@pytest.mark.parametrize("prescribed_soln", [_TestSoln(dim=1), + _TestSoln(dim=2), + _TestSoln(dim=3)]) +@pytest.mark.parametrize("order", [2, 3]) +def test_fluid_av_boundaries(ctx_factory, prescribed_soln, order): + """Check fluid boundary AV interface.""" + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + dim = prescribed_soln.dim() + + kappa = 3.0 + sigma = 5.0 + + from mirgecom.transport import SimpleTransport + transport_model = SimpleTransport(viscosity=sigma, thermal_conductivity=kappa) + + gas_model = GasModel(eos=IdealSingleGas(gas_const=1.0), + transport=transport_model) + + def _boundary_state_func(discr, btag, gas_model, state_minus, **kwargs): + actx = state_minus.array_context + bnd_discr = discr.discr_from_dd(btag) + nodes = thaw(bnd_discr.nodes(), actx) + return make_fluid_state(prescribed_soln(r=nodes, eos=gas_model.eos, + **kwargs), gas_model) + + npts_geom = 17 + a = 1.0 + b = 2.0 + mesh = _get_box_mesh(dim=dim, a=a, b=b, n=npts_geom) + from mirgecom.discretization import create_discretization_collection + discr = create_discretization_collection(actx, mesh, order=order) + nodes = thaw(discr.nodes(), actx) + cv = prescribed_soln(r=nodes, eos=gas_model.eos) + fluid_state = make_fluid_state(cv, gas_model) + + boundary_nhat = thaw(discr.normal(BTAG_ALL), actx) + + from mirgecom.boundary import ( + PrescribedFluidBoundary, + AdiabaticNoslipMovingBoundary, + IsothermalNoSlipBoundary + ) + prescribed_boundary = \ + PrescribedFluidBoundary(boundary_state_func=_boundary_state_func) + adiabatic_noslip = AdiabaticNoslipMovingBoundary() + isothermal_noslip = IsothermalNoSlipBoundary() + + fluid_boundaries = {BTAG_ALL: prescribed_boundary} + from mirgecom.navierstokes import grad_cv_operator + fluid_grad_cv = \ + grad_cv_operator(discr, gas_model, fluid_boundaries, fluid_state) + + # Put in a testing field for AV - doesn't matter what - as long as it + # is spatially-dependent in all dimensions. + av_diffusion = 0. * fluid_grad_cv + np.dot(nodes, nodes) + av_diffusion_boundary = discr.project("vol", BTAG_ALL, av_diffusion) + + # Prescribed boundaries are used for inflow/outflow-type boundaries + # where we expect to _preserve_ the soln gradient + from grudge.dof_desc import as_dofdesc + dd_bnd = as_dofdesc(BTAG_ALL) + all_faces_dd = dd_bnd.with_dtag("all_faces") + expected_av_flux_prescribed_boundary = av_diffusion_boundary@boundary_nhat + print(f"{expected_av_flux_prescribed_boundary=}") + exp_av_flux = discr.project(dd_bnd, all_faces_dd, + expected_av_flux_prescribed_boundary) + print(f"{exp_av_flux=}") + + prescribed_boundary_av_flux = \ + prescribed_boundary.av_flux(discr, BTAG_ALL, av_diffusion) + print(f"{prescribed_boundary_av_flux=}") + + bnd_flux_resid = (prescribed_boundary_av_flux - exp_av_flux) + print(f"{bnd_flux_resid=}") + assert bnd_flux_resid == 0 + + # Solid wall boundaries are expected to have 0 AV flux + wall_bnd_flux = \ + adiabatic_noslip.av_flux(discr, BTAG_ALL, av_diffusion) + print(f"adiabatic_noslip: {wall_bnd_flux=}") + assert wall_bnd_flux == 0 + + wall_bnd_flux = \ + isothermal_noslip.av_flux(discr, BTAG_ALL, av_diffusion) + print(f"isothermal_noslip: {wall_bnd_flux=}") + assert wall_bnd_flux == 0 From 7047c72734e923395df86db8e9ad923464f8443b Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Wed, 15 Jun 2022 08:07:15 -0500 Subject: [PATCH 19/21] Move toward dcoll --- examples/wave.py | 6 ++--- mirgecom/boundary.py | 3 ++- mirgecom/diffusion.py | 43 +++++++++++++++++----------------- mirgecom/gas_model.py | 5 ++-- mirgecom/inviscid.py | 5 ++-- mirgecom/io.py | 8 +++---- mirgecom/logging_quantities.py | 19 +++++++++------ mirgecom/operators.py | 12 +++++----- mirgecom/simutil.py | 2 +- mirgecom/viscous.py | 11 +++++---- mirgecom/wave.py | 15 ++++++------ 11 files changed, 70 insertions(+), 59 deletions(-) diff --git a/examples/wave.py b/examples/wave.py index f575418d0..1d15d494f 100644 --- a/examples/wave.py +++ b/examples/wave.py @@ -32,7 +32,7 @@ from grudge.eager import EagerDGDiscretization from grudge.shortcuts import make_visualizer - +from grudge.op import nodal_min from mirgecom.wave import wave_operator from mirgecom.integrators import rk4_step @@ -108,8 +108,8 @@ def main(use_profiling=False, use_logmgr=False, lazy_eval: bool = False): wave_speed = 1.0 from grudge.dt_utils import characteristic_lengthscales nodal_dt = characteristic_lengthscales(actx, discr) / wave_speed - from grudge.op import nodal_min - dt = actx.to_numpy(current_cfl * nodal_min(discr, "vol", nodal_dt))[()] + dt = actx.to_numpy(current_cfl * nodal_min(discr, "vol", + nodal_dt))[()] print("%d elements" % mesh.nelements) diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py index 42748ad01..31fba3f05 100644 --- a/mirgecom/boundary.py +++ b/mirgecom/boundary.py @@ -48,6 +48,7 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from mirgecom.fluid import make_conserved from grudge.trace_pair import TracePair +import grudge.op as op from mirgecom.viscous import viscous_facial_flux_central from mirgecom.flux import num_flux_central from mirgecom.gas_model import make_fluid_state @@ -314,7 +315,7 @@ def _boundary_quantity(self, discr, btag, quantity, local=False, **kwargs): """Get a boundary quantity on local boundary, or projected to "all_faces".""" from grudge.dof_desc import as_dofdesc btag = as_dofdesc(btag) - return quantity if local else discr.project( + return quantity if local else op.project(discr, btag, btag.with_dtag("all_faces"), quantity) def _boundary_state_pair(self, discr, btag, gas_model, state_minus, **kwargs): diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py index f41743a53..6b7286e38 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -41,6 +41,7 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.dof_desc import DOFDesc, as_dofdesc, DISCR_TAG_BASE from grudge.trace_pair import TracePair, interior_trace_pairs +import grudge.op as op def grad_flux(discr, u_tpair, *, quadrature_tag=DISCR_TAG_BASE): @@ -54,12 +55,12 @@ def grad_flux(discr, u_tpair, *, quadrature_tag=DISCR_TAG_BASE): normal_quad = thaw(discr.normal(dd_quad), actx) def to_quad(a): - return discr.project(dd, dd_quad, a) + return op.project(discr, dd, dd_quad, a) def flux(u, normal): return -u * normal - return discr.project(dd_quad, dd_allfaces_quad, flux( + return op.project(discr, dd_quad, dd_allfaces_quad, flux( to_quad(u_tpair.avg), normal_quad)) @@ -75,7 +76,7 @@ def diffusion_flux( normal_quad = thaw(discr.normal(dd_quad), actx) def to_quad(a): - return discr.project(dd, dd_quad, a) + return op.project(discr, dd, dd_quad, a) def flux(kappa, grad_u, normal): return -kappa * np.dot(grad_u, normal) @@ -87,7 +88,7 @@ def flux(kappa, grad_u, normal): to_quad(kappa_tpair.ext), to_quad(grad_u_tpair.ext), normal_quad) ) - return discr.project(dd_quad, dd_allfaces_quad, flux_tpair.avg) + return op.project(discr, dd_quad, dd_allfaces_quad, flux_tpair.avg) class DiffusionBoundary(metaclass=abc.ABCMeta): @@ -141,16 +142,16 @@ def __init__(self, value): def get_grad_flux( self, discr, dd, u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 - u_int = discr.project("vol", dd, u) + u_int = op.project(discr, "vol", dd, u) u_tpair = TracePair(dd, interior=u_int, exterior=2*self.value-u_int) return grad_flux(discr, u_tpair, quadrature_tag=quadrature_tag) def get_diffusion_flux( self, discr, dd, kappa, grad_u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 - kappa_int = discr.project("vol", dd, kappa) + kappa_int = op.project(discr, "vol", dd, kappa) kappa_tpair = TracePair(dd, interior=kappa_int, exterior=kappa_int) - grad_u_int = discr.project("vol", dd, grad_u) + grad_u_int = op.project(discr, "vol", dd, grad_u) grad_u_tpair = TracePair(dd, interior=grad_u_int, exterior=grad_u_int) return diffusion_flux( discr, kappa_tpair, grad_u_tpair, quadrature_tag=quadrature_tag) @@ -194,7 +195,7 @@ def __init__(self, value): def get_grad_flux( self, discr, dd, u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 - u_int = discr.project("vol", dd, u) + u_int = op.project(discr, "vol", dd, u) u_tpair = TracePair(dd, interior=u_int, exterior=u_int) return grad_flux(discr, u_tpair, quadrature_tag=quadrature_tag) @@ -208,10 +209,10 @@ def get_diffusion_flux( # spatially-varying kappa case (the other approach would result in a # grad_u_tpair that lives in the quadrature discretization; diffusion_flux # would need to be modified to accept such values). - kappa_int_quad = discr.project("vol", dd_quad, kappa) - value_quad = discr.project(dd, dd_quad, self.value) + kappa_int_quad = op.project(discr, "vol", dd_quad, kappa) + value_quad = op.project(discr, dd, dd_quad, self.value) flux_quad = -kappa_int_quad*value_quad - return discr.project(dd_quad, dd_allfaces_quad, flux_quad) + return op.project(discr, dd_quad, dd_allfaces_quad, flux_quad) class _DiffusionStateTag: @@ -268,10 +269,10 @@ def grad_operator(discr, boundaries, u, quadrature_tag=DISCR_TAG_BASE): dd_allfaces_quad = DOFDesc("all_faces", quadrature_tag) - return discr.inverse_mass( - discr.weak_grad(-u) - - # noqa: W504 - discr.face_mass( + return op.inverse_mass(discr, + op.weak_local_grad(discr, "vol", -u) + - 1.0 # noqa: W504 + * op.face_mass(discr, dd_allfaces_quad, sum( grad_flux(discr, u_tpair, quadrature_tag=quadrature_tag) @@ -383,13 +384,13 @@ def diffusion_operator(discr, *args, return_grad_u=False, **kwargs): grad_u = grad_operator(discr, boundaries, u, quadrature_tag=quadrature_tag) - kappa_quad = discr.project("vol", dd_quad, kappa) - grad_u_quad = discr.project("vol", dd_quad, grad_u) + kappa_quad = op.project(discr, "vol", dd_quad, kappa) + grad_u_quad = op.project(discr, "vol", dd_quad, grad_u) - diff_u = discr.inverse_mass( - discr.weak_div(dd_quad, -kappa_quad*grad_u_quad) - - # noqa: W504 - discr.face_mass( + diff_u = op.inverse_mass(discr, + op.weak_local_div(discr, dd_quad, -kappa_quad*grad_u_quad) + - 1. # noqa: W504 + * op.face_mass(discr, dd_allfaces_quad, sum( diffusion_flux( diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py index e6067b53e..a91b3b2e3 100644 --- a/mirgecom/gas_model.py +++ b/mirgecom/gas_model.py @@ -60,6 +60,7 @@ GasTransportVars ) from grudge.dof_desc import DOFDesc, as_dofdesc +import grudge.op as op from grudge.trace_pair import ( interior_trace_pairs, tracepair_with_discr_tag @@ -318,10 +319,10 @@ def project_fluid_state(discr, src, tgt, state, gas_model): Thermally consistent fluid state """ - cv_sd = discr.project(src, tgt, state.cv) + cv_sd = op.project(discr, src, tgt, state.cv) temperature_seed = None if state.is_mixture: - temperature_seed = discr.project(src, tgt, state.dv.temperature) + temperature_seed = op.project(discr, src, tgt, state.dv.temperature) return make_fluid_state(cv=cv_sd, gas_model=gas_model, temperature_seed=temperature_seed) diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py index 31ccd30c6..515aad2b2 100644 --- a/mirgecom/inviscid.py +++ b/mirgecom/inviscid.py @@ -41,6 +41,7 @@ import numpy as np from arraycontext import thaw +import grudge.op as op from mirgecom.fluid import make_conserved @@ -265,14 +266,14 @@ def inviscid_flux_on_element_boundary( from grudge.dof_desc import as_dofdesc def _interior_flux(state_pair): - return discr.project( + return op.project(discr, state_pair.dd, state_pair.dd.with_dtag("all_faces"), numerical_flux_func( state_pair, gas_model, thaw(discr.normal(state_pair.dd), state_pair.int.array_context))) def _boundary_flux(dd_bdry, boundary, state_minus): - return discr.project( + return op.project(discr, dd_bdry, dd_bdry.with_dtag("all_faces"), boundary.inviscid_divergence_flux( discr, dd_bdry, gas_model, state_minus=state_minus, diff --git a/mirgecom/io.py b/mirgecom/io.py index e6c7eb817..aac74db51 100644 --- a/mirgecom/io.py +++ b/mirgecom/io.py @@ -28,7 +28,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ - +from functools import partial +import grudge.op as op from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa @@ -54,9 +55,8 @@ def make_init_message(*, dim, order, dt, t_final, def make_status_message(*, discr, t, step, dt, cfl, dependent_vars): r"""Make simulation status and health message.""" dv = dependent_vars - from functools import partial - _min = partial(discr.nodal_min, "vol") - _max = partial(discr.nodal_max, "vol") + _min = partial(op.nodal_min, discr, "vol") + _max = partial(op.nodal_max, discr, "vol") statusmsg = ( f"Status: {step=} {t=}\n" f"------- P({_min(dv.pressure):.3g}, {_max(dv.pressure):.3g})\n" diff --git a/mirgecom/logging_quantities.py b/mirgecom/logging_quantities.py index 40dd2161f..21c9b9101 100644 --- a/mirgecom/logging_quantities.py +++ b/mirgecom/logging_quantities.py @@ -50,6 +50,8 @@ from typing import Optional, Callable import numpy as np +import grudge.op as oper + def initialize_logmgr(enable_logmgr: bool, filename: str = None, mode: str = "wu", @@ -102,18 +104,21 @@ def logmgr_add_device_memory_usage(logmgr: LogManager, queue: cl.CommandQueue): def logmgr_add_many_discretization_quantities(logmgr: LogManager, discr, dim, extract_vars_for_logging, units_for_logging): """Add default discretization quantities to the logmgr.""" - for op in ["min", "max", "L2_norm"]: + for reduction_op in ["min", "max", "L2_norm"]: for quantity in ["pressure", "temperature"]: logmgr.add_quantity(DiscretizationBasedQuantity( - discr, quantity, op, extract_vars_for_logging, units_for_logging)) + discr, quantity, reduction_op, extract_vars_for_logging, + units_for_logging)) for quantity in ["mass", "energy"]: logmgr.add_quantity(DiscretizationBasedQuantity( - discr, quantity, op, extract_vars_for_logging, units_for_logging)) + discr, quantity, reduction_op, extract_vars_for_logging, + units_for_logging)) for d in range(dim): logmgr.add_quantity(DiscretizationBasedQuantity( - discr, "momentum", op, extract_vars_for_logging, units_for_logging, + discr, "momentum", reduction_op, extract_vars_for_logging, + units_for_logging, axis=d)) @@ -257,13 +262,13 @@ def __init__(self, discr: Discretization, quantity: str, op: str, from functools import partial if op == "min": - self._discr_reduction = partial(self.discr.nodal_min, "vol") + self._discr_reduction = partial(oper.nodal_min, self.discr, "vol") self.rank_aggr = min elif op == "max": - self._discr_reduction = partial(self.discr.nodal_max, "vol") + self._discr_reduction = partial(oper.nodal_max, self.discr, "vol") self.rank_aggr = max elif op == "L2_norm": - self._discr_reduction = partial(self.discr.norm, p=2) + self._discr_reduction = partial(oper.norm, self.discr, p=2) self.rank_aggr = max else: raise ValueError(f"unknown operation {op}") diff --git a/mirgecom/operators.py b/mirgecom/operators.py index 240ee542e..b584bb5d8 100644 --- a/mirgecom/operators.py +++ b/mirgecom/operators.py @@ -56,10 +56,10 @@ def grad_operator(discr, dd_vol, dd_faces, u, flux): meshmode.dof_array.DOFArray or numpy.ndarray the dg gradient operator applied to *u* """ - return -discr.inverse_mass( + # pylint: disable=invalid-unary-operand-type + return - op.inverse_mass(discr, op.weak_local_grad(discr, dd_vol, u) - - op.face_mass(discr, dd_faces, flux) - ) + - op.face_mass(discr, dd_faces, flux)) def div_operator(discr, dd_vol, dd_faces, v, flux): @@ -87,7 +87,7 @@ def div_operator(discr, dd_vol, dd_faces, v, flux): meshmode.dof_array.DOFArray or numpy.ndarray the dg divergence operator applied to vector-valued function(s) *v*. """ - return -discr.inverse_mass( + # pylint: disable=invalid-unary-operand-type + return - op.inverse_mass(discr, op.weak_local_div(discr, dd_vol, v) - - op.face_mass(discr, dd_faces, flux) - ) + - op.face_mass(discr, dd_faces, flux)) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index 1cdd2d4cc..b9d7e498a 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -323,7 +323,7 @@ def componentwise_norms(discr, fields, order=np.inf): return map_array_container( partial(componentwise_norms, discr, order=order), fields) if len(fields) > 0: - return discr.norm(fields, order) + return op.norm(discr, fields, order) else: # FIXME: This work-around for #575 can go away after #569 return 0 diff --git a/mirgecom/viscous.py b/mirgecom/viscous.py index 8f5982ad7..1bf5d73d6 100644 --- a/mirgecom/viscous.py +++ b/mirgecom/viscous.py @@ -47,6 +47,8 @@ from meshmode.dof_array import DOFArray from arraycontext import thaw +import grudge.op as op + from mirgecom.fluid import ( velocity_gradient, species_mass_fraction_gradient, @@ -384,7 +386,6 @@ def viscous_flux_on_element_boundary( Time """ from grudge.dof_desc import as_dofdesc - from grudge.op import project dd_base = as_dofdesc("vol") @@ -392,7 +393,7 @@ def viscous_flux_on_element_boundary( # viscous fluxes across interior faces (including partition and periodic bnd) def _fvisc_divergence_flux_interior(state_pair, grad_cv_pair, grad_t_pair): - return discr.project( + return op.project(discr, state_pair.dd, state_pair.dd.with_dtag("all_faces"), numerical_flux_func( discr=discr, gas_model=gas_model, state_pair=state_pair, @@ -402,13 +403,13 @@ def _fvisc_divergence_flux_interior(state_pair, grad_cv_pair, grad_t_pair): def _fvisc_divergence_flux_boundary(dd_btag, boundary, state_minus): # Make sure we fields on the quadrature grid # restricted to the tag *btag* - return project( + return op.project( discr, dd_btag, dd_btag.with_dtag("all_faces"), boundary.viscous_divergence_flux( discr=discr, btag=dd_btag, gas_model=gas_model, state_minus=state_minus, - grad_cv_minus=project(discr, dd_base, dd_btag, grad_cv), - grad_t_minus=project(discr, dd_base, dd_btag, grad_t), + grad_cv_minus=op.project(discr, dd_base, dd_btag, grad_cv), + grad_t_minus=op.project(discr, dd_base, dd_btag, grad_t), time=time, numerical_flux_func=numerical_flux_func)) # }}} viscous flux helpers diff --git a/mirgecom/wave.py b/mirgecom/wave.py index d62b98f3f..6c89dc01c 100644 --- a/mirgecom/wave.py +++ b/mirgecom/wave.py @@ -34,6 +34,7 @@ from meshmode.dof_array import thaw from grudge.trace_pair import TracePair from grudge.eager import interior_trace_pair, cross_rank_trace_pairs +import grudge.op as op def _flux(discr, c, w_tpair): @@ -55,7 +56,7 @@ def _flux(discr, c, w_tpair): 0.5*normal*np.dot(normal, v.ext-v.int), ) - return discr.project(w_tpair.dd, "all_faces", c*flux_weak) + return op.project(discr, w_tpair.dd, "all_faces", c*flux_weak) class _WaveTag: @@ -82,19 +83,19 @@ def wave_operator(discr, c, w): u = w[0] v = w[1:] - dir_u = discr.project("vol", BTAG_ALL, u) - dir_v = discr.project("vol", BTAG_ALL, v) + dir_u = op.project(discr, "vol", BTAG_ALL, u) + dir_v = op.project(discr, "vol", BTAG_ALL, v) dir_bval = flat_obj_array(dir_u, dir_v) dir_bc = flat_obj_array(-dir_u, dir_v) return ( - discr.inverse_mass( + op.inverse_mass(discr, flat_obj_array( - -c*discr.weak_div(v), - -c*discr.weak_grad(u) + -c*op.weak_local_div(discr, "vol", v), + -c*op.weak_local_grad(discr, "vol", u) ) + # noqa: W504 - discr.face_mass( + op.face_mass(discr, _flux(discr, c=c, w_tpair=interior_trace_pair(discr, w)) + _flux(discr, c=c, w_tpair=TracePair(BTAG_ALL, interior=dir_bval, exterior=dir_bc)) From 5f6e1e37e63ffc92b98385689ca2f3ed82a4c835 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Wed, 15 Jun 2022 08:26:54 -0500 Subject: [PATCH 20/21] remove superfluous 1.0 --- mirgecom/diffusion.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py index 6b7286e38..5baa921d0 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -271,8 +271,8 @@ def grad_operator(discr, boundaries, u, quadrature_tag=DISCR_TAG_BASE): return op.inverse_mass(discr, op.weak_local_grad(discr, "vol", -u) - - 1.0 # noqa: W504 - * op.face_mass(discr, + - # noqa: W504 + op.face_mass(discr, dd_allfaces_quad, sum( grad_flux(discr, u_tpair, quadrature_tag=quadrature_tag) From c15683607ce4106f3bc438a1af3e6411df305f53 Mon Sep 17 00:00:00 2001 From: Michael Campbell Date: Wed, 15 Jun 2022 08:28:09 -0500 Subject: [PATCH 21/21] Undo production 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 6ddf8f69e..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="move-toward-dcoll-production" # 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