From 521388e86f0a98f2e099fa9cc33e7b4b235630ca Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Thu, 11 Aug 2022 12:44:57 -0500 Subject: [PATCH 01/19] add multi-volume support --- examples/doublemach-mpi.py | 22 ++-- examples/heat-source-mpi.py | 6 +- examples/hotplate-mpi.py | 23 ++-- examples/lump-mpi.py | 4 +- examples/mixture-mpi.py | 4 +- examples/poiseuille-mpi.py | 22 ++-- examples/scalar-lump-mpi.py | 4 +- examples/sod-mpi.py | 4 +- examples/vortex-mpi.py | 4 +- mirgecom/artificial_viscosity.py | 78 ++++++++----- mirgecom/boundary.py | 163 ++++++++++++++------------ mirgecom/diffusion.py | 190 ++++++++++++++++++++----------- mirgecom/discretization.py | 13 +-- mirgecom/euler.py | 25 ++-- mirgecom/filter.py | 18 +-- mirgecom/gas_model.py | 39 ++++--- mirgecom/inviscid.py | 33 +++--- mirgecom/io.py | 9 +- mirgecom/limiter.py | 16 ++- mirgecom/logging_quantities.py | 18 +-- mirgecom/navierstokes.py | 133 ++++++++++++++-------- mirgecom/operators.py | 8 +- mirgecom/simutil.py | 21 ++-- mirgecom/viscous.py | 57 ++++++---- mirgecom/wave.py | 4 +- requirements.txt | 8 +- test/test_av.py | 21 ++-- test/test_bc.py | 55 +++++---- test/test_diffusion.py | 76 ++++++------- test/test_euler.py | 16 +-- test/test_filter.py | 4 +- test/test_flux.py | 11 +- test/test_inviscid.py | 3 +- test/test_lazy.py | 16 +-- test/test_multiphysics.py | 142 +++++++++++++++++++++++ test/test_navierstokes.py | 12 +- test/test_operators.py | 15 +-- test/test_viscous.py | 22 ++-- 38 files changed, 827 insertions(+), 492 deletions(-) create mode 100644 test/test_multiphysics.py diff --git a/examples/doublemach-mpi.py b/examples/doublemach-mpi.py index 689852997..243eed84d 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 @@ -237,22 +237,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..a21279817 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 ( @@ -125,8 +125,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..181d922ce 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 @@ -221,21 +221,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..6054f3d1e 100644 --- a/examples/lump-mpi.py +++ b/examples/lump-mpi.py @@ -180,9 +180,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..07d613c04 100644 --- a/examples/mixture-mpi.py +++ b/examples/mixture-mpi.py @@ -202,9 +202,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/poiseuille-mpi.py b/examples/poiseuille-mpi.py index 2d0afbd5e..9d5506341 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 @@ -235,19 +234,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/scalar-lump-mpi.py b/examples/scalar-lump-mpi.py index 19b3f3851..3d3653cb2 100644 --- a/examples/scalar-lump-mpi.py +++ b/examples/scalar-lump-mpi.py @@ -185,9 +185,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..9645e0224 100644 --- a/examples/sod-mpi.py +++ b/examples/sod-mpi.py @@ -173,9 +173,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..77491aef6 100644 --- a/examples/vortex-mpi.py +++ b/examples/vortex-mpi.py @@ -190,9 +190,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/mirgecom/artificial_viscosity.py b/mirgecom/artificial_viscosity.py index a17000e01..e16fe9208 100644 --- a/mirgecom/artificial_viscosity.py +++ b/mirgecom/artificial_viscosity.py @@ -126,10 +126,12 @@ """ import numpy as np +from dataclasses import replace 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,10 +142,10 @@ ) from grudge.dof_desc import ( - DOFDesc, + DD_VOLUME_ALL, + DISCR_TAG_BASE, + DISCR_TAG_MODAL, as_dofdesc, - DD_VOLUME_MODAL, - DD_VOLUME ) import grudge.op as op @@ -158,9 +160,11 @@ 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, + kappa=1., s0=-6., time=0, quadrature_tag=DISCR_TAG_BASE, + volume_dd=DD_VOLUME_ALL, boundary_kwargs=None, indicator=None, divergence_numerical_flux=num_flux_central, + operator_states_quad=None, + grad_cv=None, **kwargs): r"""Compute the artificial viscosity right-hand-side. @@ -200,20 +204,25 @@ 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 + volume_dd: grudge.dof_desc.DOFDesc + The DOF descriptor of the volume on which to apply the operator. Returns ------- :class:`mirgecom.fluid.ConservedVars` The artificial viscosity operator applied to *q*. """ + boundaries = { + as_dofdesc(bdtag).domain_tag: bdry + for bdtag, bdry in boundaries.items()} + cv = fluid_state.cv actx = cv.array_context - dd_vol = DOFDesc("vol", quadrature_tag) - dd_faces = DOFDesc("all_faces", quadrature_tag) + + dd_base = volume_dd + dd_vol = dd_base.with_discr_tag(quadrature_tag) + dd_allfaces = dd_vol.trace(FACE_RESTR_ALL) from warnings import warn @@ -228,12 +237,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_base, dd_vol, u) if operator_states_quad is None: from mirgecom.gas_model import make_operator_fluid_states operator_states_quad = make_operator_fluid_states( - dcoll, fluid_state, gas_model, boundaries, quadrature_tag) + dcoll, fluid_state, gas_model, boundaries, quadrature_tag, + volume_dd=dd_base) vol_state_quad, inter_elem_bnd_states_quad, domain_bnd_states_quad = \ operator_states_quad @@ -241,21 +251,24 @@ 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, volume_dd=dd_base) 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, + volume_dd=dd_base, 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"), + dd_trace = utpair.dd + dd_allfaces = dd_trace.with_domain_tag( + replace(dd_trace.domain_tag, tag=FACE_RESTR_ALL)) + normal = actx.thaw(dcoll.normal(dd_trace)) + return op.project(dcoll, dd_trace, dd_allfaces, # 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) @@ -263,21 +276,26 @@ def central_flux_div(utpair): # 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, tag=_AVRTag)) + sum( + central_flux_div(interp_to_surf_quad(tpair=tpair)) + for tpair in interior_trace_pairs( + dcoll, r, volume_dd=dd_base, tag=_AVRTag)) # 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_vol=dd_vol, + 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, dd_allfaces, 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, volume_dd=DD_VOLUME_ALL): r"""Calculate the smoothness indicator. Parameters @@ -304,7 +322,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", volume_dd)) def indicator_prg(): """Compute the smoothness indicator for all elements.""" from arraycontext import make_loopy_program @@ -331,7 +349,7 @@ def indicator_prg(): "idof": ConcurrentDOFInameTag()}) @keyed_memoize_in(actx, (smoothness_indicator, - "highest_mode"), + "highest_mode", volume_dd), lambda grp: grp.discretization_key()) def highest_mode(grp): return actx.from_numpy( @@ -341,7 +359,9 @@ def highest_mode(grp): ) # Convert to modal solution representation - modal_map = dcoll.connection_from_dds(DD_VOLUME, DD_VOLUME_MODAL) + dd_vol = volume_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 @@ -362,7 +382,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: @@ -374,7 +394,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..3cf4a45a7 100644 --- a/mirgecom/boundary.py +++ b/mirgecom/boundary.py @@ -44,7 +44,10 @@ """ import numpy as np +from dataclasses import replace from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from meshmode.discretization.connection import FACE_RESTR_ALL +from grudge.dof_desc import VolumeDomainTag, as_dofdesc from mirgecom.fluid import make_conserved from grudge.trace_pair import TracePair import grudge.op as op @@ -67,7 +70,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 +87,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 +112,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 +126,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 +165,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 +177,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 +199,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 +213,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 +317,17 @@ 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_domain_tag( + replace(dd_bdry.domain_tag, 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 +339,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 +353,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 +393,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 +435,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 +477,20 @@ 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) + dd_vol = dd_bdry.with_domain_tag( + VolumeDomainTag(dd_bdry.domain_tag.volume_tag)) + grad_av_minus = op.project(dcoll, dd_vol, 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 +529,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 +545,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 +562,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 +602,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 +640,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 b285d0e3d..af9b7ed5f 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -36,9 +36,11 @@ import abc import numpy as np import numpy.linalg as la # noqa +from dataclasses import replace 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, DISCR_TAG_BASE, as_dofdesc from grudge.trace_pair import TracePair, interior_trace_pairs import grudge.op as op @@ -47,19 +49,20 @@ def grad_flux(dcoll, u_tpair, *, quadrature_tag=DISCR_TAG_BASE): r"""Compute the numerical flux for $\nabla u$.""" actx = u_tpair.int.array_context - dd = u_tpair.dd - dd_quad = dd.with_discr_tag(quadrature_tag) - dd_allfaces_quad = dd_quad.with_dtag("all_faces") + dd_trace = u_tpair.dd + dd_trace_quad = dd_trace.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_trace_quad.with_domain_tag( + replace(dd_trace_quad.domain_tag, tag=FACE_RESTR_ALL)) - normal_quad = actx.thaw(dcoll.normal(dd_quad)) + normal_quad = actx.thaw(dcoll.normal(dd_trace_quad)) def to_quad(a): - return op.project(dcoll, dd, dd_quad, a) + return op.project(dcoll, dd_trace, dd_trace_quad, a) def flux(u, normal): return -u * normal - return op.project(dcoll, dd_quad, dd_allfaces_quad, flux( + return op.project(dcoll, dd_trace_quad, dd_allfaces_quad, flux( to_quad(u_tpair.avg), normal_quad)) @@ -68,26 +71,27 @@ def diffusion_flux( r"""Compute the numerical flux for $\nabla \cdot (\kappa \nabla u)$.""" actx = grad_u_tpair.int[0].array_context - dd = grad_u_tpair.dd - dd_quad = dd.with_discr_tag(quadrature_tag) - dd_allfaces_quad = dd_quad.with_dtag("all_faces") + dd_trace = grad_u_tpair.dd + dd_trace_quad = dd_trace.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_trace_quad.with_domain_tag( + replace(dd_trace_quad.domain_tag, tag=FACE_RESTR_ALL)) - normal_quad = actx.thaw(dcoll.normal(dd_quad)) + normal_quad = actx.thaw(dcoll.normal(dd_trace_quad)) def to_quad(a): - return op.project(dcoll, dd, dd_quad, a) + return op.project(dcoll, dd_trace, dd_trace_quad, a) def flux(kappa, grad_u, normal): return -kappa * np.dot(grad_u, normal) - flux_tpair = TracePair(dd_quad, + flux_tpair = TracePair(dd_trace_quad, interior=flux( to_quad(kappa_tpair.int), to_quad(grad_u_tpair.int), normal_quad), exterior=flux( to_quad(kappa_tpair.ext), to_quad(grad_u_tpair.ext), normal_quad) ) - return op.project(dcoll, dd_quad, dd_allfaces_quad, flux_tpair.avg) + return op.project(dcoll, dd_trace_quad, dd_allfaces_quad, flux_tpair.avg) class DiffusionBoundary(metaclass=abc.ABCMeta): @@ -100,14 +104,16 @@ class DiffusionBoundary(metaclass=abc.ABCMeta): @abc.abstractmethod def get_grad_flux( - self, dcoll, dd, u, *, quadrature_tag=DISCR_TAG_BASE): - """Compute the flux for grad(u) on the boundary corresponding to *dd*.""" + self, dcoll, dd_vol, dd_bdry, u, *, + quadrature_tag=DISCR_TAG_BASE): + """Compute the flux for grad(u) on the boundary *dd_bdry*.""" raise NotImplementedError @abc.abstractmethod def get_diffusion_flux( - self, dcoll, dd, kappa, grad_u, *, quadrature_tag=DISCR_TAG_BASE): - """Compute the flux for diff(u) on the boundary corresponding to *dd*.""" + self, dcoll, dd_vol, dd_bdry, kappa, grad_u, *, + quadrature_tag=DISCR_TAG_BASE): + """Compute the flux for diff(u) on the boundary *dd_bdry*.""" raise NotImplementedError @@ -140,18 +146,19 @@ def __init__(self, value): self.value = value def get_grad_flux( - self, dcoll, dd, u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 - u_int = op.project(dcoll, "vol", dd, u) - u_tpair = TracePair(dd, interior=u_int, exterior=2*self.value-u_int) + self, dcoll, dd_vol, dd_bdry, u, *, + quadrature_tag=DISCR_TAG_BASE): # noqa: D102 + u_int = op.project(dcoll, dd_vol, dd_bdry, u) + u_tpair = TracePair(dd_bdry, interior=u_int, exterior=2*self.value-u_int) return grad_flux(dcoll, u_tpair, quadrature_tag=quadrature_tag) def get_diffusion_flux( - self, dcoll, dd, kappa, grad_u, *, + self, dcoll, dd_vol, dd_bdry, kappa, grad_u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 - kappa_int = op.project(dcoll, "vol", dd, kappa) - kappa_tpair = TracePair(dd, interior=kappa_int, exterior=kappa_int) - grad_u_int = op.project(dcoll, "vol", dd, grad_u) - grad_u_tpair = TracePair(dd, interior=grad_u_int, exterior=grad_u_int) + kappa_int = op.project(dcoll, dd_vol, dd_bdry, kappa) + kappa_tpair = TracePair(dd_bdry, interior=kappa_int, exterior=kappa_int) + grad_u_int = op.project(dcoll, dd_vol, dd_bdry, grad_u) + grad_u_tpair = TracePair(dd_bdry, interior=grad_u_int, exterior=grad_u_int) return diffusion_flux( dcoll, kappa_tpair, grad_u_tpair, quadrature_tag=quadrature_tag) @@ -193,25 +200,27 @@ def __init__(self, value): self.value = value def get_grad_flux( - self, dcoll, dd, u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 - u_int = op.project(dcoll, "vol", dd, u) - u_tpair = TracePair(dd, interior=u_int, exterior=u_int) + self, dcoll, dd_vol, dd_bdry, u, *, + quadrature_tag=DISCR_TAG_BASE): # noqa: D102 + u_int = op.project(dcoll, dd_vol, dd_bdry, u) + u_tpair = TracePair(dd_bdry, interior=u_int, exterior=u_int) return grad_flux(dcoll, u_tpair, quadrature_tag=quadrature_tag) def get_diffusion_flux( - self, dcoll, dd, kappa, grad_u, *, + self, dcoll, dd_vol, dd_bdry, kappa, grad_u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 - dd_quad = dd.with_discr_tag(quadrature_tag) - dd_allfaces_quad = dd_quad.with_dtag("all_faces") + dd_bdry_quad = dd_bdry.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_bdry_quad.with_domain_tag( + replace(dd_bdry_quad.domain_tag, tag=FACE_RESTR_ALL)) # Compute the flux directly instead of constructing an external grad_u value # (and the associated TracePair); this approach is simpler in the # spatially-varying kappa case (the other approach would result in a # grad_u_tpair that lives in the quadrature discretization; diffusion_flux # would need to be modified to accept such values). - kappa_int_quad = op.project(dcoll, "vol", dd_quad, kappa) - value_quad = op.project(dcoll, dd, dd_quad, self.value) + kappa_int_quad = op.project(dcoll, dd_vol, dd_bdry_quad, kappa) + value_quad = op.project(dcoll, dd_bdry, dd_bdry_quad, self.value) flux_quad = -kappa_int_quad*value_quad - return op.project(dcoll, dd_quad, dd_allfaces_quad, flux_quad) + return op.project(dcoll, dd_bdry_quad, dd_allfaces_quad, flux_quad) class _DiffusionStateTag: @@ -226,7 +235,9 @@ class _DiffusionGradTag: pass -def grad_operator(dcoll, boundaries, u, quadrature_tag=DISCR_TAG_BASE): +def grad_operator( + dcoll, boundaries, u, *, quadrature_tag=DISCR_TAG_BASE, + volume_dd=DD_VOLUME_ALL): r""" Compute the gradient of *u*. @@ -245,6 +256,8 @@ def grad_operator(dcoll, boundaries, u, quadrature_tag=DISCR_TAG_BASE): quadrature_tag: quadrature tag indicating which discretization in *dcoll* to use for overintegration + volume_dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the volume on which to apply the operator Returns ------- @@ -258,29 +271,38 @@ def grad_operator(dcoll, boundaries, u, quadrature_tag=DISCR_TAG_BASE): 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, + volume_dd=volume_dd), make_obj_array(boundaries), u) - for btag, bdry in boundaries.items(): + boundaries = { + as_dofdesc(bdtag).domain_tag: bdry + for bdtag, bdry in boundaries.items()} + + 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) + dd_vol_base = volume_dd + dd_vol_quad = dd_vol_base.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) - return op.inverse_mass(dcoll, - op.weak_local_grad(dcoll, "vol", -u) + return op.inverse_mass( + dcoll, dd_vol_base, + op.weak_local_grad(dcoll, dd_vol_base, -u) - # noqa: W504 - op.face_mass(dcoll, - dd_allfaces_quad, + op.face_mass( + dcoll, dd_allfaces_quad, sum( grad_flux(dcoll, u_tpair, quadrature_tag=quadrature_tag) for u_tpair in interior_trace_pairs( - dcoll, u, comm_tag=_DiffusionStateTag)) + dcoll, u, volume_dd=dd_vol_base, tag=_DiffusionStateTag)) + sum( - bdry.get_grad_flux( - dcoll, as_dofdesc(btag), u, quadrature_tag=quadrature_tag) - for btag, bdry in boundaries.items()) + bdry.get_grad_flux(dcoll, dd_vol_base, + dd_vol_base.with_domain_tag(bdtag), u, + quadrature_tag=quadrature_tag) + for bdtag, bdry in boundaries.items()) ) ) @@ -293,6 +315,19 @@ def _normalize_arguments(*args, **kwargs): else: pos_arg_names = ["kappa", "boundaries", "u"] + if len(args) > len(pos_arg_names): + raise TypeError( + f"diffusion_operator() takes up to {len(pos_arg_names)} positional " + f"arguments but {len(args)} were given") + + all_arg_names = [ + "alpha", "kappa", "quad_tag", "boundaries", "u", "quadrature_tag"] + for arg_name in kwargs.keys(): + if arg_name not in all_arg_names: + raise TypeError( + "diffusion_operator() got an unexpected keyword argument " + f"'{arg_name}'") + arg_dict = { arg_name: arg for arg_name, arg in zip(pos_arg_names[:len(args)], args)} @@ -325,7 +360,12 @@ def _normalize_arguments(*args, **kwargs): return kappa, boundaries, u, quadrature_tag -def diffusion_operator(dcoll, *args, return_grad_u=False, **kwargs): +def diffusion_operator( + dcoll, *args, return_grad_u=False, volume_dd=DD_VOLUME_ALL, + # Added to avoid repeated computation + # FIXME: See if there's a better way to do this + grad_u=None, + **kwargs): r""" Compute the diffusion operator. @@ -342,7 +382,7 @@ def diffusion_operator(dcoll, *args, return_grad_u=False, **kwargs): 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 @@ -352,6 +392,8 @@ def diffusion_operator(dcoll, *args, return_grad_u=False, **kwargs): quadrature_tag: quadrature tag indicating which discretization in *dcoll* to use for overintegration + volume_dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the volume on which to apply the operator Returns ------- @@ -370,38 +412,50 @@ def diffusion_operator(dcoll, *args, return_grad_u=False, **kwargs): 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, volume_dd=volume_dd), make_obj_array(boundaries), u) - for btag, bdry in boundaries.items(): + boundaries = { + as_dofdesc(bdtag).domain_tag: bdry + for bdtag, bdry in boundaries.items()} + + 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) + dd_vol_base = volume_dd + dd_vol_quad = dd_vol_base.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, + volume_dd=dd_vol_base) - 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_base, dd_vol_quad, kappa) + grad_u_quad = op.project(dcoll, dd_vol_base, dd_vol_quad, grad_u) - diff_u = op.inverse_mass(dcoll, - op.weak_local_div(dcoll, dd_quad, -kappa_quad*grad_u_quad) - - 1. # noqa: W504 - * op.face_mass(dcoll, - dd_allfaces_quad, + diff_u = op.inverse_mass( + dcoll, dd_vol_base, + op.weak_local_div(dcoll, dd_vol_quad, -kappa_quad*grad_u_quad) + - # noqa: W504 + op.face_mass( + dcoll, dd_allfaces_quad, sum( diffusion_flux( dcoll, kappa_tpair, grad_u_tpair, quadrature_tag=quadrature_tag) for kappa_tpair, grad_u_tpair in zip( - interior_trace_pairs(dcoll, kappa, comm_tag=_DiffusionKappaTag), - interior_trace_pairs(dcoll, grad_u, comm_tag=_DiffusionGradTag))) + interior_trace_pairs( + dcoll, kappa, volume_dd=dd_vol_base, tag=_DiffusionKappaTag), + interior_trace_pairs( + dcoll, grad_u, volume_dd=dd_vol_base, tag=_DiffusionGradTag)) + ) + sum( bdry.get_diffusion_flux( - dcoll, as_dofdesc(btag), kappa, grad_u, - quadrature_tag=quadrature_tag) - for btag, bdry in boundaries.items()) + dcoll, dd_vol_base, dd_vol_base.with_domain_tag(bdtag), kappa, + grad_u, quadrature_tag=quadrature_tag) + for bdtag, bdry in boundaries.items()) ) ) diff --git a/mirgecom/discretization.py b/mirgecom/discretization.py index c592f8669..40d39bfc7 100644 --- a/mirgecom/discretization.py +++ b/mirgecom/discretization.py @@ -38,11 +38,11 @@ # 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.""" 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 +51,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 5f4e8454d..f74f4f7aa 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -54,7 +54,8 @@ 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, DISCR_TAG_BASE from mirgecom.gas_model import make_operator_fluid_states from mirgecom.inviscid import ( @@ -68,7 +69,8 @@ def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, inviscid_numerical_flux_func=inviscid_facial_flux_rusanov, - quadrature_tag=None, operator_states_quad=None): + quadrature_tag=DISCR_TAG_BASE, volume_dd=DD_VOLUME_ALL, + operator_states_quad=None): r"""Compute RHS of the Euler flow equations. Returns @@ -91,7 +93,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,14 +109,17 @@ 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*. + + volume_dd: grudge.dof_desc.DOFDesc + The DOF descriptor of the volume on which to apply the operator. """ - dd_quad_vol = DOFDesc("vol", quadrature_tag) - dd_quad_faces = DOFDesc("all_faces", quadrature_tag) + dd_quad_vol = volume_dd.with_discr_tag(quadrature_tag) + dd_quad_allfaces = dd_quad_vol.trace(FACE_RESTR_ALL) if operator_states_quad is None: operator_states_quad = make_operator_fluid_states(dcoll, state, gas_model, - boundaries, quadrature_tag) + boundaries, quadrature_tag, + volume_dd=volume_dd) volume_state_quad, interior_state_pairs_quad, domain_boundary_states_quad = \ operator_states_quad @@ -125,9 +131,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, + volume_dd=volume_dd) - return -div_operator(dcoll, dd_quad_vol, dd_quad_faces, + return -div_operator(dcoll, dd_quad_vol, dd_quad_allfaces, inviscid_flux_vol, inviscid_flux_bnd) diff --git a/mirgecom/filter.py b/mirgecom/filter.py index da09d9526..cf0a1e2b0 100644 --- a/mirgecom/filter.py +++ b/mirgecom/filter.py @@ -46,11 +46,11 @@ """ 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 DISCR_TAG_MODAL, as_dofdesc -from functools import partial +from arraycontext import map_array_container from meshmode.dof_array import DOFArray @@ -205,13 +205,15 @@ def filter_modally(dcoll, dd, cutoff, mode_resp_func, field): partial(filter_modally, dcoll, dd, cutoff, mode_resp_func), field ) + dd_nodal = as_dofdesc(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 4df6481c4..41332a4b8 100644 --- a/mirgecom/gas_model.py +++ b/mirgecom/gas_model.py @@ -59,7 +59,7 @@ TransportModel, GasTransportVars ) -from grudge.dof_desc import DOFDesc, as_dofdesc +from grudge.dof_desc import DD_VOLUME_ALL, DISCR_TAG_BASE import grudge.op as op from grudge.trace_pair import ( interior_trace_pairs, @@ -384,8 +384,9 @@ class _FluidTemperatureTag: pass -def make_operator_fluid_states(dcoll, volume_state, gas_model, boundaries, - quadrature_tag=None): +def make_operator_fluid_states( + dcoll, volume_state, gas_model, boundaries, quadrature_tag=DISCR_TAG_BASE, + volume_dd=DD_VOLUME_ALL): """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 +415,15 @@ 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. + + volume_dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the volume on which to construct the states Returns ------- @@ -428,19 +432,20 @@ 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) + dd_base_vol = volume_dd + dd_quad_vol = dd_base_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_base_vol, + dd_quad_vol.with_domain_tag(bdtag), volume_state, gas_model) - for btag in boundaries + for bdtag in boundaries } # performs MPI communication of CV if needed @@ -448,7 +453,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, tag=_FluidCVTag) + for tpair in interior_trace_pairs( + dcoll, volume_state.cv, volume_dd=dd_base_vol, tag=_FluidCVTag) ] tseed_interior_pairs = None @@ -461,8 +467,9 @@ 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.temperature, - tag=_FluidTemperatureTag)] + for tpair in interior_trace_pairs( + dcoll, volume_state.temperature, volume_dd=dd_base_vol, + tag=_FluidTemperatureTag)] interior_boundary_states_quad = \ make_fluid_state_trace_pairs(cv_interior_pairs, gas_model, diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py index 6cdc7e187..87250730a 100644 --- a/mirgecom/inviscid.py +++ b/mirgecom/inviscid.py @@ -40,6 +40,8 @@ """ import numpy as np +from meshmode.discretization.connection import FACE_RESTR_ALL +from grudge.dof_desc import DD_VOLUME_ALL, DISCR_TAG_BASE import grudge.op as op from mirgecom.fluid import make_conserved @@ -224,8 +226,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, + volume_dd=DD_VOLUME_ALL): """Compute the inviscid boundary fluxes for the divergence operator. This routine encapsulates the computation of the inviscid contributions @@ -242,38 +245,42 @@ 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 + + volume_dd: grudge.dof_desc.DOFDesc + The DOF descriptor of the volume on which to compute the flux. """ - from grudge.dof_desc import as_dofdesc + dd_vol_quad = volume_dd.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): return op.project(dcoll, - dd_bdry, dd_bdry.with_dtag("all_faces"), + dd_bdry, dd_allfaces_quad, boundary.inviscid_divergence_flux( dcoll, dd_bdry, gas_model, state_minus=state_minus, numerical_flux_func=numerical_flux_func, time=time)) @@ -287,10 +294,10 @@ def _boundary_flux(dd_bdry, boundary, state_minus): # Domain boundary faces + sum( _boundary_flux( - as_dofdesc(btag).with_discr_tag(quadrature_tag), + dd_vol_quad.with_domain_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 diff --git a/mirgecom/io.py b/mirgecom/io.py index 8cb5b1829..88e64b71b 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_volume_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_volume_dd) + _max = partial(op.nodal_max, dcoll, fluid_volume_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..06d719b7f 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,23 @@ 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, volume_dd=DD_VOLUME_ALL): """Add default discretization quantities to the logmgr.""" for reduction_op in ["min", "max", "L2_norm"]: for quantity in ["pressure", "temperature"]: logmgr.add_quantity(DiscretizationBasedQuantity( dcoll, quantity, reduction_op, extract_vars_for_logging, - units_for_logging)) + units_for_logging, volume_dd=volume_dd)) for quantity in ["mass", "energy"]: logmgr.add_quantity(DiscretizationBasedQuantity( dcoll, quantity, reduction_op, extract_vars_for_logging, - units_for_logging)) + units_for_logging, volume_dd=volume_dd)) for d in range(dim): logmgr.add_quantity(DiscretizationBasedQuantity( dcoll, "momentum", reduction_op, extract_vars_for_logging, - units_for_logging, - axis=d)) + units_for_logging, axis=d, volume_dd=volume_dd)) # {{{ Package versions @@ -245,7 +245,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, volume_dd=DD_VOLUME_ALL): unit = units_logging(quantity) if name is None: @@ -262,13 +262,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, volume_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, volume_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=volume_dd) self.rank_aggr = max else: raise ValueError(f"unknown operation {op}") diff --git a/mirgecom/navierstokes.py b/mirgecom/navierstokes.py index 89e3b793e..23cc25026 100644 --- a/mirgecom/navierstokes.py +++ b/mirgecom/navierstokes.py @@ -57,14 +57,21 @@ THE SOFTWARE. """ +from dataclasses import replace 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, + DISCR_TAG_BASE, + as_dofdesc, +) import grudge.op as op @@ -98,16 +105,18 @@ 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_domain_tag( + replace(dd_trace.domain_tag, 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, + quadrature_tag=DISCR_TAG_BASE, volume_dd=DD_VOLUME_ALL, # Added to avoid repeated computation # FIXME: See if there's a better way to do this operator_states_quad=None): @@ -121,7 +130,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 +150,9 @@ def grad_cv_operator( An identifier denoting a particular quadrature discretization to use during operator evaluations. + volume_dd: grudge.dof_desc.DOFDesc + The DOF descriptor of the volume on which to apply the operator. + Returns ------- :class:`~mirgecom.fluid.ConservedVars` @@ -147,12 +160,18 @@ 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 = { + as_dofdesc(bdtag).domain_tag: bdry + for bdtag, bdry in boundaries.items()} + + dd_vol_base = volume_dd + dd_vol_quad = dd_vol_base.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) + dcoll, state, gas_model, boundaries, quadrature_tag, + volume_dd=dd_vol_base) vol_state_quad, inter_elem_bnd_states_quad, domain_bnd_states_quad = \ operator_states_quad @@ -169,18 +188,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) @@ -188,13 +205,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, + quadrature_tag=DISCR_TAG_BASE, volume_dd=DD_VOLUME_ALL, # Added to avoid repeated computation # FIXME: See if there's a better way to do this operator_states_quad=None): @@ -227,6 +244,9 @@ def grad_t_operator( An identifier denoting a particular quadrature discretization to use during operator evaluations. + volume_dd: grudge.dof_desc.DOFDesc + The DOF descriptor of the volume on which to apply the operator. + Returns ------- :class:`numpy.ndarray` @@ -234,12 +254,18 @@ 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 = { + as_dofdesc(bdtag).domain_tag: bdry + for bdtag, bdry in boundaries.items()} + + dd_vol_base = volume_dd + dd_vol_quad = dd_vol_base.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) + dcoll, state, gas_model, boundaries, quadrature_tag, + volume_dd=dd_vol_base) vol_state_quad, inter_elem_bnd_states_quad, domain_bnd_states_quad = \ operator_states_quad @@ -260,18 +286,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) @@ -279,14 +303,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, + return_gradients=False, quadrature_tag=DISCR_TAG_BASE, + volume_dd=DD_VOLUME_ALL, # Added to avoid repeated computation # FIXME: See if there's a better way to do this operator_states_quad=None, @@ -325,10 +350,18 @@ 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. + volume_dd: grudge.dof_desc.DOFDesc + The DOF descriptor of the volume on which to apply the operator. + operator_states_quad Optional iterable container providing the full fluid states (:class:`~mirgecom.gas_model.FluidState`) on the quadrature @@ -347,11 +380,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. - Returns ------- :class:`mirgecom.fluid.ConservedVars` @@ -365,9 +393,13 @@ 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 = { + as_dofdesc(bdtag).domain_tag: bdry + for bdtag, bdry in boundaries.items()} + + dd_vol_base = volume_dd + dd_vol_quad = dd_vol_base.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 @@ -378,7 +410,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) + dcoll, state, gas_model, boundaries, quadrature_tag, + volume_dd=dd_vol_base) vol_state_quad, inter_elem_bnd_states_quad, domain_bnd_states_quad = \ operator_states_quad @@ -396,7 +429,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, volume_dd=dd_vol_base, operator_states_quad=operator_states_quad) # Communicate grad(CV) and put it on the quadrature domain @@ -404,7 +437,8 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, # 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, grad_cv, tag=_NSGradCVTag) + for tpair in interior_trace_pairs( + dcoll, grad_cv, volume_dd=dd_vol_base, tag=_NSGradCVTag) ] # }}} Compute grad(CV) @@ -415,7 +449,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, volume_dd=dd_vol_base, operator_states_quad=operator_states_quad) # Create the interior face trace pairs, perform MPI exchange, interp to quad @@ -423,7 +457,8 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, # 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, grad_t, tag=_NSGradTemperatureTag) + for tpair in interior_trace_pairs( + dcoll, grad_t, volume_dd=dd_vol_base, tag=_NSGradTemperatureTag) ] # }}} compute grad(temperature) @@ -437,8 +472,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_base, dd_vol_quad, grad_cv), + grad_t=op.project(dcoll, dd_vol_base, dd_vol_quad, grad_t)) # Compute the volume contribution of the inviscid flux terms # using field values on the quadrature grid @@ -453,16 +488,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, + volume_dd=dd_vol_base) # 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, + volume_dd=dd_vol_base) ) - 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..7a35ea6b7 100644 --- a/mirgecom/operators.py +++ b/mirgecom/operators.py @@ -30,6 +30,8 @@ import grudge.op as op +from grudge.dof_desc import DISCR_TAG_BASE + def grad_operator(dcoll, dd_vol, dd_faces, u, flux): r"""Compute a DG gradient for the input *u* with flux given by *flux*. @@ -57,7 +59,8 @@ 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)) @@ -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)) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index 95627a06e..54669a03d 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_volume_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 @@ -135,7 +138,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_volume_dd, get_viscous_timestep(dcoll=dcoll, state=state)))[()] return min(t_remaining, mydt) @@ -335,7 +338,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 +347,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 +360,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 +376,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/viscous.py b/mirgecom/viscous.py index 261c0d2dd..65988999e 100644 --- a/mirgecom/viscous.py +++ b/mirgecom/viscous.py @@ -45,6 +45,8 @@ 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, DISCR_TAG_BASE import grudge.op as op @@ -335,8 +337,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, + volume_dd=DD_VOLUME_ALL): """Compute the viscous boundary fluxes for the divergence operator. This routine encapsulates the computation of the viscous contributions @@ -351,14 +354,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 +378,42 @@ 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 - """ - from grudge.dof_desc import as_dofdesc - dd_base = as_dofdesc("vol") + volume_dd: grudge.dof_desc.DOFDesc + The DOF descriptor of the volume on which to compute the flux. + """ + dd_base = volume_dd + dd_vol_quad = dd_base.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): + dd_bdry = dd_vol_quad.with_domain_tag(bdtag) return op.project( - dcoll, dd_btag, dd_btag.with_dtag("all_faces"), + dcoll, dd_bdry, dd_allfaces_quad, boundary.viscous_divergence_flux( - dcoll=dcoll, btag=dd_btag, gas_model=gas_model, + dcoll=dcoll, dd_bdry=dd_bdry, 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), + grad_cv_minus=op.project(dcoll, dd_base, dd_bdry, grad_cv), + grad_t_minus=op.project(dcoll, dd_base, dd_bdry, grad_t), time=time, numerical_flux_func=numerical_flux_func)) # }}} viscous flux helpers @@ -420,9 +425,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 +442,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, *, volume_dd=DD_VOLUME_ALL): """Routine returns the the node-local maximum stable viscous timestep. Parameters @@ -457,7 +463,8 @@ def get_viscous_timestep(dcoll, state): """ from grudge.dt_utils import characteristic_lengthscales - length_scales = characteristic_lengthscales(state.array_context, dcoll) + length_scales = characteristic_lengthscales( + state.array_context, dcoll, dd=volume_dd) nu = 0 d_alpha_max = 0 @@ -475,7 +482,7 @@ def get_viscous_timestep(dcoll, state): ) -def get_viscous_cfl(dcoll, dt, state): +def get_viscous_cfl(dcoll, dt, state, *, volume_dd=DD_VOLUME_ALL): """Calculate and return node-local CFL based on current state and timestep. Parameters @@ -498,7 +505,7 @@ def get_viscous_cfl(dcoll, dt, state): The CFL at each node. """ - return dt / get_viscous_timestep(dcoll, state=state) + return dt / get_viscous_timestep(dcoll, state=state, volume_dd=volume_dd) def get_local_max_species_diffusivity(actx, d_alpha): diff --git a/mirgecom/wave.py b/mirgecom/wave.py index 97cf60b0d..5db446b43 100644 --- a/mirgecom/wave.py +++ b/mirgecom/wave.py @@ -88,8 +88,8 @@ def wave_operator(dcoll, c, w): 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..94e801701 100644 --- a/requirements.txt +++ b/requirements.txt @@ -16,9 +16,9 @@ pyyaml --editable git+https://github.com/inducer/dagrt.git#egg=dagrt --editable git+https://github.com/inducer/leap.git#egg=leap --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/inducer/pytato.git#egg=pytato +--editable git+https://github.com/kaushikcfd/arraycontext.git#egg=arraycontext +--editable git+https://github.com/majosm/meshmode.git@production#egg=meshmode +--editable git+https://github.com/majosm/grudge.git@multi-volume-misc#egg=grudge +--editable git+https://github.com/kaushikcfd/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 22d95b604..34cb76fb5 100644 --- a/test/test_av.py +++ b/test/test_av.py @@ -205,29 +205,28 @@ def test_artificial_viscosity(ctx_factory, dim, order): nodes = actx.thaw(dcoll.nodes()) class TestBoundary: - - def cv_gradient_flux(self, disc, btag, state_minus, gas_model, **kwargs): + def cv_gradient_flux(self, disc, 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(disc.normal(btag)) + nhat = actx.thaw(disc.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, disc, btag, diffusion, **kwargs): - nhat = actx.thaw(disc.normal(btag)) - diffusion_minus = op.project(dcoll, "vol", btag, diffusion) + def av_flux(self, disc, dd_bdry, diffusion, **kwargs): + nhat = actx.thaw(disc.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(disc, btag, "all_faces", flux_weak) + return op.project(dcoll, dd_bdry, "all_faces", flux_weak) boundaries = {BTAG_ALL: TestBoundary()} @@ -402,9 +401,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) diff --git a/test/test_bc.py b/test/test_bc.py index a8b02104c..21c4d1a88 100644 --- a/test/test_bc.py +++ b/test/test_bc.py @@ -103,11 +103,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 +181,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 +316,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 +338,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,7 +348,7 @@ 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) @@ -372,7 +380,7 @@ def scalar_flux_interior(int_tpair): 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 +469,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 +529,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 +542,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 +553,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) @@ -579,11 +587,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..cf0dbb0da 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) @@ -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..bcd587f9c 100644 --- a/test/test_filter.py +++ b/test/test_filter.py @@ -230,9 +230,9 @@ 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)] 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..8aa42bbe6 100644 --- a/test/test_inviscid.py +++ b/test/test_inviscid.py @@ -37,6 +37,7 @@ ) from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from grudge.dof_desc import as_dofdesc from grudge.trace_pair import TracePair from mirgecom.fluid import make_conserved from mirgecom.eos import IdealSingleGas @@ -375,7 +376,7 @@ 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)) diff --git a/test/test_lazy.py b/test/test_lazy.py index 044376f98..9c41e37b9 100644 --- a/test/test_lazy.py +++ b/test/test_lazy.py @@ -102,15 +102,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): @@ -173,15 +173,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 +239,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..daccbd4fa --- /dev/null +++ b/test/test_multiphysics.py @@ -0,0 +1,142 @@ +__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 meshmode.mesh import BTAG_PARTITION +from grudge.dof_desc import DOFDesc, VolumeDomainTag, DISCR_TAG_BASE +from grudge.discretization import PartID +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__) + + +def get_box_mesh(dim, a, b, n): + dim_names = ["x", "y", "z"] + boundary_tag_to_face = {} + for i in range(dim): + boundary_tag_to_face["-"+str(i)] = ["-"+dim_names[i]] + boundary_tag_to_face["+"+str(i)] = ["+"+dim_names[i]] + from meshmode.mesh.generation import generate_regular_rect_mesh + return generate_regular_rect_mesh(a=(a,)*dim, b=(b,)*dim, + nelements_per_axis=(n,)*dim, boundary_tag_to_face=boundary_tag_to_face) + + +@pytest.mark.parametrize("order", [1, 2, 3, 4]) +def test_independent_volumes(actx_factory, order, visualize=False): + """Check multi-volume machinery by setting up two independent volumes.""" + actx = actx_factory() + + n = 8 + global_mesh = get_box_mesh(2, -1, 1, n) + + mgrp, = global_mesh.groups + y = global_mesh.vertices[1, mgrp.vertex_indices] + y_elem_avg = np.sum(y, axis=1)/y.shape[1] + volume_to_elements = { + "Lower": np.where(y_elem_avg < 0)[0], + "Upper": np.where(y_elem_avg > 0)[0]} + + from meshmode.mesh.processing import partition_mesh + volume_meshes = partition_mesh(global_mesh, volume_to_elements) + + dcoll = create_discretization_collection(actx, volume_meshes, order) + + dd_vol_lower = DOFDesc(VolumeDomainTag("Lower"), DISCR_TAG_BASE) + dd_vol_upper = DOFDesc(VolumeDomainTag("Upper"), DISCR_TAG_BASE) + + lower_nodes = actx.thaw(dcoll.nodes(dd=dd_vol_lower)) + upper_nodes = actx.thaw(dcoll.nodes(dd=dd_vol_upper)) + + lower_boundaries = { + dd_vol_lower.trace("-0").domain_tag: NeumannDiffusionBoundary(0.), + dd_vol_lower.trace("+0").domain_tag: NeumannDiffusionBoundary(0.), + dd_vol_lower.trace("-1").domain_tag: DirichletDiffusionBoundary(0.), + dd_vol_lower.trace(BTAG_PARTITION(PartID("Upper"))).domain_tag: + DirichletDiffusionBoundary(1.), + } + + upper_boundaries = { + dd_vol_upper.trace("-0").domain_tag: NeumannDiffusionBoundary(0.), + dd_vol_upper.trace("+0").domain_tag: NeumannDiffusionBoundary(0.), + dd_vol_upper.trace(BTAG_PARTITION(PartID("Lower"))).domain_tag: + DirichletDiffusionBoundary(0.), + dd_vol_upper.trace("+1"): DirichletDiffusionBoundary(1.), + } + + lower_u = lower_nodes[1] + 1 + upper_u = upper_nodes[1] + + u = make_obj_array([lower_u, upper_u]) + + def get_rhs(t, u): + return make_obj_array([ + diffusion_operator( + dcoll, kappa=1, boundaries=lower_boundaries, u=u[0], + volume_dd=dd_vol_lower), + diffusion_operator( + dcoll, kappa=1, boundaries=upper_boundaries, u=u[1], + volume_dd=dd_vol_upper)]) + + rhs = get_rhs(0, u) + + if visualize: + from grudge.shortcuts import make_visualizer + viz_lower = make_visualizer(dcoll, order+3, volume_dd=dd_vol_lower) + viz_upper = make_visualizer(dcoll, order+3, volume_dd=dd_vol_upper) + viz_lower.write_vtk_file( + f"multiphysics_independent_volumes_{order}_lower.vtu", [ + ("u", u[0]), + ("rhs", rhs[0]), + ]) + viz_upper.write_vtk_file( + f"multiphysics_independent_volumes_{order}_upper.vtu", [ + ("u", u[1]), + ("rhs", rhs[1]), + ]) + + linf_err_lower = actx.to_numpy(op.norm(dcoll, rhs[0], np.inf, dd=dd_vol_lower)) + linf_err_upper = actx.to_numpy(op.norm(dcoll, rhs[1], np.inf, dd=dd_vol_upper)) + + assert linf_err_lower < 1e-9 + assert linf_err_upper < 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..749fea536 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 @@ -128,14 +130,14 @@ def central_flux_interior(actx, dcoll, int_tpair): return op.project(dcoll, int_tpair.dd, dd_all_faces, 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") @@ -225,7 +227,6 @@ 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, diff --git a/test/test_viscous.py b/test/test_viscous.py index c595c9372..9b06c18c7 100644 --- a/test/test_viscous.py +++ b/test/test_viscous.py @@ -30,9 +30,12 @@ import pyopencl.clmath # noqa import logging import pytest # noqa +from dataclasses import replace 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 +164,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 @@ -173,17 +177,18 @@ def cv_flux_interior(int_tpair): dd_all_faces = int_tpair.dd.with_dtag("all_faces") return op.project(dcoll, int_tpair.dd, dd_all_faces, 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_all_faces = dd_bdry.with_domain_tag( + replace(dd_bdry.domain_tag, tag=FACE_RESTR_ALL)) + return op.project(dcoll, dd_bdry, dd_all_faces, flux_weak) for nfac in [1, 2, 4]: @@ -213,7 +218,6 @@ 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) From d7f83ce594693a44f759f6251a955bceb22dcfe8 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Mon, 15 Aug 2022 13:59:09 -0500 Subject: [PATCH 02/19] deprecate mpi_communicator argument in create_discretization_collection --- examples/autoignition-mpi.py | 6 ++---- examples/doublemach-mpi.py | 3 +-- examples/heat-source-mpi.py | 3 +-- examples/hotplate-mpi.py | 4 +--- examples/lump-mpi.py | 4 +--- examples/mixture-mpi.py | 4 +--- examples/nsmix-mpi.py | 4 +--- examples/poiseuille-mpi.py | 4 +--- examples/pulse-mpi.py | 4 +--- examples/scalar-lump-mpi.py | 4 +--- examples/sod-mpi.py | 4 +--- examples/vortex-mpi.py | 4 +--- examples/wave-mpi.py | 7 ++----- mirgecom/discretization.py | 6 ++++++ 14 files changed, 21 insertions(+), 40 deletions(-) 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 243eed84d..39c4daeb1 100644 --- a/examples/doublemach-mpi.py +++ b/examples/doublemach-mpi.py @@ -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 diff --git a/examples/heat-source-mpi.py b/examples/heat-source-mpi.py index a21279817..7218b2a29 100644 --- a/examples/heat-source-mpi.py +++ b/examples/heat-source-mpi.py @@ -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 diff --git a/examples/hotplate-mpi.py b/examples/hotplate-mpi.py index 181d922ce..001d3b968 100644 --- a/examples/hotplate-mpi.py +++ b/examples/hotplate-mpi.py @@ -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: diff --git a/examples/lump-mpi.py b/examples/lump-mpi.py index 6054f3d1e..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 diff --git a/examples/mixture-mpi.py b/examples/mixture-mpi.py index 07d613c04..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 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 9d5506341..dc3492827 100644 --- a/examples/poiseuille-mpi.py +++ b/examples/poiseuille-mpi.py @@ -167,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: 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 3d3653cb2..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 diff --git a/examples/sod-mpi.py b/examples/sod-mpi.py index 9645e0224..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 diff --git a/examples/vortex-mpi.py b/examples/vortex-mpi.py index 77491aef6..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 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/discretization.py b/mirgecom/discretization.py index 40d39bfc7..37677b7a3 100644 --- a/mirgecom/discretization.py +++ b/mirgecom/discretization.py @@ -41,6 +41,12 @@ 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 make_discretization_collection from meshmode.discretization.poly_element import ( From ed2be8865109cbf6d21d777e55235835bdff8ff7 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Fri, 26 Aug 2022 12:01:15 -0500 Subject: [PATCH 03/19] don't depend on full grudge multi-volume support yet --- requirements.txt | 8 +++--- test/test_multiphysics.py | 53 ++++++++++++++++++--------------------- 2 files changed, 29 insertions(+), 32 deletions(-) diff --git a/requirements.txt b/requirements.txt index 94e801701..9931d6448 100644 --- a/requirements.txt +++ b/requirements.txt @@ -16,9 +16,9 @@ pyyaml --editable git+https://github.com/inducer/dagrt.git#egg=dagrt --editable git+https://github.com/inducer/leap.git#egg=leap --editable git+https://github.com/inducer/modepy.git#egg=modepy ---editable git+https://github.com/kaushikcfd/arraycontext.git#egg=arraycontext ---editable git+https://github.com/majosm/meshmode.git@production#egg=meshmode ---editable git+https://github.com/majosm/grudge.git@multi-volume-misc#egg=grudge ---editable git+https://github.com/kaushikcfd/pytato.git#egg=pytato +--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/majosm/grudge.git@multiple-independent-volumes#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_multiphysics.py b/test/test_multiphysics.py index daccbd4fa..e9f0d0820 100644 --- a/test/test_multiphysics.py +++ b/test/test_multiphysics.py @@ -29,9 +29,7 @@ diffusion_operator, DirichletDiffusionBoundary, NeumannDiffusionBoundary) -from meshmode.mesh import BTAG_PARTITION from grudge.dof_desc import DOFDesc, VolumeDomainTag, DISCR_TAG_BASE -from grudge.discretization import PartID from mirgecom.discretization import create_discretization_collection from meshmode.array_context import ( # noqa pytest_generate_tests_for_pyopencl_array_context @@ -42,34 +40,35 @@ logger = logging.getLogger(__name__) -def get_box_mesh(dim, a, b, n): +@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 - return generate_regular_rect_mesh(a=(a,)*dim, b=(b,)*dim, - nelements_per_axis=(n,)*dim, boundary_tag_to_face=boundary_tag_to_face) + from meshmode.mesh.generation import generate_regular_rect_mesh -@pytest.mark.parametrize("order", [1, 2, 3, 4]) -def test_independent_volumes(actx_factory, order, visualize=False): - """Check multi-volume machinery by setting up two independent volumes.""" - actx = actx_factory() - - n = 8 - global_mesh = get_box_mesh(2, -1, 1, n) + global_mesh1 = generate_regular_rect_mesh( + a=(-1, -2), b=(1, -1), + nelements_per_axis=(n,)*dim, boundary_tag_to_face=boundary_tag_to_face) - mgrp, = global_mesh.groups - y = global_mesh.vertices[1, mgrp.vertex_indices] - y_elem_avg = np.sum(y, axis=1)/y.shape[1] - volume_to_elements = { - "Lower": np.where(y_elem_avg < 0)[0], - "Upper": np.where(y_elem_avg > 0)[0]} + global_mesh2 = generate_regular_rect_mesh( + a=(-1, 1), b=(1, 2), + nelements_per_axis=(n,)*dim, boundary_tag_to_face=boundary_tag_to_face) - from meshmode.mesh.processing import partition_mesh - volume_meshes = partition_mesh(global_mesh, volume_to_elements) + volume_meshes = { + "Lower": global_mesh1, + "Upper": global_mesh2, + } dcoll = create_discretization_collection(actx, volume_meshes, order) @@ -82,20 +81,18 @@ def test_independent_volumes(actx_factory, order, visualize=False): lower_boundaries = { dd_vol_lower.trace("-0").domain_tag: NeumannDiffusionBoundary(0.), dd_vol_lower.trace("+0").domain_tag: NeumannDiffusionBoundary(0.), - dd_vol_lower.trace("-1").domain_tag: DirichletDiffusionBoundary(0.), - dd_vol_lower.trace(BTAG_PARTITION(PartID("Upper"))).domain_tag: - DirichletDiffusionBoundary(1.), + dd_vol_lower.trace("-1").domain_tag: DirichletDiffusionBoundary(-2.), + dd_vol_lower.trace("+1").domain_tag: DirichletDiffusionBoundary(-1.), } upper_boundaries = { dd_vol_upper.trace("-0").domain_tag: NeumannDiffusionBoundary(0.), dd_vol_upper.trace("+0").domain_tag: NeumannDiffusionBoundary(0.), - dd_vol_upper.trace(BTAG_PARTITION(PartID("Lower"))).domain_tag: - DirichletDiffusionBoundary(0.), - dd_vol_upper.trace("+1"): DirichletDiffusionBoundary(1.), + dd_vol_upper.trace("-1").domain_tag: DirichletDiffusionBoundary(1.), + dd_vol_upper.trace("+1").domain_tag: DirichletDiffusionBoundary(2.), } - lower_u = lower_nodes[1] + 1 + lower_u = lower_nodes[1] upper_u = upper_nodes[1] u = make_obj_array([lower_u, upper_u]) From 191d7ca2eaf8ab90058271d7aeda090f6ec7f404 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Tue, 30 Aug 2022 11:46:58 -0500 Subject: [PATCH 04/19] change grudge branch --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 9931d6448..bf9dc1a29 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/majosm/grudge.git@multiple-independent-volumes#egg=grudge +--editable git+https://github.com/majosm/grudge.git@dof-desc-helpers#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 From 6b475f98537961f0a574ae4619397a27af8d8739 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Tue, 30 Aug 2022 10:36:07 -0500 Subject: [PATCH 05/19] use new DOFDesc methods --- mirgecom/artificial_viscosity.py | 4 +--- mirgecom/boundary.py | 9 +++------ mirgecom/diffusion.py | 10 +++------- mirgecom/navierstokes.py | 4 +--- test/test_viscous.py | 4 +--- 5 files changed, 9 insertions(+), 22 deletions(-) diff --git a/mirgecom/artificial_viscosity.py b/mirgecom/artificial_viscosity.py index e16fe9208..565fe873c 100644 --- a/mirgecom/artificial_viscosity.py +++ b/mirgecom/artificial_viscosity.py @@ -126,7 +126,6 @@ """ import numpy as np -from dataclasses import replace from pytools import memoize_in, keyed_memoize_in from functools import partial @@ -265,8 +264,7 @@ def interp_to_vol_quad(u): def central_flux_div(utpair): dd_trace = utpair.dd - dd_allfaces = dd_trace.with_domain_tag( - replace(dd_trace.domain_tag, tag=FACE_RESTR_ALL)) + dd_allfaces = dd_trace.with_boundary_tag(FACE_RESTR_ALL) normal = actx.thaw(dcoll.normal(dd_trace)) return op.project(dcoll, dd_trace, dd_allfaces, # This uses a central vector flux along nhat: diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py index 3cf4a45a7..53a7717ef 100644 --- a/mirgecom/boundary.py +++ b/mirgecom/boundary.py @@ -44,10 +44,9 @@ """ import numpy as np -from dataclasses import replace from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from meshmode.discretization.connection import FACE_RESTR_ALL -from grudge.dof_desc import VolumeDomainTag, as_dofdesc +from grudge.dof_desc import as_dofdesc from mirgecom.fluid import make_conserved from grudge.trace_pair import TracePair import grudge.op as op @@ -319,8 +318,7 @@ def __init__(self, def _boundary_quantity(self, dcoll, dd_bdry, quantity, local=False, **kwargs): """Get a boundary quantity on local boundary, or projected to "all_faces".""" - dd_allfaces = dd_bdry.with_domain_tag( - replace(dd_bdry.domain_tag, tag=FACE_RESTR_ALL)) + dd_allfaces = dd_bdry.with_boundary_tag(FACE_RESTR_ALL) return quantity if local else op.project(dcoll, dd_bdry, dd_allfaces, quantity) @@ -480,8 +478,7 @@ def _identical_grad_av(self, grad_av_minus, **kwargs): def av_flux(self, dcoll, dd_bdry, diffusion, **kwargs): """Get the diffusive fluxes for the AV operator API.""" dd_bdry = as_dofdesc(dd_bdry) - dd_vol = dd_bdry.with_domain_tag( - VolumeDomainTag(dd_bdry.domain_tag.volume_tag)) + dd_vol = dd_bdry.untrace() grad_av_minus = op.project(dcoll, dd_vol, dd_bdry, diffusion) actx = grad_av_minus.mass[0].array_context nhat = actx.thaw(dcoll.normal(dd_bdry)) diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py index af9b7ed5f..2212434a0 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -36,7 +36,6 @@ import abc import numpy as np import numpy.linalg as la # noqa -from dataclasses import replace from pytools.obj_array import make_obj_array, obj_array_vectorize_n_args from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from meshmode.discretization.connection import FACE_RESTR_ALL # noqa @@ -51,8 +50,7 @@ def grad_flux(dcoll, u_tpair, *, quadrature_tag=DISCR_TAG_BASE): dd_trace = u_tpair.dd dd_trace_quad = dd_trace.with_discr_tag(quadrature_tag) - dd_allfaces_quad = dd_trace_quad.with_domain_tag( - replace(dd_trace_quad.domain_tag, tag=FACE_RESTR_ALL)) + dd_allfaces_quad = dd_trace_quad.with_boundary_tag(FACE_RESTR_ALL) normal_quad = actx.thaw(dcoll.normal(dd_trace_quad)) @@ -73,8 +71,7 @@ def diffusion_flux( dd_trace = grad_u_tpair.dd dd_trace_quad = dd_trace.with_discr_tag(quadrature_tag) - dd_allfaces_quad = dd_trace_quad.with_domain_tag( - replace(dd_trace_quad.domain_tag, tag=FACE_RESTR_ALL)) + dd_allfaces_quad = dd_trace_quad.with_boundary_tag(FACE_RESTR_ALL) normal_quad = actx.thaw(dcoll.normal(dd_trace_quad)) @@ -210,8 +207,7 @@ def get_diffusion_flux( self, dcoll, dd_vol, dd_bdry, kappa, grad_u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 dd_bdry_quad = dd_bdry.with_discr_tag(quadrature_tag) - dd_allfaces_quad = dd_bdry_quad.with_domain_tag( - replace(dd_bdry_quad.domain_tag, tag=FACE_RESTR_ALL)) + dd_allfaces_quad = dd_bdry_quad.with_boundary_tag(FACE_RESTR_ALL) # Compute the flux directly instead of constructing an external grad_u value # (and the associated TracePair); this approach is simpler in the # spatially-varying kappa case (the other approach would result in a diff --git a/mirgecom/navierstokes.py b/mirgecom/navierstokes.py index 23cc25026..98c924950 100644 --- a/mirgecom/navierstokes.py +++ b/mirgecom/navierstokes.py @@ -57,7 +57,6 @@ THE SOFTWARE. """ -from dataclasses import replace from functools import partial from meshmode.discretization.connection import FACE_RESTR_ALL @@ -106,8 +105,7 @@ def _gradient_flux_interior(dcoll, numerical_flux_func, tpair): from arraycontext import outer actx = tpair.int.array_context dd_trace = tpair.dd - dd_allfaces = dd_trace.with_domain_tag( - replace(dd_trace.domain_tag, tag=FACE_RESTR_ALL)) + 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_trace, dd_allfaces, flux) diff --git a/test/test_viscous.py b/test/test_viscous.py index 9b06c18c7..198164ab1 100644 --- a/test/test_viscous.py +++ b/test/test_viscous.py @@ -30,7 +30,6 @@ import pyopencl.clmath # noqa import logging import pytest # noqa -from dataclasses import replace from pytools.obj_array import make_obj_array from meshmode.discretization.connection import FACE_RESTR_ALL @@ -186,8 +185,7 @@ def cv_flux_boundary(dd_bdry): 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 = dd_bdry.with_domain_tag( - replace(dd_bdry.domain_tag, tag=FACE_RESTR_ALL)) + dd_all_faces = dd_bdry.with_boundary_tag(FACE_RESTR_ALL) return op.project(dcoll, dd_bdry, dd_all_faces, flux_weak) for nfac in [1, 2, 4]: From e861f112461ff9137acdc1b94a2bf6190433e128 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Tue, 30 Aug 2022 11:47:15 -0500 Subject: [PATCH 06/19] clean up some DOFDesc stuff --- mirgecom/boundary.py | 3 +-- mirgecom/diffusion.py | 28 ++++++++++++---------------- 2 files changed, 13 insertions(+), 18 deletions(-) diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py index 53a7717ef..94c4e934b 100644 --- a/mirgecom/boundary.py +++ b/mirgecom/boundary.py @@ -478,8 +478,7 @@ def _identical_grad_av(self, grad_av_minus, **kwargs): def av_flux(self, dcoll, dd_bdry, diffusion, **kwargs): """Get the diffusive fluxes for the AV operator API.""" dd_bdry = as_dofdesc(dd_bdry) - dd_vol = dd_bdry.untrace() - grad_av_minus = op.project(dcoll, dd_vol, dd_bdry, diffusion) + 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(dd_bdry)) grad_av_plus = self._bnd_grad_av_func( diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py index 2212434a0..4162f4333 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -101,15 +101,13 @@ class DiffusionBoundary(metaclass=abc.ABCMeta): @abc.abstractmethod def get_grad_flux( - self, dcoll, dd_vol, dd_bdry, u, *, - quadrature_tag=DISCR_TAG_BASE): + self, dcoll, dd_bdry, u, *, quadrature_tag=DISCR_TAG_BASE): """Compute the flux for grad(u) on the boundary *dd_bdry*.""" raise NotImplementedError @abc.abstractmethod def get_diffusion_flux( - self, dcoll, dd_vol, dd_bdry, kappa, grad_u, *, - quadrature_tag=DISCR_TAG_BASE): + self, dcoll, dd_bdry, kappa, grad_u, *, quadrature_tag=DISCR_TAG_BASE): """Compute the flux for diff(u) on the boundary *dd_bdry*.""" raise NotImplementedError @@ -143,15 +141,15 @@ def __init__(self, value): self.value = value def get_grad_flux( - self, dcoll, dd_vol, dd_bdry, u, *, - quadrature_tag=DISCR_TAG_BASE): # noqa: D102 - u_int = op.project(dcoll, dd_vol, dd_bdry, u) + self, dcoll, dd_bdry, u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 + u_int = op.project(dcoll, dd_bdry.untrace(), dd_bdry, u) u_tpair = TracePair(dd_bdry, interior=u_int, exterior=2*self.value-u_int) return grad_flux(dcoll, u_tpair, quadrature_tag=quadrature_tag) def get_diffusion_flux( - self, dcoll, dd_vol, dd_bdry, kappa, grad_u, *, + self, dcoll, dd_bdry, kappa, grad_u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 + dd_vol = dd_bdry.untrace() kappa_int = op.project(dcoll, dd_vol, dd_bdry, kappa) kappa_tpair = TracePair(dd_bdry, interior=kappa_int, exterior=kappa_int) grad_u_int = op.project(dcoll, dd_vol, dd_bdry, grad_u) @@ -197,14 +195,13 @@ def __init__(self, value): self.value = value def get_grad_flux( - self, dcoll, dd_vol, dd_bdry, u, *, - quadrature_tag=DISCR_TAG_BASE): # noqa: D102 - u_int = op.project(dcoll, dd_vol, dd_bdry, u) + self, dcoll, dd_bdry, u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 + u_int = op.project(dcoll, dd_bdry.untrace(), dd_bdry, u) u_tpair = TracePair(dd_bdry, interior=u_int, exterior=u_int) return grad_flux(dcoll, u_tpair, quadrature_tag=quadrature_tag) def get_diffusion_flux( - self, dcoll, dd_vol, dd_bdry, kappa, grad_u, *, + self, dcoll, dd_bdry, kappa, grad_u, *, quadrature_tag=DISCR_TAG_BASE): # noqa: D102 dd_bdry_quad = dd_bdry.with_discr_tag(quadrature_tag) dd_allfaces_quad = dd_bdry_quad.with_boundary_tag(FACE_RESTR_ALL) @@ -213,7 +210,7 @@ def get_diffusion_flux( # spatially-varying kappa case (the other approach would result in a # grad_u_tpair that lives in the quadrature discretization; diffusion_flux # would need to be modified to accept such values). - kappa_int_quad = op.project(dcoll, dd_vol, dd_bdry_quad, kappa) + kappa_int_quad = op.project(dcoll, dd_bdry.untrace(), dd_bdry_quad, kappa) value_quad = op.project(dcoll, dd_bdry, dd_bdry_quad, self.value) flux_quad = -kappa_int_quad*value_quad return op.project(dcoll, dd_bdry_quad, dd_allfaces_quad, flux_quad) @@ -295,8 +292,7 @@ def grad_operator( for u_tpair in interior_trace_pairs( dcoll, u, volume_dd=dd_vol_base, tag=_DiffusionStateTag)) + sum( - bdry.get_grad_flux(dcoll, dd_vol_base, - dd_vol_base.with_domain_tag(bdtag), u, + bdry.get_grad_flux(dcoll, dd_vol_base.with_domain_tag(bdtag), u, quadrature_tag=quadrature_tag) for bdtag, bdry in boundaries.items()) ) @@ -449,7 +445,7 @@ def diffusion_operator( ) + sum( bdry.get_diffusion_flux( - dcoll, dd_vol_base, dd_vol_base.with_domain_tag(bdtag), kappa, + dcoll, dd_vol_base.with_domain_tag(bdtag), kappa, grad_u, quadrature_tag=quadrature_tag) for bdtag, bdry in boundaries.items()) ) From 85e8dcbc1ac27a87232e4e21653dc054a2472426 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Wed, 31 Aug 2022 14:51:52 -0500 Subject: [PATCH 07/19] tweak indepedent volumes test --- test/test_multiphysics.py | 74 ++++++++++++++++++--------------------- 1 file changed, 35 insertions(+), 39 deletions(-) diff --git a/test/test_multiphysics.py b/test/test_multiphysics.py index e9f0d0820..f440b69f3 100644 --- a/test/test_multiphysics.py +++ b/test/test_multiphysics.py @@ -57,77 +57,73 @@ def test_independent_volumes(actx_factory, order, visualize=False): from meshmode.mesh.generation import generate_regular_rect_mesh - global_mesh1 = generate_regular_rect_mesh( - a=(-1, -2), b=(1, -1), - nelements_per_axis=(n,)*dim, boundary_tag_to_face=boundary_tag_to_face) - - global_mesh2 = generate_regular_rect_mesh( - a=(-1, 1), b=(1, 2), + 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 = { - "Lower": global_mesh1, - "Upper": global_mesh2, + "vol1": mesh, + "vol2": mesh, } dcoll = create_discretization_collection(actx, volume_meshes, order) - dd_vol_lower = DOFDesc(VolumeDomainTag("Lower"), DISCR_TAG_BASE) - dd_vol_upper = DOFDesc(VolumeDomainTag("Upper"), DISCR_TAG_BASE) + dd_vol1 = DOFDesc(VolumeDomainTag("vol1"), DISCR_TAG_BASE) + dd_vol2 = DOFDesc(VolumeDomainTag("vol2"), DISCR_TAG_BASE) - lower_nodes = actx.thaw(dcoll.nodes(dd=dd_vol_lower)) - upper_nodes = actx.thaw(dcoll.nodes(dd=dd_vol_upper)) + nodes1 = actx.thaw(dcoll.nodes(dd=dd_vol1)) + nodes2 = actx.thaw(dcoll.nodes(dd=dd_vol2)) - lower_boundaries = { - dd_vol_lower.trace("-0").domain_tag: NeumannDiffusionBoundary(0.), - dd_vol_lower.trace("+0").domain_tag: NeumannDiffusionBoundary(0.), - dd_vol_lower.trace("-1").domain_tag: DirichletDiffusionBoundary(-2.), - dd_vol_lower.trace("+1").domain_tag: DirichletDiffusionBoundary(-1.), + 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.), } - upper_boundaries = { - dd_vol_upper.trace("-0").domain_tag: NeumannDiffusionBoundary(0.), - dd_vol_upper.trace("+0").domain_tag: NeumannDiffusionBoundary(0.), - dd_vol_upper.trace("-1").domain_tag: DirichletDiffusionBoundary(1.), - dd_vol_upper.trace("+1").domain_tag: DirichletDiffusionBoundary(2.), + 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.), } - lower_u = lower_nodes[1] - upper_u = upper_nodes[1] + u1 = nodes1[0] + u2 = nodes2[1] - u = make_obj_array([lower_u, upper_u]) + u = make_obj_array([u1, u2]) def get_rhs(t, u): return make_obj_array([ diffusion_operator( - dcoll, kappa=1, boundaries=lower_boundaries, u=u[0], - volume_dd=dd_vol_lower), + dcoll, kappa=1, boundaries=boundaries1, u=u[0], + volume_dd=dd_vol1), diffusion_operator( - dcoll, kappa=1, boundaries=upper_boundaries, u=u[1], - volume_dd=dd_vol_upper)]) + dcoll, kappa=1, boundaries=boundaries2, u=u[1], + volume_dd=dd_vol2)]) rhs = get_rhs(0, u) if visualize: from grudge.shortcuts import make_visualizer - viz_lower = make_visualizer(dcoll, order+3, volume_dd=dd_vol_lower) - viz_upper = make_visualizer(dcoll, order+3, volume_dd=dd_vol_upper) - viz_lower.write_vtk_file( - f"multiphysics_independent_volumes_{order}_lower.vtu", [ + 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]), ]) - viz_upper.write_vtk_file( - f"multiphysics_independent_volumes_{order}_upper.vtu", [ + viz2.write_vtk_file( + f"multiphysics_independent_volumes_{order}_2.vtu", [ ("u", u[1]), ("rhs", rhs[1]), ]) - linf_err_lower = actx.to_numpy(op.norm(dcoll, rhs[0], np.inf, dd=dd_vol_lower)) - linf_err_upper = actx.to_numpy(op.norm(dcoll, rhs[1], np.inf, dd=dd_vol_upper)) + 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_err_lower < 1e-9 - assert linf_err_upper < 1e-9 + assert linf_err1 < 1e-9 + assert linf_err2 < 1e-9 if __name__ == "__main__": From 07405f93c46449de986320c7240e8f8298acf5d4 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Wed, 31 Aug 2022 14:48:56 -0500 Subject: [PATCH 08/19] add multiple independent volume example --- examples/multiple-volumes-mpi-lazy.py | 1 + examples/multiple-volumes-mpi.py | 418 ++++++++++++++++++++++++++ mirgecom/logging_quantities.py | 11 +- 3 files changed, 427 insertions(+), 3 deletions(-) create mode 120000 examples/multiple-volumes-mpi-lazy.py create mode 100644 examples/multiple-volumes-mpi.py 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..7d30d9476 --- /dev/null +++ b/examples/multiple-volumes-mpi.py @@ -0,0 +1,418 @@ +"""Demonstrate multiple independent volumes.""" + +__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 +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.euler import extract_vars_for_logging +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, + volume_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, + volume_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/mirgecom/logging_quantities.py b/mirgecom/logging_quantities.py index 06d719b7f..e22bc0913 100644 --- a/mirgecom/logging_quantities.py +++ b/mirgecom/logging_quantities.py @@ -105,20 +105,25 @@ 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, volume_dd=DD_VOLUME_ALL): """Add default discretization quantities to the logmgr.""" + if volume_dd != DD_VOLUME_ALL: + suffix = f"_{volume_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, volume_dd=volume_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, volume_dd=volume_dd)) for d in range(dim): logmgr.add_quantity(DiscretizationBasedQuantity( - dcoll, "momentum", reduction_op, extract_vars_for_logging, + dcoll, "momentum"+suffix, reduction_op, extract_vars_for_logging, units_for_logging, axis=d, volume_dd=volume_dd)) From a6940d2d2b27f4a553b67d76153845e4e9b0272f Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Wed, 7 Sep 2022 10:29:23 -0500 Subject: [PATCH 09/19] add missing btag -> bdtag promotion in euler_operator --- mirgecom/euler.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/mirgecom/euler.py b/mirgecom/euler.py index f99a87934..b8dcab15f 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -117,6 +117,10 @@ def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, comm_tag: Hashable Tag for distributed communication """ + boundaries = { + as_dofdesc(bdtag).domain_tag: bdry + for bdtag, bdry in boundaries.items()} + dd_quad_vol = volume_dd.with_discr_tag(quadrature_tag) dd_quad_allfaces = dd_quad_vol.trace(FACE_RESTR_ALL) From fd3c2c417c691592d64b9f694a07769bf2f04ba0 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Thu, 1 Sep 2022 16:34:12 -0500 Subject: [PATCH 10/19] unify dd naming --- examples/multiple-volumes-mpi.py | 4 +- mirgecom/artificial_viscosity.py | 58 ++++++++++++++------------ mirgecom/diffusion.py | 70 +++++++++++++++++++------------ mirgecom/euler.py | 32 ++++++++++---- mirgecom/filter.py | 21 ++++++---- mirgecom/gas_model.py | 32 +++++++++----- mirgecom/inviscid.py | 30 ++++++++++---- mirgecom/io.py | 6 +-- mirgecom/logging_quantities.py | 20 ++++----- mirgecom/navierstokes.py | 71 ++++++++++++++++++++------------ mirgecom/operators.py | 12 +++--- mirgecom/simutil.py | 7 +++- mirgecom/viscous.py | 59 ++++++++++++++++++-------- test/test_av.py | 7 ++-- test/test_bc.py | 21 +++++----- test/test_euler.py | 2 +- test/test_filter.py | 9 ++-- test/test_inviscid.py | 5 ++- test/test_lazy.py | 7 ++-- test/test_multiphysics.py | 6 +-- test/test_operators.py | 12 +++--- test/test_viscous.py | 12 +++--- 22 files changed, 307 insertions(+), 196 deletions(-) diff --git a/examples/multiple-volumes-mpi.py b/examples/multiple-volumes-mpi.py index 7d30d9476..a1c763736 100644 --- a/examples/multiple-volumes-mpi.py +++ b/examples/multiple-volumes-mpi.py @@ -187,7 +187,7 @@ def units(quantity): for i in range(nvolumes): logmgr_add_many_discretization_quantities( logmgr, dcoll, dim, partial(extract_vars, i), units, - volume_dd=volume_dds[i]) + dd=volume_dds[i]) vis_timer = IntervalTimer("t_vis", "Time spent visualizing") logmgr.add_quantity(vis_timer) @@ -349,7 +349,7 @@ def my_rhs(t, state): dcoll, state=fluid_state, time=t, boundaries={dd.trace(BTAG_ALL).domain_tag: wall}, gas_model=gas_model, quadrature_tag=quadrature_tag, - volume_dd=dd, comm_tag=dd) + dd=dd, comm_tag=dd) for dd, fluid_state in zip(volume_dds, fluid_states)]) current_dt = my_get_timestep( diff --git a/mirgecom/artificial_viscosity.py b/mirgecom/artificial_viscosity.py index 48e9ea457..dab7e02bc 100644 --- a/mirgecom/artificial_viscosity.py +++ b/mirgecom/artificial_viscosity.py @@ -142,6 +142,7 @@ from grudge.dof_desc import ( DD_VOLUME_ALL, + VolumeDomainTag, DISCR_TAG_BASE, DISCR_TAG_MODAL, as_dofdesc, @@ -156,9 +157,8 @@ class _AVRTag: def av_laplacian_operator(dcoll, boundaries, fluid_state, alpha, gas_model=None, kappa=1., s0=-6., time=0, quadrature_tag=DISCR_TAG_BASE, - volume_dd=DD_VOLUME_ALL, boundary_kwargs=None, - indicator=None, divergence_numerical_flux=num_flux_central, - comm_tag=None, + 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): @@ -201,8 +201,9 @@ def av_laplacian_operator(dcoll, boundaries, fluid_state, alpha, gas_model=None, An optional identifier denoting a particular quadrature discretization to use during operator evaluations. - volume_dd: grudge.dof_desc.DOFDesc - The DOF descriptor of the volume on which to apply the operator. + 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 @@ -219,9 +220,14 @@ def av_laplacian_operator(dcoll, boundaries, fluid_state, alpha, gas_model=None, cv = fluid_state.cv actx = cv.array_context - dd_base = volume_dd - dd_vol = dd_base.with_discr_tag(quadrature_tag) - dd_allfaces = dd_vol.trace(FACE_RESTR_ALL) + 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 @@ -236,13 +242,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, dd_base, 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, - volume_dd=dd_base, comm_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 @@ -250,27 +256,28 @@ 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, volume_dd=dd_base) + 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, - volume_dd=dd_base, + 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_trace = utpair.dd - dd_allfaces = dd_trace.with_boundary_tag(FACE_RESTR_ALL) - normal = actx.thaw(dcoll.normal(dd_trace)) - return op.project(dcoll, dd_trace, dd_allfaces, + 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 = ( @@ -278,23 +285,23 @@ def central_flux_div(utpair): sum( central_flux_div(interp_to_surf_quad(tpair=tpair)) for tpair in interior_trace_pairs( - dcoll, r, volume_dd=dd_base, comm_tag=(_AVRTag, comm_tag))) + dcoll, r, volume_dd=dd_vol, comm_tag=(_AVRTag, comm_tag))) # Contributions from boundary fluxes + sum( bdry.av_flux( dcoll, - dd_vol=dd_vol, 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_allfaces, 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, volume_dd=DD_VOLUME_ALL): +def smoothness_indicator(dcoll, u, kappa=1.0, s0=-6.0, dd=DD_VOLUME_ALL): r"""Calculate the smoothness indicator. Parameters @@ -321,7 +328,7 @@ def smoothness_indicator(dcoll, u, kappa=1.0, s0=-6.0, volume_dd=DD_VOLUME_ALL): actx = u.array_context - @memoize_in(actx, (smoothness_indicator, "smooth_comp_knl", volume_dd)) + @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 @@ -347,8 +354,7 @@ def indicator_prg(): return lp.tag_inames(t_unit, {"iel": ConcurrentElementInameTag(), "idof": ConcurrentDOFInameTag()}) - @keyed_memoize_in(actx, (smoothness_indicator, - "highest_mode", volume_dd), + @keyed_memoize_in(actx, (smoothness_indicator, "highest_mode", dd), lambda grp: grp.discretization_key()) def highest_mode(grp): return actx.from_numpy( @@ -358,7 +364,7 @@ def highest_mode(grp): ) # Convert to modal solution representation - dd_vol = volume_dd + 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) diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py index ca2d67f24..390154a95 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -40,7 +40,12 @@ from pytools.obj_array import make_obj_array, obj_array_vectorize_n_args from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from meshmode.discretization.connection import FACE_RESTR_ALL # noqa -from grudge.dof_desc import DD_VOLUME_ALL, DISCR_TAG_BASE, as_dofdesc +from grudge.dof_desc import ( + DD_VOLUME_ALL, + VolumeDomainTag, + DISCR_TAG_BASE, + as_dofdesc, +) from grudge.trace_pair import ( TracePair, interior_trace_pairs, @@ -206,8 +211,8 @@ class _DiffusionGradTag: def grad_operator( - dcoll, boundaries, u, *, quadrature_tag=DISCR_TAG_BASE, - volume_dd=DD_VOLUME_ALL, comm_tag=None): + dcoll, boundaries, u, *, quadrature_tag=DISCR_TAG_BASE, dd=DD_VOLUME_ALL, + comm_tag=None): r""" Compute the gradient of *u*. @@ -226,8 +231,9 @@ def grad_operator( quadrature_tag: quadrature tag indicating which discretization in *dcoll* to use for overintegration - volume_dd: grudge.dof_desc.DOFDesc - the DOF descriptor of the volume on which to apply the operator + 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 @@ -244,7 +250,7 @@ def grad_operator( return obj_array_vectorize_n_args( lambda boundaries_i, u_i: grad_operator( dcoll, boundaries_i, u_i, quadrature_tag=quadrature_tag, - volume_dd=volume_dd), + dd=dd), make_obj_array(boundaries), u) actx = u.array_context @@ -258,8 +264,13 @@ def grad_operator( raise TypeError(f"Unrecognized boundary type for tag {bdtag}. " "Must be an instance of DiffusionBoundary.") - dd_vol_base = volume_dd - dd_vol_quad = dd_vol_base.with_discr_tag(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) @@ -274,21 +285,21 @@ def interior_flux(u_tpair): def boundary_flux(bdtag, bdry): dd_bdry_quad = dd_vol_quad.with_domain_tag(bdtag) - u_minus_quad = op.project(dcoll, dd_vol_base, dd_bdry_quad, u) + u_minus_quad = op.project(dcoll, dd_vol, dd_bdry_quad, u) return op.project( dcoll, dd_bdry_quad, dd_allfaces_quad, bdry.get_grad_flux(dcoll, dd_bdry_quad, u_minus_quad)) return op.inverse_mass( - dcoll, dd_vol_base, - op.weak_local_grad(dcoll, dd_vol_base, -u) + dcoll, dd_vol, + op.weak_local_grad(dcoll, dd_vol, -u) - # noqa: W504 op.face_mass( dcoll, dd_allfaces_quad, sum( interior_flux(u_tpair) for u_tpair in interior_trace_pairs( - dcoll, u, volume_dd=dd_vol_base, + dcoll, u, volume_dd=dd_vol, comm_tag=(_DiffusionStateTag, comm_tag))) + sum( boundary_flux(bdtag, bdry) @@ -299,7 +310,7 @@ def boundary_flux(bdtag, bdry): def diffusion_operator( dcoll, kappa, boundaries, u, *, return_grad_u=False, - quadrature_tag=DISCR_TAG_BASE, volume_dd=DD_VOLUME_ALL, 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): @@ -329,8 +340,9 @@ def diffusion_operator( quadrature_tag: quadrature tag indicating which discretization in *dcoll* to use for overintegration - volume_dd: grudge.dof_desc.DOFDesc - the DOF descriptor of the volume on which to apply the operator + 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 @@ -349,7 +361,7 @@ 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, volume_dd=volume_dd), + quadrature_tag=quadrature_tag, dd=dd), make_obj_array(boundaries), u) actx = u.array_context @@ -363,17 +375,21 @@ def diffusion_operator( raise TypeError(f"Unrecognized boundary type for tag {bdtag}. " "Must be an instance of DiffusionBoundary.") - dd_vol_base = volume_dd - dd_vol_quad = dd_vol_base.with_discr_tag(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) if grad_u is None: grad_u = grad_operator( - dcoll, boundaries, u, quadrature_tag=quadrature_tag, - volume_dd=dd_vol_base) + dcoll, boundaries, u, quadrature_tag=quadrature_tag, dd=dd_vol) - kappa_quad = op.project(dcoll, dd_vol_base, dd_vol_quad, kappa) - grad_u_quad = op.project(dcoll, dd_vol_base, dd_vol_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) @@ -388,15 +404,15 @@ def interior_flux(kappa_tpair, grad_u_tpair): def boundary_flux(bdtag, bdry): dd_bdry_quad = dd_vol_quad.with_domain_tag(bdtag) - kappa_minus_quad = op.project(dcoll, dd_vol_base, dd_bdry_quad, kappa) - grad_u_minus_quad = op.project(dcoll, dd_vol_base, dd_bdry_quad, grad_u) + 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_quad, dd_allfaces_quad, bdry.get_diffusion_flux( dcoll, dd_bdry_quad, kappa_minus_quad, grad_u_minus_quad)) diff_u = op.inverse_mass( - dcoll, dd_vol_base, + dcoll, dd_vol, op.weak_local_div(dcoll, dd_vol_quad, -kappa_quad*grad_u_quad) - # noqa: W504 op.face_mass( @@ -405,10 +421,10 @@ def boundary_flux(bdtag, bdry): interior_flux(kappa_tpair, grad_u_tpair) for kappa_tpair, grad_u_tpair in zip( interior_trace_pairs( - dcoll, kappa, volume_dd=dd_vol_base, + dcoll, kappa, volume_dd=dd_vol, comm_tag=(_DiffusionKappaTag, comm_tag)), interior_trace_pairs( - dcoll, grad_u, volume_dd=dd_vol_base, + dcoll, grad_u, volume_dd=dd_vol, comm_tag=(_DiffusionGradTag, comm_tag))) ) + sum( diff --git a/mirgecom/euler.py b/mirgecom/euler.py index b8dcab15f..34395ad8f 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -55,7 +55,12 @@ import numpy as np # noqa from meshmode.discretization.connection import FACE_RESTR_ALL -from grudge.dof_desc import DD_VOLUME_ALL, DISCR_TAG_BASE +from grudge.dof_desc import ( + DD_VOLUME_ALL, + VolumeDomainTag, + DISCR_TAG_BASE, + as_dofdesc, +) from mirgecom.gas_model import make_operator_fluid_states from mirgecom.inviscid import ( @@ -69,7 +74,7 @@ def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, inviscid_numerical_flux_func=inviscid_facial_flux_rusanov, - quadrature_tag=DISCR_TAG_BASE, volume_dd=DD_VOLUME_ALL, + quadrature_tag=DISCR_TAG_BASE, dd=DD_VOLUME_ALL, comm_tag=None, operator_states_quad=None): r"""Compute RHS of the Euler flow equations. @@ -111,23 +116,32 @@ def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, An optional identifier denoting a particular quadrature discretization to use during operator evaluations. - volume_dd: grudge.dof_desc.DOFDesc - The DOF descriptor of the volume on which to apply the operator. + 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 """ boundaries = { as_dofdesc(bdtag).domain_tag: bdry for bdtag, bdry in boundaries.items()} - dd_quad_vol = volume_dd.with_discr_tag(quadrature_tag) - dd_quad_allfaces = dd_quad_vol.trace(FACE_RESTR_ALL) + 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, - volume_dd=volume_dd, comm_tag=comm_tag) + dd=dd_vol, comm_tag=comm_tag) volume_state_quad, interior_state_pairs_quad, domain_boundary_states_quad = \ operator_states_quad @@ -140,9 +154,9 @@ def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, 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, - volume_dd=volume_dd) + dd=dd_vol) - return -div_operator(dcoll, dd_quad_vol, dd_quad_allfaces, + 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 cf0a1e2b0..f2f12235f 100644 --- a/mirgecom/filter.py +++ b/mirgecom/filter.py @@ -48,7 +48,11 @@ import numpy as np from functools import partial -from grudge.dof_desc import DISCR_TAG_MODAL, as_dofdesc +from grudge.dof_desc import ( + DD_VOLUME_ALL, + DISCR_TAG_BASE, + DISCR_TAG_MODAL, +) from arraycontext import map_array_container @@ -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,10 +206,13 @@ 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 ) - dd_nodal = as_dofdesc(dd) + 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) diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py index 9d622a984..8d8f0f48a 100644 --- a/mirgecom/gas_model.py +++ b/mirgecom/gas_model.py @@ -59,7 +59,11 @@ TransportModel, GasTransportVars ) -from grudge.dof_desc import DD_VOLUME_ALL, DISCR_TAG_BASE +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, @@ -386,7 +390,7 @@ class _FluidTemperatureTag: def make_operator_fluid_states( dcoll, volume_state, gas_model, boundaries, quadrature_tag=DISCR_TAG_BASE, - volume_dd=DD_VOLUME_ALL, comm_tag=None): + 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 @@ -422,8 +426,9 @@ def make_operator_fluid_states( An identifier denoting a particular quadrature discretization to use during operator evaluations. - volume_dd: grudge.dof_desc.DOFDesc - the DOF descriptor of the volume on which to construct the states + 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 @@ -438,15 +443,20 @@ def make_operator_fluid_states( boundary domain tags in *boundaries*, all on the quadrature grid (if specified). """ - dd_base_vol = volume_dd - dd_quad_vol = dd_base_vol.with_discr_tag(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) # 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 = { - bdtag: project_fluid_state(dcoll, dd_base_vol, - dd_quad_vol.with_domain_tag(bdtag), + bdtag: project_fluid_state(dcoll, dd_vol, + dd_vol_quad.with_domain_tag(bdtag), volume_state, gas_model) for bdtag in boundaries } @@ -457,7 +467,7 @@ def make_operator_fluid_states( # discretization (if any) interp_to_surf_quad(tpair=tpair) for tpair in interior_trace_pairs( - dcoll, volume_state.cv, volume_dd=dd_base_vol, + dcoll, volume_state.cv, volume_dd=dd_vol, comm_tag=(_FluidCVTag, comm_tag)) ] @@ -472,7 +482,7 @@ def make_operator_fluid_states( # discretization (if any) interp_to_surf_quad(tpair=tpair) for tpair in interior_trace_pairs( - dcoll, volume_state.temperature, volume_dd=dd_base_vol, + dcoll, volume_state.temperature, volume_dd=dd_vol, comm_tag=(_FluidTemperatureTag, comm_tag))] interior_boundary_states_quad = \ @@ -481,7 +491,7 @@ def make_operator_fluid_states( # 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 87250730a..ecd2a72f0 100644 --- a/mirgecom/inviscid.py +++ b/mirgecom/inviscid.py @@ -41,7 +41,11 @@ import numpy as np from meshmode.discretization.connection import FACE_RESTR_ALL -from grudge.dof_desc import DD_VOLUME_ALL, DISCR_TAG_BASE +from grudge.dof_desc import ( + DD_VOLUME_ALL, + VolumeDomainTag, + DISCR_TAG_BASE, +) import grudge.op as op from mirgecom.fluid import make_conserved @@ -228,7 +232,7 @@ def inviscid_flux_on_element_boundary( dcoll, gas_model, boundaries, interior_state_pairs, domain_boundary_states, quadrature_tag=DISCR_TAG_BASE, numerical_flux_func=inviscid_facial_flux_rusanov, time=0.0, - volume_dd=DD_VOLUME_ALL): + dd=DD_VOLUME_ALL): """Compute the inviscid boundary fluxes for the divergence operator. This routine encapsulates the computation of the inviscid contributions @@ -265,10 +269,17 @@ def inviscid_flux_on_element_boundary( time: float Time - volume_dd: grudge.dof_desc.DOFDesc - The DOF descriptor of the volume on which to compute the flux. + dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the discretization on which the fluid lives. Must be + a volume on the base discretization. """ - dd_vol_quad = volume_dd.with_discr_tag(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) def _interior_flux(state_pair): @@ -278,11 +289,12 @@ def _interior_flux(state_pair): 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_allfaces_quad, + 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 @@ -294,7 +306,7 @@ def _boundary_flux(dd_bdry, boundary, state_minus): # Domain boundary faces + sum( _boundary_flux( - dd_vol_quad.with_domain_tag(bdtag), + bdtag, boundary, domain_boundary_states[bdtag]) for bdtag, boundary in boundaries.items()) diff --git a/mirgecom/io.py b/mirgecom/io.py index 88e64b71b..0adc69b9c 100644 --- a/mirgecom/io.py +++ b/mirgecom/io.py @@ -55,11 +55,11 @@ def make_init_message(*, dim, order, dt, t_final, def make_status_message( - *, dcoll, t, step, dt, cfl, dependent_vars, fluid_volume_dd=DD_VOLUME_ALL): + *, 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, fluid_volume_dd) - _max = partial(op.nodal_max, dcoll, fluid_volume_dd) + _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/logging_quantities.py b/mirgecom/logging_quantities.py index e22bc0913..2ea3974f8 100644 --- a/mirgecom/logging_quantities.py +++ b/mirgecom/logging_quantities.py @@ -103,10 +103,10 @@ 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, volume_dd=DD_VOLUME_ALL): + extract_vars_for_logging, units_for_logging, dd=DD_VOLUME_ALL): """Add default discretization quantities to the logmgr.""" - if volume_dd != DD_VOLUME_ALL: - suffix = f"_{volume_dd.domain_tag.tag}" + if dd != DD_VOLUME_ALL: + suffix = f"_{dd.domain_tag.tag}" else: suffix = "" @@ -114,17 +114,17 @@ def logmgr_add_many_discretization_quantities(logmgr: LogManager, dcoll, dim, for quantity in ["pressure"+suffix, "temperature"+suffix]: logmgr.add_quantity(DiscretizationBasedQuantity( dcoll, quantity, reduction_op, extract_vars_for_logging, - units_for_logging, volume_dd=volume_dd)) + units_for_logging, dd=dd)) for quantity in ["mass"+suffix, "energy"+suffix]: logmgr.add_quantity(DiscretizationBasedQuantity( dcoll, quantity, reduction_op, extract_vars_for_logging, - units_for_logging, volume_dd=volume_dd)) + units_for_logging, dd=dd)) for d in range(dim): logmgr.add_quantity(DiscretizationBasedQuantity( dcoll, "momentum"+suffix, reduction_op, extract_vars_for_logging, - units_for_logging, axis=d, volume_dd=volume_dd)) + units_for_logging, axis=d, dd=dd)) # {{{ Package versions @@ -250,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, volume_dd=DD_VOLUME_ALL): + axis: Optional[int] = None, dd=DD_VOLUME_ALL): unit = units_logging(quantity) if name is None: @@ -267,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, volume_dd) + 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, volume_dd) + 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, dd=volume_dd) + 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 ab576e8f0..deb9fb55c 100644 --- a/mirgecom/navierstokes.py +++ b/mirgecom/navierstokes.py @@ -68,6 +68,7 @@ ) from grudge.dof_desc import ( DD_VOLUME_ALL, + VolumeDomainTag, DISCR_TAG_BASE, as_dofdesc, ) @@ -114,7 +115,7 @@ def _gradient_flux_interior(dcoll, numerical_flux_func, tpair): def grad_cv_operator( dcoll, gas_model, boundaries, state, *, time=0.0, numerical_flux_func=num_flux_central, - quadrature_tag=DISCR_TAG_BASE, volume_dd=DD_VOLUME_ALL, 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): @@ -148,8 +149,9 @@ def grad_cv_operator( An identifier denoting a particular quadrature discretization to use during operator evaluations. - volume_dd: grudge.dof_desc.DOFDesc - The DOF descriptor of the volume on which to apply the operator. + 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 @@ -165,14 +167,19 @@ def grad_cv_operator( as_dofdesc(bdtag).domain_tag: bdry for bdtag, bdry in boundaries.items()} - dd_vol_base = volume_dd - dd_vol_quad = dd_vol_base.with_discr_tag(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) if operator_states_quad is None: operator_states_quad = make_operator_fluid_states( dcoll, state, gas_model, boundaries, quadrature_tag, - volume_dd=dd_vol_base, comm_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 @@ -212,7 +219,7 @@ def grad_cv_operator( def grad_t_operator( dcoll, gas_model, boundaries, state, *, time=0.0, numerical_flux_func=num_flux_central, - quadrature_tag=DISCR_TAG_BASE, volume_dd=DD_VOLUME_ALL, 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): @@ -245,8 +252,9 @@ def grad_t_operator( An identifier denoting a particular quadrature discretization to use during operator evaluations. - volume_dd: grudge.dof_desc.DOFDesc - The DOF descriptor of the volume on which to apply the operator. + 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 @@ -262,14 +270,19 @@ def grad_t_operator( as_dofdesc(bdtag).domain_tag: bdry for bdtag, bdry in boundaries.items()} - dd_vol_base = volume_dd - dd_vol_quad = dd_vol_base.with_discr_tag(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) if operator_states_quad is None: operator_states_quad = make_operator_fluid_states( dcoll, state, gas_model, boundaries, quadrature_tag, - volume_dd=dd_vol_base, comm_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 @@ -315,7 +328,7 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, gradient_numerical_flux_func=num_flux_central, viscous_numerical_flux_func=viscous_facial_flux_central, return_gradients=False, quadrature_tag=DISCR_TAG_BASE, - volume_dd=DD_VOLUME_ALL, comm_tag=None, + 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, @@ -363,8 +376,9 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, An identifier denoting a particular quadrature discretization to use during operator evaluations. - volume_dd: grudge.dof_desc.DOFDesc - The DOF descriptor of the volume on which to apply the operator. + 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 @@ -404,8 +418,13 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, as_dofdesc(bdtag).domain_tag: bdry for bdtag, bdry in boundaries.items()} - dd_vol_base = volume_dd - dd_vol_quad = dd_vol_base.with_discr_tag(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) # Make model-consistent fluid state data (i.e. CV *and* DV) for: @@ -418,7 +437,7 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, if operator_states_quad is None: operator_states_quad = make_operator_fluid_states( dcoll, state, gas_model, boundaries, quadrature_tag, - volume_dd=dd_vol_base, comm_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 @@ -436,7 +455,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, volume_dd=dd_vol_base, + quadrature_tag=quadrature_tag, dd=dd_vol, operator_states_quad=operator_states_quad) # Communicate grad(CV) and put it on the quadrature domain @@ -445,7 +464,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, volume_dd=dd_vol_base, comm_tag=(_NSGradCVTag, comm_tag)) + dcoll, grad_cv, volume_dd=dd_vol, comm_tag=(_NSGradCVTag, comm_tag)) ] # }}} Compute grad(CV) @@ -456,7 +475,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, volume_dd=dd_vol_base, + 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 @@ -465,7 +484,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_t, volume_dd=dd_vol_base, + dcoll, grad_t, volume_dd=dd_vol, comm_tag=(_NSGradTemperatureTag, comm_tag)) ] @@ -480,8 +499,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_vol_base, dd_vol_quad, grad_cv), - grad_t=op.project(dcoll, dd_vol_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 @@ -497,14 +516,14 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, 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, - volume_dd=dd_vol_base) + 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, - volume_dd=dd_vol_base) + dd=dd_vol) ) ns_rhs = div_operator(dcoll, dd_vol_quad, dd_allfaces_quad, vol_term, bnd_term) diff --git a/mirgecom/operators.py b/mirgecom/operators.py index 7a35ea6b7..3984f1656 100644 --- a/mirgecom/operators.py +++ b/mirgecom/operators.py @@ -33,7 +33,7 @@ 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 @@ -43,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 @@ -62,10 +62,10 @@ def grad_operator(dcoll, dd_vol, dd_faces, u, flux): 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 @@ -75,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 @@ -94,4 +94,4 @@ def div_operator(dcoll, dd_vol, dd_faces, v, flux): 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 34373c267..e2e3943f5 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -94,7 +94,7 @@ def check_step(step, interval): def get_sim_timestep( dcoll, state, t, dt, cfl, t_final, constant_cfl=False, - fluid_volume_dd=DD_VOLUME_ALL): + 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 @@ -125,6 +125,9 @@ def get_sim_timestep( 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 ------- @@ -138,7 +141,7 @@ def get_sim_timestep( from grudge.op import nodal_min mydt = state.array_context.to_numpy( cfl * nodal_min( - dcoll, fluid_volume_dd, + dcoll, fluid_dd, get_viscous_timestep(dcoll=dcoll, state=state)))[()] return min(t_remaining, mydt) diff --git a/mirgecom/viscous.py b/mirgecom/viscous.py index 65988999e..9336faf0a 100644 --- a/mirgecom/viscous.py +++ b/mirgecom/viscous.py @@ -46,7 +46,11 @@ 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, DISCR_TAG_BASE +from grudge.dof_desc import ( + DD_VOLUME_ALL, + VolumeDomainTag, + DISCR_TAG_BASE, +) import grudge.op as op @@ -339,7 +343,7 @@ def viscous_flux_on_element_boundary( domain_boundary_states, grad_cv, interior_grad_cv_pairs, grad_t, interior_grad_t_pairs, quadrature_tag=DISCR_TAG_BASE, numerical_flux_func=viscous_facial_flux_central, time=0.0, - volume_dd=DD_VOLUME_ALL): + dd=DD_VOLUME_ALL): """Compute the viscous boundary fluxes for the divergence operator. This routine encapsulates the computation of the viscous contributions @@ -387,11 +391,17 @@ def viscous_flux_on_element_boundary( time: float Time - volume_dd: grudge.dof_desc.DOFDesc - The DOF descriptor of the volume on which to compute the flux. + dd: grudge.dof_desc.DOFDesc + the DOF descriptor of the discretization on which the fluid lives. Must be + a volume on the base discretization. """ - dd_base = volume_dd - dd_vol_quad = dd_base.with_discr_tag(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) # {{{ - Viscous flux helpers - @@ -405,15 +415,15 @@ def _fvisc_divergence_flux_interior(state_pair, grad_cv_pair, grad_t_pair): grad_cv_pair=grad_cv_pair, grad_t_pair=grad_t_pair)) # viscous part of bcs applied here - def _fvisc_divergence_flux_boundary(bdtag, boundary, state_minus): - dd_bdry = dd_vol_quad.with_domain_tag(bdtag) + 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_bdry, dd_allfaces_quad, + dcoll, dd_bdry_quad, dd_allfaces_quad, boundary.viscous_divergence_flux( - dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model, - state_minus=state_minus, - grad_cv_minus=op.project(dcoll, dd_base, dd_bdry, grad_cv), - grad_t_minus=op.project(dcoll, dd_base, dd_bdry, 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 @@ -442,7 +452,7 @@ def _fvisc_divergence_flux_boundary(bdtag, boundary, state_minus): return bnd_term -def get_viscous_timestep(dcoll, state, *, volume_dd=DD_VOLUME_ALL): +def get_viscous_timestep(dcoll, state, *, dd=DD_VOLUME_ALL): """Routine returns the the node-local maximum stable viscous timestep. Parameters @@ -455,16 +465,26 @@ def get_viscous_timestep(dcoll, state, *, volume_dd=DD_VOLUME_ALL): 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, dd=volume_dd) + state.array_context, dcoll, dd=dd) nu = 0 d_alpha_max = 0 @@ -482,7 +502,7 @@ def get_viscous_timestep(dcoll, state, *, volume_dd=DD_VOLUME_ALL): ) -def get_viscous_cfl(dcoll, dt, state, *, volume_dd=DD_VOLUME_ALL): +def get_viscous_cfl(dcoll, dt, state, *, dd=DD_VOLUME_ALL): """Calculate and return node-local CFL based on current state and timestep. Parameters @@ -499,13 +519,18 @@ def get_viscous_cfl(dcoll, dt, state, *, volume_dd=DD_VOLUME_ALL): 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, volume_dd=volume_dd) + return dt / get_viscous_timestep(dcoll, state=state, dd=dd) def get_local_max_species_diffusivity(actx, d_alpha): diff --git a/test/test_av.py b/test/test_av.py index 216b110b2..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, @@ -443,11 +444,11 @@ def _boundary_state_func(dcoll, dd_bdry, 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 21c4d1a88..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 @@ -355,7 +356,7 @@ 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") + 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) @@ -368,13 +369,13 @@ 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=}") @@ -564,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 @@ -577,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) diff --git a/test/test_euler.py b/test/test_euler.py index cf0dbb0da..cefbad856 100644 --- a/test/test_euler.py +++ b/test/test_euler.py @@ -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 diff --git a/test/test_filter.py b/test/test_filter.py index bcd587f9c..edb1143e1 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) @@ -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}") @@ -237,8 +235,7 @@ def polyfn(coeff): # , x_vec): 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_inviscid.py b/test/test_inviscid.py index 8aa42bbe6..c78b6d79b 100644 --- a/test/test_inviscid.py +++ b/test/test_inviscid.py @@ -37,6 +37,7 @@ ) 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 @@ -334,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): @@ -383,7 +384,7 @@ def inf_norm(data): 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 9c41e37b9..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 @@ -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) diff --git a/test/test_multiphysics.py b/test/test_multiphysics.py index f440b69f3..78742397e 100644 --- a/test/test_multiphysics.py +++ b/test/test_multiphysics.py @@ -96,11 +96,9 @@ def test_independent_volumes(actx_factory, order, visualize=False): def get_rhs(t, u): return make_obj_array([ diffusion_operator( - dcoll, kappa=1, boundaries=boundaries1, u=u[0], - volume_dd=dd_vol1), + dcoll, kappa=1, boundaries=boundaries1, u=u[0], dd=dd_vol1), diffusion_operator( - dcoll, kappa=1, boundaries=boundaries2, u=u[1], - volume_dd=dd_vol2)]) + dcoll, kappa=1, boundaries=boundaries2, u=u[1], dd=dd_vol2)]) rhs = get_rhs(0, u) diff --git a/test/test_operators.py b/test/test_operators.py index 749fea536..ab6d6e865 100644 --- a/test/test_operators.py +++ b/test/test_operators.py @@ -126,8 +126,8 @@ 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, dd_bdry): @@ -140,8 +140,8 @@ def central_flux_boundary(actx, dcoll, soln_func, dd_bdry): 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]) @@ -228,8 +228,8 @@ def sym_eval(expr, x_vec): from mirgecom.operators import grad_operator 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 198164ab1..1b156b3e9 100644 --- a/test/test_viscous.py +++ b/test/test_viscous.py @@ -173,8 +173,8 @@ 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(dd_bdry): boundary_discr = dcoll.discr_from_dd(dd_bdry) @@ -185,8 +185,8 @@ def cv_flux_boundary(dd_bdry): 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 = dd_bdry.with_boundary_tag(FACE_RESTR_ALL) - return op.project(dcoll, dd_bdry, 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]: @@ -217,8 +217,8 @@ def inf_norm(x): cv_int_tpairs, boundaries) 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) + 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 From ef3863eb4a73f55e9e48ddb802e4e87c5ac77d7a Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Tue, 6 Sep 2022 16:23:26 -0500 Subject: [PATCH 11/19] add missing dd argument to get_inviscid_timestep --- mirgecom/inviscid.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py index ecd2a72f0..6501ea178 100644 --- a/mirgecom/inviscid.py +++ b/mirgecom/inviscid.py @@ -315,7 +315,7 @@ def _boundary_flux(bdtag, boundary, state_minus_quad): 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. @@ -330,15 +330,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 ) From 4aa83f1c51899c43b47b2417398d8de68bd39534 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Wed, 7 Sep 2022 10:26:08 -0500 Subject: [PATCH 12/19] remove deprecated += --- test/test_filter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_filter.py b/test/test_filter.py index edb1143e1..6e3484b65 100644 --- a/test/test_filter.py +++ b/test/test_filter.py @@ -205,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 From 8c7de8a2b3a08a13c84a23fd05db720defdbe75f Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Wed, 7 Sep 2022 13:03:11 -0500 Subject: [PATCH 13/19] add missing btag -> bdtag promotion in make_operator_fluid_states --- mirgecom/gas_model.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py index 8d8f0f48a..555d65903 100644 --- a/mirgecom/gas_model.py +++ b/mirgecom/gas_model.py @@ -63,6 +63,7 @@ DD_VOLUME_ALL, VolumeDomainTag, DISCR_TAG_BASE, + as_dofdesc, ) import grudge.op as op from grudge.trace_pair import ( @@ -443,6 +444,10 @@ def make_operator_fluid_states( boundary domain tags in *boundaries*, all on the quadrature grid (if specified). """ + boundaries = { + as_dofdesc(bdtag).domain_tag: bdry + for bdtag, bdry in boundaries.items()} + if not isinstance(dd.domain_tag, VolumeDomainTag): raise TypeError("dd must represent a volume") if dd.discretization_tag != DISCR_TAG_BASE: From 1d235eef48b622fbab1f3e314a2bdde523fc5c97 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Wed, 7 Sep 2022 13:20:13 -0500 Subject: [PATCH 14/19] wrap btag -> bdtag promotion into a function --- mirgecom/artificial_viscosity.py | 7 +++---- mirgecom/diffusion.py | 10 +++------- mirgecom/euler.py | 6 ++---- mirgecom/gas_model.py | 6 ++---- mirgecom/inviscid.py | 3 +++ mirgecom/navierstokes.py | 14 ++++---------- mirgecom/utils.py | 13 +++++++++++++ mirgecom/viscous.py | 4 ++++ 8 files changed, 34 insertions(+), 29 deletions(-) diff --git a/mirgecom/artificial_viscosity.py b/mirgecom/artificial_viscosity.py index dab7e02bc..346bc6c6e 100644 --- a/mirgecom/artificial_viscosity.py +++ b/mirgecom/artificial_viscosity.py @@ -145,9 +145,10 @@ VolumeDomainTag, DISCR_TAG_BASE, DISCR_TAG_MODAL, - as_dofdesc, ) +from mirgecom.utils import normalize_boundaries + import grudge.op as op @@ -213,9 +214,7 @@ def av_laplacian_operator(dcoll, boundaries, fluid_state, alpha, gas_model=None, :class:`mirgecom.fluid.ConservedVars` The artificial viscosity operator applied to *q*. """ - boundaries = { - as_dofdesc(bdtag).domain_tag: bdry - for bdtag, bdry in boundaries.items()} + boundaries = normalize_boundaries(boundaries) cv = fluid_state.cv actx = cv.array_context diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py index 390154a95..6726e85e7 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -44,7 +44,6 @@ DD_VOLUME_ALL, VolumeDomainTag, DISCR_TAG_BASE, - as_dofdesc, ) from grudge.trace_pair import ( TracePair, @@ -52,6 +51,7 @@ tracepair_with_discr_tag, ) import grudge.op as op +from mirgecom.utils import normalize_boundaries def grad_facial_flux(u_tpair, normal): @@ -255,9 +255,7 @@ def grad_operator( actx = u.array_context - boundaries = { - as_dofdesc(bdtag).domain_tag: bdry - for bdtag, bdry in boundaries.items()} + boundaries = normalize_boundaries(boundaries) for bdtag, bdry in boundaries.items(): if not isinstance(bdry, DiffusionBoundary): @@ -366,9 +364,7 @@ def diffusion_operator( actx = u.array_context - boundaries = { - as_dofdesc(bdtag).domain_tag: bdry - for bdtag, bdry in boundaries.items()} + boundaries = normalize_boundaries(boundaries) for bdtag, bdry in boundaries.items(): if not isinstance(bdry, DiffusionBoundary): diff --git a/mirgecom/euler.py b/mirgecom/euler.py index 34395ad8f..e583367a5 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -59,7 +59,6 @@ DD_VOLUME_ALL, VolumeDomainTag, DISCR_TAG_BASE, - as_dofdesc, ) from mirgecom.gas_model import make_operator_fluid_states @@ -70,6 +69,7 @@ ) from mirgecom.operators import div_operator +from mirgecom.utils import normalize_boundaries def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, @@ -125,9 +125,7 @@ def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, Tag for distributed communication """ - boundaries = { - as_dofdesc(bdtag).domain_tag: bdry - for bdtag, bdry in boundaries.items()} + boundaries = normalize_boundaries(boundaries) if not isinstance(dd.domain_tag, VolumeDomainTag): raise TypeError("dd must represent a volume") diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py index 555d65903..f11e7aa14 100644 --- a/mirgecom/gas_model.py +++ b/mirgecom/gas_model.py @@ -63,13 +63,13 @@ DD_VOLUME_ALL, VolumeDomainTag, DISCR_TAG_BASE, - as_dofdesc, ) 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) @@ -444,9 +444,7 @@ def make_operator_fluid_states( boundary domain tags in *boundaries*, all on the quadrature grid (if specified). """ - boundaries = { - as_dofdesc(bdtag).domain_tag: bdry - for bdtag, bdry in boundaries.items()} + boundaries = normalize_boundaries(boundaries) if not isinstance(dd.domain_tag, VolumeDomainTag): raise TypeError("dd must represent a volume") diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py index 6501ea178..d70228a87 100644 --- a/mirgecom/inviscid.py +++ b/mirgecom/inviscid.py @@ -48,6 +48,7 @@ ) import grudge.op as op from mirgecom.fluid import make_conserved +from mirgecom.utils import normalize_boundaries def inviscid_flux(state): @@ -273,6 +274,8 @@ def inviscid_flux_on_element_boundary( the DOF descriptor of the discretization on which the fluid lives. Must be a volume on the base discretization. """ + 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: diff --git a/mirgecom/navierstokes.py b/mirgecom/navierstokes.py index deb9fb55c..91a01259a 100644 --- a/mirgecom/navierstokes.py +++ b/mirgecom/navierstokes.py @@ -70,7 +70,6 @@ DD_VOLUME_ALL, VolumeDomainTag, DISCR_TAG_BASE, - as_dofdesc, ) import grudge.op as op @@ -91,6 +90,7 @@ div_operator, grad_operator ) from mirgecom.gas_model import make_operator_fluid_states +from mirgecom.utils import normalize_boundaries class _NSGradCVTag: @@ -163,9 +163,7 @@ def grad_cv_operator( CV object with vector components representing the gradient of the fluid conserved variables. """ - boundaries = { - as_dofdesc(bdtag).domain_tag: bdry - for bdtag, bdry in boundaries.items()} + boundaries = normalize_boundaries(boundaries) if not isinstance(dd.domain_tag, VolumeDomainTag): raise TypeError("dd must represent a volume") @@ -266,9 +264,7 @@ def grad_t_operator( Array of :class:`~meshmode.dof_array.DOFArray` representing the gradient of the fluid temperature. """ - boundaries = { - as_dofdesc(bdtag).domain_tag: bdry - for bdtag, bdry in boundaries.items()} + boundaries = normalize_boundaries(boundaries) if not isinstance(dd.domain_tag, VolumeDomainTag): raise TypeError("dd must represent a volume") @@ -414,9 +410,7 @@ 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.") - boundaries = { - as_dofdesc(bdtag).domain_tag: bdry - for bdtag, bdry in boundaries.items()} + boundaries = normalize_boundaries(boundaries) if not isinstance(dd.domain_tag, VolumeDomainTag): raise TypeError("dd must represent a volume") 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 9336faf0a..1790fd41f 100644 --- a/mirgecom/viscous.py +++ b/mirgecom/viscous.py @@ -60,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 @@ -395,6 +397,8 @@ def viscous_flux_on_element_boundary( the DOF descriptor of the discretization on which the fluid lives. Must be a volume on the base discretization. """ + 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: From 3dfca1bf3d46a5c0f2fbc0d1ff340abe90266b6f Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Fri, 16 Sep 2022 10:01:23 -0500 Subject: [PATCH 15/19] fix doc build issue --- doc/conf.py | 5 +++++ 1 file changed, 5 insertions(+) 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.*"), +] From bbfbd77d1b5d9cd38f4b24b831df0ebc04f60d3a Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Fri, 16 Sep 2022 12:11:41 -0500 Subject: [PATCH 16/19] clarify what multiple-volumes example does --- examples/multiple-volumes-mpi.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/examples/multiple-volumes-mpi.py b/examples/multiple-volumes-mpi.py index a1c763736..80a59e40e 100644 --- a/examples/multiple-volumes-mpi.py +++ b/examples/multiple-volumes-mpi.py @@ -1,4 +1,9 @@ -"""Demonstrate multiple independent volumes.""" +""" +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 From 17090ab1f6e9bdebe4b1915e0a4874e2e109306b Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Fri, 16 Sep 2022 12:29:34 -0500 Subject: [PATCH 17/19] clarify solution setup in test_independent_volumes --- test/test_multiphysics.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/test_multiphysics.py b/test/test_multiphysics.py index 78742397e..7805b2481 100644 --- a/test/test_multiphysics.py +++ b/test/test_multiphysics.py @@ -74,6 +74,9 @@ def test_independent_volumes(actx_factory, order, visualize=False): 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.), From ab0ecc69123a67813865d433be6d1fd1976ffce9 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Fri, 16 Sep 2022 12:24:37 -0500 Subject: [PATCH 18/19] group euler imports --- examples/multiple-volumes-mpi.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/multiple-volumes-mpi.py b/examples/multiple-volumes-mpi.py index 80a59e40e..67d466d80 100644 --- a/examples/multiple-volumes-mpi.py +++ b/examples/multiple-volumes-mpi.py @@ -41,7 +41,10 @@ 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 +from mirgecom.euler import ( + euler_operator, + extract_vars_for_logging +) from mirgecom.simutil import ( get_sim_timestep, generate_and_distribute_mesh @@ -61,7 +64,6 @@ make_fluid_state ) from logpyle import IntervalTimer, set_dt -from mirgecom.euler import extract_vars_for_logging from mirgecom.logging_quantities import ( initialize_logmgr, logmgr_add_many_discretization_quantities, From 0280f0dca167cab7c7d0ea6a8b14341ead73ef0f Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Mon, 19 Sep 2022 16:53:56 -0500 Subject: [PATCH 19/19] set grudge branch --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index bf9dc1a29..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/majosm/grudge.git@dof-desc-helpers#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