diff --git a/doc/conf.py b/doc/conf.py index b26fbdd74..196570247 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -91,3 +91,8 @@ rst_prolog = """ .. |mirgecom| replace:: *MIRGE-Com* """ + +# FIXME: Remove when grudge#280 gets merged +nitpick_ignore_regex = [ + ("py:class", r".*BoundaryDomainTag.*"), +] diff --git a/examples/autoignition-mpi.py b/examples/autoignition-mpi.py index bcddb2b99..4eb18b4dd 100644 --- a/examples/autoignition-mpi.py +++ b/examples/autoignition-mpi.py @@ -170,8 +170,7 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, generate_mesh) local_nelements = local_mesh.nelements - dcoll = create_discretization_collection(actx, local_mesh, order=order, - mpi_communicator=comm) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) ones = dcoll.zeros(actx) + 1.0 @@ -320,8 +319,7 @@ def get_fluid_state(cv, tseed): else: rst_cv = restart_data["cv"] old_dcoll = \ - create_discretization_collection(actx, local_mesh, order=rst_order, - mpi_communicator=comm) + create_discretization_collection(actx, local_mesh, order=rst_order) from meshmode.discretization.connection import make_same_mesh_connection connection = make_same_mesh_connection(actx, dcoll.discr_from_dd("vol"), old_dcoll.discr_from_dd("vol")) diff --git a/examples/doublemach-mpi.py b/examples/doublemach-mpi.py index 689852997..39c4daeb1 100644 --- a/examples/doublemach-mpi.py +++ b/examples/doublemach-mpi.py @@ -30,7 +30,7 @@ from functools import partial from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.dof_desc import DTAG_BOUNDARY +from grudge.dof_desc import BoundaryDomainTag from grudge.shortcuts import make_visualizer @@ -192,8 +192,7 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, from mirgecom.discretization import create_discretization_collection order = 3 - dcoll = create_discretization_collection(actx, local_mesh, order=order, - mpi_communicator=comm) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) from grudge.dof_desc import DISCR_TAG_QUAD @@ -237,22 +236,22 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, eos = IdealSingleGas() gas_model = GasModel(eos=eos, transport=transport_model) - def _boundary_state(dcoll, btag, gas_model, state_minus, **kwargs): + def _boundary_state(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) 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(), + BoundaryDomainTag("ic1"): + PrescribedFluidBoundary(boundary_state_func=_boundary_state), + BoundaryDomainTag("ic2"): + PrescribedFluidBoundary(boundary_state_func=_boundary_state), + BoundaryDomainTag("ic3"): + PrescribedFluidBoundary(boundary_state_func=_boundary_state), + BoundaryDomainTag("wall"): AdiabaticNoslipMovingBoundary(), + BoundaryDomainTag("out"): AdiabaticNoslipMovingBoundary(), } if rst_filename: diff --git a/examples/heat-source-mpi.py b/examples/heat-source-mpi.py index 893d7e759..7218b2a29 100644 --- a/examples/heat-source-mpi.py +++ b/examples/heat-source-mpi.py @@ -30,7 +30,7 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa import grudge.op as op from grudge.shortcuts import make_visualizer -from grudge.dof_desc import DTAG_BOUNDARY +from grudge.dof_desc import BoundaryDomainTag from mirgecom.discretization import create_discretization_collection from mirgecom.integrators import rk4_step from mirgecom.diffusion import ( @@ -111,8 +111,7 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, order = 3 - dcoll = create_discretization_collection(actx, local_mesh, order=order, - mpi_communicator=comm) + dcoll = create_discretization_collection(actx, local_mesh, order=order) if dim == 2: # no deep meaning here, just a fudge factor @@ -125,8 +124,8 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, nodes = actx.thaw(dcoll.nodes()) boundaries = { - DTAG_BOUNDARY("dirichlet"): DirichletDiffusionBoundary(0.), - DTAG_BOUNDARY("neumann"): NeumannDiffusionBoundary(0.) + BoundaryDomainTag("dirichlet"): DirichletDiffusionBoundary(0.), + BoundaryDomainTag("neumann"): NeumannDiffusionBoundary(0.) } u = dcoll.zeros(actx) diff --git a/examples/hotplate-mpi.py b/examples/hotplate-mpi.py index a8082c6b3..001d3b968 100644 --- a/examples/hotplate-mpi.py +++ b/examples/hotplate-mpi.py @@ -31,7 +31,7 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.shortcuts import make_visualizer -from grudge.dof_desc import DTAG_BOUNDARY +from grudge.dof_desc import BoundaryDomainTag from mirgecom.discretization import create_discretization_collection from mirgecom.fluid import make_conserved @@ -168,9 +168,7 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, local_nelements = local_mesh.nelements order = 1 - dcoll = create_discretization_collection( - actx, local_mesh, order=order, mpi_communicator=comm - ) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) if logmgr: @@ -221,21 +219,22 @@ def tramp_2d(x_vec, eos, cv=None, **kwargs): exact = initializer(x_vec=nodes, eos=gas_model.eos) - def _boundary_state(dcoll, btag, gas_model, state_minus, **kwargs): + def _boundary_state(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, **kwargs), gas_model) - boundaries = {DTAG_BOUNDARY("-1"): - PrescribedFluidBoundary(boundary_state_func=_boundary_state), - DTAG_BOUNDARY("+1"): - PrescribedFluidBoundary(boundary_state_func=_boundary_state), - DTAG_BOUNDARY("-2"): IsothermalNoSlipBoundary( - wall_temperature=bottom_boundary_temperature), - DTAG_BOUNDARY("+2"): IsothermalNoSlipBoundary( - wall_temperature=top_boundary_temperature)} + boundaries = { + BoundaryDomainTag("-1"): PrescribedFluidBoundary( + boundary_state_func=_boundary_state), + BoundaryDomainTag("+1"): PrescribedFluidBoundary( + boundary_state_func=_boundary_state), + BoundaryDomainTag("-2"): IsothermalNoSlipBoundary( + wall_temperature=bottom_boundary_temperature), + BoundaryDomainTag("+2"): IsothermalNoSlipBoundary( + wall_temperature=top_boundary_temperature)} if rst_filename: current_t = restart_data["t"] diff --git a/examples/lump-mpi.py b/examples/lump-mpi.py index 30cc53293..b26d9a58f 100644 --- a/examples/lump-mpi.py +++ b/examples/lump-mpi.py @@ -146,9 +146,7 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, local_nelements = local_mesh.nelements order = 3 - dcoll = create_discretization_collection( - actx, local_mesh, order=order, mpi_communicator=comm - ) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) vis_timer = None @@ -180,9 +178,9 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, from mirgecom.gas_model import GasModel, make_fluid_state gas_model = GasModel(eos=eos) - def boundary_solution(dcoll, btag, gas_model, state_minus, **kwargs): + def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, **kwargs), gas_model) diff --git a/examples/mixture-mpi.py b/examples/mixture-mpi.py index 8dd9e98ed..08fc74979 100644 --- a/examples/mixture-mpi.py +++ b/examples/mixture-mpi.py @@ -146,9 +146,7 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, local_nelements = local_mesh.nelements order = 3 - dcoll = create_discretization_collection( - actx, local_mesh, order=order, mpi_communicator=comm - ) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) vis_timer = None @@ -202,9 +200,9 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, initializer = MixtureInitializer(dim=dim, nspecies=nspecies, massfractions=y0s, velocity=velocity) - def boundary_solution(dcoll, btag, gas_model, state_minus, **kwargs): + def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, **kwargs), gas_model, diff --git a/examples/multiple-volumes-mpi-lazy.py b/examples/multiple-volumes-mpi-lazy.py new file mode 120000 index 000000000..ab0b3b489 --- /dev/null +++ b/examples/multiple-volumes-mpi-lazy.py @@ -0,0 +1 @@ +multiple-volumes-mpi.py \ No newline at end of file diff --git a/examples/multiple-volumes-mpi.py b/examples/multiple-volumes-mpi.py new file mode 100644 index 000000000..67d466d80 --- /dev/null +++ b/examples/multiple-volumes-mpi.py @@ -0,0 +1,425 @@ +""" +Demonstrate multiple non-interacting volumes. + +Runs several acoustic pulse simulations with different pulse amplitudes +simultaneously. +""" + +__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 +from mirgecom.mpi import mpi_entry_point +import numpy as np +from functools import partial +import pyopencl as cl +from pytools.obj_array import make_obj_array + +from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from grudge.shortcuts import make_visualizer +from grudge.dof_desc import VolumeDomainTag, DISCR_TAG_BASE, DISCR_TAG_QUAD, DOFDesc + +from mirgecom.discretization import create_discretization_collection +from mirgecom.euler import ( + euler_operator, + extract_vars_for_logging +) +from mirgecom.simutil import ( + get_sim_timestep, + generate_and_distribute_mesh +) +from mirgecom.io import make_init_message + +from mirgecom.integrators import rk4_step +from mirgecom.steppers import advance_state +from mirgecom.boundary import AdiabaticSlipBoundary +from mirgecom.initializers import ( + Lump, + AcousticPulse +) +from mirgecom.eos import IdealSingleGas +from mirgecom.gas_model import ( + GasModel, + make_fluid_state +) +from logpyle import IntervalTimer, set_dt +from mirgecom.logging_quantities import ( + initialize_logmgr, + logmgr_add_many_discretization_quantities, + logmgr_add_cl_device_info, + logmgr_add_device_memory_usage, + set_sim_state +) + +logger = logging.getLogger(__name__) + + +class MyRuntimeError(RuntimeError): + """Simple exception to kill the simulation.""" + + pass + + +@mpi_entry_point +def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, + use_overintegration=False, lazy=False, use_leap=False, use_profiling=False, + casename=None, rst_filename=None): + """Drive the example.""" + cl_ctx = ctx_factory() + + if casename is None: + casename = "mirgecom" + + from mpi4py import MPI + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + num_parts = comm.Get_size() + + from mirgecom.simutil import global_reduce as _global_reduce + global_reduce = partial(_global_reduce, comm=comm) + + 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) + + from mirgecom.simutil import get_reasonable_memory_pool + alloc = get_reasonable_memory_pool(cl_ctx, queue) + + if lazy: + actx = actx_class(comm, queue, mpi_base_tag=12000, allocator=alloc) + else: + actx = actx_class(comm, queue, allocator=alloc, force_device_scalars=True) + + # timestepping control + current_step = 0 + if use_leap: + from leap.rk import RK4MethodBuilder + timestepper = RK4MethodBuilder("state") + else: + timestepper = rk4_step + t_final = 0.1 + current_cfl = 1.0 + current_dt = .005 + current_t = 0 + constant_cfl = False + + # some i/o frequencies + nstatus = 1 + nrestart = 5 + nviz = 10 + nhealth = 1 + + dim = 2 + + # Run simulations with several different pulse amplitudes simultaneously + pulse_amplitudes = [0.01, 0.1, 1.0] + nvolumes = len(pulse_amplitudes) + + 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_prototype_mesh = restart_data["local_prototype_mesh"] + global_prototype_nelements = restart_data["global_prototype_nelements"] + assert restart_data["num_parts"] == num_parts + else: # generate the grids from scratch + from meshmode.mesh.generation import generate_regular_rect_mesh + generate_mesh = partial(generate_regular_rect_mesh, + a=(-1,)*dim, b=(1,)*dim, nelements_per_axis=(16,)*dim) + local_prototype_mesh, global_prototype_nelements = \ + generate_and_distribute_mesh(comm, generate_mesh) + + volume_to_local_mesh = {i: local_prototype_mesh for i in range(nvolumes)} + + local_nelements = local_prototype_mesh.nelements * nvolumes + global_nelements = global_prototype_nelements * nvolumes + + order = 3 + dcoll = create_discretization_collection(actx, volume_to_local_mesh, order=order) + + volume_dds = [ + DOFDesc(VolumeDomainTag(i), DISCR_TAG_BASE) + for i in range(nvolumes)] + + if use_overintegration: + quadrature_tag = DISCR_TAG_QUAD + else: + quadrature_tag = None + + vis_timer = None + + if logmgr: + logmgr_add_cl_device_info(logmgr, queue) + logmgr_add_device_memory_usage(logmgr, queue) + + def extract_vars(i, dim, cvs, eos): + name_to_field = extract_vars_for_logging(dim, cvs[i], eos) + return { + name + f"_{i}": field + for name, field in name_to_field.items()} + + def units(quantity): + return "" + + for i in range(nvolumes): + logmgr_add_many_discretization_quantities( + logmgr, dcoll, dim, partial(extract_vars, i), units, + dd=volume_dds[i]) + + vis_timer = IntervalTimer("t_vis", "Time spent visualizing") + logmgr.add_quantity(vis_timer) + + logmgr.add_watches([ + ("step.max", "step = {value}, "), + ("t_sim.max", "sim time: {value:1.6e} s\n"), + ("t_step.max", "------- step walltime: {value:6g} s, "), + ("t_log.max", "log walltime: {value:6g} s\n") + ]) + + for i in range(nvolumes): + logmgr.add_watches([ + (f"min_pressure_{i}", "------- P (vol. " + str(i) + + ") (min, max) (Pa) = ({value:1.9e}, "), + (f"max_pressure_{i}", "{value:1.9e})\n"), + ]) + + eos = IdealSingleGas() + gas_model = GasModel(eos=eos) + wall = AdiabaticSlipBoundary() + if rst_filename: + current_t = restart_data["t"] + current_step = restart_data["step"] + current_cvs = restart_data["cvs"] + 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 + def init(nodes, pulse_amplitude): + vel = np.zeros(shape=(dim,)) + orig = np.zeros(shape=(dim,)) + background = Lump( + dim=dim, center=orig, velocity=vel, rhoamp=0.0)(nodes) + return AcousticPulse( + dim=dim, + amplitude=pulse_amplitude, + width=0.1, + center=orig)(x_vec=nodes, cv=background, eos=eos) + current_cvs = make_obj_array([ + init(actx.thaw(dcoll.nodes(dd)), pulse_amplitude) + for dd, pulse_amplitude in zip(volume_dds, pulse_amplitudes)]) + + current_fluid_states = [make_fluid_state(cv, gas_model) for cv in current_cvs] + + visualizers = [make_visualizer(dcoll, volume_dd=dd) for dd in volume_dds] + + initname = "multiple-volumes" + 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_get_timestep(step, t, dt, fluid_states): + return min([ + get_sim_timestep( + dcoll, fluid_state, t, dt, current_cfl, t_final, constant_cfl) + for fluid_state in fluid_states]) + + def my_write_viz(step, t, cvs, dvs=None): + if dvs is None: + dvs = [eos.dependent_vars(cv) for cv in cvs] + for i in range(nvolumes): + viz_fields = [ + ("cv", cvs[i]), + ("dv", dvs[i])] + from mirgecom.simutil import write_visfile + write_visfile( + dcoll, viz_fields, visualizers[i], vizname=casename + f"-{i}", + step=step, t=t, overwrite=True, vis_timer=vis_timer, comm=comm) + + def my_write_restart(step, t, cvs): + rst_fname = rst_pattern.format(cname=casename, step=step, rank=rank) + if rst_fname != rst_filename: + rst_data = { + "local_prototype_mesh": local_prototype_mesh, + "cvs": cvs, + "t": t, + "step": step, + "order": order, + "global_nelements": global_nelements, + "num_parts": num_parts + } + from mirgecom.restart import write_restart_file + write_restart_file(actx, rst_data, rst_fname, comm) + + def my_health_check(pressures): + health_error = False + for dd, pressure in zip(volume_dds, pressures): + from mirgecom.simutil import check_naninf_local, check_range_local + if check_naninf_local(dcoll, dd, pressure) \ + or check_range_local(dcoll, dd, pressure, 1e-2, 10): + health_error = True + logger.info(f"{rank=}: Invalid pressure data found.") + break + return health_error + + def my_pre_step(step, t, dt, state): + cvs = state + fluid_states = [make_fluid_state(cv, gas_model) for cv in cvs] + dvs = [fluid_state.dv for fluid_state in fluid_states] + + 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: + pressures = [dv.pressure for dv in dvs] + health_errors = global_reduce(my_health_check(pressures), op="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, cvs=cvs) + + if do_viz: + my_write_viz(step=step, t=t, cvs=cvs, dvs=dvs) + + except MyRuntimeError: + if rank == 0: + logger.info("Errors detected; attempting graceful exit.") + my_write_viz(step=step, t=t, cvs=cvs) + my_write_restart(step=step, t=t, cvs=cvs) + raise + + dt = my_get_timestep(step=step, t=t, dt=dt, fluid_states=fluid_states) + + return cvs, 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): + cvs = state + fluid_states = [make_fluid_state(cv, gas_model) for cv in cvs] + return make_obj_array([ + euler_operator( + dcoll, state=fluid_state, time=t, + boundaries={dd.trace(BTAG_ALL).domain_tag: wall}, + gas_model=gas_model, quadrature_tag=quadrature_tag, + dd=dd, comm_tag=dd) + for dd, fluid_state in zip(volume_dds, fluid_states)]) + + current_dt = my_get_timestep( + current_step, current_t, current_dt, current_fluid_states) + + current_step, current_t, current_cvs = \ + advance_state(rhs=my_rhs, timestepper=timestepper, + pre_step_callback=my_pre_step, + post_step_callback=my_post_step, dt=current_dt, + state=current_cvs, t=current_t, t_final=t_final) + + # Dump the final data + if rank == 0: + logger.info("Checkpointing final state ...") + final_fluid_states = [make_fluid_state(cv, gas_model) for cv in current_cvs] + final_dvs = [fluid_state.dv for fluid_state in final_fluid_states] + + my_write_viz(step=current_step, t=current_t, cvs=current_cvs, dvs=final_dvs) + my_write_restart(step=current_step, t=current_t, cvs=current_cvs) + + 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 = "multiple-volumes" + 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(actx_class, use_logmgr=args.log, use_overintegration=args.overintegration, + use_leap=args.leap, use_profiling=args.profiling, lazy=lazy, + casename=casename, rst_filename=rst_filename) + +# vim: foldmethod=marker diff --git a/examples/nsmix-mpi.py b/examples/nsmix-mpi.py index 9c160ac83..cdf81a07a 100644 --- a/examples/nsmix-mpi.py +++ b/examples/nsmix-mpi.py @@ -153,9 +153,7 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, local_nelements = local_mesh.nelements order = 1 - dcoll = create_discretization_collection( - actx, local_mesh, order=order, mpi_communicator=comm - ) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) ones = dcoll.zeros(actx) + 1.0 diff --git a/examples/poiseuille-mpi.py b/examples/poiseuille-mpi.py index 2d0afbd5e..dc3492827 100644 --- a/examples/poiseuille-mpi.py +++ b/examples/poiseuille-mpi.py @@ -32,8 +32,7 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.shortcuts import make_visualizer -from grudge.dof_desc import DTAG_BOUNDARY -from grudge.dof_desc import DISCR_TAG_QUAD +from grudge.dof_desc import BoundaryDomainTag, DISCR_TAG_QUAD from mirgecom.discretization import create_discretization_collection from mirgecom.fluid import make_conserved @@ -168,9 +167,7 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, local_nelements = local_mesh.nelements order = 2 - dcoll = create_discretization_collection( - actx, local_mesh, order=order, mpi_communicator=comm - ) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) if use_overintegration: @@ -235,19 +232,22 @@ def poiseuille_2d(x_vec, eos, cv=None, **kwargs): transport=SimpleTransport(viscosity=mu)) exact = initializer(x_vec=nodes, eos=gas_model.eos) - def _boundary_solution(dcoll, btag, gas_model, state_minus, **kwargs): + def _boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, cv=state_minus.cv, **kwargs), gas_model) - boundaries = {DTAG_BOUNDARY("-1"): - PrescribedFluidBoundary(boundary_state_func=_boundary_solution), - DTAG_BOUNDARY("+1"): - PrescribedFluidBoundary(boundary_state_func=_boundary_solution), - DTAG_BOUNDARY("-2"): AdiabaticNoslipMovingBoundary(), - DTAG_BOUNDARY("+2"): AdiabaticNoslipMovingBoundary()} + boundaries = { + BoundaryDomainTag("-1"): + PrescribedFluidBoundary( + boundary_state_func=_boundary_solution), + BoundaryDomainTag("+1"): + PrescribedFluidBoundary( + boundary_state_func=_boundary_solution), + BoundaryDomainTag("-2"): AdiabaticNoslipMovingBoundary(), + BoundaryDomainTag("+2"): AdiabaticNoslipMovingBoundary()} if rst_filename: current_t = restart_data["t"] diff --git a/examples/pulse-mpi.py b/examples/pulse-mpi.py index 3b0ca775a..13ef986c8 100644 --- a/examples/pulse-mpi.py +++ b/examples/pulse-mpi.py @@ -152,9 +152,7 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, local_nelements = local_mesh.nelements order = 1 - dcoll = create_discretization_collection( - actx, local_mesh, order=order, mpi_communicator=comm - ) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) if use_overintegration: diff --git a/examples/scalar-lump-mpi.py b/examples/scalar-lump-mpi.py index 19b3f3851..6b0ee947c 100644 --- a/examples/scalar-lump-mpi.py +++ b/examples/scalar-lump-mpi.py @@ -145,9 +145,7 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, local_nelements = local_mesh.nelements order = 3 - dcoll = create_discretization_collection( - actx, local_mesh, order=order, mpi_communicator=comm - ) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) vis_timer = None @@ -185,9 +183,9 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, spec_y0s=spec_y0s, spec_amplitudes=spec_amplitudes) - def boundary_solution(dcoll, btag, gas_model, state_minus, **kwargs): + def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, **kwargs), gas_model) diff --git a/examples/sod-mpi.py b/examples/sod-mpi.py index 5e8ecc226..a8de0e204 100644 --- a/examples/sod-mpi.py +++ b/examples/sod-mpi.py @@ -144,9 +144,7 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, local_nelements = local_mesh.nelements order = 1 - dcoll = create_discretization_collection( - actx, local_mesh, order=order, mpi_communicator=comm - ) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) vis_timer = None @@ -173,9 +171,9 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, eos = IdealSingleGas() gas_model = GasModel(eos=eos) - def boundary_solution(dcoll, btag, gas_model, state_minus, **kwargs): + def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, **kwargs), gas_model) diff --git a/examples/vortex-mpi.py b/examples/vortex-mpi.py index 87b9585e7..47b69d933 100644 --- a/examples/vortex-mpi.py +++ b/examples/vortex-mpi.py @@ -149,9 +149,7 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, local_nelements = local_mesh.nelements order = 3 - dcoll = create_discretization_collection( - actx, local_mesh, order=order, mpi_communicator=comm - ) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) vis_timer = None @@ -190,9 +188,9 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, initializer = Vortex2D(center=orig, velocity=vel) gas_model = GasModel(eos=eos) - def boundary_solution(dcoll, btag, gas_model, state_minus, **kwargs): + def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, **kwargs), gas_model) diff --git a/examples/wave-mpi.py b/examples/wave-mpi.py index 093235ec5..515071fc8 100644 --- a/examples/wave-mpi.py +++ b/examples/wave-mpi.py @@ -130,9 +130,7 @@ def main(actx_class, snapshot_pattern="wave-mpi-{step:04d}-{rank:04d}.pkl", order = 3 - dcoll = create_discretization_collection( - actx, local_mesh, order=order, mpi_communicator=comm - ) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) current_cfl = 0.485 @@ -161,8 +159,7 @@ def main(actx_class, snapshot_pattern="wave-mpi-{step:04d}-{rank:04d}.pkl", old_order = restart_data["order"] if old_order != order: old_dcoll = create_discretization_collection( - actx, local_mesh, order=old_order, mpi_communicator=comm - ) + actx, local_mesh, order=old_order) from meshmode.discretization.connection import make_same_mesh_connection connection = make_same_mesh_connection(actx, dcoll.discr_from_dd("vol"), old_dcoll.discr_from_dd("vol")) diff --git a/mirgecom/artificial_viscosity.py b/mirgecom/artificial_viscosity.py index 06e49af06..346bc6c6e 100644 --- a/mirgecom/artificial_viscosity.py +++ b/mirgecom/artificial_viscosity.py @@ -130,6 +130,7 @@ from pytools import memoize_in, keyed_memoize_in from functools import partial from meshmode.dof_array import DOFArray +from meshmode.discretization.connection import FACE_RESTR_ALL from mirgecom.flux import num_flux_central from mirgecom.operators import div_operator @@ -140,12 +141,14 @@ ) from grudge.dof_desc import ( - DOFDesc, - as_dofdesc, - DD_VOLUME_MODAL, - DD_VOLUME + DD_VOLUME_ALL, + VolumeDomainTag, + DISCR_TAG_BASE, + DISCR_TAG_MODAL, ) +from mirgecom.utils import normalize_boundaries + import grudge.op as op @@ -154,10 +157,12 @@ class _AVRTag: def av_laplacian_operator(dcoll, 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, - comm_tag=None, **kwargs): + kappa=1., s0=-6., time=0, quadrature_tag=DISCR_TAG_BASE, + dd=DD_VOLUME_ALL, boundary_kwargs=None, indicator=None, + divergence_numerical_flux=num_flux_central, comm_tag=None, + operator_states_quad=None, + grad_cv=None, + **kwargs): r"""Compute the artificial viscosity right-hand-side. Computes the the right-hand-side term for artificial viscosity. @@ -196,10 +201,10 @@ def av_laplacian_operator(dcoll, boundaries, fluid_state, alpha, gas_model=None, 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 + dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the discretization on which *fluid_state* lives. + Must be a volume on the base discretization. comm_tag: Hashable Tag for distributed communication @@ -209,10 +214,19 @@ def av_laplacian_operator(dcoll, boundaries, fluid_state, alpha, gas_model=None, :class:`mirgecom.fluid.ConservedVars` The artificial viscosity operator applied to *q*. """ + boundaries = normalize_boundaries(boundaries) + cv = fluid_state.cv actx = cv.array_context - dd_vol = DOFDesc("vol", quadrature_tag) - dd_faces = DOFDesc("all_faces", quadrature_tag) + + if not isinstance(dd.domain_tag, VolumeDomainTag): + raise TypeError("dd must represent a volume") + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + + dd_vol = dd + dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) from warnings import warn @@ -227,12 +241,13 @@ def av_laplacian_operator(dcoll, boundaries, fluid_state, alpha, gas_model=None, interp_to_surf_quad = partial(tracepair_with_discr_tag, dcoll, quadrature_tag) def interp_to_vol_quad(u): - return op.project(dcoll, "vol", dd_vol, u) + return op.project(dcoll, dd_vol, dd_vol_quad, u) if operator_states_quad is None: from mirgecom.gas_model import make_operator_fluid_states operator_states_quad = make_operator_fluid_states( - dcoll, fluid_state, gas_model, boundaries, quadrature_tag, comm_tag) + dcoll, fluid_state, gas_model, boundaries, quadrature_tag, + dd=dd_vol, comm_tag=comm_tag) vol_state_quad, inter_elem_bnd_states_quad, domain_bnd_states_quad = \ operator_states_quad @@ -240,45 +255,52 @@ def interp_to_vol_quad(u): # Get smoothness indicator based on mass component if indicator is None: indicator = smoothness_indicator(dcoll, fluid_state.mass_density, - kappa=kappa, s0=s0) + kappa=kappa, s0=s0, dd=dd_vol) if grad_cv is None: from mirgecom.navierstokes import grad_cv_operator grad_cv = grad_cv_operator(dcoll, gas_model, boundaries, fluid_state, time=time, quadrature_tag=quadrature_tag, - operator_states_quad=operator_states_quad, - comm_tag=comm_tag) + dd=dd_vol, + comm_tag=comm_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 = actx.thaw(dcoll.normal(dd)) - return op.project(dcoll, dd, dd.with_dtag("all_faces"), + def central_flux_div(utpair_quad): + dd_trace_quad = utpair_quad.dd + dd_allfaces_quad = dd_trace_quad.with_boundary_tag(FACE_RESTR_ALL) + normal_quad = actx.thaw(dcoll.normal(dd_trace_quad)) + return op.project(dcoll, dd_trace_quad, dd_allfaces_quad, # 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) + divergence_numerical_flux( + utpair_quad.int, utpair_quad.ext)@normal_quad) # 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(dcoll, r, - comm_tag=(_AVRTag, comm_tag))) + sum( + central_flux_div(interp_to_surf_quad(tpair=tpair)) + for tpair in interior_trace_pairs( + dcoll, r, volume_dd=dd_vol, comm_tag=(_AVRTag, comm_tag))) # Contributions from boundary fluxes - + sum(boundaries[btag].av_flux( - dcoll, - btag=as_dofdesc(btag).with_discr_tag(quadrature_tag), - diffusion=r) for btag in boundaries) + + sum( + bdry.av_flux( + dcoll, + dd_bdry=dd_vol.with_domain_tag(bdtag), + diffusion=r) + for bdtag, bdry in boundaries.items()) ) # Return the AV RHS term - return -div_operator(dcoll, dd_vol, dd_faces, interp_to_vol_quad(r), r_bnd) + return -div_operator( + dcoll, dd_vol_quad, dd_allfaces_quad, interp_to_vol_quad(r), r_bnd) -def smoothness_indicator(dcoll, u, kappa=1.0, s0=-6.0): +def smoothness_indicator(dcoll, u, kappa=1.0, s0=-6.0, dd=DD_VOLUME_ALL): r"""Calculate the smoothness indicator. Parameters @@ -305,7 +327,7 @@ def smoothness_indicator(dcoll, u, kappa=1.0, s0=-6.0): actx = u.array_context - @memoize_in(actx, (smoothness_indicator, "smooth_comp_knl")) + @memoize_in(actx, (smoothness_indicator, "smooth_comp_knl", dd)) def indicator_prg(): """Compute the smoothness indicator for all elements.""" from arraycontext import make_loopy_program @@ -331,8 +353,7 @@ def indicator_prg(): return lp.tag_inames(t_unit, {"iel": ConcurrentElementInameTag(), "idof": ConcurrentDOFInameTag()}) - @keyed_memoize_in(actx, (smoothness_indicator, - "highest_mode"), + @keyed_memoize_in(actx, (smoothness_indicator, "highest_mode", dd), lambda grp: grp.discretization_key()) def highest_mode(grp): return actx.from_numpy( @@ -342,7 +363,9 @@ def highest_mode(grp): ) # Convert to modal solution representation - modal_map = dcoll.connection_from_dds(DD_VOLUME, DD_VOLUME_MODAL) + dd_vol = dd + dd_modal = dd_vol.with_discr_tag(DISCR_TAG_MODAL) + modal_map = dcoll.connection_from_dds(dd_vol, dd_modal) uhat = modal_map(u) # Compute smoothness indicator value @@ -363,7 +386,7 @@ def highest_mode(grp): ))) .reshape(-1, 1)), uhat[grp.index].shape)) - for grp in dcoll.discr_from_dd("vol").groups + for grp in dcoll.discr_from_dd(dd_vol).groups ) ) else: @@ -375,7 +398,7 @@ def highest_mode(grp): vec=uhat[grp.index], modes_active_flag=highest_mode(grp) )["result"] - for grp in dcoll.discr_from_dd("vol").groups + for grp in dcoll.discr_from_dd(dd_vol).groups ) ) diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py index 34c789824..94c4e934b 100644 --- a/mirgecom/boundary.py +++ b/mirgecom/boundary.py @@ -45,6 +45,8 @@ import numpy as np from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from meshmode.discretization.connection import FACE_RESTR_ALL +from grudge.dof_desc import as_dofdesc from mirgecom.fluid import make_conserved from grudge.trace_pair import TracePair import grudge.op as op @@ -67,7 +69,7 @@ class FluidBoundary(metaclass=ABCMeta): """ @abstractmethod - def inviscid_divergence_flux(self, dcoll, btag, gas_model, state_minus, + def inviscid_divergence_flux(self, dcoll, dd_bdry, gas_model, state_minus, numerical_flux_func, **kwargs): """Get the inviscid boundary flux for the divergence operator. @@ -84,11 +86,12 @@ def inviscid_divergence_flux(self, dcoll, btag, gas_model, state_minus, Fluid state object with the conserved state, and dependent quantities for the (-) side of the boundary specified by - *btag*. + *dd_bdry*. - btag: + dd_bdry: - Boundary tag indicating which domain boundary to process + Boundary DOF descriptor (or object convertible to one) indicating which + domain boundary to process gas_model: :class:`~mirgecom.gas_model.GasModel` @@ -108,7 +111,7 @@ def inviscid_divergence_flux(self, dcoll, btag, gas_model, state_minus, """ @abstractmethod - def viscous_divergence_flux(self, dcoll, btag, gas_model, state_minus, + def viscous_divergence_flux(self, dcoll, dd_bdry, gas_model, state_minus, grad_cv_minus, grad_t_minus, numerical_flux_func, **kwargs): """Get the viscous boundary flux for the divergence operator. @@ -122,25 +125,26 @@ def viscous_divergence_flux(self, dcoll, btag, gas_model, state_minus, A discretization collection encapsulating the DG elements - btag: + dd_bdry: - Boundary tag indicating which domain boundary to process + Boundary DOF descriptor (or object convertible to one) indicating which + domain boundary to process state_minus: :class:`~mirgecom.gas_model.FluidState` Fluid state object with the conserved state, and dependent quantities for the (-) side of the boundary specified - by *btag*. + by *dd_bdry*. grad_cv_minus: :class:`~mirgecom.fluid.ConservedVars` The gradient of the conserved quantities on the (-) side - of the boundary specified by *btag*. + of the boundary specified by *dd_bdry*. grad_t_minus: numpy.ndarray The gradient of the fluid temperature on the (-) side - of the boundary specified by *btag*. + of the boundary specified by *dd_bdry*. gas_model: :class:`~mirgecom.gas_model.GasModel` @@ -160,7 +164,7 @@ def viscous_divergence_flux(self, dcoll, btag, gas_model, state_minus, """ @abstractmethod - def cv_gradient_flux(self, dcoll, btag, gas_model, state_minus, **kwargs): + def cv_gradient_flux(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): """Get the boundary flux for the gradient of the fluid conserved variables. This routine returns the facial flux used by the gradient operator to @@ -172,15 +176,16 @@ def cv_gradient_flux(self, dcoll, btag, gas_model, state_minus, **kwargs): A discretization collection encapsulating the DG elements - btag: + dd_bdry: - Boundary tag indicating which domain boundary to process + Boundary DOF descriptor (or object convertible to one) indicating which + domain boundary to process state_minus: :class:`~mirgecom.gas_model.FluidState` Fluid state object with the conserved state, and dependent quantities for the (-) side of the boundary specified by - *btag*. + *dd_bdry*. gas_model: :class:`~mirgecom.gas_model.GasModel` @@ -193,7 +198,7 @@ def cv_gradient_flux(self, dcoll, btag, gas_model, state_minus, **kwargs): """ @abstractmethod - def temperature_gradient_flux(self, dcoll, btag, gas_model, state_minus, + def temperature_gradient_flux(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): """Get the boundary flux for the gradient of the fluid temperature. @@ -207,15 +212,16 @@ def temperature_gradient_flux(self, dcoll, btag, gas_model, state_minus, A discretization collection encapsulating the DG elements - btag: + dd_bdry: - Boundary tag indicating which domain boundary to process + Boundary DOF descriptor (or object convertible to one) indicating which + domain boundary to process state_minus: :class:`~mirgecom.gas_model.FluidState` Fluid state object with the conserved state, and dependent quantities for the (-) side of the boundary specified by - *btag*. + *dd_bdry*. gas_model: :class:`~mirgecom.gas_model.GasModel` @@ -310,17 +316,16 @@ def __init__(self, if not self._bnd_grad_temperature_func: self._bnd_grad_temperature_func = self._identical_grad_temperature - def _boundary_quantity(self, dcoll, btag, quantity, local=False, **kwargs): + def _boundary_quantity(self, dcoll, dd_bdry, 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) + dd_allfaces = dd_bdry.with_boundary_tag(FACE_RESTR_ALL) return quantity if local else op.project(dcoll, - btag, btag.with_dtag("all_faces"), quantity) + dd_bdry, dd_allfaces, quantity) - def _boundary_state_pair(self, dcoll, btag, gas_model, state_minus, **kwargs): - return TracePair(btag, + def _boundary_state_pair(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): + return TracePair(dd_bdry, interior=state_minus, - exterior=self._bnd_state_func(dcoll=dcoll, btag=btag, + exterior=self._bnd_state_func(dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model, state_minus=state_minus, **kwargs)) @@ -332,9 +337,9 @@ def _boundary_state_pair(self, dcoll, btag, gas_model, state_minus, **kwargs): # {{{ Default boundary helpers # Returns temperature(+) for boundaries that prescribe CV(+) - def _temperature_for_prescribed_state(self, dcoll, btag, + def _temperature_for_prescribed_state(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): - boundary_state = self._bnd_state_func(dcoll=dcoll, btag=btag, + boundary_state = self._bnd_state_func(dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model, state_minus=state_minus, **kwargs) @@ -346,38 +351,38 @@ def _identical_state(self, state_minus, **kwargs): def _identical_grad_cv(self, grad_cv_minus, **kwargs): return grad_cv_minus - def _identical_grad_temperature(self, grad_t_minus, **kwargs): + def _identical_grad_temperature(self, dcoll, dd_bdry, grad_t_minus, **kwargs): return grad_t_minus # Returns the flux to be used by the gradient operator when computing the # gradient of the fluid solution on boundaries that prescribe CV(+). - def _gradient_flux_for_prescribed_cv(self, dcoll, btag, gas_model, state_minus, - **kwargs): + def _gradient_flux_for_prescribed_cv(self, dcoll, dd_bdry, gas_model, + state_minus, **kwargs): # Use prescribed external state and gradient numerical flux function - boundary_state = self._bnd_state_func(dcoll=dcoll, btag=btag, + boundary_state = self._bnd_state_func(dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model, state_minus=state_minus, **kwargs) - cv_pair = TracePair(btag, + cv_pair = TracePair(dd_bdry, interior=state_minus.cv, exterior=boundary_state.cv) actx = state_minus.array_context - nhat = actx.thaw(dcoll.normal(btag)) + nhat = actx.thaw(dcoll.normal(dd_bdry)) from arraycontext import outer return outer(self._grad_num_flux_func(cv_pair.int, cv_pair.ext), nhat) # Returns the flux to be used by the gradient operator when computing the # gradient of fluid temperature using prescribed fluid temperature(+). - def _gradient_flux_for_prescribed_temperature(self, dcoll, btag, gas_model, + def _gradient_flux_for_prescribed_temperature(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): # Feed a boundary temperature to numerical flux for grad op actx = state_minus.array_context - nhat = actx.thaw(dcoll.normal(btag)) - bnd_tpair = TracePair(btag, + nhat = actx.thaw(dcoll.normal(dd_bdry)) + bnd_tpair = TracePair(dd_bdry, interior=state_minus.temperature, exterior=self._bnd_temperature_func( - dcoll=dcoll, btag=btag, gas_model=gas_model, + dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model, state_minus=state_minus, **kwargs)) from arraycontext import outer return outer(self._grad_num_flux_func(bnd_tpair.int, bnd_tpair.ext), nhat) @@ -386,39 +391,39 @@ def _gradient_flux_for_prescribed_temperature(self, dcoll, btag, gas_model, # divergence of inviscid fluid transport flux using the boundary's # prescribed CV(+). def _inviscid_flux_for_prescribed_state( - self, dcoll, btag, gas_model, state_minus, + self, dcoll, dd_bdry, gas_model, state_minus, numerical_flux_func=inviscid_facial_flux_rusanov, **kwargs): # Use a prescribed boundary state and the numerical flux function - boundary_state_pair = self._boundary_state_pair(dcoll=dcoll, btag=btag, + boundary_state_pair = self._boundary_state_pair(dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model, state_minus=state_minus, **kwargs) - normal = state_minus.array_context.thaw(dcoll.normal(btag)) + normal = state_minus.array_context.thaw(dcoll.normal(dd_bdry)) return numerical_flux_func(boundary_state_pair, gas_model, normal) # Returns the flux to be used by the divergence operator when computing the # divergence of viscous fluid transport flux using the boundary's # prescribed CV(+). def _viscous_flux_for_prescribed_state( - self, dcoll, btag, gas_model, state_minus, grad_cv_minus, grad_t_minus, - numerical_flux_func=viscous_facial_flux_central, **kwargs): + self, dcoll, dd_bdry, gas_model, state_minus, grad_cv_minus, + grad_t_minus, numerical_flux_func=viscous_facial_flux_central, **kwargs): state_pair = self._boundary_state_pair( - dcoll=dcoll, btag=btag, gas_model=gas_model, state_minus=state_minus, - **kwargs) + dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model, + state_minus=state_minus, **kwargs) grad_cv_pair = \ - TracePair(btag, interior=grad_cv_minus, + TracePair(dd_bdry, interior=grad_cv_minus, exterior=self._bnd_grad_cv_func( - dcoll=dcoll, btag=btag, gas_model=gas_model, + dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model, state_minus=state_minus, grad_cv_minus=grad_cv_minus, grad_t_minus=grad_t_minus)) grad_t_pair = \ TracePair( - btag, interior=grad_t_minus, + dd_bdry, interior=grad_t_minus, exterior=self._bnd_grad_temperature_func( - dcoll=dcoll, btag=btag, gas_model=gas_model, + dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model, state_minus=state_minus, grad_cv_minus=grad_cv_minus, grad_t_minus=grad_t_minus)) @@ -428,32 +433,37 @@ def _viscous_flux_for_prescribed_state( # }}} Default boundary helpers - def inviscid_divergence_flux(self, dcoll, btag, gas_model, state_minus, + def inviscid_divergence_flux(self, dcoll, dd_bdry, gas_model, state_minus, numerical_flux_func=inviscid_facial_flux_rusanov, **kwargs): """Get the inviscid boundary flux for the divergence operator.""" - return self._inviscid_flux_func(dcoll, btag, gas_model, state_minus, + dd_bdry = as_dofdesc(dd_bdry) + return self._inviscid_flux_func(dcoll, dd_bdry, gas_model, state_minus, numerical_flux_func=numerical_flux_func, **kwargs) - def cv_gradient_flux(self, dcoll, btag, gas_model, state_minus, **kwargs): - """Get the cv flux for *btag* for use in the gradient operator.""" + def cv_gradient_flux(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): + """Get the flux for *dd_bdry* for use in grad(CV).""" + dd_bdry = as_dofdesc(dd_bdry) return self._cv_gradient_flux_func( - dcoll=dcoll, btag=btag, gas_model=gas_model, state_minus=state_minus, - **kwargs) + dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model, + state_minus=state_minus, **kwargs) - def temperature_gradient_flux(self, dcoll, btag, gas_model, state_minus, + def temperature_gradient_flux(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): - """Get the "temperature flux" for *btag* for use in the gradient operator.""" - return self._temperature_grad_flux_func(dcoll, btag, gas_model, state_minus, - **kwargs) + """Get the flux for *dd_bdry* for use in grad(T).""" + dd_bdry = as_dofdesc(dd_bdry) + return self._temperature_grad_flux_func(dcoll, dd_bdry, gas_model, + state_minus, **kwargs) - def viscous_divergence_flux(self, dcoll, btag, gas_model, state_minus, + def viscous_divergence_flux(self, dcoll, dd_bdry, gas_model, state_minus, grad_cv_minus, grad_t_minus, numerical_flux_func=viscous_facial_flux_central, **kwargs): - """Get the viscous flux for *btag* for use in the divergence operator.""" - return self._viscous_flux_func(dcoll=dcoll, btag=btag, gas_model=gas_model, + """Get the viscous flux for *dd_bdry* for use in the divergence operator.""" + dd_bdry = as_dofdesc(dd_bdry) + return self._viscous_flux_func(dcoll=dcoll, dd_bdry=dd_bdry, + gas_model=gas_model, state_minus=state_minus, grad_cv_minus=grad_cv_minus, grad_t_minus=grad_t_minus, @@ -465,17 +475,18 @@ def viscous_divergence_flux(self, dcoll, btag, gas_model, state_minus, def _identical_grad_av(self, grad_av_minus, **kwargs): return grad_av_minus - def av_flux(self, dcoll, btag, diffusion, **kwargs): + def av_flux(self, dcoll, dd_bdry, diffusion, **kwargs): """Get the diffusive fluxes for the AV operator API.""" - grad_av_minus = op.project(dcoll, "vol", btag, diffusion) + dd_bdry = as_dofdesc(dd_bdry) + grad_av_minus = op.project(dcoll, dd_bdry.untrace(), dd_bdry, diffusion) actx = grad_av_minus.mass[0].array_context - nhat = actx.thaw(dcoll.normal(btag)) + nhat = actx.thaw(dcoll.normal(dd_bdry)) grad_av_plus = self._bnd_grad_av_func( - dcoll=dcoll, btag=btag, grad_av_minus=grad_av_minus, **kwargs) - bnd_grad_pair = TracePair(btag, interior=grad_av_minus, + dcoll=dcoll, dd_bdry=dd_bdry, grad_av_minus=grad_av_minus, **kwargs) + bnd_grad_pair = TracePair(dd_bdry, 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(dcoll, btag, num_flux, **kwargs) + return self._boundary_quantity(dcoll, dd_bdry, num_flux, **kwargs) # }}} @@ -514,7 +525,7 @@ def __init__(self): boundary_grad_av_func=self.adiabatic_slip_grad_av ) - def adiabatic_slip_state(self, dcoll, btag, gas_model, state_minus, **kwargs): + def adiabatic_slip_state(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): """Get the exterior solution on the boundary. The exterior solution is set such that there will be vanishing @@ -530,7 +541,7 @@ def adiabatic_slip_state(self, dcoll, btag, gas_model, state_minus, **kwargs): actx = state_minus.array_context # Grab a unit normal to the boundary - nhat = actx.thaw(dcoll.normal(btag)) + nhat = actx.thaw(dcoll.normal(dd_bdry)) # Subtract out the 2*wall-normal component # of velocity from the velocity at the wall to @@ -547,12 +558,12 @@ def adiabatic_slip_state(self, dcoll, 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, dcoll, btag, grad_av_minus, **kwargs): + def adiabatic_slip_grad_av(self, dcoll, dd_bdry, 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 = actx.thaw(dcoll.normal(btag)) + nhat = actx.thaw(dcoll.normal(dd_bdry)) # Subtract 2*wall-normal component of q # to enforce q=0 on the wall @@ -587,7 +598,8 @@ def __init__(self, wall_velocity=None, dim=2): raise ValueError(f"Specified wall velocity must be {dim}-vector.") self._wall_velocity = wall_velocity - def adiabatic_noslip_state(self, dcoll, btag, gas_model, state_minus, **kwargs): + def adiabatic_noslip_state( + self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): """Get the exterior solution on the boundary. Sets the external state s.t. $v^+ = -v^-$, giving vanishing contact velocity @@ -624,7 +636,8 @@ def __init__(self, wall_temperature=300): boundary_temperature_func=self.temperature_bc ) - def isothermal_noslip_state(self, dcoll, btag, gas_model, state_minus, **kwargs): + def isothermal_noslip_state( + self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): r"""Get the interior and exterior solution (*state_minus*) on the boundary. Sets the external state s.t. $v^+ = -v^-$, giving vanishing contact velocity diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py index 2a3ab6fc0..6726e85e7 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -39,13 +39,19 @@ import numpy.linalg as la # noqa from pytools.obj_array import make_obj_array, obj_array_vectorize_n_args from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.dof_desc import DOFDesc, as_dofdesc, DISCR_TAG_BASE +from meshmode.discretization.connection import FACE_RESTR_ALL # noqa +from grudge.dof_desc import ( + DD_VOLUME_ALL, + VolumeDomainTag, + DISCR_TAG_BASE, +) from grudge.trace_pair import ( TracePair, interior_trace_pairs, tracepair_with_discr_tag, ) import grudge.op as op +from mirgecom.utils import normalize_boundaries def grad_facial_flux(u_tpair, normal): @@ -205,7 +211,8 @@ class _DiffusionGradTag: def grad_operator( - dcoll, boundaries, u, quadrature_tag=DISCR_TAG_BASE, comm_tag=None): + dcoll, boundaries, u, *, quadrature_tag=DISCR_TAG_BASE, dd=DD_VOLUME_ALL, + comm_tag=None): r""" Compute the gradient of *u*. @@ -224,6 +231,9 @@ def grad_operator( quadrature_tag: quadrature tag indicating which discretization in *dcoll* to use for overintegration + dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the discretization on which *u* lives. Must be a volume + on the base discretization. comm_tag: Hashable Tag for distributed communication @@ -239,17 +249,27 @@ def grad_operator( raise TypeError("boundaries must be the same length as u") return obj_array_vectorize_n_args( lambda boundaries_i, u_i: grad_operator( - dcoll, boundaries_i, u_i, quadrature_tag=quadrature_tag), + dcoll, boundaries_i, u_i, quadrature_tag=quadrature_tag, + dd=dd), make_obj_array(boundaries), u) actx = u.array_context - for btag, bdry in boundaries.items(): + boundaries = normalize_boundaries(boundaries) + + for bdtag, bdry in boundaries.items(): if not isinstance(bdry, DiffusionBoundary): - raise TypeError(f"Unrecognized boundary type for tag {btag}. " + raise TypeError(f"Unrecognized boundary type for tag {bdtag}. " "Must be an instance of DiffusionBoundary.") - dd_allfaces_quad = DOFDesc("all_faces", quadrature_tag) + if not isinstance(dd.domain_tag, VolumeDomainTag): + raise TypeError("dd must represent a volume") + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + + dd_vol = dd + dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) interp_to_surf_quad = partial(tracepair_with_discr_tag, dcoll, quadrature_tag) @@ -261,33 +281,37 @@ def interior_flux(u_tpair): dcoll, dd_trace_quad, dd_allfaces_quad, grad_facial_flux(u_tpair_quad, normal_quad)) - def boundary_flux(dd_bdry, bdry): - u_minus = op.project(dcoll, "vol", dd_bdry, u) + def boundary_flux(bdtag, bdry): + dd_bdry_quad = dd_vol_quad.with_domain_tag(bdtag) + u_minus_quad = op.project(dcoll, dd_vol, dd_bdry_quad, u) return op.project( - dcoll, dd_bdry, dd_allfaces_quad, - bdry.get_grad_flux(dcoll, dd_bdry, u_minus)) + dcoll, dd_bdry_quad, dd_allfaces_quad, + bdry.get_grad_flux(dcoll, dd_bdry_quad, u_minus_quad)) - return op.inverse_mass(dcoll, - op.weak_local_grad(dcoll, "vol", -u) + return op.inverse_mass( + dcoll, dd_vol, + op.weak_local_grad(dcoll, dd_vol, -u) - # noqa: W504 - op.face_mass(dcoll, - dd_allfaces_quad, + op.face_mass( + dcoll, dd_allfaces_quad, sum( interior_flux(u_tpair) for u_tpair in interior_trace_pairs( - dcoll, u, comm_tag=(_DiffusionStateTag, comm_tag))) + dcoll, u, volume_dd=dd_vol, + comm_tag=(_DiffusionStateTag, comm_tag))) + sum( - boundary_flux( - as_dofdesc(btag).with_discr_tag(quadrature_tag), - bdry) - for btag, bdry in boundaries.items()) + boundary_flux(bdtag, bdry) + for bdtag, bdry in boundaries.items()) ) ) def diffusion_operator( dcoll, kappa, boundaries, u, *, return_grad_u=False, - quadrature_tag=DISCR_TAG_BASE, comm_tag=None): + quadrature_tag=DISCR_TAG_BASE, dd=DD_VOLUME_ALL, comm_tag=None, + # Added to avoid repeated computation + # FIXME: See if there's a better way to do this + grad_u=None): r""" Compute the diffusion operator. @@ -304,7 +328,7 @@ def diffusion_operator( kappa: numbers.Number or meshmode.dof_array.DOFArray the conductivity value(s) boundaries: - dictionary (or list of dictionaries) mapping boundary tags to + dictionary (or list of dictionaries) mapping boundary domain tags to :class:`DiffusionBoundary` instances u: meshmode.dof_array.DOFArray or numpy.ndarray the DOF array (or object array of DOF arrays) to which the operator should be @@ -314,6 +338,9 @@ def diffusion_operator( quadrature_tag: quadrature tag indicating which discretization in *dcoll* to use for overintegration + dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the discretization on which *u* lives. Must be a volume + on the base discretization. comm_tag: Hashable Tag for distributed communication @@ -332,23 +359,33 @@ def diffusion_operator( return obj_array_vectorize_n_args( lambda boundaries_i, u_i: diffusion_operator( dcoll, kappa, boundaries_i, u_i, return_grad_u=return_grad_u, - quadrature_tag=quadrature_tag), + quadrature_tag=quadrature_tag, dd=dd), make_obj_array(boundaries), u) actx = u.array_context - for btag, bdry in boundaries.items(): + boundaries = normalize_boundaries(boundaries) + + for bdtag, bdry in boundaries.items(): if not isinstance(bdry, DiffusionBoundary): - raise TypeError(f"Unrecognized boundary type for tag {btag}. " + raise TypeError(f"Unrecognized boundary type for tag {bdtag}. " "Must be an instance of DiffusionBoundary.") - dd_quad = DOFDesc("vol", quadrature_tag) - dd_allfaces_quad = DOFDesc("all_faces", quadrature_tag) + if not isinstance(dd.domain_tag, VolumeDomainTag): + raise TypeError("dd must represent a volume") + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + + dd_vol = dd + dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) - grad_u = grad_operator(dcoll, boundaries, u, quadrature_tag=quadrature_tag) + if grad_u is None: + grad_u = grad_operator( + dcoll, boundaries, u, quadrature_tag=quadrature_tag, dd=dd_vol) - kappa_quad = op.project(dcoll, "vol", dd_quad, kappa) - grad_u_quad = op.project(dcoll, "vol", dd_quad, grad_u) + kappa_quad = op.project(dcoll, dd_vol, dd_vol_quad, kappa) + grad_u_quad = op.project(dcoll, dd_vol, dd_vol_quad, grad_u) interp_to_surf_quad = partial(tracepair_with_discr_tag, dcoll, quadrature_tag) @@ -361,31 +398,34 @@ def interior_flux(kappa_tpair, grad_u_tpair): dcoll, dd_trace_quad, dd_allfaces_quad, diffusion_facial_flux(kappa_tpair_quad, grad_u_tpair_quad, normal_quad)) - def boundary_flux(dd_bdry, bdry): - kappa_minus = op.project(dcoll, "vol", dd_bdry, kappa) - grad_u_minus = op.project(dcoll, "vol", dd_bdry, grad_u) + def boundary_flux(bdtag, bdry): + dd_bdry_quad = dd_vol_quad.with_domain_tag(bdtag) + kappa_minus_quad = op.project(dcoll, dd_vol, dd_bdry_quad, kappa) + grad_u_minus_quad = op.project(dcoll, dd_vol, dd_bdry_quad, grad_u) return op.project( - dcoll, dd_bdry, dd_allfaces_quad, + dcoll, dd_bdry_quad, dd_allfaces_quad, bdry.get_diffusion_flux( - dcoll, dd_bdry, kappa_minus, grad_u_minus)) + dcoll, dd_bdry_quad, kappa_minus_quad, grad_u_minus_quad)) - diff_u = op.inverse_mass(dcoll, - op.weak_local_div(dcoll, dd_quad, -kappa_quad*grad_u_quad) + diff_u = op.inverse_mass( + dcoll, dd_vol, + op.weak_local_div(dcoll, dd_vol_quad, -kappa_quad*grad_u_quad) - # noqa: W504 - op.face_mass(dcoll, - dd_allfaces_quad, + op.face_mass( + dcoll, dd_allfaces_quad, sum( interior_flux(kappa_tpair, grad_u_tpair) for kappa_tpair, grad_u_tpair in zip( interior_trace_pairs( - dcoll, kappa, comm_tag=(_DiffusionKappaTag, comm_tag)), + dcoll, kappa, volume_dd=dd_vol, + comm_tag=(_DiffusionKappaTag, comm_tag)), interior_trace_pairs( - dcoll, grad_u, comm_tag=(_DiffusionGradTag, comm_tag)))) + dcoll, grad_u, volume_dd=dd_vol, + comm_tag=(_DiffusionGradTag, comm_tag))) + ) + sum( - boundary_flux( - as_dofdesc(btag).with_discr_tag(quadrature_tag), - bdry) - for btag, bdry in boundaries.items()) + boundary_flux(bdtag, bdry) + for bdtag, bdry in boundaries.items()) ) ) diff --git a/mirgecom/discretization.py b/mirgecom/discretization.py index c592f8669..37677b7a3 100644 --- a/mirgecom/discretization.py +++ b/mirgecom/discretization.py @@ -38,11 +38,17 @@ # we can replace it more easily when we refactor the drivers and # examples to use discretization collections, and change it centrally # when we want to change it. -def create_discretization_collection(actx, mesh, order, *, mpi_communicator=None, - quadrature_order=-1): +def create_discretization_collection(actx, volume_meshes, order, *, + mpi_communicator=None, quadrature_order=-1): """Create and return a grudge DG discretization collection.""" + if mpi_communicator is not None: + from warnings import warn + warn( + "mpi_communicator argument is deprecated and will disappear in Q4 2022.", + DeprecationWarning, stacklevel=2) + from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD - from grudge.discretization import DiscretizationCollection + from grudge.discretization import make_discretization_collection from meshmode.discretization.poly_element import ( QuadratureSimplexGroupFactory, PolynomialRecursiveNodesGroupFactory @@ -51,12 +57,11 @@ def create_discretization_collection(actx, mesh, order, *, mpi_communicator=None if quadrature_order < 0: quadrature_order = 2*order+1 - return DiscretizationCollection( - actx, mesh, + return make_discretization_collection( + actx, volume_meshes, discr_tag_to_group_factory={ DISCR_TAG_BASE: PolynomialRecursiveNodesGroupFactory(order=order, family="lgl"), DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(quadrature_order), - }, - mpi_communicator=mpi_communicator + } ) diff --git a/mirgecom/euler.py b/mirgecom/euler.py index 5426fcccc..e583367a5 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -54,7 +54,12 @@ import numpy as np # noqa -from grudge.dof_desc import DOFDesc +from meshmode.discretization.connection import FACE_RESTR_ALL +from grudge.dof_desc import ( + DD_VOLUME_ALL, + VolumeDomainTag, + DISCR_TAG_BASE, +) from mirgecom.gas_model import make_operator_fluid_states from mirgecom.inviscid import ( @@ -64,11 +69,14 @@ ) from mirgecom.operators import div_operator +from mirgecom.utils import normalize_boundaries def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, inviscid_numerical_flux_func=inviscid_facial_flux_rusanov, - quadrature_tag=None, comm_tag=None, operator_states_quad=None): + quadrature_tag=DISCR_TAG_BASE, dd=DD_VOLUME_ALL, + comm_tag=None, + operator_states_quad=None): r"""Compute RHS of the Euler flow equations. Returns @@ -91,7 +99,8 @@ def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, boundaries - Dictionary of boundary functions, one for each valid btag + Dictionary of boundary functions, one for each valid + :class:`~grudge.dof_desc.BoundaryDomainTag` time @@ -106,17 +115,31 @@ def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, An optional identifier denoting a particular quadrature discretization to use during operator evaluations. - The default value is *None*. + + dd: grudge.dof_desc.DOFDesc + + the DOF descriptor of the discretization on which *state* lives. Must be a + volume on the base discretization. comm_tag: Hashable + Tag for distributed communication """ - dd_quad_vol = DOFDesc("vol", quadrature_tag) - dd_quad_faces = DOFDesc("all_faces", quadrature_tag) + boundaries = normalize_boundaries(boundaries) + + if not isinstance(dd.domain_tag, VolumeDomainTag): + raise TypeError("dd must represent a volume") + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + + dd_vol = dd + dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) if operator_states_quad is None: operator_states_quad = make_operator_fluid_states( - dcoll, state, gas_model, boundaries, quadrature_tag, comm_tag) + dcoll, state, gas_model, boundaries, quadrature_tag, + dd=dd_vol, comm_tag=comm_tag) volume_state_quad, interior_state_pairs_quad, domain_boundary_states_quad = \ operator_states_quad @@ -128,9 +151,10 @@ def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, inviscid_flux_bnd = inviscid_flux_on_element_boundary( dcoll, gas_model, boundaries, interior_state_pairs_quad, domain_boundary_states_quad, quadrature_tag=quadrature_tag, - numerical_flux_func=inviscid_numerical_flux_func, time=time) + numerical_flux_func=inviscid_numerical_flux_func, time=time, + dd=dd_vol) - return -div_operator(dcoll, dd_quad_vol, dd_quad_faces, + return -div_operator(dcoll, dd_vol_quad, dd_allfaces_quad, inviscid_flux_vol, inviscid_flux_bnd) diff --git a/mirgecom/filter.py b/mirgecom/filter.py index da09d9526..f2f12235f 100644 --- a/mirgecom/filter.py +++ b/mirgecom/filter.py @@ -46,11 +46,15 @@ """ import numpy as np -import grudge.dof_desc as dof_desc +from functools import partial -from arraycontext import map_array_container +from grudge.dof_desc import ( + DD_VOLUME_ALL, + DISCR_TAG_BASE, + DISCR_TAG_MODAL, +) -from functools import partial +from arraycontext import map_array_container from meshmode.dof_array import DOFArray @@ -167,7 +171,7 @@ def apply_spectral_filter(actx, modal_field, discr, cutoff, ) -def filter_modally(dcoll, dd, cutoff, mode_resp_func, field): +def filter_modally(dcoll, cutoff, mode_resp_func, field, *, dd=DD_VOLUME_ALL): """Stand-alone procedural interface to spectral filtering. For each element group in the discretization, and restriction, @@ -185,15 +189,15 @@ def filter_modally(dcoll, dd, cutoff, mode_resp_func, field): ---------- dcoll: :class:`grudge.discretization.DiscretizationCollection` Grudge discretization with boundaries object - dd: :class:`grudge.dof_desc.DOFDesc` or as accepted by - :func:`grudge.dof_desc.as_dofdesc` - Describe the type of DOF vector on which to operate. cutoff: int Mode below which *field* will not be filtered mode_resp_func: Modal response function returns a filter coefficient for input mode id field: :class:`mirgecom.fluid.ConservedVars` An array container containing the relevant field(s) to filter. + dd: grudge.dof_desc.DOFDesc + Describe the type of DOF vector on which to operate. Must be on the base + discretization. Returns ------- @@ -202,16 +206,21 @@ def filter_modally(dcoll, dd, cutoff, mode_resp_func, field): """ if not isinstance(field, DOFArray): return map_array_container( - partial(filter_modally, dcoll, dd, cutoff, mode_resp_func), field + partial(filter_modally, dcoll, cutoff, mode_resp_func, dd=dd), field ) + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + + dd_nodal = dd + dd_modal = dd_nodal.with_discr_tag(DISCR_TAG_MODAL) + + discr = dcoll.discr_from_dd(dd_nodal) + actx = field.array_context - dd = dof_desc.as_dofdesc(dd) - dd_modal = dof_desc.DD_VOLUME_MODAL - discr = dcoll.discr_from_dd(dd) - modal_map = dcoll.connection_from_dds(dd, dd_modal) - nodal_map = dcoll.connection_from_dds(dd_modal, dd) + modal_map = dcoll.connection_from_dds(dd_nodal, dd_modal) + nodal_map = dcoll.connection_from_dds(dd_modal, dd_nodal) field = modal_map(field) field = apply_spectral_filter(actx, field, discr, cutoff, mode_resp_func) diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py index 18d2fbbff..f11e7aa14 100644 --- a/mirgecom/gas_model.py +++ b/mirgecom/gas_model.py @@ -59,12 +59,17 @@ TransportModel, GasTransportVars ) -from grudge.dof_desc import DOFDesc, as_dofdesc +from grudge.dof_desc import ( + DD_VOLUME_ALL, + VolumeDomainTag, + DISCR_TAG_BASE, +) import grudge.op as op from grudge.trace_pair import ( interior_trace_pairs, tracepair_with_discr_tag ) +from mirgecom.utils import normalize_boundaries @dataclass(frozen=True) @@ -384,8 +389,9 @@ class _FluidTemperatureTag: pass -def make_operator_fluid_states(dcoll, volume_state, gas_model, boundaries, - quadrature_tag=None, comm_tag=None): +def make_operator_fluid_states( + dcoll, volume_state, gas_model, boundaries, quadrature_tag=DISCR_TAG_BASE, + dd=DD_VOLUME_ALL, comm_tag=None): """Prepare gas model-consistent fluid states for use in fluid operators. This routine prepares a model-consistent fluid state for each of the volume and @@ -414,12 +420,16 @@ def make_operator_fluid_states(dcoll, volume_state, gas_model, boundaries, The physical model constructs for the gas_model boundaries - Dictionary of boundary functions, one for each valid btag + Dictionary of boundary functions, one for each valid + :class:`~grudge.dof_desc.BoundaryDomainTag`. quadrature_tag - An optional identifier denoting a particular quadrature - discretization to use during operator evaluations. - The default value is *None*. + An identifier denoting a particular quadrature discretization to use during + operator evaluations. + + dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the discretization on which *volume_state* lives. Must + be a volume on the base discretization. comm_tag: Hashable Tag for distributed communication @@ -431,19 +441,27 @@ def make_operator_fluid_states(dcoll, volume_state, gas_model, boundaries, Thermally consistent fluid state for the volume, fluid state trace pairs for the internal boundaries, and a dictionary of fluid states keyed by - valid btag in boundaries, all on the quadrature grid (if specified). + boundary domain tags in *boundaries*, all on the quadrature grid (if + specified). """ - dd_base_vol = DOFDesc("vol") - dd_quad_vol = DOFDesc("vol", quadrature_tag) + boundaries = normalize_boundaries(boundaries) + + if not isinstance(dd.domain_tag, VolumeDomainTag): + raise TypeError("dd must represent a volume") + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + + dd_vol = dd + dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag) # project pair to the quadrature discretization and update dd to quad interp_to_surf_quad = partial(tracepair_with_discr_tag, dcoll, quadrature_tag) domain_boundary_states_quad = { - btag: project_fluid_state(dcoll, dd_base_vol, - as_dofdesc(btag).with_discr_tag(quadrature_tag), + bdtag: project_fluid_state(dcoll, dd_vol, + dd_vol_quad.with_domain_tag(bdtag), volume_state, gas_model) - for btag in boundaries + for bdtag in boundaries } # performs MPI communication of CV if needed @@ -451,7 +469,8 @@ def make_operator_fluid_states(dcoll, volume_state, gas_model, boundaries, # Get the interior trace pairs onto the surface quadrature # discretization (if any) interp_to_surf_quad(tpair=tpair) - for tpair in interior_trace_pairs(dcoll, volume_state.cv, + for tpair in interior_trace_pairs( + dcoll, volume_state.cv, volume_dd=dd_vol, comm_tag=(_FluidCVTag, comm_tag)) ] @@ -466,7 +485,7 @@ def make_operator_fluid_states(dcoll, volume_state, gas_model, boundaries, # discretization (if any) interp_to_surf_quad(tpair=tpair) for tpair in interior_trace_pairs( - dcoll, volume_state.temperature, + dcoll, volume_state.temperature, volume_dd=dd_vol, comm_tag=(_FluidTemperatureTag, comm_tag))] interior_boundary_states_quad = \ @@ -475,7 +494,7 @@ def make_operator_fluid_states(dcoll, volume_state, gas_model, boundaries, # Interpolate the fluid state to the volume quadrature grid # (this includes the conserved and dependent quantities) - volume_state_quad = project_fluid_state(dcoll, dd_base_vol, dd_quad_vol, + volume_state_quad = project_fluid_state(dcoll, dd_vol, dd_vol_quad, volume_state, gas_model) return \ diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py index 6cdc7e187..d70228a87 100644 --- a/mirgecom/inviscid.py +++ b/mirgecom/inviscid.py @@ -40,8 +40,15 @@ """ import numpy as np +from meshmode.discretization.connection import FACE_RESTR_ALL +from grudge.dof_desc import ( + DD_VOLUME_ALL, + VolumeDomainTag, + DISCR_TAG_BASE, +) import grudge.op as op from mirgecom.fluid import make_conserved +from mirgecom.utils import normalize_boundaries def inviscid_flux(state): @@ -224,8 +231,9 @@ def inviscid_facial_flux_hll(state_pair, gas_model, normal): def inviscid_flux_on_element_boundary( dcoll, gas_model, boundaries, interior_state_pairs, - domain_boundary_states, quadrature_tag=None, - numerical_flux_func=inviscid_facial_flux_rusanov, time=0.0): + domain_boundary_states, quadrature_tag=DISCR_TAG_BASE, + numerical_flux_func=inviscid_facial_flux_rusanov, time=0.0, + dd=DD_VOLUME_ALL): """Compute the inviscid boundary fluxes for the divergence operator. This routine encapsulates the computation of the inviscid contributions @@ -242,40 +250,54 @@ def inviscid_flux_on_element_boundary( The physical model constructs for the gas_model boundaries - Dictionary of boundary functions, one for each valid btag + Dictionary of boundary functions, one for each valid + :class:`~grudge.dof_desc.BoundaryDomainTag` interior_state_pairs A :class:`~mirgecom.gas_model.FluidState` TracePair for each internal face. domain_boundary_states A dictionary of boundary-restricted :class:`~mirgecom.gas_model.FluidState`, - keyed by btags in *boundaries*. + keyed by boundary domain tags in *boundaries*. quadrature_tag - An optional identifier denoting a particular quadrature - discretization to use during operator evaluations. - The default value is *None*. + An identifier denoting a particular quadrature discretization to use during + operator evaluations. numerical_flux_func The numerical flux function to use in computing the boundary flux. time: float Time + + dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the discretization on which the fluid lives. Must be + a volume on the base discretization. """ - from grudge.dof_desc import as_dofdesc + boundaries = normalize_boundaries(boundaries) + + if not isinstance(dd.domain_tag, VolumeDomainTag): + raise TypeError("dd must represent a volume") + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + + dd_vol = dd + dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) def _interior_flux(state_pair): return op.project(dcoll, - state_pair.dd, state_pair.dd.with_dtag("all_faces"), + state_pair.dd, dd_allfaces_quad, numerical_flux_func( state_pair, gas_model, state_pair.int.array_context.thaw(dcoll.normal(state_pair.dd)))) - def _boundary_flux(dd_bdry, boundary, state_minus): + def _boundary_flux(bdtag, boundary, state_minus_quad): + dd_bdry_quad = dd_vol_quad.with_domain_tag(bdtag) return op.project(dcoll, - dd_bdry, dd_bdry.with_dtag("all_faces"), + dd_bdry_quad, dd_allfaces_quad, boundary.inviscid_divergence_flux( - dcoll, dd_bdry, gas_model, state_minus=state_minus, + dcoll, dd_bdry_quad, gas_model, state_minus=state_minus_quad, numerical_flux_func=numerical_flux_func, time=time)) # Compute interface contributions @@ -287,16 +309,16 @@ def _boundary_flux(dd_bdry, boundary, state_minus): # Domain boundary faces + sum( _boundary_flux( - as_dofdesc(btag).with_discr_tag(quadrature_tag), + bdtag, boundary, - domain_boundary_states[btag]) - for btag, boundary in boundaries.items()) + domain_boundary_states[bdtag]) + for bdtag, boundary in boundaries.items()) ) return inviscid_flux_bnd -def get_inviscid_timestep(dcoll, state): +def get_inviscid_timestep(dcoll, state, dd=DD_VOLUME_ALL): """Return node-local stable timestep estimate for an inviscid fluid. The maximum stable timestep is computed from the acoustic wavespeed. @@ -311,15 +333,25 @@ def get_inviscid_timestep(dcoll, state): Full fluid conserved and thermal state + dd: grudge.dof_desc.DOFDesc + + the DOF descriptor of the discretization on which *state* lives. Must be + a volume on the base discretization. + Returns ------- class:`~meshmode.dof_array.DOFArray` The maximum stable timestep at each node. """ + if not isinstance(dd.domain_tag, VolumeDomainTag): + raise TypeError("dd must represent a volume") + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + from grudge.dt_utils import characteristic_lengthscales return ( - characteristic_lengthscales(state.array_context, dcoll) + characteristic_lengthscales(state.array_context, dcoll, dd=dd) / state.wavespeed ) diff --git a/mirgecom/io.py b/mirgecom/io.py index 8cb5b1829..0adc69b9c 100644 --- a/mirgecom/io.py +++ b/mirgecom/io.py @@ -32,6 +32,8 @@ import grudge.op as op from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from grudge.dof_desc import DD_VOLUME_ALL + def make_init_message(*, dim, order, dt, t_final, nstatus, nviz, cfl, constant_cfl, @@ -52,11 +54,12 @@ def make_init_message(*, dim, order, dt, t_final, ) -def make_status_message(*, dcoll, t, step, dt, cfl, dependent_vars): +def make_status_message( + *, dcoll, t, step, dt, cfl, dependent_vars, fluid_dd=DD_VOLUME_ALL): r"""Make simulation status and health message.""" dv = dependent_vars - _min = partial(op.nodal_min, dcoll, "vol") - _max = partial(op.nodal_max, dcoll, "vol") + _min = partial(op.nodal_min, dcoll, fluid_dd) + _max = partial(op.nodal_max, dcoll, fluid_dd) statusmsg = ( f"Status: {step=} {t=}\n" f"------- P({_min(dv.pressure):.3g}, {_max(dv.pressure):.3g})\n" diff --git a/mirgecom/limiter.py b/mirgecom/limiter.py index 3625928e2..f3822f7c2 100644 --- a/mirgecom/limiter.py +++ b/mirgecom/limiter.py @@ -32,11 +32,13 @@ from pytools import memoize_in from grudge.discretization import DiscretizationCollection +from grudge.dof_desc import DD_VOLUME_ALL import grudge.op as op def bound_preserving_limiter(dcoll: DiscretizationCollection, field, - mmin=0.0, mmax=None, modify_average=False): + mmin=0.0, mmax=None, modify_average=False, + dd=DD_VOLUME_ALL): r"""Implement a slope limiter for bound-preserving properties. The implementation is summarized in [Zhang_2011]_, Sec. 2.3, Eq. 2.9, @@ -76,6 +78,8 @@ def bound_preserving_limiter(dcoll: DiscretizationCollection, field, Optional float with the target upper bound. Default to None. modify_average: bool Flag to avoid modification the cell average. Defaults to False. + dd: grudge.dof_desc.DOFDesc + The DOF descriptor corresponding to *field*. Returns ------- @@ -84,21 +88,21 @@ def bound_preserving_limiter(dcoll: DiscretizationCollection, field, """ actx = field.array_context - @memoize_in(dcoll, (bound_preserving_limiter, "cell_volume")) + @memoize_in(dcoll, (bound_preserving_limiter, "cell_volume", dd)) def cell_volumes(dcoll): - return op.elementwise_integral(dcoll, dcoll.zeros(actx) + 1.0) + return op.elementwise_integral(dcoll, dd, dcoll.zeros(actx) + 1.0) cell_size = cell_volumes(dcoll) # Compute cell averages of the state - cell_avgs = 1.0/cell_size*op.elementwise_integral(dcoll, field) + cell_avgs = 1.0/cell_size*op.elementwise_integral(dcoll, dd, field) # Bound cell average in case it doesn't respect the realizability if modify_average: cell_avgs = actx.np.where(actx.np.greater(cell_avgs, mmin), cell_avgs, mmin) # Compute elementwise max/mins of the field - mmin_i = op.elementwise_min(dcoll, field) + mmin_i = op.elementwise_min(dcoll, dd, field) # Linear scaling of polynomial coefficients _theta = actx.np.minimum( @@ -110,7 +114,7 @@ def cell_volumes(dcoll): if mmax is not None: cell_avgs = actx.np.where(actx.np.greater(cell_avgs, mmax), mmax, cell_avgs) - mmax_i = op.elementwise_max(dcoll, field) + mmax_i = op.elementwise_max(dcoll, dd, field) _theta = actx.np.minimum( _theta, actx.np.where(actx.np.greater(mmax_i, mmax), diff --git a/mirgecom/logging_quantities.py b/mirgecom/logging_quantities.py index 1ff9b6813..2ea3974f8 100644 --- a/mirgecom/logging_quantities.py +++ b/mirgecom/logging_quantities.py @@ -50,6 +50,7 @@ from typing import Optional, Callable import numpy as np +from grudge.dof_desc import DD_VOLUME_ALL import grudge.op as oper @@ -102,24 +103,28 @@ def logmgr_add_device_memory_usage(logmgr: LogManager, queue: cl.CommandQueue): def logmgr_add_many_discretization_quantities(logmgr: LogManager, dcoll, dim, - extract_vars_for_logging, units_for_logging): + extract_vars_for_logging, units_for_logging, dd=DD_VOLUME_ALL): """Add default discretization quantities to the logmgr.""" + if dd != DD_VOLUME_ALL: + suffix = f"_{dd.domain_tag.tag}" + else: + suffix = "" + for reduction_op in ["min", "max", "L2_norm"]: - for quantity in ["pressure", "temperature"]: + for quantity in ["pressure"+suffix, "temperature"+suffix]: logmgr.add_quantity(DiscretizationBasedQuantity( dcoll, quantity, reduction_op, extract_vars_for_logging, - units_for_logging)) + units_for_logging, dd=dd)) - for quantity in ["mass", "energy"]: + for quantity in ["mass"+suffix, "energy"+suffix]: logmgr.add_quantity(DiscretizationBasedQuantity( dcoll, quantity, reduction_op, extract_vars_for_logging, - units_for_logging)) + units_for_logging, dd=dd)) for d in range(dim): logmgr.add_quantity(DiscretizationBasedQuantity( - dcoll, "momentum", reduction_op, extract_vars_for_logging, - units_for_logging, - axis=d)) + dcoll, "momentum"+suffix, reduction_op, extract_vars_for_logging, + units_for_logging, axis=d, dd=dd)) # {{{ Package versions @@ -245,7 +250,7 @@ class DiscretizationBasedQuantity(PostLogQuantity, StateConsumer): def __init__(self, dcoll: DiscretizationCollection, quantity: str, op: str, extract_vars_for_logging, units_logging, name: str = None, - axis: Optional[int] = None): + axis: Optional[int] = None, dd=DD_VOLUME_ALL): unit = units_logging(quantity) if name is None: @@ -262,13 +267,13 @@ def __init__(self, dcoll: DiscretizationCollection, quantity: str, op: str, from functools import partial if op == "min": - self._discr_reduction = partial(oper.nodal_min, self.dcoll, "vol") + self._discr_reduction = partial(oper.nodal_min, self.dcoll, dd) self.rank_aggr = min elif op == "max": - self._discr_reduction = partial(oper.nodal_max, self.dcoll, "vol") + self._discr_reduction = partial(oper.nodal_max, self.dcoll, dd) self.rank_aggr = max elif op == "L2_norm": - self._discr_reduction = partial(oper.norm, self.dcoll, p=2) + self._discr_reduction = partial(oper.norm, self.dcoll, p=2, dd=dd) self.rank_aggr = max else: raise ValueError(f"unknown operation {op}") diff --git a/mirgecom/navierstokes.py b/mirgecom/navierstokes.py index 01bbe3da6..91a01259a 100644 --- a/mirgecom/navierstokes.py +++ b/mirgecom/navierstokes.py @@ -59,12 +59,18 @@ from functools import partial +from meshmode.discretization.connection import FACE_RESTR_ALL + from grudge.trace_pair import ( TracePair, interior_trace_pairs, tracepair_with_discr_tag ) -from grudge.dof_desc import DOFDesc, as_dofdesc, DISCR_TAG_BASE +from grudge.dof_desc import ( + DD_VOLUME_ALL, + VolumeDomainTag, + DISCR_TAG_BASE, +) import grudge.op as op @@ -84,6 +90,7 @@ div_operator, grad_operator ) from mirgecom.gas_model import make_operator_fluid_states +from mirgecom.utils import normalize_boundaries class _NSGradCVTag: @@ -98,16 +105,17 @@ def _gradient_flux_interior(dcoll, numerical_flux_func, tpair): """Compute interior face flux for gradient operator.""" from arraycontext import outer actx = tpair.int.array_context - dd = tpair.dd - normal = actx.thaw(dcoll.normal(dd)) + dd_trace = tpair.dd + dd_allfaces = dd_trace.with_boundary_tag(FACE_RESTR_ALL) + normal = actx.thaw(dcoll.normal(dd_trace)) flux = outer(numerical_flux_func(tpair.int, tpair.ext), normal) - return op.project(dcoll, dd, dd.with_dtag("all_faces"), flux) + return op.project(dcoll, dd_trace, dd_allfaces, flux) def grad_cv_operator( dcoll, gas_model, boundaries, state, *, time=0.0, numerical_flux_func=num_flux_central, - quadrature_tag=DISCR_TAG_BASE, comm_tag=None, + quadrature_tag=DISCR_TAG_BASE, dd=DD_VOLUME_ALL, comm_tag=None, # Added to avoid repeated computation # FIXME: See if there's a better way to do this operator_states_quad=None): @@ -121,7 +129,8 @@ def grad_cv_operator( quantities. boundaries - Dictionary of boundary functions keyed by btags + Dictionary of boundary functions, one for each valid + :class:`~grudge.dof_desc.BoundaryDomainTag` time Time @@ -140,6 +149,10 @@ def grad_cv_operator( An identifier denoting a particular quadrature discretization to use during operator evaluations. + dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the discretization on which *state* lives. Must be a + volume on the base discretization. + comm_tag: Hashable Tag for distributed communication @@ -150,13 +163,21 @@ def grad_cv_operator( CV object with vector components representing the gradient of the fluid conserved variables. """ - dd_vol_quad = DOFDesc("vol", quadrature_tag) - dd_faces_quad = DOFDesc("all_faces", quadrature_tag) + boundaries = normalize_boundaries(boundaries) + + if not isinstance(dd.domain_tag, VolumeDomainTag): + raise TypeError("dd must represent a volume") + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + + dd_vol = dd + dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) if operator_states_quad is None: operator_states_quad = make_operator_fluid_states( dcoll, state, gas_model, boundaries, quadrature_tag, - comm_tag) + dd=dd_vol, comm_tag=comm_tag) vol_state_quad, inter_elem_bnd_states_quad, domain_bnd_states_quad = \ operator_states_quad @@ -173,18 +194,16 @@ def grad_cv_operator( # Domain boundaries sum(op.project( - dcoll, as_dofdesc(btag).with_discr_tag(quadrature_tag), - as_dofdesc(btag).with_discr_tag(quadrature_tag).with_dtag("all_faces"), + dcoll, dd_vol_quad.with_domain_tag(bdtag), + dd_allfaces_quad, bdry.cv_gradient_flux( dcoll, - # Make sure we get the state on the quadrature grid - # restricted to the tag *btag* - as_dofdesc(btag).with_discr_tag(quadrature_tag), + dd_vol_quad.with_domain_tag(bdtag), gas_model=gas_model, - state_minus=domain_bnd_states_quad[btag], + state_minus=domain_bnd_states_quad[bdtag], time=time, numerical_flux_func=numerical_flux_func)) - for btag, bdry in boundaries.items()) + for bdtag, bdry in boundaries.items()) # Interior boundaries + sum(get_interior_flux(tpair) for tpair in cv_interior_pairs) @@ -192,13 +211,13 @@ def grad_cv_operator( # [Bassi_1997]_ eqn 15 (s = grad_q) return grad_operator( - dcoll, dd_vol_quad, dd_faces_quad, vol_state_quad.cv, cv_flux_bnd) + dcoll, dd_vol_quad, dd_allfaces_quad, vol_state_quad.cv, cv_flux_bnd) def grad_t_operator( dcoll, gas_model, boundaries, state, *, time=0.0, numerical_flux_func=num_flux_central, - quadrature_tag=DISCR_TAG_BASE, comm_tag=None, + quadrature_tag=DISCR_TAG_BASE, dd=DD_VOLUME_ALL, comm_tag=None, # Added to avoid repeated computation # FIXME: See if there's a better way to do this operator_states_quad=None): @@ -231,6 +250,10 @@ def grad_t_operator( An identifier denoting a particular quadrature discretization to use during operator evaluations. + dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the discretization on which *state* lives. Must be a + volume on the base discretization. + comm_tag: Hashable Tag for distributed communication @@ -241,12 +264,21 @@ def grad_t_operator( Array of :class:`~meshmode.dof_array.DOFArray` representing the gradient of the fluid temperature. """ - dd_vol_quad = DOFDesc("vol", quadrature_tag) - dd_faces_quad = DOFDesc("all_faces", quadrature_tag) + boundaries = normalize_boundaries(boundaries) + + if not isinstance(dd.domain_tag, VolumeDomainTag): + raise TypeError("dd must represent a volume") + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + + dd_vol = dd + dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) if operator_states_quad is None: operator_states_quad = make_operator_fluid_states( - dcoll, state, gas_model, boundaries, quadrature_tag, comm_tag) + dcoll, state, gas_model, boundaries, quadrature_tag, + dd=dd_vol, comm_tag=comm_tag) vol_state_quad, inter_elem_bnd_states_quad, domain_bnd_states_quad = \ operator_states_quad @@ -267,18 +299,16 @@ def grad_t_operator( # Domain boundaries sum(op.project( - dcoll, as_dofdesc(btag).with_discr_tag(quadrature_tag), - as_dofdesc(btag).with_discr_tag(quadrature_tag).with_dtag("all_faces"), + dcoll, dd_vol_quad.with_domain_tag(bdtag), + dd_allfaces_quad, bdry.temperature_gradient_flux( dcoll, - # Make sure we get the state on the quadrature grid - # restricted to the tag *btag* - as_dofdesc(btag).with_discr_tag(quadrature_tag), + dd_vol_quad.with_domain_tag(bdtag), gas_model=gas_model, - state_minus=domain_bnd_states_quad[btag], + state_minus=domain_bnd_states_quad[bdtag], time=time, numerical_flux_func=numerical_flux_func)) - for btag, bdry in boundaries.items()) + for bdtag, bdry in boundaries.items()) # Interior boundaries + sum(get_interior_flux(tpair) for tpair in t_interior_pairs) @@ -286,15 +316,15 @@ def grad_t_operator( # Fluxes in-hand, compute the gradient of temperature return grad_operator( - dcoll, dd_vol_quad, dd_faces_quad, vol_state_quad.temperature, t_flux_bnd) + dcoll, dd_vol_quad, dd_allfaces_quad, vol_state_quad.temperature, t_flux_bnd) def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, inviscid_numerical_flux_func=inviscid_facial_flux_rusanov, gradient_numerical_flux_func=num_flux_central, viscous_numerical_flux_func=viscous_facial_flux_central, - quadrature_tag=DISCR_TAG_BASE, return_gradients=False, - comm_tag=None, + return_gradients=False, quadrature_tag=DISCR_TAG_BASE, + dd=DD_VOLUME_ALL, comm_tag=None, # Added to avoid repeated computation # FIXME: See if there's a better way to do this operator_states_quad=None, @@ -333,10 +363,22 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, Optional callable function to return the numerical flux to be used when computing gradients in the Navier-Stokes operator. + return_gradients + Optional boolean (defaults to false) indicating whether to return + $\nabla(\text{CV})$ and $\nabla(T)$ along with the RHS for the Navier-Stokes + equations. Useful for debugging and visualization. + quadrature_tag An identifier denoting a particular quadrature discretization to use during operator evaluations. + dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the discretization on which *state* lives. Must be a + volume on the base discretization. + + comm_tag: Hashable + Tag for distributed communication + operator_states_quad Optional iterable container providing the full fluid states (:class:`~mirgecom.gas_model.FluidState`) on the quadrature @@ -355,14 +397,6 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, provided, the operator will calculate it with :func:`~mirgecom.navierstokes.grad_t_operator`. - return_gradients - Optional boolean (defaults to false) indicating whether to return - $\nabla(\text{CV})$ and $\nabla(T)$ along with the RHS for the Navier-Stokes - equations. Useful for debugging and visualization. - - comm_tag: Hashable - Tag for distributed communication - Returns ------- :class:`mirgecom.fluid.ConservedVars` @@ -376,9 +410,16 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, if not state.is_viscous: raise ValueError("Navier-Stokes operator expects viscous gas model.") - dd_base = as_dofdesc("vol") - dd_vol_quad = DOFDesc("vol", quadrature_tag) - dd_faces_quad = DOFDesc("all_faces", quadrature_tag) + boundaries = normalize_boundaries(boundaries) + + if not isinstance(dd.domain_tag, VolumeDomainTag): + raise TypeError("dd must represent a volume") + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + + dd_vol = dd + dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) # Make model-consistent fluid state data (i.e. CV *and* DV) for: # - Volume: vol_state_quad @@ -389,7 +430,8 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, # otherwise they stay on the interpolatory/base domain. if operator_states_quad is None: operator_states_quad = make_operator_fluid_states( - dcoll, state, gas_model, boundaries, quadrature_tag, comm_tag) + dcoll, state, gas_model, boundaries, quadrature_tag, + dd=dd_vol, comm_tag=comm_tag) vol_state_quad, inter_elem_bnd_states_quad, domain_bnd_states_quad = \ operator_states_quad @@ -407,7 +449,7 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, grad_cv = grad_cv_operator( dcoll, gas_model, boundaries, state, time=time, numerical_flux_func=gradient_numerical_flux_func, - quadrature_tag=quadrature_tag, + quadrature_tag=quadrature_tag, dd=dd_vol, operator_states_quad=operator_states_quad) # Communicate grad(CV) and put it on the quadrature domain @@ -416,7 +458,7 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, # discretization (if any) interp_to_surf_quad(tpair=tpair) for tpair in interior_trace_pairs( - dcoll, grad_cv, comm_tag=(_NSGradCVTag, comm_tag)) + dcoll, grad_cv, volume_dd=dd_vol, comm_tag=(_NSGradCVTag, comm_tag)) ] # }}} Compute grad(CV) @@ -427,7 +469,7 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, grad_t = grad_t_operator( dcoll, gas_model, boundaries, state, time=time, numerical_flux_func=gradient_numerical_flux_func, - quadrature_tag=quadrature_tag, + quadrature_tag=quadrature_tag, dd=dd_vol, operator_states_quad=operator_states_quad) # Create the interior face trace pairs, perform MPI exchange, interp to quad @@ -436,7 +478,8 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, # discretization (if any) interp_to_surf_quad(tpair=tpair) for tpair in interior_trace_pairs( - dcoll, grad_t, comm_tag=(_NSGradTemperatureTag, comm_tag)) + dcoll, grad_t, volume_dd=dd_vol, + comm_tag=(_NSGradTemperatureTag, comm_tag)) ] # }}} compute grad(temperature) @@ -450,8 +493,8 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, # using field values on the quadrature grid viscous_flux(state=vol_state_quad, # Interpolate gradients to the quadrature grid - grad_cv=op.project(dcoll, dd_base, dd_vol_quad, grad_cv), - grad_t=op.project(dcoll, dd_base, dd_vol_quad, grad_t)) + grad_cv=op.project(dcoll, dd_vol, dd_vol_quad, grad_cv), + grad_t=op.project(dcoll, dd_vol, dd_vol_quad, grad_t)) # Compute the volume contribution of the inviscid flux terms # using field values on the quadrature grid @@ -466,16 +509,18 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, dcoll, gas_model, boundaries, inter_elem_bnd_states_quad, domain_bnd_states_quad, grad_cv, grad_cv_interior_pairs, grad_t, grad_t_interior_pairs, quadrature_tag=quadrature_tag, - numerical_flux_func=viscous_numerical_flux_func, time=time) + numerical_flux_func=viscous_numerical_flux_func, time=time, + dd=dd_vol) # All surface contributions from the inviscid fluxes - inviscid_flux_on_element_boundary( dcoll, gas_model, boundaries, inter_elem_bnd_states_quad, domain_bnd_states_quad, quadrature_tag=quadrature_tag, - numerical_flux_func=inviscid_numerical_flux_func, time=time) + numerical_flux_func=inviscid_numerical_flux_func, time=time, + dd=dd_vol) ) - ns_rhs = div_operator(dcoll, dd_vol_quad, dd_faces_quad, vol_term, bnd_term) + ns_rhs = div_operator(dcoll, dd_vol_quad, dd_allfaces_quad, vol_term, bnd_term) if return_gradients: return ns_rhs, grad_cv, grad_t return ns_rhs diff --git a/mirgecom/operators.py b/mirgecom/operators.py index 23a913334..3984f1656 100644 --- a/mirgecom/operators.py +++ b/mirgecom/operators.py @@ -30,8 +30,10 @@ import grudge.op as op +from grudge.dof_desc import DISCR_TAG_BASE -def grad_operator(dcoll, dd_vol, dd_faces, u, flux): + +def grad_operator(dcoll, dd_vol, dd_allfaces, u, flux): r"""Compute a DG gradient for the input *u* with flux given by *flux*. Parameters @@ -41,7 +43,7 @@ def grad_operator(dcoll, dd_vol, dd_faces, u, flux): dd_vol: grudge.dof_desc.DOFDesc the degree-of-freedom tag associated with the volume discretization. This determines the type of quadrature to be used. - dd_faces: grudge.dof_desc.DOFDesc + dd_allfaces: grudge.dof_desc.DOFDesc the degree-of-freedom tag associated with the surface discretization. This determines the type of quadrature to be used. u: meshmode.dof_array.DOFArray or numpy.ndarray @@ -57,12 +59,13 @@ def grad_operator(dcoll, dd_vol, dd_faces, u, flux): the dg gradient operator applied to *u* """ # pylint: disable=invalid-unary-operand-type - return - op.inverse_mass(dcoll, + return -op.inverse_mass( + dcoll, dd_vol.with_discr_tag(DISCR_TAG_BASE), op.weak_local_grad(dcoll, dd_vol, u) - - op.face_mass(dcoll, dd_faces, flux)) + - op.face_mass(dcoll, dd_allfaces, flux)) -def div_operator(dcoll, dd_vol, dd_faces, v, flux): +def div_operator(dcoll, dd_vol, dd_allfaces, v, flux): r"""Compute DG divergence of vector-valued func *v* with flux given by *flux*. Parameters @@ -72,7 +75,7 @@ def div_operator(dcoll, dd_vol, dd_faces, v, flux): dd_vol: grudge.dof_desc.DOFDesc the degree-of-freedom tag associated with the volume discretization. This determines the type of quadrature to be used. - dd_faces: grudge.dof_desc.DOFDesc + dd_allfaces: grudge.dof_desc.DOFDesc the degree-of-freedom tag associated with the surface discretization. This determines the type of quadrature to be used. v: numpy.ndarray @@ -88,6 +91,7 @@ def div_operator(dcoll, dd_vol, dd_faces, v, flux): the dg divergence operator applied to vector-valued function(s) *v*. """ # pylint: disable=invalid-unary-operand-type - return - op.inverse_mass(dcoll, + return -op.inverse_mass( + dcoll, dd_vol.with_discr_tag(DISCR_TAG_BASE), op.weak_local_div(dcoll, dd_vol, v) - - op.face_mass(dcoll, dd_faces, flux)) + - op.face_mass(dcoll, dd_allfaces, flux)) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index 952a22c78..e2e3943f5 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -66,6 +66,7 @@ from typing import List, Dict from grudge.discretization import DiscretizationCollection +from grudge.dof_desc import DD_VOLUME_ALL logger = logging.getLogger(__name__) @@ -91,7 +92,9 @@ def check_step(step, interval): return False -def get_sim_timestep(dcoll, state, t, dt, cfl, t_final, constant_cfl=False): +def get_sim_timestep( + dcoll, state, t, dt, cfl, t_final, constant_cfl=False, + fluid_dd=DD_VOLUME_ALL): """Return the maximum stable timestep for a typical fluid simulation. This routine returns *dt*, the users defined constant timestep, or *max_dt*, the @@ -122,6 +125,9 @@ def get_sim_timestep(dcoll, state, t, dt, cfl, t_final, constant_cfl=False): The current CFL number constant_cfl: bool True if running constant CFL mode + fluid_dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the discretization on which *state* lives. Must be a + volume on the base discretization. Returns ------- @@ -135,7 +141,7 @@ def get_sim_timestep(dcoll, state, t, dt, cfl, t_final, constant_cfl=False): from grudge.op import nodal_min mydt = state.array_context.to_numpy( cfl * nodal_min( - dcoll, "vol", + dcoll, fluid_dd, get_viscous_timestep(dcoll=dcoll, state=state)))[()] return min(t_remaining, mydt) @@ -335,7 +341,7 @@ def check_naninf_local(dcoll: DiscretizationCollection, dd: str, return not np.isfinite(s) -def compare_fluid_solutions(dcoll, red_state, blue_state): +def compare_fluid_solutions(dcoll, red_state, blue_state, *, dd=DD_VOLUME_ALL): """Return inf norm of (*red_state* - *blue_state*) for each component. .. note:: @@ -344,12 +350,12 @@ def compare_fluid_solutions(dcoll, red_state, blue_state): actx = red_state.array_context resid = red_state - blue_state resid_errs = actx.to_numpy( - flatten(componentwise_norms(dcoll, resid, order=np.inf), actx)) + flatten(componentwise_norms(dcoll, resid, order=np.inf, dd=dd), actx)) return resid_errs.tolist() -def componentwise_norms(dcoll, fields, order=np.inf): +def componentwise_norms(dcoll, fields, order=np.inf, *, dd=DD_VOLUME_ALL): """Return the *order*-norm for each component of *fields*. .. note:: @@ -357,15 +363,15 @@ def componentwise_norms(dcoll, fields, order=np.inf): """ if not isinstance(fields, DOFArray): return map_array_container( - partial(componentwise_norms, dcoll, order=order), fields) + partial(componentwise_norms, dcoll, order=order, dd=dd), fields) if len(fields) > 0: - return op.norm(dcoll, fields, order) + return op.norm(dcoll, fields, order, dd=dd) else: # FIXME: This work-around for #575 can go away after #569 return 0 -def max_component_norm(dcoll, fields, order=np.inf): +def max_component_norm(dcoll, fields, order=np.inf, *, dd=DD_VOLUME_ALL): """Return the max *order*-norm over the components of *fields*. .. note:: @@ -373,7 +379,7 @@ def max_component_norm(dcoll, fields, order=np.inf): """ actx = fields.array_context return max(actx.to_numpy(flatten( - componentwise_norms(dcoll, fields, order), actx))) + componentwise_norms(dcoll, fields, order, dd=dd), actx))) def generate_and_distribute_mesh(comm, generate_mesh): diff --git a/mirgecom/utils.py b/mirgecom/utils.py index 6e72a1e2b..c33fcb3c9 100644 --- a/mirgecom/utils.py +++ b/mirgecom/utils.py @@ -3,6 +3,7 @@ .. autoclass:: StatisticsAccumulator .. autofunction:: asdict_shallow .. autofunction:: force_evaluation +.. autofunction:: normalize_boundaries """ __copyright__ = """ @@ -117,3 +118,15 @@ def force_evaluation(actx, x): if actx is None: return x return actx.freeze_thaw(x) + + +def normalize_boundaries(boundaries): + """ + Normalize the keys of *boundaries*. + + Promotes boundary tags to :class:`grudge.dof_desc.BoundaryDomainTag`. + """ + from grudge.dof_desc import as_dofdesc + return { + as_dofdesc(key).domain_tag: bdry + for key, bdry in boundaries.items()} diff --git a/mirgecom/viscous.py b/mirgecom/viscous.py index 261c0d2dd..1790fd41f 100644 --- a/mirgecom/viscous.py +++ b/mirgecom/viscous.py @@ -45,6 +45,12 @@ import numpy as np from meshmode.dof_array import DOFArray +from meshmode.discretization.connection import FACE_RESTR_ALL +from grudge.dof_desc import ( + DD_VOLUME_ALL, + VolumeDomainTag, + DISCR_TAG_BASE, +) import grudge.op as op @@ -54,6 +60,8 @@ make_conserved ) +from mirgecom.utils import normalize_boundaries + # low level routine works with numpy arrays and can be tested without # a full grid + fluid state, etc @@ -335,8 +343,9 @@ def viscous_facial_flux_central(dcoll, state_pair, grad_cv_pair, grad_t_pair, def viscous_flux_on_element_boundary( dcoll, gas_model, boundaries, interior_state_pairs, domain_boundary_states, grad_cv, interior_grad_cv_pairs, - grad_t, interior_grad_t_pairs, quadrature_tag=None, - numerical_flux_func=viscous_facial_flux_central, time=0.0): + grad_t, interior_grad_t_pairs, quadrature_tag=DISCR_TAG_BASE, + numerical_flux_func=viscous_facial_flux_central, time=0.0, + dd=DD_VOLUME_ALL): """Compute the viscous boundary fluxes for the divergence operator. This routine encapsulates the computation of the viscous contributions @@ -351,14 +360,15 @@ def viscous_flux_on_element_boundary( The physical model constructs for the gas model boundaries - Dictionary of boundary functions, one for each valid btag + Dictionary of boundary functions, one for each valid + :class:`~grudge.dof_desc.BoundaryDomainTag` interior_state_pairs Trace pairs of :class:`~mirgecom.gas_model.FluidState` for the interior faces domain_boundary_states A dictionary of boundary-restricted :class:`~mirgecom.gas_model.FluidState`, - keyed by btags in *boundaries*. + keyed by boundary domain tags in *boundaries*. grad_cv: :class:`~mirgecom.fluid.ConservedVars` The gradient of the fluid conserved quantities. @@ -374,41 +384,50 @@ def viscous_flux_on_element_boundary( Trace pairs for the temperature gradient on interior faces quadrature_tag - An optional identifier denoting a particular quadrature - discretization to use during operator evaluations. - The default value is *None*. + An identifier denoting a particular quadrature discretization to use during + operator evaluations. numerical_flux_func The numerical flux function to use in computing the boundary flux. time: float Time + + dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the discretization on which the fluid lives. Must be + a volume on the base discretization. """ - from grudge.dof_desc import as_dofdesc + boundaries = normalize_boundaries(boundaries) - dd_base = as_dofdesc("vol") + if not isinstance(dd.domain_tag, VolumeDomainTag): + raise TypeError("dd must represent a volume") + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + + dd_vol = dd + dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) # {{{ - Viscous flux helpers - # viscous fluxes across interior faces (including partition and periodic bnd) def _fvisc_divergence_flux_interior(state_pair, grad_cv_pair, grad_t_pair): return op.project(dcoll, - state_pair.dd, state_pair.dd.with_dtag("all_faces"), + state_pair.dd, dd_allfaces_quad, numerical_flux_func( dcoll=dcoll, gas_model=gas_model, state_pair=state_pair, grad_cv_pair=grad_cv_pair, grad_t_pair=grad_t_pair)) # viscous part of bcs applied here - def _fvisc_divergence_flux_boundary(dd_btag, boundary, state_minus): - # Make sure we fields on the quadrature grid - # restricted to the tag *btag* + def _fvisc_divergence_flux_boundary(bdtag, boundary, state_minus_quad): + dd_bdry_quad = dd_vol_quad.with_domain_tag(bdtag) return op.project( - dcoll, dd_btag, dd_btag.with_dtag("all_faces"), + dcoll, dd_bdry_quad, dd_allfaces_quad, boundary.viscous_divergence_flux( - dcoll=dcoll, btag=dd_btag, gas_model=gas_model, - state_minus=state_minus, - grad_cv_minus=op.project(dcoll, dd_base, dd_btag, grad_cv), - grad_t_minus=op.project(dcoll, dd_base, dd_btag, grad_t), + dcoll=dcoll, dd_bdry=dd_bdry_quad, gas_model=gas_model, + state_minus=state_minus_quad, + grad_cv_minus=op.project(dcoll, dd_vol, dd_bdry_quad, grad_cv), + grad_t_minus=op.project(dcoll, dd_vol, dd_bdry_quad, grad_t), time=time, numerical_flux_func=numerical_flux_func)) # }}} viscous flux helpers @@ -420,9 +439,10 @@ def _fvisc_divergence_flux_boundary(dd_btag, boundary, state_minus): ( # Domain boundary contributions for the viscous terms sum(_fvisc_divergence_flux_boundary( - as_dofdesc(btag).with_discr_tag(quadrature_tag), - boundary, domain_boundary_states[btag]) - for btag, boundary in boundaries.items()) + bdtag, + boundary, + domain_boundary_states[bdtag]) + for bdtag, boundary in boundaries.items()) # Interior interface contributions for the viscous terms + sum( @@ -436,7 +456,7 @@ def _fvisc_divergence_flux_boundary(dd_btag, boundary, state_minus): return bnd_term -def get_viscous_timestep(dcoll, state): +def get_viscous_timestep(dcoll, state, *, dd=DD_VOLUME_ALL): """Routine returns the the node-local maximum stable viscous timestep. Parameters @@ -449,15 +469,26 @@ def get_viscous_timestep(dcoll, state): Full fluid conserved and thermal state + dd: grudge.dof_desc.DOFDesc + + the DOF descriptor of the discretization on which *state* lives. Must be + a volume on the base discretization. + Returns ------- :class:`~meshmode.dof_array.DOFArray` The maximum stable timestep at each node. """ + if not isinstance(dd.domain_tag, VolumeDomainTag): + raise TypeError("dd must represent a volume") + if dd.discretization_tag != DISCR_TAG_BASE: + raise ValueError("dd must belong to the base discretization") + from grudge.dt_utils import characteristic_lengthscales - length_scales = characteristic_lengthscales(state.array_context, dcoll) + length_scales = characteristic_lengthscales( + state.array_context, dcoll, dd=dd) nu = 0 d_alpha_max = 0 @@ -475,7 +506,7 @@ def get_viscous_timestep(dcoll, state): ) -def get_viscous_cfl(dcoll, dt, state): +def get_viscous_cfl(dcoll, dt, state, *, dd=DD_VOLUME_ALL): """Calculate and return node-local CFL based on current state and timestep. Parameters @@ -492,13 +523,18 @@ def get_viscous_cfl(dcoll, dt, state): The full fluid conserved and thermal state + dd: grudge.dof_desc.DOFDesc + + the DOF descriptor of the discretization on which *state* lives. Must be + a volume on the base discretization. + Returns ------- :class:`~meshmode.dof_array.DOFArray` The CFL at each node. """ - return dt / get_viscous_timestep(dcoll, state=state) + return dt / get_viscous_timestep(dcoll, state=state, dd=dd) def get_local_max_species_diffusivity(actx, d_alpha): diff --git a/mirgecom/wave.py b/mirgecom/wave.py index fe0470b2a..a41d0cb29 100644 --- a/mirgecom/wave.py +++ b/mirgecom/wave.py @@ -90,8 +90,8 @@ def wave_operator(dcoll, c, w, *, comm_tag=None): return ( op.inverse_mass(dcoll, flat_obj_array( - -c*op.weak_local_div(dcoll, "vol", v), - -c*op.weak_local_grad(dcoll, "vol", u) + -c*op.weak_local_div(dcoll, v), + -c*op.weak_local_grad(dcoll, u) ) + # noqa: W504 op.face_mass(dcoll, diff --git a/requirements.txt b/requirements.txt index 4dc2b033d..dbf618ba8 100644 --- a/requirements.txt +++ b/requirements.txt @@ -18,7 +18,7 @@ pyyaml --editable git+https://github.com/inducer/modepy.git#egg=modepy --editable git+https://github.com/inducer/arraycontext.git#egg=arraycontext --editable git+https://github.com/inducer/meshmode.git#egg=meshmode ---editable git+https://github.com/inducer/grudge.git#egg=grudge +--editable git+https://github.com/majosm/grudge.git@multi-main#egg=grudge --editable git+https://github.com/inducer/pytato.git#egg=pytato --editable git+https://github.com/ecisneros8/pyrometheus.git#egg=pyrometheus --editable git+https://github.com/illinois-ceesd/logpyle.git#egg=logpyle diff --git a/test/test_av.py b/test/test_av.py index ae1eea928..20451eeb2 100644 --- a/test/test_av.py +++ b/test/test_av.py @@ -35,6 +35,7 @@ as pytest_generate_tests ) from meshmode.mesh import BTAG_ALL +from meshmode.discretization.connection import FACE_RESTR_ALL import grudge.op as op from mirgecom.artificial_viscosity import ( av_laplacian_operator, @@ -205,29 +206,28 @@ def test_artificial_viscosity(ctx_factory, dim, order): nodes = actx.thaw(dcoll.nodes()) class TestBoundary: - - def cv_gradient_flux(self, dcoll, btag, state_minus, gas_model, **kwargs): + def cv_gradient_flux(self, dcoll, dd_bdry, state_minus, gas_model, **kwargs): cv_int = state_minus.cv from grudge.trace_pair import TracePair - bnd_pair = TracePair(btag, + bnd_pair = TracePair(dd_bdry, interior=cv_int, exterior=cv_int) - nhat = actx.thaw(dcoll.normal(btag)) + nhat = actx.thaw(dcoll.normal(dd_bdry)) 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, dcoll, btag, diffusion, **kwargs): - nhat = actx.thaw(dcoll.normal(btag)) - diffusion_minus = op.project(dcoll, "vol", btag, diffusion) + def av_flux(self, dcoll, dd_bdry, diffusion, **kwargs): + nhat = actx.thaw(dcoll.normal(dd_bdry)) + diffusion_minus = op.project(dcoll, "vol", dd_bdry, diffusion) diffusion_plus = diffusion_minus from grudge.trace_pair import TracePair - bnd_grad_pair = TracePair(btag, interior=diffusion_minus, + bnd_grad_pair = TracePair(dd_bdry, 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 op.project(dcoll, btag, "all_faces", flux_weak) + return op.project(dcoll, dd_bdry, "all_faces", flux_weak) boundaries = {BTAG_ALL: TestBoundary()} @@ -402,9 +402,9 @@ def test_fluid_av_boundaries(ctx_factory, prescribed_soln, order): gas_model = GasModel(eos=IdealSingleGas(gas_const=1.0), transport=transport_model) - def _boundary_state_func(dcoll, btag, gas_model, state_minus, **kwargs): + def _boundary_state_func(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(prescribed_soln(r=nodes, eos=gas_model.eos, **kwargs), gas_model) @@ -444,11 +444,11 @@ def _boundary_state_func(dcoll, btag, gas_model, state_minus, **kwargs): # 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") + dd_bdry = as_dofdesc(BTAG_ALL) + dd_allfaces = dd_bdry.with_boundary_tag(FACE_RESTR_ALL) expected_av_flux_prescribed_boundary = av_diffusion_boundary@boundary_nhat print(f"{expected_av_flux_prescribed_boundary=}") - exp_av_flux = op.project(dcoll, dd_bnd, all_faces_dd, + exp_av_flux = op.project(dcoll, dd_bdry, dd_allfaces, expected_av_flux_prescribed_boundary) print(f"{exp_av_flux=}") diff --git a/test/test_bc.py b/test/test_bc.py index a8b02104c..b4d9db488 100644 --- a/test/test_bc.py +++ b/test/test_bc.py @@ -30,6 +30,7 @@ import pytest from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from meshmode.discretization.connection import FACE_RESTR_ALL from mirgecom.initializers import Lump from mirgecom.boundary import AdiabaticSlipBoundary from mirgecom.eos import IdealSingleGas @@ -103,11 +104,14 @@ def bnd_norm(vec): return actx.to_numpy(op.norm(dcoll, vec, p=np.inf, dd=BTAG_ALL)) state_plus = \ - wall.adiabatic_slip_state(dcoll, btag=BTAG_ALL, gas_model=gas_model, - state_minus=state_minus) + wall.adiabatic_slip_state( + dcoll, dd_bdry=BTAG_ALL, gas_model=gas_model, + state_minus=state_minus) - bnd_pair = TracePair(BTAG_ALL, interior=state_minus.cv, - exterior=state_plus.cv) + bnd_pair = TracePair( + as_dofdesc(BTAG_ALL), + interior=state_minus.cv, + exterior=state_plus.cv) # check that mass and energy are preserved mass_resid = bnd_pair.int.mass - bnd_pair.ext.mass @@ -178,14 +182,18 @@ def bnd_norm(vec): state=fluid_state, gas_model=gas_model) - bnd_soln = wall.adiabatic_slip_state(dcoll, btag=BTAG_ALL, + bnd_soln = wall.adiabatic_slip_state(dcoll, dd_bdry=BTAG_ALL, gas_model=gas_model, state_minus=interior_soln) - bnd_pair = TracePair(BTAG_ALL, interior=interior_soln.cv, - exterior=bnd_soln.cv) - state_pair = TracePair(BTAG_ALL, interior=interior_soln, - exterior=bnd_soln) + bnd_pair = TracePair( + as_dofdesc(BTAG_ALL), + interior=interior_soln.cv, + exterior=bnd_soln.cv) + state_pair = TracePair( + as_dofdesc(BTAG_ALL), + interior=interior_soln, + exterior=bnd_soln) # Check the total velocity component normal # to each surface. It should be zero. The @@ -309,10 +317,11 @@ def scalar_flux_interior(int_tpair): print(f"{cv_flux_int=}") wall_state = wall.isothermal_noslip_state( - dcoll, btag=BTAG_ALL, gas_model=gas_model, state_minus=state_minus) + dcoll, dd_bdry=BTAG_ALL, gas_model=gas_model, + state_minus=state_minus) print(f"{wall_state=}") - cv_grad_flux_wall = wall.cv_gradient_flux(dcoll, btag=BTAG_ALL, + cv_grad_flux_wall = wall.cv_gradient_flux(dcoll, dd_bdry=BTAG_ALL, gas_model=gas_model, state_minus=state_minus) @@ -330,7 +339,7 @@ def scalar_flux_interior(int_tpair): t_int_tpair = interior_trace_pair(dcoll, temper) t_flux_int = scalar_flux_interior(t_int_tpair) - t_flux_bc = wall.temperature_gradient_flux(dcoll, btag=BTAG_ALL, + t_flux_bc = wall.temperature_gradient_flux(dcoll, dd_bdry=BTAG_ALL, gas_model=gas_model, state_minus=state_minus) @@ -340,14 +349,14 @@ def scalar_flux_interior(int_tpair): t_flux_bnd = t_flux_bc + t_flux_int - i_flux_bc = wall.inviscid_divergence_flux(dcoll, btag=BTAG_ALL, + i_flux_bc = wall.inviscid_divergence_flux(dcoll, dd_bdry=BTAG_ALL, gas_model=gas_model, state_minus=state_minus) nhat = actx.thaw(dcoll.normal(state_pair.dd)) bnd_flux = flux_func(state_pair, gas_model, nhat) dd = state_pair.dd - dd_allfaces = dd.with_dtag("all_faces") + dd_allfaces = dd.with_boundary_tag(FACE_RESTR_ALL) i_flux_int = op.project(dcoll, dd, dd_allfaces, bnd_flux) bc_dd = as_dofdesc(BTAG_ALL) i_flux_bc = op.project(dcoll, bc_dd, dd_allfaces, i_flux_bc) @@ -360,19 +369,19 @@ def scalar_flux_interior(int_tpair): from mirgecom.operators import grad_operator dd_vol = as_dofdesc("vol") - dd_faces = as_dofdesc("all_faces") + dd_allfaces = as_dofdesc("all_faces") grad_cv_minus = \ op.project(dcoll, "vol", BTAG_ALL, - grad_operator(dcoll, dd_vol, dd_faces, + grad_operator(dcoll, dd_vol, dd_allfaces, uniform_state.cv, cv_flux_bnd)) grad_t_minus = op.project(dcoll, "vol", BTAG_ALL, - grad_operator(dcoll, dd_vol, dd_faces, + grad_operator(dcoll, dd_vol, dd_allfaces, temper, t_flux_bnd)) print(f"{grad_cv_minus=}") print(f"{grad_t_minus=}") - v_flux_bc = wall.viscous_divergence_flux(dcoll, btag=BTAG_ALL, + v_flux_bc = wall.viscous_divergence_flux(dcoll, dd_bdry=BTAG_ALL, gas_model=gas_model, state_minus=state_minus, grad_cv_minus=grad_cv_minus, @@ -461,9 +470,9 @@ def test_prescribed(actx_factory, prescribed_soln, flux_func): # (*note): Most people will never change these as they are used internally # to compute a DG gradient of Q and temperature. - def _boundary_state_func(dcoll, btag, gas_model, state_minus, **kwargs): + def _boundary_state_func(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(prescribed_soln(r=nodes, eos=gas_model.eos, **kwargs), gas_model) @@ -521,7 +530,7 @@ def scalar_flux_interior(int_tpair): cv_int_tpair = interior_trace_pair(dcoll, cv) cv_flux_int = scalar_flux_interior(cv_int_tpair) - cv_flux_bc = domain_boundary.cv_gradient_flux(dcoll, btag=BTAG_ALL, + cv_flux_bc = domain_boundary.cv_gradient_flux(dcoll, dd_bdry=BTAG_ALL, gas_model=gas_model, state_minus=state_minus) @@ -534,7 +543,7 @@ def scalar_flux_interior(int_tpair): t_int_tpair = interior_trace_pair(dcoll, temper) t_flux_int = scalar_flux_interior(t_int_tpair) t_flux_bc = \ - domain_boundary.temperature_gradient_flux(dcoll, btag=BTAG_ALL, + domain_boundary.temperature_gradient_flux(dcoll, dd_bdry=BTAG_ALL, gas_model=gas_model, state_minus=state_minus) @@ -545,7 +554,7 @@ def scalar_flux_interior(int_tpair): t_flux_bnd = t_flux_bc + t_flux_int i_flux_bc = \ - domain_boundary.inviscid_divergence_flux(dcoll, btag=BTAG_ALL, + domain_boundary.inviscid_divergence_flux(dcoll, dd_bdry=BTAG_ALL, gas_model=gas_model, state_minus=state_minus) @@ -556,9 +565,9 @@ def scalar_flux_interior(int_tpair): nhat = actx.thaw(dcoll.normal(state_pair.dd)) bnd_flux = flux_func(state_pair, gas_model, nhat) dd = state_pair.dd - dd_allfaces = dd.with_dtag("all_faces") - bc_dd = as_dofdesc(BTAG_ALL) - i_flux_bc = op.project(dcoll, bc_dd, dd_allfaces, i_flux_bc) + dd_allfaces = dd.with_boundary_tag(FACE_RESTR_ALL) + dd_bdry = as_dofdesc(BTAG_ALL) + i_flux_bc = op.project(dcoll, dd_bdry, dd_allfaces, i_flux_bc) i_flux_int = op.project(dcoll, dd, dd_allfaces, bnd_flux) i_flux_bnd = i_flux_bc + i_flux_int @@ -569,9 +578,9 @@ def scalar_flux_interior(int_tpair): from mirgecom.operators import grad_operator dd_vol = as_dofdesc("vol") - dd_faces = as_dofdesc("all_faces") - grad_cv = grad_operator(dcoll, dd_vol, dd_faces, cv, cv_flux_bnd) - grad_t = grad_operator(dcoll, dd_vol, dd_faces, temper, t_flux_bnd) + dd_allfaces = as_dofdesc("all_faces") + grad_cv = grad_operator(dcoll, dd_vol, dd_allfaces, cv, cv_flux_bnd) + grad_t = grad_operator(dcoll, dd_vol, dd_allfaces, temper, t_flux_bnd) grad_cv_minus = op.project(dcoll, "vol", BTAG_ALL, grad_cv) grad_t_minus = op.project(dcoll, "vol", BTAG_ALL, grad_t) @@ -579,11 +588,10 @@ def scalar_flux_interior(int_tpair): print(f"{grad_t_minus=}") v_flux_bc = \ - domain_boundary.viscous_divergence_flux(dcoll=dcoll, btag=BTAG_ALL, - gas_model=gas_model, - state_minus=state_minus, - grad_cv_minus=grad_cv_minus, - grad_t_minus=grad_t_minus) + domain_boundary.viscous_divergence_flux( + dcoll=dcoll, dd_bdry=BTAG_ALL, gas_model=gas_model, + state_minus=state_minus, grad_cv_minus=grad_cv_minus, + grad_t_minus=grad_t_minus) print(f"{v_flux_bc=}") bc_soln = \ domain_boundary._boundary_state_pair(dcoll, BTAG_ALL, gas_model, diff --git a/test/test_diffusion.py b/test/test_diffusion.py index a210b415e..862f835eb 100644 --- a/test/test_diffusion.py +++ b/test/test_diffusion.py @@ -38,7 +38,7 @@ DirichletDiffusionBoundary, NeumannDiffusionBoundary) from meshmode.dof_array import DOFArray -from grudge.dof_desc import DTAG_BOUNDARY, DISCR_TAG_BASE, DISCR_TAG_QUAD +from grudge.dof_desc import BoundaryDomainTag, DISCR_TAG_BASE, DISCR_TAG_QUAD from mirgecom.discretization import create_discretization_collection from meshmode.array_context import ( # noqa pytest_generate_tests_for_pyopencl_array_context @@ -133,14 +133,14 @@ def get_boundaries(self, dcoll, actx, t): boundaries = {} for i in range(self.dim-1): - lower_btag = DTAG_BOUNDARY("-"+str(i)) - upper_btag = DTAG_BOUNDARY("+"+str(i)) - boundaries[lower_btag] = NeumannDiffusionBoundary(0.) - boundaries[upper_btag] = NeumannDiffusionBoundary(0.) - lower_btag = DTAG_BOUNDARY("-"+str(self.dim-1)) - upper_btag = DTAG_BOUNDARY("+"+str(self.dim-1)) - boundaries[lower_btag] = DirichletDiffusionBoundary(0.) - boundaries[upper_btag] = DirichletDiffusionBoundary(0.) + lower_bdtag = BoundaryDomainTag("-"+str(i)) + upper_bdtag = BoundaryDomainTag("+"+str(i)) + boundaries[lower_bdtag] = NeumannDiffusionBoundary(0.) + boundaries[upper_bdtag] = NeumannDiffusionBoundary(0.) + lower_bdtag = BoundaryDomainTag("-"+str(self.dim-1)) + upper_bdtag = BoundaryDomainTag("+"+str(self.dim-1)) + boundaries[lower_bdtag] = DirichletDiffusionBoundary(0.) + boundaries[upper_bdtag] = DirichletDiffusionBoundary(0.) return boundaries @@ -179,18 +179,18 @@ def get_boundaries(self, dcoll, actx, t): boundaries = {} for i in range(self.dim-1): - lower_btag = DTAG_BOUNDARY("-"+str(i)) - upper_btag = DTAG_BOUNDARY("+"+str(i)) - upper_grad_u = op.project(dcoll, "vol", upper_btag, exact_grad_u) - normal = actx.thaw(dcoll.normal(upper_btag)) + lower_bdtag = BoundaryDomainTag("-"+str(i)) + upper_bdtag = BoundaryDomainTag("+"+str(i)) + upper_grad_u = op.project(dcoll, "vol", upper_bdtag, exact_grad_u) + normal = actx.thaw(dcoll.normal(upper_bdtag)) upper_grad_u_dot_n = np.dot(upper_grad_u, normal) - boundaries[lower_btag] = NeumannDiffusionBoundary(0.) - boundaries[upper_btag] = NeumannDiffusionBoundary(upper_grad_u_dot_n) - lower_btag = DTAG_BOUNDARY("-"+str(self.dim-1)) - upper_btag = DTAG_BOUNDARY("+"+str(self.dim-1)) - upper_u = op.project(dcoll, "vol", upper_btag, exact_u) - boundaries[lower_btag] = DirichletDiffusionBoundary(0.) - boundaries[upper_btag] = DirichletDiffusionBoundary(upper_u) + boundaries[lower_bdtag] = NeumannDiffusionBoundary(0.) + boundaries[upper_bdtag] = NeumannDiffusionBoundary(upper_grad_u_dot_n) + lower_bdtag = BoundaryDomainTag("-"+str(self.dim-1)) + upper_bdtag = BoundaryDomainTag("+"+str(self.dim-1)) + upper_u = op.project(dcoll, "vol", upper_bdtag, exact_u) + boundaries[lower_bdtag] = DirichletDiffusionBoundary(0.) + boundaries[upper_bdtag] = DirichletDiffusionBoundary(upper_u) return boundaries @@ -227,14 +227,14 @@ def get_boundaries(self, dcoll, actx, t): boundaries = {} for i in range(self.dim-1): - lower_btag = DTAG_BOUNDARY("-"+str(i)) - upper_btag = DTAG_BOUNDARY("+"+str(i)) - boundaries[lower_btag] = NeumannDiffusionBoundary(0.) - boundaries[upper_btag] = NeumannDiffusionBoundary(0.) - lower_btag = DTAG_BOUNDARY("-"+str(self.dim-1)) - upper_btag = DTAG_BOUNDARY("+"+str(self.dim-1)) - boundaries[lower_btag] = DirichletDiffusionBoundary(0.) - boundaries[upper_btag] = DirichletDiffusionBoundary(0.) + lower_bdtag = BoundaryDomainTag("-"+str(i)) + upper_bdtag = BoundaryDomainTag("+"+str(i)) + boundaries[lower_bdtag] = NeumannDiffusionBoundary(0.) + boundaries[upper_bdtag] = NeumannDiffusionBoundary(0.) + lower_bdtag = BoundaryDomainTag("-"+str(self.dim-1)) + upper_bdtag = BoundaryDomainTag("+"+str(self.dim-1)) + boundaries[lower_bdtag] = DirichletDiffusionBoundary(0.) + boundaries[upper_bdtag] = DirichletDiffusionBoundary(0.) return boundaries @@ -265,14 +265,14 @@ def get_boundaries(self, dcoll, actx, t): boundaries = {} for i in range(self.dim-1): - lower_btag = DTAG_BOUNDARY("-"+str(i)) - upper_btag = DTAG_BOUNDARY("+"+str(i)) - boundaries[lower_btag] = NeumannDiffusionBoundary(0.) - boundaries[upper_btag] = NeumannDiffusionBoundary(0.) - lower_btag = DTAG_BOUNDARY("-"+str(self.dim-1)) - upper_btag = DTAG_BOUNDARY("+"+str(self.dim-1)) - boundaries[lower_btag] = DirichletDiffusionBoundary(0.) - boundaries[upper_btag] = DirichletDiffusionBoundary(0.) + lower_bdtag = BoundaryDomainTag("-"+str(i)) + upper_bdtag = BoundaryDomainTag("+"+str(i)) + boundaries[lower_bdtag] = NeumannDiffusionBoundary(0.) + boundaries[upper_bdtag] = NeumannDiffusionBoundary(0.) + lower_bdtag = BoundaryDomainTag("-"+str(self.dim-1)) + upper_bdtag = BoundaryDomainTag("+"+str(self.dim-1)) + boundaries[lower_bdtag] = DirichletDiffusionBoundary(0.) + boundaries[upper_bdtag] = DirichletDiffusionBoundary(0.) return boundaries @@ -411,8 +411,8 @@ def test_diffusion_discontinuous_kappa(actx_factory, order, visualize=False): kappa = kappa_lower * lower_mask + kappa_upper * upper_mask boundaries = { - DTAG_BOUNDARY("-0"): DirichletDiffusionBoundary(0.), - DTAG_BOUNDARY("+0"): DirichletDiffusionBoundary(1.), + BoundaryDomainTag("-0"): DirichletDiffusionBoundary(0.), + BoundaryDomainTag("+0"): DirichletDiffusionBoundary(1.), } flux = -kappa_lower*kappa_upper/(kappa_lower + kappa_upper) diff --git a/test/test_euler.py b/test/test_euler.py index 3498e3514..cefbad856 100644 --- a/test/test_euler.py +++ b/test/test_euler.py @@ -270,9 +270,9 @@ def test_vortex_rhs(actx_factory, order, use_overintegration, numerical_flux_fun gas_model = GasModel(eos=IdealSingleGas()) fluid_state = make_fluid_state(vortex_soln, gas_model) - def _vortex_boundary(dcoll, btag, gas_model, state_minus, **kwargs): + def _vortex_boundary(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(vortex(x_vec=nodes, **kwargs), gas_model) @@ -350,9 +350,9 @@ def test_lump_rhs(actx_factory, dim, order, use_overintegration, gas_model = GasModel(eos=IdealSingleGas()) fluid_state = make_fluid_state(lump_soln, gas_model) - def _lump_boundary(dcoll, btag, gas_model, state_minus, **kwargs): + def _lump_boundary(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(lump(x_vec=nodes, cv=state_minus, **kwargs), gas_model) @@ -445,9 +445,9 @@ def test_multilump_rhs(actx_factory, dim, order, v0, use_overintegration, gas_model = GasModel(eos=IdealSingleGas()) fluid_state = make_fluid_state(lump_soln, gas_model) - def _my_boundary(dcoll, btag, gas_model, state_minus, **kwargs): + def _my_boundary(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(lump(x_vec=nodes, **kwargs), gas_model) @@ -603,7 +603,7 @@ def rhs(t, q): write_soln(state=cv) cv = rk4_step(cv, t, dt, rhs) - cv = filter_modally(dcoll, "vol", cutoff, frfunc, cv) + cv = filter_modally(dcoll, cutoff, frfunc, cv) fluid_state = make_fluid_state(cv, gas_model) t += dt @@ -664,9 +664,9 @@ def test_isentropic_vortex(actx_factory, order, use_overintegration, initializer = Vortex2D(center=orig, velocity=vel) casename = "Vortex" - def _vortex_boundary(dcoll, btag, state_minus, gas_model, **kwargs): + def _vortex_boundary(dcoll, dd_bdry, state_minus, gas_model, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(initializer(x_vec=nodes, **kwargs), gas_model) diff --git a/test/test_filter.py b/test/test_filter.py index 661c664dd..6e3484b65 100644 --- a/test/test_filter.py +++ b/test/test_filter.py @@ -188,8 +188,7 @@ def test_filter_function(actx_factory, dim, order, do_viz=False): uniform_soln = initr(t=0, x_vec=nodes) from mirgecom.filter import filter_modally - filtered_soln = filter_modally(dcoll, "vol", cutoff, - frfunc, uniform_soln) + filtered_soln = filter_modally(dcoll, cutoff, frfunc, uniform_soln) soln_resid = uniform_soln - filtered_soln from mirgecom.simutil import componentwise_norms max_errors = componentwise_norms(dcoll, soln_resid, np.inf) @@ -206,7 +205,7 @@ def polyfn(coeff): # , x_vec): r = nodes[0] result = 0 for n, a in enumerate(coeff): - result += a * r ** n + result = result + a * r ** n return make_obj_array([result]) # Any order {cutoff} and below fields should be unharmed @@ -214,8 +213,7 @@ def polyfn(coeff): # , x_vec): field_order = int(cutoff) coeff = [1.0 / (i + 1) for i in range(field_order + 1)] field = polyfn(coeff=coeff) - filtered_field = filter_modally(dcoll, "vol", cutoff, - frfunc, field) + filtered_field = filter_modally(dcoll, cutoff, frfunc, field) soln_resid = field - filtered_field max_errors = [actx.to_numpy(op.norm(dcoll, v, np.inf)) for v in soln_resid] logger.info(f"Field = {field}") @@ -230,15 +228,14 @@ def polyfn(coeff): # , x_vec): from grudge.shortcuts import make_visualizer vis = make_visualizer(dcoll, order) - from grudge.dof_desc import DD_VOLUME_MODAL, DD_VOLUME + from grudge.dof_desc import DD_VOLUME_ALL, DD_VOLUME_ALL_MODAL - modal_map = dcoll.connection_from_dds(DD_VOLUME, DD_VOLUME_MODAL) + modal_map = dcoll.connection_from_dds(DD_VOLUME_ALL, DD_VOLUME_ALL_MODAL) for field_order in range(cutoff+1, cutoff+4): coeff = [1.0 / (i + 1) for i in range(field_order+1)] field = polyfn(coeff=coeff) - filtered_field = filter_modally(dcoll, "vol", cutoff, - frfunc, field) + filtered_field = filter_modally(dcoll, cutoff, frfunc, field) unfiltered_spectrum = modal_map(field) filtered_spectrum = modal_map(filtered_field) diff --git a/test/test_flux.py b/test/test_flux.py index 12bc6f2b5..a72dec470 100644 --- a/test/test_flux.py +++ b/test/test_flux.py @@ -34,6 +34,7 @@ from pytools.obj_array import make_obj_array from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from meshmode.dof_array import DOFArray +from grudge.dof_desc import as_dofdesc from grudge.trace_pair import TracePair from mirgecom.fluid import make_conserved from mirgecom.eos import IdealSingleGas @@ -129,7 +130,10 @@ def test_lfr_flux(actx_factory, nspecies, dim, norm_dir, vel_mag): mag = np.linalg.norm(normal) normal = norm_dir*normal/mag - state_pair = TracePair("vol", interior=fluid_state_int, exterior=fluid_state_ext) + state_pair = TracePair( + as_dofdesc("vol"), + interior=fluid_state_int, + exterior=fluid_state_ext) # code passes in fluxes in the direction of the surface normal, # so we will too @@ -275,7 +279,10 @@ def test_hll_flux(actx_factory, nspecies, dim, norm_dir, vel_mag): mag = np.linalg.norm(normal) normal = norm_dir*normal/mag - state_pair = TracePair("vol", interior=fluid_state_int, exterior=fluid_state_ext) + state_pair = TracePair( + as_dofdesc("vol"), + interior=fluid_state_int, + exterior=fluid_state_ext) from mirgecom.inviscid import inviscid_facial_flux_hll flux_bnd = inviscid_facial_flux_hll(state_pair, gas_model, normal) diff --git a/test/test_inviscid.py b/test/test_inviscid.py index 83dca1cb2..c78b6d79b 100644 --- a/test/test_inviscid.py +++ b/test/test_inviscid.py @@ -37,6 +37,8 @@ ) from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from meshmode.discretization.connection import FACE_RESTR_ALL +from grudge.dof_desc import as_dofdesc from grudge.trace_pair import TracePair from mirgecom.fluid import make_conserved from mirgecom.eos import IdealSingleGas @@ -333,7 +335,7 @@ def test_facial_flux(actx_factory, nspecies, order, dim, num_flux): nhat = actx.thaw(dcoll.normal(interior_state_pair.dd)) bnd_flux = num_flux(interior_state_pair, gas_model, nhat) dd = interior_state_pair.dd - dd_allfaces = dd.with_dtag("all_faces") + dd_allfaces = dd.with_boundary_tag(FACE_RESTR_ALL) interior_face_flux = op.project(dcoll, dd, dd_allfaces, bnd_flux) def inf_norm(data): @@ -375,14 +377,14 @@ def inf_norm(data): momentum=dir_mom, species_mass=dir_mf) dir_bval = make_conserved(dim, mass=dir_mass, energy=dir_e, momentum=dir_mom, species_mass=dir_mf) - state_tpair = TracePair(BTAG_ALL, + state_tpair = TracePair(as_dofdesc(BTAG_ALL), interior=make_fluid_state(dir_bval, gas_model), exterior=make_fluid_state(dir_bc, gas_model)) nhat = actx.thaw(dcoll.normal(state_tpair.dd)) bnd_flux = num_flux(state_tpair, gas_model, nhat) dd = state_tpair.dd - dd_allfaces = dd.with_dtag("all_faces") + dd_allfaces = dd.with_boundary_tag(FACE_RESTR_ALL) boundary_flux = op.project(dcoll, dd, dd_allfaces, bnd_flux) assert inf_norm(boundary_flux.mass) < tolerance diff --git a/test/test_lazy.py b/test/test_lazy.py index 044376f98..413a7fe6c 100644 --- a/test/test_lazy.py +++ b/test/test_lazy.py @@ -33,6 +33,7 @@ PyOpenCLArrayContext, PytatoPyOpenCLArrayContext ) +from meshmode.discretization.connection import FACE_RESTR_ALL from mirgecom.discretization import create_discretization_collection import grudge.op as op from meshmode.array_context import ( # noqa @@ -102,15 +103,15 @@ def _isclose(dcoll, x, y, rel_tol=1e-9, abs_tol=0, return_operands=False): # cl_ctx = ctx_factory() # actx, dcoll = _op_test_fixture(cl_ctx) # -# from grudge.dof_desc import DTAG_BOUNDARY +# from grudge.dof_desc import BoundaryDomainTag # from mirgecom.diffusion import ( # _gradient_operator, # DirichletDiffusionBoundary, # NeumannDiffusionBoundary) # # boundaries = { -# DTAG_BOUNDARY("x"): DirichletDiffusionBoundary(0), -# DTAG_BOUNDARY("y"): NeumannDiffusionBoundary(0) +# BoundaryDomainTag("x"): DirichletDiffusionBoundary(0), +# BoundaryDomainTag("y"): NeumannDiffusionBoundary(0) # } # # def op(u): @@ -132,11 +133,11 @@ def test_lazy_op_divergence(op_test_data, order): from mirgecom.operators import div_operator dd_vol = as_dofdesc("vol") - dd_faces = as_dofdesc("all_faces") + dd_allfaces = as_dofdesc("all_faces") def get_flux(u_tpair): dd = u_tpair.dd - dd_allfaces = dd.with_dtag("all_faces") + dd_allfaces = dd.with_boundary_tag(FACE_RESTR_ALL) normal = dcoll.normal(dd) actx = u_tpair.int[0].array_context normal = actx.thaw(normal) @@ -144,7 +145,7 @@ def get_flux(u_tpair): return op.project(dcoll, dd, dd_allfaces, flux) def div_op(u): - return div_operator(dcoll, dd_vol, dd_faces, + return div_operator(dcoll, dd_vol, dd_allfaces, u, get_flux(interior_trace_pair(dcoll, u))) lazy_op = lazy_actx.compile(div_op) @@ -173,15 +174,15 @@ def test_lazy_op_diffusion(op_test_data, order): eager_actx, lazy_actx, get_dcoll = op_test_data dcoll = get_dcoll(order) - from grudge.dof_desc import DTAG_BOUNDARY + from grudge.dof_desc import BoundaryDomainTag from mirgecom.diffusion import ( diffusion_operator, DirichletDiffusionBoundary, NeumannDiffusionBoundary) boundaries = { - DTAG_BOUNDARY("x"): DirichletDiffusionBoundary(0), - DTAG_BOUNDARY("y"): NeumannDiffusionBoundary(0) + BoundaryDomainTag("x"): DirichletDiffusionBoundary(0), + BoundaryDomainTag("y"): NeumannDiffusionBoundary(0) } def diffusion_op(kappa, u): @@ -239,9 +240,9 @@ def _get_scalar_lump(): dim=2, nspecies=3, velocity=np.ones(2), spec_y0s=np.ones(3), spec_amplitudes=np.ones(3)) - def _my_boundary(dcoll, btag, gas_model, state_minus, **kwargs): + def _my_boundary(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(init(x_vec=nodes, eos=gas_model.eos, **kwargs), gas_model) diff --git a/test/test_multiphysics.py b/test/test_multiphysics.py new file mode 100644 index 000000000..7805b2481 --- /dev/null +++ b/test/test_multiphysics.py @@ -0,0 +1,136 @@ +__copyright__ = """Copyright (C) 2022 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 +import pyopencl.array as cla # noqa +import pyopencl.clmath as clmath # noqa +from pytools.obj_array import make_obj_array +import grudge.op as op +from mirgecom.diffusion import ( + diffusion_operator, + DirichletDiffusionBoundary, + NeumannDiffusionBoundary) +from grudge.dof_desc import DOFDesc, VolumeDomainTag, DISCR_TAG_BASE +from mirgecom.discretization import create_discretization_collection +from meshmode.array_context import ( # noqa + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests) +import pytest + +import logging +logger = logging.getLogger(__name__) + + +@pytest.mark.parametrize("order", [1, 2, 3]) +def test_independent_volumes(actx_factory, order, visualize=False): + """Check multi-volume machinery by setting up two independent volumes.""" + actx = actx_factory() + + n = 8 + + dim = 2 + + dim_names = ["x", "y", "z"] + boundary_tag_to_face = {} + for i in range(dim): + boundary_tag_to_face["-"+str(i)] = ["-"+dim_names[i]] + boundary_tag_to_face["+"+str(i)] = ["+"+dim_names[i]] + + from meshmode.mesh.generation import generate_regular_rect_mesh + + mesh = generate_regular_rect_mesh( + a=(-1,)*dim, b=(1,)*dim, + nelements_per_axis=(n,)*dim, boundary_tag_to_face=boundary_tag_to_face) + + volume_meshes = { + "vol1": mesh, + "vol2": mesh, + } + + dcoll = create_discretization_collection(actx, volume_meshes, order) + + dd_vol1 = DOFDesc(VolumeDomainTag("vol1"), DISCR_TAG_BASE) + dd_vol2 = DOFDesc(VolumeDomainTag("vol2"), DISCR_TAG_BASE) + + nodes1 = actx.thaw(dcoll.nodes(dd=dd_vol1)) + nodes2 = actx.thaw(dcoll.nodes(dd=dd_vol2)) + + # Set solution to x for volume 1 + # Set solution to y for volume 2 + + boundaries1 = { + dd_vol1.trace("-0").domain_tag: DirichletDiffusionBoundary(-1.), + dd_vol1.trace("+0").domain_tag: DirichletDiffusionBoundary(1.), + dd_vol1.trace("-1").domain_tag: NeumannDiffusionBoundary(0.), + dd_vol1.trace("+1").domain_tag: NeumannDiffusionBoundary(0.), + } + + boundaries2 = { + dd_vol2.trace("-0").domain_tag: NeumannDiffusionBoundary(0.), + dd_vol2.trace("+0").domain_tag: NeumannDiffusionBoundary(0.), + dd_vol2.trace("-1").domain_tag: DirichletDiffusionBoundary(-1.), + dd_vol2.trace("+1").domain_tag: DirichletDiffusionBoundary(1.), + } + + u1 = nodes1[0] + u2 = nodes2[1] + + u = make_obj_array([u1, u2]) + + def get_rhs(t, u): + return make_obj_array([ + diffusion_operator( + dcoll, kappa=1, boundaries=boundaries1, u=u[0], dd=dd_vol1), + diffusion_operator( + dcoll, kappa=1, boundaries=boundaries2, u=u[1], dd=dd_vol2)]) + + rhs = get_rhs(0, u) + + if visualize: + from grudge.shortcuts import make_visualizer + viz1 = make_visualizer(dcoll, order+3, volume_dd=dd_vol1) + viz2 = make_visualizer(dcoll, order+3, volume_dd=dd_vol2) + viz1.write_vtk_file( + f"multiphysics_independent_volumes_{order}_1.vtu", [ + ("u", u[0]), + ("rhs", rhs[0]), + ]) + viz2.write_vtk_file( + f"multiphysics_independent_volumes_{order}_2.vtu", [ + ("u", u[1]), + ("rhs", rhs[1]), + ]) + + linf_err1 = actx.to_numpy(op.norm(dcoll, rhs[0], np.inf, dd=dd_vol1)) + linf_err2 = actx.to_numpy(op.norm(dcoll, rhs[1], np.inf, dd=dd_vol2)) + + assert linf_err1 < 1e-9 + assert linf_err2 < 1e-9 + + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from pytest import main + main([__file__]) diff --git a/test/test_navierstokes.py b/test/test_navierstokes.py index 86cb95dfc..d46654ea0 100644 --- a/test/test_navierstokes.py +++ b/test/test_navierstokes.py @@ -294,10 +294,10 @@ def get_boundaries(self): """Get the boundary condition dictionary: prescribed exact by default.""" from mirgecom.gas_model import make_fluid_state - def _boundary_state_func(dcoll, btag, gas_model, state_minus, time=0, + def _boundary_state_func(dcoll, dd_bdry, gas_model, state_minus, time=0, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(self.get_solution(x=nodes, t=time), gas_model) @@ -676,10 +676,10 @@ def test_shear_flow(actx_factory, dim, flow_direction, order): from pytools.convergence import EOCRecorder eoc_energy = EOCRecorder() - def _boundary_state_func(dcoll, btag, gas_model, state_minus, time=0, + def _boundary_state_func(dcoll, dd_bdry, gas_model, state_minus, time=0, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) boundary_cv = exact_soln(x=nodes) return make_fluid_state(boundary_cv, gas_model) @@ -932,10 +932,10 @@ def test_roy_mms(actx_factory, order, dim, u_0, v_0, w_0, a_r, a_p, a_u, logger.info(f"{source_norms=}") logger.info(f"{source_eval=}") - def _boundary_state_func(dcoll, btag, gas_model, state_minus, time=0, + def _boundary_state_func(dcoll, dd_bdry, gas_model, state_minus, time=0, **kwargs): actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(btag) + bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) boundary_cv = evaluate(sym_cv, x=nodes, t=time) return make_fluid_state(boundary_cv, gas_model) diff --git a/test/test_operators.py b/test/test_operators.py index bd28375ff..ab6d6e865 100644 --- a/test/test_operators.py +++ b/test/test_operators.py @@ -35,6 +35,7 @@ import pymbolic as pmbl # noqa import pymbolic.primitives as prim from meshmode.mesh import BTAG_ALL +from grudge.dof_desc import as_dofdesc from mirgecom.flux import num_flux_central from mirgecom.fluid import ( make_conserved @@ -50,8 +51,9 @@ def _elbnd_flux(dcoll, compute_interior_flux, compute_boundary_flux, int_tpair, boundaries): - return (compute_interior_flux(int_tpair) - + sum(compute_boundary_flux(btag) for btag in boundaries)) + return ( + compute_interior_flux(int_tpair) + + sum(compute_boundary_flux(as_dofdesc(bdtag)) for bdtag in boundaries)) # Box grid generator widget lifted from @majosm and slightly bent @@ -124,22 +126,22 @@ def central_flux_interior(actx, dcoll, int_tpair): normal = actx.thaw(dcoll.normal(int_tpair.dd)) from arraycontext import outer flux_weak = outer(num_flux_central(int_tpair.int, int_tpair.ext), normal) - dd_all_faces = int_tpair.dd.with_dtag("all_faces") - return op.project(dcoll, int_tpair.dd, dd_all_faces, flux_weak) + dd_allfaces = int_tpair.dd.with_dtag("all_faces") + return op.project(dcoll, int_tpair.dd, dd_allfaces, flux_weak) -def central_flux_boundary(actx, dcoll, soln_func, btag): +def central_flux_boundary(actx, dcoll, soln_func, dd_bdry): """Compute a central flux for boundary faces.""" - boundary_discr = dcoll.discr_from_dd(btag) + boundary_discr = dcoll.discr_from_dd(dd_bdry) bnd_nodes = actx.thaw(boundary_discr.nodes()) soln_bnd = soln_func(x_vec=bnd_nodes) - bnd_nhat = actx.thaw(dcoll.normal(btag)) + bnd_nhat = actx.thaw(dcoll.normal(dd_bdry)) from grudge.trace_pair import TracePair - bnd_tpair = TracePair(btag, interior=soln_bnd, exterior=soln_bnd) + bnd_tpair = TracePair(dd_bdry, interior=soln_bnd, exterior=soln_bnd) from arraycontext import outer flux_weak = outer(num_flux_central(bnd_tpair.int, bnd_tpair.ext), bnd_nhat) - dd_all_faces = bnd_tpair.dd.with_dtag("all_faces") - return op.project(dcoll, bnd_tpair.dd, dd_all_faces, flux_weak) + dd_allfaces = bnd_tpair.dd.with_dtag("all_faces") + return op.project(dcoll, bnd_tpair.dd, dd_allfaces, flux_weak) @pytest.mark.parametrize("dim", [1, 2, 3]) @@ -225,10 +227,9 @@ def sym_eval(expr, x_vec): test_data_int_tpair, boundaries) from mirgecom.operators import grad_operator - from grudge.dof_desc import as_dofdesc dd_vol = as_dofdesc("vol") - dd_faces = as_dofdesc("all_faces") - test_grad = grad_operator(dcoll, dd_vol, dd_faces, + dd_allfaces = as_dofdesc("all_faces") + test_grad = grad_operator(dcoll, dd_vol, dd_allfaces, test_data, test_data_flux_bnd) print(f"{test_grad=}") diff --git a/test/test_viscous.py b/test/test_viscous.py index c595c9372..1b156b3e9 100644 --- a/test/test_viscous.py +++ b/test/test_viscous.py @@ -32,7 +32,9 @@ import pytest # noqa from pytools.obj_array import make_obj_array +from meshmode.discretization.connection import FACE_RESTR_ALL from meshmode.mesh import BTAG_ALL +from grudge.dof_desc import as_dofdesc import grudge.op as op from grudge.trace_pair import interior_trace_pairs from mirgecom.discretization import create_discretization_collection @@ -161,8 +163,9 @@ def test_poiseuille_fluxes(actx_factory, order, kappa): def _elbnd_flux(dcoll, compute_interior_flux, compute_boundary_flux, int_tpairs, boundaries): - return ((sum(compute_interior_flux(int_tpair) for int_tpair in int_tpairs)) - + sum(compute_boundary_flux(btag) for btag in boundaries)) + return ( + sum(compute_interior_flux(int_tpair) for int_tpair in int_tpairs) + + sum(compute_boundary_flux(as_dofdesc(bdtag)) for bdtag in boundaries)) from mirgecom.flux import num_flux_central @@ -170,20 +173,20 @@ def cv_flux_interior(int_tpair): normal = actx.thaw(dcoll.normal(int_tpair.dd)) from arraycontext import outer flux_weak = outer(num_flux_central(int_tpair.int, int_tpair.ext), normal) - dd_all_faces = int_tpair.dd.with_dtag("all_faces") - return op.project(dcoll, int_tpair.dd, dd_all_faces, flux_weak) + dd_allfaces = int_tpair.dd.with_boundary_tag(FACE_RESTR_ALL) + return op.project(dcoll, int_tpair.dd, dd_allfaces, flux_weak) - def cv_flux_boundary(btag): - boundary_discr = dcoll.discr_from_dd(btag) + def cv_flux_boundary(dd_bdry): + boundary_discr = dcoll.discr_from_dd(dd_bdry) bnd_nodes = actx.thaw(boundary_discr.nodes()) cv_bnd = initializer(x_vec=bnd_nodes, eos=eos) - bnd_nhat = actx.thaw(dcoll.normal(btag)) + bnd_nhat = actx.thaw(dcoll.normal(dd_bdry)) from grudge.trace_pair import TracePair - bnd_tpair = TracePair(btag, interior=cv_bnd, exterior=cv_bnd) + bnd_tpair = TracePair(dd_bdry, interior=cv_bnd, exterior=cv_bnd) from arraycontext import outer flux_weak = outer(num_flux_central(bnd_tpair.int, bnd_tpair.ext), bnd_nhat) - dd_all_faces = bnd_tpair.dd.with_dtag("all_faces") - return op.project(dcoll, bnd_tpair.dd, dd_all_faces, flux_weak) + dd_allfaces = dd_bdry.with_boundary_tag(FACE_RESTR_ALL) + return op.project(dcoll, dd_bdry, dd_allfaces, flux_weak) for nfac in [1, 2, 4]: @@ -213,10 +216,9 @@ def inf_norm(x): cv_flux_bnd = _elbnd_flux(dcoll, cv_flux_interior, cv_flux_boundary, cv_int_tpairs, boundaries) from mirgecom.operators import grad_operator - from grudge.dof_desc import as_dofdesc dd_vol = as_dofdesc("vol") - dd_faces = as_dofdesc("all_faces") - grad_cv = grad_operator(dcoll, dd_vol, dd_faces, cv, cv_flux_bnd) + dd_allfaces = as_dofdesc("all_faces") + grad_cv = grad_operator(dcoll, dd_vol, dd_allfaces, cv, cv_flux_bnd) xp_grad_cv = initializer.exact_grad(x_vec=nodes, eos=eos, cv_exact=cv) xp_grad_v = 1/cv.mass * xp_grad_cv.momentum