Fix thermal interface for orthotropic cases (#1009)

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,20 +262,38 @@ 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)
 
     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,

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

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,15 +396,13 @@ 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", [
                         ("cv", fluid_state.cv),
                         ("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: