diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py index f8e411387..392ce95ad 100644 --- a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py +++ b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py @@ -4,7 +4,7 @@ Couples a fluid subdomain governed by the compressible Navier-Stokes equations (:mod:`mirgecom.navierstokes`) with a wall subdomain governed by the heat equation (:mod:`mirgecom.diffusion`) through temperature and heat flux. This -radiation can optionally include a sink term representing emitted radiation. +coupling can optionally include a sink term representing emitted radiation. In the non-radiating case, coupling enforces continuity of temperature and heat flux .. math:: @@ -234,9 +234,15 @@ def viscous_divergence_flux( lengthscales_minus = project_from_base( dcoll, dd_bdry, self._lengthscales_minus) - tau = ( - self._penalty_amount * state_bc.thermal_conductivity - / lengthscales_minus) + if isinstance(state_bc.thermal_conductivity, np.ndarray): + # orthotropic materials + actx = state_minus.array_context + normal = actx.thaw(dcoll.normal(dd_bdry)) + kappa_bc = np.dot(normal, state_bc.thermal_conductivity*normal) + else: + kappa_bc = state_bc.thermal_conductivity + + tau = self._penalty_amount * kappa_bc / lengthscales_minus t_minus = state_minus.temperature t_plus = self.temperature_plus( @@ -256,12 +262,28 @@ def __init__( self._t_plus = t_plus self._grad_t_plus = grad_t_plus - def kappa_plus(self, dcoll, dd_bdry): + def _normal_kappa_plus(self, dcoll, dd_bdry): + # project kappa plus in case of overintegration + if isinstance(self._kappa_plus, np.ndarray): + # orthotropic materials + actx = self._t_plus.array_context + normal = actx.thaw(dcoll.normal(dd_bdry)) + kappa_plus = project_from_base(dcoll, dd_bdry, self._kappa_plus) + return np.dot(normal, kappa_plus*normal) return project_from_base(dcoll, dd_bdry, self._kappa_plus) + def _normal_kappa_minus(self, dcoll, dd_bdry, kappa): + # state minus is already in the quadrature domain + if isinstance(kappa, np.ndarray): + # orthotropic materials + actx = self._t_plus.array_context + normal = actx.thaw(dcoll.normal(dd_bdry)) + return np.dot(normal, kappa*normal) + return kappa + def kappa_bc(self, dcoll, dd_bdry, kappa_minus): - kappa_plus = project_from_base(dcoll, dd_bdry, self._kappa_plus) - return harmonic_mean(kappa_minus, kappa_plus) + return harmonic_mean(kappa_minus, + project_from_base(dcoll, dd_bdry, self._kappa_plus)) def temperature_plus(self, dcoll, dd_bdry): return project_from_base(dcoll, dd_bdry, self._t_plus) @@ -269,7 +291,9 @@ def temperature_plus(self, dcoll, dd_bdry): def temperature_bc(self, dcoll, dd_bdry, kappa_minus, 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_plus = self._normal_kappa_plus(dcoll, dd_bdry) + kappa_minus = self._normal_kappa_minus(dcoll, dd_bdry, + kappa_minus + t_minus*0.0) kappa_sum = actx.np.where( actx.np.greater(kappa_minus + kappa_plus, 0*kappa_minus), kappa_minus + kappa_plus, @@ -366,9 +390,7 @@ def state_bc( cv_minus = state_minus.cv - kappa_minus = ( - # Make sure it has an array context - state_minus.tv.thermal_conductivity + 0*state_minus.mass_density) + kappa_minus = state_minus.tv.thermal_conductivity # Grab a unit normal to the boundary nhat = actx.thaw(dcoll.normal(dd_bdry)) @@ -425,9 +447,7 @@ def temperature_plus( return self._thermally_coupled.temperature_plus(dcoll, dd_bdry) def temperature_bc(self, dcoll, dd_bdry, state_minus, **kwargs): # noqa: D102 - kappa_minus = ( - # Make sure it has an array context - state_minus.tv.thermal_conductivity + 0*state_minus.mass_density) + kappa_minus = state_minus.tv.thermal_conductivity return self._thermally_coupled.temperature_bc( dcoll, dd_bdry, kappa_minus, state_minus.temperature) @@ -506,9 +526,7 @@ def state_bc( self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): # noqa: D102 dd_bdry = as_dofdesc(dd_bdry) - kappa_minus = ( - # Make sure it has an array context - state_minus.tv.thermal_conductivity + 0*state_minus.mass_density) + kappa_minus = state_minus.tv.thermal_conductivity mom_bc = self._no_slip.momentum_bc(state_minus.momentum_density) @@ -544,9 +562,7 @@ def temperature_plus( return self._thermally_coupled.temperature_plus(dcoll, dd_bdry) def temperature_bc(self, dcoll, dd_bdry, state_minus, **kwargs): # noqa: D102 - kappa_minus = ( - # Make sure it has an array context - state_minus.tv.thermal_conductivity + 0*state_minus.mass_density) + kappa_minus = state_minus.tv.thermal_conductivity return self._thermally_coupled.temperature_bc( dcoll, dd_bdry, kappa_minus, state_minus.temperature) @@ -599,6 +615,7 @@ def get_grad_flux( normal = actx.thaw(dcoll.normal(dd_bdry)) kappa_plus = project_from_base(dcoll, dd_bdry, self.kappa_plus) + kappa_tpair = TracePair( dd_bdry, interior=kappa_minus, exterior=kappa_plus) @@ -619,6 +636,7 @@ def get_diffusion_flux( normal = actx.thaw(dcoll.normal(dd_bdry)) kappa_plus = project_from_base(dcoll, dd_bdry, self.kappa_plus) + kappa_tpair = TracePair( dd_bdry, interior=kappa_minus, exterior=kappa_plus) @@ -1596,7 +1614,8 @@ def basic_coupled_ns_heat_operator( quadrature_tag=DISCR_TAG_BASE, limiter_func=None, return_gradients=False, - use_esdg=False): + use_esdg=False, + inviscid_terms_on=True): r""" Simple implementation of a thermally-coupled fluid/wall operator. @@ -1776,7 +1795,7 @@ def basic_coupled_ns_heat_operator( my_ns_operator = partial(ns_operator, viscous_numerical_flux_func=viscous_facial_flux_harmonic, - use_esdg=use_esdg) + use_esdg=use_esdg, inviscid_terms_on=inviscid_terms_on) ns_result = my_ns_operator( dcoll, gas_model, fluid_state, fluid_all_boundaries, time=time, quadrature_tag=quadrature_tag, dd=fluid_dd, diff --git a/test/test_diffusion.py b/test/test_diffusion.py index f43bf32f6..08e4d7e10 100644 --- a/test/test_diffusion.py +++ b/test/test_diffusion.py @@ -731,8 +731,7 @@ def test_orthotropic_diffusion(actx_factory): dcoll = create_discretization_collection(actx, mesh, order) nodes = actx.thaw(dcoll.nodes()) - zeros = actx.np.zeros_like(nodes[0]) - kappa = make_obj_array([kappa0 + zeros, kappa1 + zeros]) + kappa = make_obj_array([kappa0, kappa1]) u = 30.0*nodes[0] + 60.0*nodes[1] # exercise Neumann BC @@ -753,8 +752,8 @@ def test_orthotropic_diffusion(actx_factory): assert err_grad_y < 1.e-9 diff_flux = diffusion_flux(kappa, grad_u) - flux_x = -(kappa0 + zeros)*grad_u[0] - flux_y = -(kappa1 + zeros)*grad_u[1] + flux_x = -kappa0*grad_u[0] + flux_y = -kappa1*grad_u[1] err_flux_x = actx.to_numpy(op.norm(dcoll, diff_flux[0] - flux_x, np.inf)) err_flux_y = actx.to_numpy(op.norm(dcoll, diff_flux[1] - flux_y, np.inf)) assert err_flux_x < 1.e-9 @@ -785,8 +784,8 @@ def make_dirichlet_bc(btag): assert err_grad_y < 1.e-9 diff_flux = diffusion_flux(kappa, grad_u) - flux_x = -(kappa0 + zeros)*grad_u[0] - flux_y = -(kappa1 + zeros)*grad_u[1] + flux_x = -kappa0*grad_u[0] + flux_y = -kappa1*grad_u[1] err_flux_x = actx.to_numpy(op.norm(dcoll, diff_flux[0] - flux_x, np.inf)) err_flux_y = actx.to_numpy(op.norm(dcoll, diff_flux[1] - flux_y, np.inf)) assert err_flux_x < 1.e-9 diff --git a/test/test_multiphysics.py b/test/test_multiphysics.py index 12ad33caa..48ea892cb 100644 --- a/test/test_multiphysics.py +++ b/test/test_multiphysics.py @@ -26,11 +26,7 @@ import pyopencl.array as cla # noqa import pyopencl.clmath as clmath # noqa from pytools.obj_array import make_obj_array -import pymbolic as pmbl import grudge.op as op -from mirgecom.symbolic import ( - grad as sym_grad, - evaluate) from mirgecom.simutil import ( max_component_norm, get_box_mesh @@ -223,7 +219,6 @@ def test_thermally_coupled_fluid_wall( wall_density = 10*fluid_density wall_heat_capacity = fluid_heat_capacity wall_kappa = 10*fluid_kappa - base_wall_temp = 600 fluid_boundaries = { @@ -246,6 +241,7 @@ def test_thermally_coupled_fluid_wall( interface_flux = ( -fluid_kappa * wall_kappa / (fluid_kappa + wall_kappa) * (base_fluid_temp - base_wall_temp)) + fluid_alpha = fluid_kappa/(fluid_density * fluid_heat_capacity) wall_alpha = wall_kappa/(wall_density * wall_heat_capacity) @@ -262,7 +258,7 @@ def perturb_func(alpha, x, t): # This perturbation function is nonzero at the interface, so the two alphas # need to be the same (otherwise the perturbations will decay at different # rates and a discontinuity will form) - assert abs(fluid_alpha - wall_alpha) < 1e-12 + assert abs(fluid_alpha - np.max(actx.to_numpy(wall_alpha))) < 1e-12 fluid_perturb_func = partial(perturb_func, fluid_alpha) wall_perturb_func = partial(perturb_func, wall_alpha) @@ -297,14 +293,6 @@ def wall_func(x, t): ("temp", wall_temp), ]) - # Add a source term to the momentum equations to cancel out the pressure term - sym_fluid_temp = fluid_func(pmbl.var("x"), pmbl.var("t")) - sym_fluid_pressure = fluid_density * r * sym_fluid_temp - sym_momentum_source = sym_grad(2, sym_fluid_pressure) - - def momentum_source_func(x, t): - return evaluate(sym_momentum_source, x=x, t=t) - def get_rhs(t, state): fluid_state = make_fluid_state(cv=state[0], gas_model=gas_model) wall_temperature = state[1] @@ -315,10 +303,8 @@ def get_rhs(t, state): fluid_boundaries, wall_boundaries, fluid_state, wall_kappa, wall_temperature, time=t, - quadrature_tag=quadrature_tag) - fluid_rhs = replace( - fluid_rhs, - momentum=fluid_rhs.momentum + momentum_source_func(fluid_nodes, t)) + quadrature_tag=quadrature_tag, + inviscid_terms_on=False) wall_rhs = wall_energy_rhs / (wall_density * wall_heat_capacity) return make_obj_array([fluid_rhs, wall_rhs]) @@ -410,7 +396,6 @@ def cv_from_temp(temp): expected_fluid_temp = fluid_func(fluid_nodes, t) expected_wall_temp = wall_func(wall_nodes, t) rhs = get_rhs(t, state) - momentum_source = momentum_source_func(fluid_nodes, t) viz_fluid.write_vtk_file( "thermally_coupled_accuracy_" f"{viz_suffix}_fluid_{step}.vtu", [ @@ -418,7 +403,6 @@ def cv_from_temp(temp): ("dv", fluid_state.dv), ("expected_temp", expected_fluid_temp), ("rhs", rhs[0]), - ("momentum_source", momentum_source), ]) viz_wall.write_vtk_file( "thermally_coupled_accuracy_" @@ -558,10 +542,98 @@ def test_thermally_coupled_fluid_wall_with_radiation( fluid_rhs = fluid_rhs.energy solid_rhs = wall_energy_rhs/(wall_cp*wall_rho) + # FIXME check convergence instead? assert actx.to_numpy(op.norm(dcoll, fluid_rhs, np.inf)) < 1e-4 assert actx.to_numpy(op.norm(dcoll, solid_rhs, np.inf)) < 1e-4 +@pytest.mark.parametrize("use_overintegration", [False, True]) +@pytest.mark.parametrize("use_noslip", [False, True]) +@pytest.mark.parametrize("use_radiation", [False, True]) +def test_orthotropic_flux( + actx_factory, use_overintegration, use_radiation, use_noslip, + visualize=False): + """Check the RHS shape for orthotropic kappa cases.""" + actx = actx_factory() + + dim = 2 + order = 3 + nelems = 48 + + global_mesh = get_box_mesh(dim, -2, 1, nelems) + + mgrp, = global_mesh.groups + x = global_mesh.vertices[0, mgrp.vertex_indices] + x_elem_avg = np.sum(x, axis=1)/x.shape[0] + volume_to_elements = { + "Fluid": np.where(x_elem_avg > 0)[0], + "Solid": np.where(x_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=order, quadrature_order=2*order+1) + + quadrature_tag = DISCR_TAG_QUAD if use_overintegration else None + + dd_vol_fluid = DOFDesc(VolumeDomainTag("Fluid"), DISCR_TAG_BASE) + dd_vol_solid = DOFDesc(VolumeDomainTag("Solid"), DISCR_TAG_BASE) + + fluid_nodes = actx.thaw(dcoll.nodes(dd=dd_vol_fluid)) + solid_nodes = actx.thaw(dcoll.nodes(dd=dd_vol_solid)) + + eos = IdealSingleGas() + transport = SimpleTransport(viscosity=0.0, thermal_conductivity=10.0) + gas_model = GasModel(eos=eos, transport=transport) + + # Fluid cv + + mom_x = 0.0*fluid_nodes[0] + mom_y = 0.0*fluid_nodes[0] + momentum = make_obj_array([mom_x, mom_y]) + + temperature = 2.0 + 0.0*fluid_nodes[0] + mass = 1.0 + 0.0*fluid_nodes[0] + energy = mass*eos.get_internal_energy(temperature) + + fluid_cv = make_conserved(dim=dim, mass=mass, energy=energy, + momentum=momentum) + + fluid_state = make_fluid_state(cv=fluid_cv, gas_model=gas_model) + + # Made-up wall material + wall_kappa = make_obj_array([1.0, 1.0]) + wall_emissivity = 1.0 + wall_temperature = 1.0 + 0.0*solid_nodes[0] + + fluid_boundaries = { + dd_vol_fluid.trace("-2").domain_tag: AdiabaticNoslipWallBoundary(), + dd_vol_fluid.trace("+2").domain_tag: AdiabaticNoslipWallBoundary(), + dd_vol_fluid.trace("+1").domain_tag: + IsothermalWallBoundary(wall_temperature=2.0), + } + + solid_boundaries = { + dd_vol_solid.trace("-2").domain_tag: NeumannDiffusionBoundary(0.), + dd_vol_solid.trace("+2").domain_tag: NeumannDiffusionBoundary(0.), + dd_vol_solid.trace("-1").domain_tag: DirichletDiffusionBoundary(1.0), + } + + fluid_rhs, wall_energy_rhs = coupled_ns_heat_operator( + dcoll, gas_model, dd_vol_fluid, dd_vol_solid, fluid_boundaries, + solid_boundaries, fluid_state, wall_kappa, wall_temperature, + time=0.0, quadrature_tag=quadrature_tag, interface_noslip=use_noslip, + interface_radiation=use_radiation, sigma=2.0, + ambient_temperature=0.0, wall_emissivity=wall_emissivity, + inviscid_terms_on=False + ) + + from meshmode.dof_array import DOFArray + assert isinstance(fluid_rhs.energy, DOFArray) + assert isinstance(wall_energy_rhs, DOFArray) + + if __name__ == "__main__": import sys if len(sys.argv) > 1: