diff --git a/examples/thermally-coupled-mpi.py b/examples/thermally-coupled-mpi.py index f9d2cf2c5..a51eab839 100644 --- a/examples/thermally-coupled-mpi.py +++ b/examples/thermally-coupled-mpi.py @@ -41,9 +41,7 @@ DISCR_TAG_QUAD, DOFDesc, ) -from mirgecom.diffusion import ( - NeumannDiffusionBoundary, -) +from mirgecom.diffusion import NeumannDiffusionBoundary from mirgecom.discretization import create_discretization_collection from mirgecom.simutil import ( get_sim_timestep, @@ -60,7 +58,7 @@ from mirgecom.fluid import make_conserved from mirgecom.gas_model import ( GasModel, - make_fluid_state + make_fluid_state, ) from logpyle import IntervalTimer, set_dt from mirgecom.euler import extract_vars_for_logging @@ -71,9 +69,8 @@ logmgr_add_device_memory_usage, set_sim_state ) - from mirgecom.multiphysics.thermally_coupled_fluid_wall import ( - coupled_ns_heat_operator, + basic_coupled_ns_heat_operator as coupled_ns_heat_operator, ) logger = logging.getLogger(__name__) @@ -499,13 +496,13 @@ def my_post_step(step, t, dt, state): def my_rhs(t, state, return_gradients=False): fluid_state = make_fluid_state(cv=state[0], gas_model=gas_model) wall_temperature = state[1] + ns_heat_result = coupled_ns_heat_operator( dcoll, gas_model, dd_vol_fluid, dd_vol_wall, fluid_boundaries, wall_boundaries, - fluid_state, - wall_kappa, wall_temperature, + fluid_state, wall_kappa, wall_temperature, time=t, return_gradients=return_gradients, quadrature_tag=quadrature_tag) diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py index cfcc8d37e..8be93f6ab 100644 --- a/mirgecom/boundary.py +++ b/mirgecom/boundary.py @@ -19,6 +19,7 @@ .. autoclass:: RiemannInflowBoundary .. autoclass:: RiemannOutflowBoundary .. autoclass:: PressureOutflowBoundary +.. autoclass:: IsothermalSlipWallBoundary .. autoclass:: IsothermalWallBoundary .. autoclass:: AdiabaticSlipBoundary .. autoclass:: AdiabaticNoslipWallBoundary @@ -61,6 +62,7 @@ from mirgecom.flux import num_flux_central from mirgecom.gas_model import make_fluid_state, replace_fluid_state from pytools.obj_array import make_obj_array +from mirgecom.utils import project_from_base from mirgecom.inviscid import inviscid_facial_flux_rusanov @@ -1045,7 +1047,7 @@ def __init__(self): def state_plus( self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): - """Return state with zero normal-component velocity for an adiabatic wall.""" + """Return state with reflected normal-component velocity.""" actx = state_minus.array_context # Grab a unit normal to the boundary @@ -1057,7 +1059,7 @@ def state_plus( def state_bc( self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): - """Return state with zero normal-component velocity for an adiabatic wall.""" + """Return state with zero normal-component velocity.""" actx = state_minus.array_context # Grab a unit normal to the boundary @@ -1613,6 +1615,102 @@ def outflow_state(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): smoothness_beta=state_minus.smoothness_beta) +class IsothermalSlipWallBoundary(MengaldoBoundaryCondition): + r"""Isothermal viscous slip wall boundary. + + This class implements an isothermal slip wall consistent with the prescriptions + in [Mengaldo_2014]_. + + .. automethod:: __init__ + + .. automethod:: state_bc + .. automethod:: temperature_bc + .. automethod:: grad_cv_bc + .. automethod:: grad_temperature_bc + .. automethod:: state_plus + """ + + def __init__(self, wall_temperature=300): + """Initialize the boundary condition object.""" + self._wall_temp = wall_temperature + self._slip = _SlipBoundaryComponent() + self._impermeable = _ImpermeableBoundaryComponent() + + def temperature_bc(self, dcoll, dd_bdry, state_minus, **kwargs): + """Get temperature value used in grad(T).""" + wall_temp = project_from_base(dcoll, dd_bdry, self._wall_temp) + return 0*state_minus.temperature + wall_temp + + def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): + """Return BC fluid state.""" + dd_bdry = as_dofdesc(dd_bdry) + + nhat = state_minus.array_context.thaw(dcoll.normal(dd_bdry)) + + mom_bc = self._slip.momentum_bc(state_minus.momentum_density, nhat) + t_bc = self.temperature_bc(dcoll, dd_bdry, state_minus, **kwargs) + + internal_energy_bc = gas_model.eos.get_internal_energy( + temperature=t_bc, + species_mass_fractions=state_minus.species_mass_fractions) + + total_energy_bc = ( + state_minus.mass_density * internal_energy_bc + + 0.5*np.dot(mom_bc, mom_bc)/state_minus.mass_density) + + return replace_fluid_state( + state_minus, gas_model, + energy=total_energy_bc, + momentum=mom_bc) + + def grad_cv_bc(self, dcoll, dd_bdry, gas_model, state_minus, grad_cv_minus, + normal, **kwargs): + """ + Return grad(CV) used in the boundary calculation of viscous flux. + + Specify the velocity gradients on the external state to ensure zero + energy and momentum flux due to shear stresses. + + Gradients of species mass fractions are set to zero in the normal direction + to ensure zero flux of species across the boundary. + """ + dd_bdry = as_dofdesc(dd_bdry) + normal = state_minus.array_context.thaw(dcoll.normal(dd_bdry)) + state_bc = self.state_bc( + dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model, + state_minus=state_minus, **kwargs) + + grad_v_bc = self._slip.grad_velocity_bc( + state_minus, state_bc, grad_cv_minus, normal) + + grad_mom_bc = ( + state_bc.mass_density * grad_v_bc + + np.outer(state_bc.velocity, grad_cv_minus.mass)) + + grad_species_mass_bc = self._impermeable.grad_species_mass_bc( + state_minus, grad_cv_minus, normal) + + return grad_cv_minus.replace( + momentum=grad_mom_bc, + species_mass=grad_species_mass_bc) + + def grad_temperature_bc(self, dcoll, dd_bdry, grad_t_minus, **kwargs): + """Return BC on grad(temperature).""" + # Mengaldo Eqns (50-51) + return grad_t_minus + + def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): + """Return state with reflected normal-component velocity.""" + actx = state_minus.array_context + + # Grab a unit normal to the boundary + nhat = actx.thaw(dcoll.normal(dd_bdry)) + + # set the normal momentum to 0 + mom_plus = self._slip.momentum_plus(state_minus.momentum_density, nhat) + return replace_fluid_state(state_minus, gas_model, momentum=mom_plus) + + class IsothermalWallBoundary(MengaldoBoundaryCondition): r"""Isothermal viscous wall boundary. @@ -1636,7 +1734,8 @@ def __init__(self, wall_temperature=300): def temperature_bc(self, dcoll, dd_bdry, state_minus, **kwargs): """Get temperature value used in grad(T).""" - return 0*state_minus.temperature + self._wall_temp + wall_temp = project_from_base(dcoll, dd_bdry, self._wall_temp) + return 0*state_minus.temperature + wall_temp def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): """Return BC fluid state.""" diff --git a/mirgecom/multiphysics/__init__.py b/mirgecom/multiphysics/__init__.py index bd412d3bd..3af04f8a0 100644 --- a/mirgecom/multiphysics/__init__.py +++ b/mirgecom/multiphysics/__init__.py @@ -26,4 +26,38 @@ __doc__ = """ .. automodule:: mirgecom.multiphysics.thermally_coupled_fluid_wall + +.. autofunction:: make_interface_boundaries """ + + +def make_interface_boundaries(bdry_factories, inter_volume_tpairs): + """ + Create volume-pairwise interface boundaries from inter-volume data. + + Return a :class:`dict` mapping a (directional) pair of + :class:`grudge.dof_desc.DOFDesc` *(other_vol_dd, self_vol_dd)* representing + neighboring volumes to a :class:`dict` of boundary objects. Specifically, + *interface_boundaries[other_vol_dd, self_vol_dd]* maps each interface boundary + :class:`~grudge.dof_desc.DOFDesc` to a boundary object. + + Parameters + ---------- + + bdry_factories + + Mapping from directional volume :class:`~grudge.dof_desc.DOFDesc` pair to + a function that takes a :class:`grudge.trace_pair.TracePair` and returns + an interface boundary object. + + inter_volume_tpairs + + Mapping from directional volume :class:`~grudge.dof_desc.DOFDesc` pair to + a :class:`list` of :class:`grudge.trace_pair.TracePair` (as is returned by + :func:`grudge.trace_pair.inter_volume_trace_pairs`). + """ + return { + vol_dd_pair: { + tpair.dd: bdry_factories[vol_dd_pair](tpair) + for tpair in tpairs} + for vol_dd_pair, tpairs in inter_volume_tpairs.items()} diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py index b91839011..d5720a706 100644 --- a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py +++ b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py @@ -13,9 +13,19 @@ at the interface. -.. autofunction:: get_interface_boundaries -.. autofunction:: coupled_grad_t_operator -.. autofunction:: coupled_ns_heat_operator +Boundary Setup Functions +^^^^^^^^^^^^^^^^^^^^^^^^ + +.. autofunction:: add_interface_boundaries_no_grad +.. autofunction:: add_interface_boundaries + +Basic Coupled Operators +^^^^^^^^^^^^^^^^^^^^^^^ + +.. autofunction:: basic_coupled_ns_heat_operator + +Boundary Conditions +^^^^^^^^^^^^^^^^^^^ .. autoclass:: InterfaceFluidBoundary .. autoclass:: InterfaceFluidSlipBoundary @@ -47,10 +57,13 @@ THE SOFTWARE. """ -from dataclasses import replace +from dataclasses import dataclass, replace import numpy as np from abc import abstractmethod +from functools import partial +from arraycontext import dataclass_array_container +from meshmode.dof_array import DOFArray from grudge.trace_pair import ( TracePair, inter_volume_trace_pairs @@ -85,13 +98,15 @@ grad_operator as wall_grad_t_operator, diffusion_operator, ) +from mirgecom.multiphysics import make_interface_boundaries +from mirgecom.utils import project_from_base -class _TemperatureInterVolTag: +class _ThermalDataNoGradInterVolTag: pass -class _KappaInterVolTag: +class _ThermalDataInterVolTag: pass @@ -119,20 +134,11 @@ class _WallOperatorTag: pass -def _project_from_base(dcoll, dd_bdry, field): - """Project *field* from *DISCR_TAG_BASE* to the same discr. as *dd_bdry*.""" - if dd_bdry.discretization_tag is not DISCR_TAG_BASE: - dd_bdry_base = dd_bdry.with_discr_tag(DISCR_TAG_BASE) - return op.project(dcoll, dd_bdry_base, dd_bdry, field) - else: - return field - - # FIXME: Interior penalty should probably use an average of the lengthscales on # both sides of the interface class InterfaceFluidBoundary(MengaldoBoundaryCondition): r""" - Abstract interface for the fluid side of the fluid-wall interface boundary. + Abstract interface for the fluid side of the fluid-wall interface. Extends :class:`~mirgecom.boundary.MengaldoBoundaryCondition` to include an interior penalty on the heat flux: @@ -215,7 +221,7 @@ def viscous_divergence_flux( state_minus=state_minus, numerical_flux_func=numerical_flux_func, grad_cv_minus=grad_cv_minus, grad_t_minus=grad_t_minus, **kwargs) - lengthscales_minus = _project_from_base( + lengthscales_minus = project_from_base( dcoll, dd_bdry, self._lengthscales_minus) tau = ( @@ -235,42 +241,36 @@ def viscous_divergence_flux( class _ThermallyCoupledHarmonicMeanBoundaryComponent: def __init__( - self, kappa_plus, t_plus, grad_t_plus=None, - use_kappa_weighted_t_bc=False): + self, kappa_plus, t_plus, grad_t_plus=None): self._kappa_plus = kappa_plus self._t_plus = t_plus self._grad_t_plus = grad_t_plus - self._use_kappa_weighted_t_bc = use_kappa_weighted_t_bc def kappa_plus(self, dcoll, dd_bdry): - return _project_from_base(dcoll, dd_bdry, self._kappa_plus) + return project_from_base(dcoll, dd_bdry, self._kappa_plus) def kappa_bc(self, dcoll, dd_bdry, kappa_minus): - kappa_plus = _project_from_base(dcoll, dd_bdry, self._kappa_plus) + kappa_plus = project_from_base(dcoll, dd_bdry, self._kappa_plus) return harmonic_mean(kappa_minus, kappa_plus) def temperature_plus(self, dcoll, dd_bdry): - return _project_from_base(dcoll, dd_bdry, self._t_plus) + return project_from_base(dcoll, dd_bdry, self._t_plus) def temperature_bc(self, dcoll, dd_bdry, kappa_minus, t_minus): - t_plus = _project_from_base(dcoll, dd_bdry, self._t_plus) - if self._use_kappa_weighted_t_bc: - actx = t_minus.array_context - kappa_plus = _project_from_base(dcoll, dd_bdry, self._kappa_plus) - kappa_sum = actx.np.where( - actx.np.greater(kappa_minus + kappa_plus, 0*kappa_minus), - kappa_minus + kappa_plus, - 0*kappa_minus + 1) - return (t_minus * kappa_minus + t_plus * kappa_plus)/kappa_sum - else: - return (t_minus + t_plus)/2 - - def grad_temperature_bc( - self, dcoll, dd_bdry, grad_t_minus): + t_plus = project_from_base(dcoll, dd_bdry, self._t_plus) + actx = t_minus.array_context + kappa_plus = project_from_base(dcoll, dd_bdry, self._kappa_plus) + kappa_sum = actx.np.where( + actx.np.greater(kappa_minus + kappa_plus, 0*kappa_minus), + kappa_minus + kappa_plus, + 0*kappa_minus + 1) + return (t_minus * kappa_minus + t_plus * kappa_plus)/kappa_sum + + def grad_temperature_bc(self, dcoll, dd_bdry, grad_t_minus): if self._grad_t_plus is None: raise ValueError( "Boundary does not have external temperature gradient data.") - grad_t_plus = _project_from_base(dcoll, dd_bdry, self._grad_t_plus) + grad_t_plus = project_from_base(dcoll, dd_bdry, self._grad_t_plus) return (grad_t_plus + grad_t_minus)/2 @@ -295,8 +295,7 @@ class InterfaceFluidSlipBoundary(InterfaceFluidBoundary): def __init__( self, kappa_plus, t_plus, grad_t_plus=None, - heat_flux_penalty_amount=None, lengthscales_minus=None, - use_kappa_weighted_grad_t_flux=False): + heat_flux_penalty_amount=None, lengthscales_minus=None): r""" Initialize InterfaceFluidSlipBoundary. @@ -325,12 +324,6 @@ def __init__( lengthscales_minus: :class:`meshmode.dof_array.DOFArray` or None Characteristic mesh spacing $h^-$. - - use_kappa_weighted_grad_t_flux: bool - - Indicates whether the temperature gradient flux at the interface should - be computed using a simple average of temperatures or by weighting the - temperature from each side by its respective thermal conductivity. """ InterfaceFluidBoundary.__init__( self, @@ -340,8 +333,7 @@ def __init__( self._thermally_coupled = _ThermallyCoupledHarmonicMeanBoundaryComponent( kappa_plus=kappa_plus, t_plus=t_plus, - grad_t_plus=grad_t_plus, - use_kappa_weighted_t_bc=use_kappa_weighted_grad_t_flux) + grad_t_plus=grad_t_plus) self._slip = _SlipBoundaryComponent() self._impermeable = _ImpermeableBoundaryComponent() @@ -450,8 +442,7 @@ class InterfaceFluidNoslipBoundary(InterfaceFluidBoundary): def __init__( self, kappa_plus, t_plus, grad_t_plus=None, - heat_flux_penalty_amount=None, lengthscales_minus=None, - use_kappa_weighted_grad_t_flux=False): + heat_flux_penalty_amount=None, lengthscales_minus=None): r""" Initialize InterfaceFluidNoslipBoundary. @@ -480,12 +471,6 @@ def __init__( lengthscales_minus: :class:`meshmode.dof_array.DOFArray` or None Characteristic mesh spacing $h^-$. - - use_kappa_weighted_grad_t_flux: bool - - Indicates whether the temperature gradient flux at the interface should - be computed using a simple average of temperatures or by weighting the - temperature from each side by its respective thermal conductivity. """ InterfaceFluidBoundary.__init__( self, @@ -495,8 +480,7 @@ def __init__( self._thermally_coupled = _ThermallyCoupledHarmonicMeanBoundaryComponent( kappa_plus=kappa_plus, t_plus=t_plus, - grad_t_plus=grad_t_plus, - use_kappa_weighted_t_bc=use_kappa_weighted_grad_t_flux) + grad_t_plus=grad_t_plus) self._no_slip = _NoSlipBoundaryComponent() self._impermeable = _ImpermeableBoundaryComponent() @@ -577,7 +561,7 @@ def __init__(self, kappa_plus, u_plus, grad_u_plus=None): r""" Initialize InterfaceWallBoundary. - Argument *grad_u_plus*, is only required if the boundary will be used to + Argument *grad_u_plus* is only required if the boundary will be used to compute the heat flux. Parameters @@ -586,11 +570,11 @@ def __init__(self, kappa_plus, u_plus, grad_u_plus=None): Thermal conductivity from the fluid side. - t_plus: :class:`meshmode.dof_array.DOFArray` + u_plus: :class:`meshmode.dof_array.DOFArray` Temperature from the fluid side. - grad_t_plus: :class:`meshmode.dof_array.DOFArray` or None + grad_u_plus: :class:`meshmode.dof_array.DOFArray` or None Temperature gradient from the fluid side. """ @@ -604,11 +588,11 @@ def get_grad_flux( actx = u_minus.array_context normal = actx.thaw(dcoll.normal(dd_bdry)) - kappa_plus = _project_from_base(dcoll, dd_bdry, self.kappa_plus) + kappa_plus = project_from_base(dcoll, dd_bdry, self.kappa_plus) kappa_tpair = TracePair( dd_bdry, interior=kappa_minus, exterior=kappa_plus) - u_plus = _project_from_base(dcoll, dd_bdry, self.u_plus) + u_plus = project_from_base(dcoll, dd_bdry, self.u_plus) u_tpair = TracePair(dd_bdry, interior=u_minus, exterior=u_plus) return numerical_flux_func(kappa_tpair, u_tpair, normal) @@ -624,14 +608,14 @@ def get_diffusion_flux( actx = u_minus.array_context normal = actx.thaw(dcoll.normal(dd_bdry)) - kappa_plus = _project_from_base(dcoll, dd_bdry, self.kappa_plus) + kappa_plus = project_from_base(dcoll, dd_bdry, self.kappa_plus) kappa_tpair = TracePair( dd_bdry, interior=kappa_minus, exterior=kappa_plus) - u_plus = _project_from_base(dcoll, dd_bdry, self.u_plus) + u_plus = project_from_base(dcoll, dd_bdry, self.u_plus) u_tpair = TracePair(dd_bdry, interior=u_minus, exterior=u_plus) - grad_u_plus = _project_from_base(dcoll, dd_bdry, self.grad_u_plus) + grad_u_plus = project_from_base(dcoll, dd_bdry, self.grad_u_plus) grad_u_tpair = TracePair( dd_bdry, interior=grad_u_minus, exterior=grad_u_plus) @@ -643,80 +627,200 @@ def get_diffusion_flux( penalty_amount=penalty_amount) -def _kappa_inter_volume_trace_pairs( +@dataclass_array_container +@dataclass(frozen=True) +class _ThermalDataNoGrad: + kappa: DOFArray + temperature: DOFArray + + +@dataclass_array_container +@dataclass(frozen=True) +class _ThermalData(_ThermalDataNoGrad): + grad_temperature: np.ndarray + + +def _make_thermal_data(kappa, temperature, grad_temperature=None): + if not isinstance(kappa, DOFArray): + kappa = kappa * (0*temperature + 1) + + if grad_temperature is not None: + thermal_data = _ThermalData(kappa, temperature, grad_temperature) + else: + thermal_data = _ThermalDataNoGrad(kappa, temperature) + + return thermal_data + + +def _get_interface_trace_pairs_no_grad( dcoll, - gas_model, fluid_dd, wall_dd, - fluid_state, wall_kappa): - """Exchange thermal conductivity across the fluid-wall interface.""" - actx = fluid_state.array_context - fluid_kappa = fluid_state.thermal_conductivity - - # Promote constant-valued kappas to DOFArrays - from meshmode.dof_array import DOFArray - if not isinstance(fluid_kappa, DOFArray): - fluid_kappa = fluid_kappa * (dcoll.zeros(actx, dd=fluid_dd) + 1) - if not isinstance(wall_kappa, DOFArray): - wall_kappa = wall_kappa * (dcoll.zeros(actx, dd=wall_dd) + 1) - - pairwise_kappa = { - (fluid_dd, wall_dd): (fluid_kappa, wall_kappa)} + fluid_kappa, wall_kappa, + fluid_temperature, wall_temperature, + *, + comm_tag=None): + pairwise_thermal_data = { + (fluid_dd, wall_dd): ( + _make_thermal_data(fluid_kappa, fluid_temperature), + _make_thermal_data(wall_kappa, wall_temperature))} return inter_volume_trace_pairs( - dcoll, pairwise_kappa, comm_tag=_KappaInterVolTag) + dcoll, pairwise_thermal_data, + comm_tag=(_ThermalDataNoGradInterVolTag, comm_tag)) -def _temperature_inter_volume_trace_pairs( +def _get_interface_trace_pairs( dcoll, - gas_model, fluid_dd, wall_dd, - fluid_state, wall_temperature): - """Exchange temperature across the fluid-wall interface.""" - pairwise_temperature = { - (fluid_dd, wall_dd): - (fluid_state.temperature, wall_temperature)} + fluid_kappa, wall_kappa, + fluid_temperature, wall_temperature, + fluid_grad_temperature, wall_grad_temperature, + *, + comm_tag=None): + pairwise_thermal_data = { + (fluid_dd, wall_dd): ( + _make_thermal_data( + fluid_kappa, + fluid_temperature, + fluid_grad_temperature), + _make_thermal_data( + wall_kappa, + wall_temperature, + wall_grad_temperature))} + return inter_volume_trace_pairs( - dcoll, pairwise_temperature, comm_tag=_TemperatureInterVolTag) + dcoll, pairwise_thermal_data, + comm_tag=(_ThermalDataInterVolTag, comm_tag)) -def _grad_temperature_inter_volume_trace_pairs( +def _get_interface_boundaries_no_grad( dcoll, - gas_model, fluid_dd, wall_dd, - fluid_grad_temperature, wall_grad_temperature): - """Exchange temperature gradient across the fluid-wall interface.""" - pairwise_grad_temperature = { - (fluid_dd, wall_dd): - (fluid_grad_temperature, wall_grad_temperature)} - return inter_volume_trace_pairs( - dcoll, pairwise_grad_temperature, comm_tag=_GradTemperatureInterVolTag) + fluid_kappa, wall_kappa, + fluid_temperature, wall_temperature, + *, + interface_noslip=True, + quadrature_tag=DISCR_TAG_BASE, + comm_tag=None): + interface_tpairs = _get_interface_trace_pairs_no_grad( + dcoll, + fluid_dd, wall_dd, + fluid_kappa, wall_kappa, + fluid_temperature, wall_temperature, + comm_tag=comm_tag) + if interface_noslip: + fluid_bc_class = InterfaceFluidNoslipBoundary + else: + fluid_bc_class = InterfaceFluidSlipBoundary + + def make_fluid_boundary(interface_tpair): + return fluid_bc_class( + interface_tpair.ext.kappa, + interface_tpair.ext.temperature) + + def make_wall_boundary(interface_tpair): + return InterfaceWallBoundary( + interface_tpair.ext.kappa, + interface_tpair.ext.temperature) + + bdry_factories = { + (wall_dd, fluid_dd): make_fluid_boundary, + (fluid_dd, wall_dd): make_wall_boundary} -def get_interface_boundaries( + interface_boundaries = make_interface_boundaries( + bdry_factories, interface_tpairs) + + fluid_interface_boundaries = interface_boundaries[wall_dd, fluid_dd] + wall_interface_boundaries = interface_boundaries[fluid_dd, wall_dd] + + return fluid_interface_boundaries, wall_interface_boundaries + + +def _get_interface_boundaries( + dcoll, + fluid_dd, wall_dd, + fluid_kappa, wall_kappa, + fluid_temperature, wall_temperature, + fluid_grad_temperature, wall_grad_temperature, + *, + interface_noslip=True, + wall_penalty_amount=None, + quadrature_tag=DISCR_TAG_BASE, + comm_tag=None): + if wall_penalty_amount is None: + # FIXME: After verifying the form of the penalty term, figure out what value + # makes sense to use as a default here + wall_penalty_amount = 0.05 + + interface_tpairs = _get_interface_trace_pairs( + dcoll, + fluid_dd, wall_dd, + fluid_kappa, wall_kappa, + fluid_temperature, wall_temperature, + fluid_grad_temperature, wall_grad_temperature, + comm_tag=comm_tag) + + if interface_noslip: + fluid_bc_class = InterfaceFluidNoslipBoundary + else: + fluid_bc_class = InterfaceFluidSlipBoundary + + # Diffusion operator passes lengthscales_minus into the boundary flux + # functions, but NS doesn't; thus we need to pass lengthscales into + # the fluid boundary condition constructor + from grudge.dt_utils import characteristic_lengthscales + fluid_lengthscales = ( + characteristic_lengthscales( + fluid_temperature.array_context, dcoll, fluid_dd) + * (0*fluid_temperature+1)) + + def make_fluid_boundary(interface_tpair): + return fluid_bc_class( + interface_tpair.ext.kappa, + interface_tpair.ext.temperature, + interface_tpair.ext.grad_temperature, + wall_penalty_amount, + lengthscales_minus=op.project(dcoll, + fluid_dd, interface_tpair.dd, fluid_lengthscales)) + + def make_wall_boundary(interface_tpair): + return InterfaceWallBoundary( + interface_tpair.ext.kappa, + interface_tpair.ext.temperature, + interface_tpair.ext.grad_temperature) + + bdry_factories = { + (wall_dd, fluid_dd): make_fluid_boundary, + (fluid_dd, wall_dd): make_wall_boundary} + + interface_boundaries = make_interface_boundaries( + bdry_factories, interface_tpairs) + + fluid_interface_boundaries = interface_boundaries[wall_dd, fluid_dd] + wall_interface_boundaries = interface_boundaries[fluid_dd, wall_dd] + + return fluid_interface_boundaries, wall_interface_boundaries + + +def add_interface_boundaries_no_grad( dcoll, gas_model, fluid_dd, wall_dd, fluid_state, wall_kappa, wall_temperature, - fluid_grad_temperature=None, wall_grad_temperature=None, + fluid_boundaries, wall_boundaries, *, interface_noslip=True, - use_kappa_weighted_grad_flux_in_fluid=False, wall_penalty_amount=None, quadrature_tag=DISCR_TAG_BASE, - # Added to avoid repeated computation - # FIXME: See if there's a better way to do this - _kappa_inter_vol_tpairs=None, - _temperature_inter_vol_tpairs=None, - _grad_temperature_inter_vol_tpairs=None): + comm_tag=None): """ - Get the fluid-wall interface boundaries. + Include the fluid-wall interface boundaries (without temperature gradient). - Return a tuple `(fluid_interface_boundaries, wall_interface_boundaries)` in - which each of the two entries is a mapping from each interface boundary's - :class:`grudge.dof_desc.BoundaryDomainTag` to a boundary condition object - compatible with that subdomain's operators. The map contains one entry for - the collection of faces whose opposite face reside on the current MPI rank - and one-per-rank for each collection of faces whose opposite face resides on - a different rank. + Return a tuple `(fluid_all_boundaries, wall_all_boundaries)` that adds boundaries + to *fluid_boundaries* and *wall_boundaries* that represent the volume interfaces. + One entry is added for the collection of faces whose opposite face reside on the + current MPI rank and one-per-rank for each collection of faces whose opposite + face resides on a different rank. Parameters ---------- @@ -751,28 +855,21 @@ def get_interface_boundaries( Temperature for the wall volume. - fluid_grad_temperature: numpy.ndarray or None + fluid_boundaries - Temperature gradient for the fluid volume. Only needed if boundaries will - be used to compute viscous fluxes. + Dictionary of boundary functions, one for each valid non-interface + :class:`~grudge.dof_desc.BoundaryDomainTag` on the fluid subdomain. - wall_grad_temperature: numpy.ndarray or None + wall_boundaries - Temperature gradient for the wall volume. Only needed if boundaries will - be used to compute diffusion fluxes. + Dictionary of boundary functions, one for each valid non-interface + :class:`~grudge.dof_desc.BoundaryDomainTag` on the wall subdomain. interface_noslip: bool If `True`, interface boundaries on the fluid side will be treated as no-slip walls. If `False` they will be treated as slip walls. - use_kappa_weighted_grad_flux_in_fluid: bool - - Indicates whether the temperature gradient flux on the fluid side of the - interface should be computed using a simple average of temperatures or by - weighting the temperature from each side by its respective thermal - conductivity. - wall_penalty_amount: float Coefficient $c$ for the interior penalty on the heat flux. See @@ -783,112 +880,143 @@ def get_interface_boundaries( An identifier denoting a particular quadrature discretization to use during operator evaluations. + + comm_tag: Hashable + Tag for distributed communication """ - if interface_noslip: - fluid_bc_class = InterfaceFluidNoslipBoundary - else: - fluid_bc_class = InterfaceFluidSlipBoundary + fluid_interface_boundaries_no_grad, wall_interface_boundaries_no_grad = \ + _get_interface_boundaries_no_grad( + dcoll, + fluid_dd, wall_dd, + fluid_state.tv.thermal_conductivity, wall_kappa, + fluid_state.temperature, wall_temperature, + interface_noslip=interface_noslip, + quadrature_tag=quadrature_tag, + comm_tag=comm_tag) - assert ( - (fluid_grad_temperature is None) == (wall_grad_temperature is None)), ( - "Expected both fluid_grad_temperature and wall_grad_temperature or neither") + fluid_all_boundaries_no_grad = {} + fluid_all_boundaries_no_grad.update(fluid_boundaries) + fluid_all_boundaries_no_grad.update(fluid_interface_boundaries_no_grad) - include_gradient = fluid_grad_temperature is not None + wall_all_boundaries_no_grad = {} + wall_all_boundaries_no_grad.update(wall_boundaries) + wall_all_boundaries_no_grad.update(wall_interface_boundaries_no_grad) - # Exchange thermal conductivity, temperature, and (optionally) temperature - # gradient + return fluid_all_boundaries_no_grad, wall_all_boundaries_no_grad - if _kappa_inter_vol_tpairs is None: - kappa_inter_vol_tpairs = _kappa_inter_volume_trace_pairs( - dcoll, - gas_model, - fluid_dd, wall_dd, - fluid_state, wall_kappa) - else: - kappa_inter_vol_tpairs = _kappa_inter_vol_tpairs - if _temperature_inter_vol_tpairs is None: - temperature_inter_vol_tpairs = _temperature_inter_volume_trace_pairs( +def add_interface_boundaries( + dcoll, + gas_model, + fluid_dd, wall_dd, + fluid_state, wall_kappa, wall_temperature, + fluid_grad_temperature, wall_grad_temperature, + fluid_boundaries, wall_boundaries, + *, + interface_noslip=True, + wall_penalty_amount=None, + quadrature_tag=DISCR_TAG_BASE, + comm_tag=None): + """ + Include the fluid-wall interface boundaries. + + Return a tuple `(fluid_all_boundaries, wall_all_boundaries)` that adds boundaries + to *fluid_boundaries* and *wall_boundaries* that represent the volume interfaces. + One entry is added for the collection of faces whose opposite face reside on the + current MPI rank and one-per-rank for each collection of faces whose opposite + face resides on a different rank. + + Parameters + ---------- + + dcoll: class:`~grudge.discretization.DiscretizationCollection` + + A discretization collection encapsulating the DG elements + + gas_model: :class:`~mirgecom.gas_model.GasModel` + + Physical gas model including equation of state, transport, + and kinetic properties as required by fluid state + + fluid_dd: :class:`grudge.dof_desc.DOFDesc` + + DOF descriptor for the fluid volume. + + wall_dd: :class:`grudge.dof_desc.DOFDesc` + + DOF descriptor for the wall volume. + + fluid_state: :class:`~mirgecom.gas_model.FluidState` + + Fluid state object with the conserved state and dependent + quantities for the fluid volume. + + wall_kappa: float or :class:`meshmode.dof_array.DOFArray` + + Thermal conductivity for the wall volume. + + wall_temperature: :class:`meshmode.dof_array.DOFArray` + + Temperature for the wall volume. + + fluid_grad_temperature: numpy.ndarray + + Temperature gradient for the fluid volume. + + wall_grad_temperature: numpy.ndarray + + Temperature gradient for the wall volume. + + fluid_boundaries + + Dictionary of boundary functions, one for each valid non-interface + :class:`~grudge.dof_desc.BoundaryDomainTag` on the fluid subdomain. + + wall_boundaries + + Dictionary of boundary functions, one for each valid non-interface + :class:`~grudge.dof_desc.BoundaryDomainTag` on the wall subdomain. + + interface_noslip: bool + + If `True`, interface boundaries on the fluid side will be treated as + no-slip walls. If `False` they will be treated as slip walls. + + wall_penalty_amount: float + + Coefficient $c$ for the interior penalty on the heat flux. See + :class:`~mirgecom.multiphysics.thermally_coupled_fluid_wall.InterfaceFluidBoundary` + for details. + + quadrature_tag + + An identifier denoting a particular quadrature discretization to use during + operator evaluations. + + comm_tag: Hashable + Tag for distributed communication + """ + fluid_interface_boundaries, wall_interface_boundaries = \ + _get_interface_boundaries( dcoll, - gas_model, fluid_dd, wall_dd, - fluid_state, wall_temperature) - else: - temperature_inter_vol_tpairs = _temperature_inter_vol_tpairs - - if include_gradient: - if _grad_temperature_inter_vol_tpairs is None: - grad_temperature_inter_vol_tpairs = \ - _grad_temperature_inter_volume_trace_pairs( - dcoll, - gas_model, - fluid_dd, wall_dd, - fluid_grad_temperature, wall_grad_temperature) - else: - grad_temperature_inter_vol_tpairs = _grad_temperature_inter_vol_tpairs - else: - grad_temperature_inter_vol_tpairs = None - - # Set up the interface boundaries - - if include_gradient: - - # Diffusion operator passes lengthscales_minus into the boundary flux - # functions, but NS doesn't; thus we need to pass lengthscales into - # the fluid boundary condition constructor - from grudge.dt_utils import characteristic_lengthscales - fluid_lengthscales = ( - characteristic_lengthscales( - fluid_state.array_context, dcoll, fluid_dd) - * (0*fluid_state.temperature+1)) - - # Construct interface boundaries with temperature gradient - - fluid_interface_boundaries = { - kappa_tpair.dd.domain_tag: fluid_bc_class( - kappa_tpair.ext, - temperature_tpair.ext, - grad_temperature_tpair.ext, - wall_penalty_amount, - lengthscales_minus=op.project(dcoll, - fluid_dd, temperature_tpair.dd, fluid_lengthscales), - use_kappa_weighted_grad_t_flux=use_kappa_weighted_grad_flux_in_fluid) - for kappa_tpair, temperature_tpair, grad_temperature_tpair in zip( - kappa_inter_vol_tpairs[wall_dd, fluid_dd], - temperature_inter_vol_tpairs[wall_dd, fluid_dd], - grad_temperature_inter_vol_tpairs[wall_dd, fluid_dd])} - - wall_interface_boundaries = { - kappa_tpair.dd.domain_tag: InterfaceWallBoundary( - kappa_tpair.ext, - temperature_tpair.ext, - grad_temperature_tpair.ext) - for kappa_tpair, temperature_tpair, grad_temperature_tpair in zip( - kappa_inter_vol_tpairs[fluid_dd, wall_dd], - temperature_inter_vol_tpairs[fluid_dd, wall_dd], - grad_temperature_inter_vol_tpairs[fluid_dd, wall_dd])} - else: + fluid_state.tv.thermal_conductivity, wall_kappa, + fluid_state.temperature, wall_temperature, + fluid_grad_temperature, wall_grad_temperature, + interface_noslip=interface_noslip, + wall_penalty_amount=wall_penalty_amount, + quadrature_tag=quadrature_tag, + comm_tag=comm_tag) - # Construct interface boundaries without temperature gradient - - fluid_interface_boundaries = { - kappa_tpair.dd.domain_tag: fluid_bc_class( - kappa_tpair.ext, - temperature_tpair.ext, - use_kappa_weighted_grad_t_flux=use_kappa_weighted_grad_flux_in_fluid) - for kappa_tpair, temperature_tpair in zip( - kappa_inter_vol_tpairs[wall_dd, fluid_dd], - temperature_inter_vol_tpairs[wall_dd, fluid_dd])} - - wall_interface_boundaries = { - kappa_tpair.dd.domain_tag: InterfaceWallBoundary( - kappa_tpair.ext, - temperature_tpair.ext) - for kappa_tpair, temperature_tpair in zip( - kappa_inter_vol_tpairs[fluid_dd, wall_dd], - temperature_inter_vol_tpairs[fluid_dd, wall_dd])} + fluid_all_boundaries = {} + fluid_all_boundaries.update(fluid_boundaries) + fluid_all_boundaries.update(fluid_interface_boundaries) - return fluid_interface_boundaries, wall_interface_boundaries + wall_all_boundaries = {} + wall_all_boundaries.update(wall_boundaries) + wall_all_boundaries.update(wall_interface_boundaries) + + return fluid_all_boundaries, wall_all_boundaries def coupled_grad_t_operator( @@ -900,19 +1028,21 @@ def coupled_grad_t_operator( *, time=0., interface_noslip=True, - use_kappa_weighted_grad_flux_in_fluid=False, + use_kappa_weighted_grad_flux_in_fluid=None, quadrature_tag=DISCR_TAG_BASE, fluid_numerical_flux_func=num_flux_central, # Added to avoid repeated computation # FIXME: See if there's a better way to do this - _kappa_inter_vol_tpairs=None, - _temperature_inter_vol_tpairs=None, _fluid_operator_states_quad=None, - _fluid_interface_boundaries_no_grad=None, - _wall_interface_boundaries_no_grad=None): + _fluid_all_boundaries_no_grad=None, + _wall_all_boundaries_no_grad=None): r""" Compute $\nabla T$ on the fluid and wall subdomains. + Deprecated; set up interface boundaries explicitly via + :func:`add_interface_boundaries_no_grad` and include them when calling the + individual operators instead. + Parameters ---------- @@ -990,6 +1120,25 @@ def coupled_grad_t_operator( The tuple `(fluid_grad_temperature, wall_grad_temperature)`. """ + from warnings import warn + warn( + "coupled_grad_t_operator is deprecated and will disappear in Q3 2023. " + "Set up interface boundaries explicitly via " + ":func:`add_interface_boundaries_no_grad` and include them when calling the " + "individual operators instead.", DeprecationWarning, stacklevel=2) + + if use_kappa_weighted_grad_flux_in_fluid is None: + warn( + "Default value of use_kappa_weighted_grad_flux_in_fluid has changed " + "from False to True as False is no longer allowed. Explicitly set it " + "to True to suppress this warning.", UserWarning, stacklevel=2) + use_kappa_weighted_grad_flux_in_fluid = True + elif not use_kappa_weighted_grad_flux_in_fluid: + warn( + "Setting use_kappa_weighted_grad_flux_in_fluid to True; False is no " + "longer allowed. Explicitly set it to True to suppress this warning.", + UserWarning, stacklevel=2) + fluid_boundaries = { as_dofdesc(bdtag).domain_tag: bdry for bdtag, bdry in fluid_boundaries.items()} @@ -997,42 +1146,30 @@ def coupled_grad_t_operator( as_dofdesc(bdtag).domain_tag: bdry for bdtag, bdry in wall_boundaries.items()} - # Construct boundaries for the fluid-wall interface; no temperature gradient + # Include boundaries for the fluid-wall interface; no temperature gradient # yet because that's what we're trying to compute assert ( - (_fluid_interface_boundaries_no_grad is None) - == (_wall_interface_boundaries_no_grad is None)), ( - "Expected both _fluid_interface_boundaries_no_grad and " - "_wall_interface_boundaries_no_grad or neither") + (_fluid_all_boundaries_no_grad is None) + == (_wall_all_boundaries_no_grad is None)), ( + "Expected both _fluid_all_boundaries_no_grad and " + "_wall_all_boundaries_no_grad or neither") - if _fluid_interface_boundaries_no_grad is None: + if _fluid_all_boundaries_no_grad is None: # Note: We don't need to supply wall_penalty_amount here since we're only # using these to compute the temperature gradient - fluid_interface_boundaries_no_grad, wall_interface_boundaries_no_grad = \ - get_interface_boundaries( + fluid_all_boundaries_no_grad, wall_all_boundaries_no_grad = \ + add_interface_boundaries_no_grad( dcoll, gas_model, fluid_dd, wall_dd, fluid_state, wall_kappa, wall_temperature, + fluid_boundaries, wall_boundaries, interface_noslip=interface_noslip, - use_kappa_weighted_grad_flux_in_fluid=( - use_kappa_weighted_grad_flux_in_fluid), - _kappa_inter_vol_tpairs=_kappa_inter_vol_tpairs, - _temperature_inter_vol_tpairs=_temperature_inter_vol_tpairs) + quadrature_tag=quadrature_tag) else: - fluid_interface_boundaries_no_grad = _fluid_interface_boundaries_no_grad - wall_interface_boundaries_no_grad = _wall_interface_boundaries_no_grad - - # Augment the domain boundaries with the interface boundaries - - fluid_all_boundaries_no_grad = {} - fluid_all_boundaries_no_grad.update(fluid_boundaries) - fluid_all_boundaries_no_grad.update(fluid_interface_boundaries_no_grad) - - wall_all_boundaries_no_grad = {} - wall_all_boundaries_no_grad.update(wall_boundaries) - wall_all_boundaries_no_grad.update(wall_interface_boundaries_no_grad) + fluid_all_boundaries_no_grad = _fluid_all_boundaries_no_grad + wall_all_boundaries_no_grad = _wall_all_boundaries_no_grad # Compute the subdomain gradient operators using the augmented boundaries @@ -1057,7 +1194,7 @@ def coupled_ns_heat_operator( *, time=0., interface_noslip=True, - use_kappa_weighted_grad_flux_in_fluid=False, + use_kappa_weighted_grad_flux_in_fluid=None, wall_penalty_amount=None, quadrature_tag=DISCR_TAG_BASE, limiter_func=None, @@ -1072,6 +1209,10 @@ def coupled_ns_heat_operator( fluid-wall interface that are needed to enforce continuity of temperature and heat flux. + Deprecated; set up interface boundaries explicitly via + :func:`add_interface_boundaries` and include them when calling the individual + operators instead. + Parameters ---------- @@ -1126,6 +1267,12 @@ def coupled_ns_heat_operator( If `True`, interface boundaries on the fluid side will be treated as no-slip walls. If `False` they will be treated as slip walls. + wall_penalty_amount: float + + Coefficient $c$ for the interior penalty on the heat flux. See + :class:`~mirgecom.multiphysics.thermally_coupled_fluid_wall.InterfaceFluidBoundary` + for details. + use_kappa_weighted_grad_flux_in_fluid: bool Indicates whether the temperature gradient flux on the fluid side of the @@ -1133,12 +1280,6 @@ def coupled_ns_heat_operator( weighting the temperature from each side by its respective thermal conductivity. - wall_penalty_amount: float - - Coefficient $c$ for the interior penalty on the heat flux. See - :class:`~mirgecom.multiphysics.thermally_coupled_fluid_wall.InterfaceFluidBoundary` - for details. - quadrature_tag: An identifier denoting a particular quadrature discretization to use during @@ -1174,6 +1315,25 @@ def coupled_ns_heat_operator( The tuple `(fluid_rhs, wall_rhs)`. """ + from warnings import warn + warn( + "coupled_ns_heat_operator is deprecated and will disappear in Q3 2023. " + "Set up interface boundaries explicitly via " + ":func:`add_interface_boundaries` and include them when calling the " + "individual operators instead.", DeprecationWarning, stacklevel=2) + + if use_kappa_weighted_grad_flux_in_fluid is None: + warn( + "Default value of use_kappa_weighted_grad_flux_in_fluid has changed " + "from False to True as False is no longer allowed. Explicitly set it " + "to True to suppress this warning.", UserWarning, stacklevel=2) + use_kappa_weighted_grad_flux_in_fluid = True + elif not use_kappa_weighted_grad_flux_in_fluid: + warn( + "Setting use_kappa_weighted_grad_flux_in_fluid to True; False is no " + "longer allowed. Explicitly set it to True to suppress this warning.", + UserWarning, stacklevel=2) + if wall_penalty_amount is None: # FIXME: After verifying the form of the penalty term, figure out what value # makes sense to use as a default here @@ -1186,53 +1346,25 @@ def coupled_ns_heat_operator( as_dofdesc(bdtag).domain_tag: bdry for bdtag, bdry in wall_boundaries.items()} - # Pre-exchange kappa and temperature since we will need them in multiple steps - - kappa_inter_vol_tpairs = _kappa_inter_volume_trace_pairs( - dcoll, - gas_model, - fluid_dd, wall_dd, - fluid_state, wall_kappa) - - # FIXME: Maybe better to project CV and recompute T instead? - temperature_inter_vol_tpairs = _temperature_inter_volume_trace_pairs( - dcoll, - gas_model, - fluid_dd, wall_dd, - fluid_state, wall_temperature) - - # Construct boundaries for the fluid-wall interface; no temperature gradient + # Include boundaries for the fluid-wall interface; no temperature gradient # yet because we need to compute it - - fluid_interface_boundaries_no_grad, wall_interface_boundaries_no_grad = \ - get_interface_boundaries( - dcoll=dcoll, - gas_model=gas_model, - fluid_dd=fluid_dd, wall_dd=wall_dd, - fluid_state=fluid_state, wall_kappa=wall_kappa, - wall_temperature=wall_temperature, + fluid_all_boundaries_no_grad, wall_all_boundaries_no_grad = \ + add_interface_boundaries_no_grad( + dcoll, + gas_model, + fluid_dd, wall_dd, + fluid_state, wall_kappa, wall_temperature, + fluid_boundaries, wall_boundaries, interface_noslip=interface_noslip, - use_kappa_weighted_grad_flux_in_fluid=( - use_kappa_weighted_grad_flux_in_fluid), - _kappa_inter_vol_tpairs=kappa_inter_vol_tpairs, - _temperature_inter_vol_tpairs=temperature_inter_vol_tpairs) - - # Augment the domain boundaries with the interface boundaries (fluid only; - # needed for make_operator_fluid_states) - - fluid_all_boundaries_no_grad = {} - fluid_all_boundaries_no_grad.update(fluid_boundaries) - fluid_all_boundaries_no_grad.update(fluid_interface_boundaries_no_grad) + quadrature_tag=quadrature_tag) # Get the operator fluid states - fluid_operator_states_quad = make_operator_fluid_states( dcoll, fluid_state, gas_model, fluid_all_boundaries_no_grad, quadrature_tag, dd=fluid_dd, comm_tag=_FluidOpStatesTag, limiter_func=limiter_func) # Compute the temperature gradient for both subdomains - fluid_grad_temperature, wall_grad_temperature = coupled_grad_t_operator( dcoll, gas_model, @@ -1241,60 +1373,215 @@ def coupled_ns_heat_operator( fluid_state, wall_kappa, wall_temperature, time=time, interface_noslip=interface_noslip, - use_kappa_weighted_grad_flux_in_fluid=( - use_kappa_weighted_grad_flux_in_fluid), quadrature_tag=quadrature_tag, fluid_numerical_flux_func=fluid_gradient_numerical_flux_func, - _kappa_inter_vol_tpairs=kappa_inter_vol_tpairs, - _temperature_inter_vol_tpairs=temperature_inter_vol_tpairs, _fluid_operator_states_quad=fluid_operator_states_quad, - _fluid_interface_boundaries_no_grad=fluid_interface_boundaries_no_grad, - _wall_interface_boundaries_no_grad=wall_interface_boundaries_no_grad) + _fluid_all_boundaries_no_grad=fluid_all_boundaries_no_grad, + _wall_all_boundaries_no_grad=wall_all_boundaries_no_grad) - # Construct boundaries for the fluid-wall interface, now with the temperature + # Include boundaries for the fluid-wall interface, now with the temperature # gradient - - fluid_interface_boundaries, wall_interface_boundaries = \ - get_interface_boundaries( - dcoll=dcoll, - gas_model=gas_model, - fluid_dd=fluid_dd, wall_dd=wall_dd, - fluid_state=fluid_state, wall_kappa=wall_kappa, - wall_temperature=wall_temperature, - fluid_grad_temperature=fluid_grad_temperature, - wall_grad_temperature=wall_grad_temperature, + fluid_all_boundaries, wall_all_boundaries = \ + add_interface_boundaries( + dcoll, + gas_model, + fluid_dd, wall_dd, + fluid_state, wall_kappa, wall_temperature, + fluid_grad_temperature, wall_grad_temperature, + fluid_boundaries, wall_boundaries, interface_noslip=interface_noslip, - use_kappa_weighted_grad_flux_in_fluid=( - use_kappa_weighted_grad_flux_in_fluid), wall_penalty_amount=wall_penalty_amount, - _kappa_inter_vol_tpairs=kappa_inter_vol_tpairs, - _temperature_inter_vol_tpairs=temperature_inter_vol_tpairs) - - # Augment the domain boundaries with the interface boundaries - - fluid_all_boundaries = {} - fluid_all_boundaries.update(fluid_boundaries) - fluid_all_boundaries.update(fluid_interface_boundaries) - - wall_all_boundaries = {} - wall_all_boundaries.update(wall_boundaries) - wall_all_boundaries.update(wall_interface_boundaries) + quadrature_tag=quadrature_tag) # Compute the subdomain NS/diffusion operators using the augmented boundaries - ns_result = ns_operator( + my_ns_operator = partial(ns_operator, + inviscid_numerical_flux_func=inviscid_numerical_flux_func, + viscous_numerical_flux_func=viscous_numerical_flux_func) + ns_result = my_ns_operator( dcoll, gas_model, fluid_state, fluid_all_boundaries, time=time, quadrature_tag=quadrature_tag, dd=fluid_dd, - inviscid_numerical_flux_func=inviscid_numerical_flux_func, - viscous_numerical_flux_func=viscous_numerical_flux_func, return_gradients=return_gradients, operator_states_quad=fluid_operator_states_quad, grad_t=fluid_grad_temperature, comm_tag=_FluidOperatorTag) + diffusion_result = diffusion_operator( + dcoll, wall_kappa, wall_all_boundaries, wall_temperature, + penalty_amount=wall_penalty_amount, quadrature_tag=quadrature_tag, + return_grad_u=return_gradients, dd=wall_dd, grad_u=wall_grad_temperature, + comm_tag=_WallOperatorTag) + if return_gradients: fluid_rhs, fluid_grad_cv, fluid_grad_temperature = ns_result + wall_rhs, wall_grad_temperature = diffusion_result + return ( + fluid_rhs, wall_rhs, fluid_grad_cv, fluid_grad_temperature, + wall_grad_temperature) else: - fluid_rhs = ns_result + return ns_result, diffusion_result + + +def basic_coupled_ns_heat_operator( + dcoll, + gas_model, + fluid_dd, wall_dd, + fluid_boundaries, wall_boundaries, + fluid_state, wall_kappa, wall_temperature, + *, + time=0., + interface_noslip=True, + wall_penalty_amount=None, + quadrature_tag=DISCR_TAG_BASE, + limiter_func=None, + return_gradients=False): + r""" + Simple implementation of a thermally-coupled fluid/wall operator. + + Computes the RHS for a two-volume domain coupled by temperature and heat flux, + by augmenting *fluid_boundaries* and *wall_boundaries* with the boundaries for + the fluid-wall interface and calling the respective NS/diffusion operators. + + Parameters + ---------- + + dcoll: class:`~grudge.discretization.DiscretizationCollection` + + A discretization collection encapsulating the DG elements + + gas_model: :class:`~mirgecom.gas_model.GasModel` + + Physical gas model including equation of state, transport, + and kinetic properties as required by fluid state + + fluid_dd: :class:`grudge.dof_desc.DOFDesc` + + DOF descriptor for the fluid volume. + + wall_dd: :class:`grudge.dof_desc.DOFDesc` + + DOF descriptor for the wall volume. + + fluid_boundaries: + + Dictionary of boundary objects for the fluid subdomain, one for each + :class:`~grudge.dof_desc.BoundaryDomainTag` that represents a domain + boundary. + + wall_boundaries: + + Dictionary of boundary objects for the wall subdomain, one for each + :class:`~grudge.dof_desc.BoundaryDomainTag` that represents a domain + boundary. + + fluid_state: :class:`~mirgecom.gas_model.FluidState` + + Fluid state object with the conserved state and dependent + quantities for the fluid volume. + + wall_kappa: float or :class:`meshmode.dof_array.DOFArray` + + Thermal conductivity for the wall volume. + + wall_temperature: :class:`meshmode.dof_array.DOFArray` + + Temperature for the wall volume. + + time: + + Time + + interface_noslip: bool + + If `True`, interface boundaries on the fluid side will be treated as + no-slip walls. If `False` they will be treated as slip walls. + + wall_penalty_amount: float + + Coefficient $c$ for the interior penalty on the heat flux. See + :class:`~mirgecom.multiphysics.thermally_coupled_fluid_wall.InterfaceFluidBoundary` + for details. Not used if *interface_radiation* is `True`. + + quadrature_tag: + + An identifier denoting a particular quadrature discretization to use during + operator evaluations. + + limiter_func: + + Callable function to be passed to + :func:`~mirgecom.gas_model.make_operator_fluid_states` + that filters or limits the produced fluid states. This is used to keep + species mass fractions in physical and realizable states, for example. + + Returns + ------- + + The tuple `(fluid_rhs, wall_rhs)`. + """ + if wall_penalty_amount is None: + # FIXME: After verifying the form of the penalty term, figure out what value + # makes sense to use as a default here + wall_penalty_amount = 0.05 + + fluid_boundaries = { + as_dofdesc(bdtag).domain_tag: bdry + for bdtag, bdry in fluid_boundaries.items()} + wall_boundaries = { + as_dofdesc(bdtag).domain_tag: bdry + for bdtag, bdry in wall_boundaries.items()} + + # Include boundaries for the fluid-wall interface; no temperature gradient + # yet because we need to compute it + fluid_all_boundaries_no_grad, wall_all_boundaries_no_grad = \ + add_interface_boundaries_no_grad( + dcoll, + gas_model, + fluid_dd, wall_dd, + fluid_state, wall_kappa, wall_temperature, + fluid_boundaries, wall_boundaries, + interface_noslip=interface_noslip, + quadrature_tag=quadrature_tag) + + # Get the operator fluid states + fluid_operator_states_quad = make_operator_fluid_states( + dcoll, fluid_state, gas_model, fluid_all_boundaries_no_grad, + quadrature_tag, dd=fluid_dd, comm_tag=_FluidOpStatesTag, + limiter_func=limiter_func) + + # Compute the temperature gradient for both subdomains + fluid_grad_temperature = fluid_grad_t_operator( + dcoll, gas_model, fluid_all_boundaries_no_grad, fluid_state, + time=time, quadrature_tag=quadrature_tag, + dd=fluid_dd, operator_states_quad=fluid_operator_states_quad, + comm_tag=_FluidGradTag) + wall_grad_temperature = wall_grad_t_operator( + dcoll, wall_kappa, wall_all_boundaries_no_grad, wall_temperature, + quadrature_tag=quadrature_tag, dd=wall_dd, comm_tag=_WallGradTag) + + # Include boundaries for the fluid-wall interface, now with the temperature + # gradient + fluid_all_boundaries, wall_all_boundaries = \ + add_interface_boundaries( + dcoll, + gas_model, + fluid_dd, wall_dd, + fluid_state, wall_kappa, wall_temperature, + fluid_grad_temperature, wall_grad_temperature, + fluid_boundaries, wall_boundaries, + interface_noslip=interface_noslip, + wall_penalty_amount=wall_penalty_amount, + quadrature_tag=quadrature_tag) + + # Compute the subdomain NS/diffusion operators using the augmented boundaries + + my_ns_operator = partial(ns_operator, + viscous_numerical_flux_func=viscous_facial_flux_harmonic) + ns_result = my_ns_operator( + dcoll, gas_model, fluid_state, fluid_all_boundaries, + time=time, quadrature_tag=quadrature_tag, dd=fluid_dd, + return_gradients=return_gradients, + operator_states_quad=fluid_operator_states_quad, + grad_t=fluid_grad_temperature, comm_tag=_FluidOperatorTag) diffusion_result = diffusion_operator( dcoll, wall_kappa, wall_all_boundaries, wall_temperature, @@ -1303,13 +1590,10 @@ def coupled_ns_heat_operator( comm_tag=_WallOperatorTag) if return_gradients: + fluid_rhs, fluid_grad_cv, fluid_grad_temperature = ns_result wall_rhs, wall_grad_temperature = diffusion_result - else: - wall_rhs = diffusion_result - - if return_gradients: return ( fluid_rhs, wall_rhs, fluid_grad_cv, fluid_grad_temperature, wall_grad_temperature) else: - return fluid_rhs, wall_rhs + return ns_result, diffusion_result diff --git a/mirgecom/utils.py b/mirgecom/utils.py index 8bc94e3b4..83123b436 100644 --- a/mirgecom/utils.py +++ b/mirgecom/utils.py @@ -4,6 +4,7 @@ .. autofunction:: asdict_shallow .. autofunction:: force_evaluation .. autofunction:: normalize_boundaries +.. autofunction:: project_from_base """ __copyright__ = """ @@ -130,3 +131,17 @@ def normalize_boundaries(boundaries): return { as_dofdesc(key).domain_tag: bdry for key, bdry in boundaries.items()} + + +def project_from_base(dcoll, tgt_dd, field): + """Project *field* from *DISCR_TAG_BASE* to the same discr. as *tgt_dd*.""" + from grudge.dof_desc import DISCR_TAG_BASE, as_dofdesc + from grudge.op import project + + tgt_dd = as_dofdesc(tgt_dd) + + if tgt_dd.discretization_tag is not DISCR_TAG_BASE: + tgt_dd_base = tgt_dd.with_discr_tag(DISCR_TAG_BASE) + return project(dcoll, tgt_dd_base, tgt_dd, field) + else: + return field diff --git a/test/test_multiphysics.py b/test/test_multiphysics.py index c9e7b8483..5dcaf958f 100644 --- a/test/test_multiphysics.py +++ b/test/test_multiphysics.py @@ -51,7 +51,7 @@ IsothermalWallBoundary, ) from mirgecom.multiphysics.thermally_coupled_fluid_wall import ( - coupled_ns_heat_operator + basic_coupled_ns_heat_operator as coupled_ns_heat_operator, ) from meshmode.array_context import ( # noqa pytest_generate_tests_for_pyopencl_array_context