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,7 +262,7 @@ def __init__(
         self._t_plus = t_plus
         self._grad_t_plus = grad_t_plus
 
-    def proj_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
@@ -266,7 +272,7 @@ def proj_kappa_plus(self, dcoll, dd_bdry):
             return np.dot(normal, kappa_plus*normal)
         return project_from_base(dcoll, dd_bdry, self._kappa_plus)
 
-    def proj_kappa_minus(self, dcoll, dd_bdry, kappa):
+    def _normal_kappa_minus(self, dcoll, dd_bdry, kappa):
         # state minus is already in the quadrature domain
         if isinstance(kappa, np.ndarray):
             # orthotropic materials
@@ -276,18 +282,18 @@ def proj_kappa_minus(self, dcoll, dd_bdry, kappa):
         return kappa
 
     def kappa_bc(self, dcoll, dd_bdry, kappa_minus):
-        return harmonic_mean(self.proj_kappa_minus(dcoll, dd_bdry, kappa_minus),
-                             self.proj_kappa_plus(dcoll, dd_bdry))
+        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 = self.proj_kappa_plus(dcoll, dd_bdry)
-        kappa_minus = self.proj_kappa_minus(dcoll, dd_bdry,
-                                            kappa_minus + t_minus*0.0)
+        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,

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

@@ -224,13 +224,7 @@ def test_thermally_coupled_fluid_wall(
         # Made-up wall material
         wall_density = 10*fluid_density
         wall_heat_capacity = fluid_heat_capacity
-        if orthotropic_kappa:
-            wall_kappa = make_obj_array([10*fluid_kappa, 10*fluid_kappa])
-            _wall_kappa = wall_kappa[0]  # dummy variable to simplify the test
-        else:
-            wall_kappa = 10*fluid_kappa
-            _wall_kappa = wall_kappa  # dummy variable to simplify the test
-
+        wall_kappa = 10*fluid_kappa
         base_wall_temp = 600
 
         fluid_boundaries = {
@@ -248,20 +242,20 @@ def test_thermally_coupled_fluid_wall(
         }
 
         interface_temp = (
-            (fluid_kappa * base_fluid_temp + _wall_kappa * base_wall_temp)
-            / (fluid_kappa + _wall_kappa))
+            (fluid_kappa * base_fluid_temp + wall_kappa * base_wall_temp)
+            / (fluid_kappa + wall_kappa))
         interface_flux = (
-            -fluid_kappa * _wall_kappa / (fluid_kappa + _wall_kappa)
+            -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)
+        wall_alpha = wall_kappa/(wall_density * wall_heat_capacity)
 
         def steady_func(kappa, x, t):
             return interface_temp - interface_flux/kappa * x[1]
 
         fluid_steady_func = partial(steady_func, fluid_kappa)
-        wall_steady_func = partial(steady_func, _wall_kappa)
+        wall_steady_func = partial(steady_func, wall_kappa)
 
         def perturb_func(alpha, x, t):
             w = 1.5 * np.pi
@@ -556,10 +550,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: