From 89218e5a21a50e6861701b67204bad995a64234c Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Sun, 3 Mar 2024 08:53:29 -0600
Subject: [PATCH 01/15] Force kappa be a DOFArray early on

---
 mirgecom/diffusion.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py
index ab9759084..ec6abd34b 100644
--- a/mirgecom/diffusion.py
+++ b/mirgecom/diffusion.py
@@ -695,6 +695,9 @@ def grad_operator(
     dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag)
     dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL)
 
+    # Make sure it has an array context
+    kappa = kappa + dcoll.zeros(array_context=actx, dd=dd_vol)
+
     if kappa_tpairs is None:
         kappa_tpairs = interior_trace_pairs(
             dcoll, kappa, volume_dd=dd_vol, comm_tag=(_DiffusionKappaTag, comm_tag))
@@ -829,6 +832,9 @@ def diffusion_operator(
     dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag)
     dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL)
 
+    # Make sure it has an array context
+    kappa = kappa + dcoll.zeros(array_context=actx, dd=dd_vol)
+
     kappa_tpairs = interior_trace_pairs(
         dcoll, kappa, volume_dd=dd_vol, comm_tag=(_DiffusionKappaTag, comm_tag))
 

From f604edf371fefec6277f4a52c6b2fd16f8d6711d Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Sun, 3 Mar 2024 10:00:25 -0600
Subject: [PATCH 02/15] Fix normal projection case for orthotropic
 condutivities

---
 .../thermally_coupled_fluid_wall.py           | 65 ++++++++++++++-----
 1 file changed, 47 insertions(+), 18 deletions(-)

diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
index f8e411387..4aa463e16 100644
--- a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
+++ b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
@@ -256,12 +256,28 @@ def __init__(
         self._t_plus = t_plus
         self._grad_t_plus = grad_t_plus
 
-    def kappa_plus(self, dcoll, dd_bdry):
+    def proj_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 proj_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(self.proj_kappa_minus(dcoll, dd_bdry, kappa_minus),
+                             self.proj_kappa_plus(dcoll, dd_bdry))
 
     def temperature_plus(self, dcoll, dd_bdry):
         return project_from_base(dcoll, dd_bdry, self._t_plus)
@@ -269,7 +285,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.proj_kappa_plus(dcoll, dd_bdry)
+        kappa_minus = self.proj_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 +384,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 +441,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 +520,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 +556,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 +609,13 @@ def get_grad_flux(
         normal = actx.thaw(dcoll.normal(dd_bdry))
 
         kappa_plus = project_from_base(dcoll, dd_bdry, self.kappa_plus)
+
+#        # orthotropic materials
+#        if isinstance(kappa_minus, np.ndarray):
+#            kappa_minus = np.dot(normal, kappa_minus*normal)
+#        if isinstance(kappa_plus, np.ndarray):
+#            kappa_plus = np.dot(normal, kappa_plus*normal)
+
         kappa_tpair = TracePair(
             dd_bdry, interior=kappa_minus, exterior=kappa_plus)
 
@@ -619,6 +636,13 @@ def get_diffusion_flux(
         normal = actx.thaw(dcoll.normal(dd_bdry))
 
         kappa_plus = project_from_base(dcoll, dd_bdry, self.kappa_plus)
+
+#        # orthotropic materials
+#        if isinstance(kappa_minus, np.ndarray):
+#            kappa_minus = np.dot(normal, kappa_minus*normal)
+#        if isinstance(kappa_plus, np.ndarray):
+#            kappa_plus = np.dot(normal, kappa_plus*normal)
+
         kappa_tpair = TracePair(
             dd_bdry, interior=kappa_minus, exterior=kappa_plus)
 
@@ -700,6 +724,10 @@ def get_grad_flux(
         actx = u_minus.array_context
         normal = actx.thaw(dcoll.normal(dd_bdry))
 
+#        # orthotropic materials
+#        if isinstance(kappa_minus, np.ndarray):
+#            kappa_minus = np.dot(normal, kappa_minus*normal)
+
         kappa_tpair = TracePair(
             dd_bdry, interior=kappa_minus, exterior=kappa_minus)
         u_tpair = TracePair(dd_bdry, interior=u_minus, exterior=u_minus)
@@ -1596,7 +1624,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 +1805,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,

From 333659b2a7d2cb536ab1199894a1a940d3070a68 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Sun, 3 Mar 2024 10:00:43 -0600
Subject: [PATCH 03/15] Add orthotropic case to test_multiphysics

---
 test/test_multiphysics.py | 56 +++++++++++++++++++--------------------
 1 file changed, 27 insertions(+), 29 deletions(-)

diff --git a/test/test_multiphysics.py b/test/test_multiphysics.py
index 12ad33caa..b00fe4e37 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
@@ -65,6 +61,10 @@
 logger = logging.getLogger(__name__)
 
 
+import os  # noqa
+os.environ["CUDA_VISIBLE_DEVICES"] = "-1"  # noqa
+
+
 @pytest.mark.parametrize("order", [1, 2, 3])
 def test_independent_volumes(actx_factory, order, visualize=False):
     """Check multi-volume machinery by setting up two independent volumes."""
@@ -154,8 +154,10 @@ def get_rhs(t, u):
 
 @pytest.mark.parametrize("order", [2, 3])
 @pytest.mark.parametrize("use_overintegration", [False, True])
+@pytest.mark.parametrize("orthotropic_kappa", [False, True])
 def test_thermally_coupled_fluid_wall(
-        actx_factory, order, use_overintegration, visualize=False):
+        actx_factory, order, use_overintegration, orthotropic_kappa,
+        visualize=False):
     """Check the thermally-coupled fluid/wall interface."""
     actx = actx_factory()
 
@@ -222,7 +224,12 @@ def test_thermally_coupled_fluid_wall(
         # Made-up wall material
         wall_density = 10*fluid_density
         wall_heat_capacity = fluid_heat_capacity
-        wall_kappa = 10*fluid_kappa
+        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
 
         base_wall_temp = 600
 
@@ -241,19 +248,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
@@ -262,7 +270,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 +305,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 +315,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 +408,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 +415,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_"
@@ -461,10 +457,12 @@ def cv_from_temp(temp):
         or eoc_rec_wall.max_error() < 1e-11)
 
 
-@pytest.mark.parametrize("order", [1, 3])
+@pytest.mark.parametrize("order", [2, 3])
 @pytest.mark.parametrize("use_overintegration", [False, True])
+@pytest.mark.parametrize("orthotropic_kappa", [False, True])
 def test_thermally_coupled_fluid_wall_with_radiation(
-        actx_factory, order, use_overintegration, visualize=False):
+        actx_factory, order, use_overintegration, orthotropic_kappa,
+        visualize=False):
     """Check the thermally-coupled fluid/wall interface with radiation.
 
     Analytic solution prescribed as initial condition, then the RHS is assessed
@@ -523,7 +521,7 @@ def test_thermally_coupled_fluid_wall_with_radiation(
     # Made-up wall material
     wall_rho = 1.0
     wall_cp = 1000.0
-    wall_kappa = 1.0
+    wall_kappa = make_obj_array([1.0, 1.0]) if orthotropic_kappa else 1.0
     wall_emissivity = 1.0
 
     base_fluid_temp = 2.0

From 9a148d3e90e164b5919c919eebdd2a4d6f0d26ce Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Sun, 3 Mar 2024 10:06:17 -0600
Subject: [PATCH 04/15] Remove useless normal projection in diffusion operators

---
 .../multiphysics/thermally_coupled_fluid_wall.py | 16 ----------------
 1 file changed, 16 deletions(-)

diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
index 4aa463e16..0636c152e 100644
--- a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
+++ b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
@@ -610,12 +610,6 @@ def get_grad_flux(
 
         kappa_plus = project_from_base(dcoll, dd_bdry, self.kappa_plus)
 
-#        # orthotropic materials
-#        if isinstance(kappa_minus, np.ndarray):
-#            kappa_minus = np.dot(normal, kappa_minus*normal)
-#        if isinstance(kappa_plus, np.ndarray):
-#            kappa_plus = np.dot(normal, kappa_plus*normal)
-
         kappa_tpair = TracePair(
             dd_bdry, interior=kappa_minus, exterior=kappa_plus)
 
@@ -637,12 +631,6 @@ def get_diffusion_flux(
 
         kappa_plus = project_from_base(dcoll, dd_bdry, self.kappa_plus)
 
-#        # orthotropic materials
-#        if isinstance(kappa_minus, np.ndarray):
-#            kappa_minus = np.dot(normal, kappa_minus*normal)
-#        if isinstance(kappa_plus, np.ndarray):
-#            kappa_plus = np.dot(normal, kappa_plus*normal)
-
         kappa_tpair = TracePair(
             dd_bdry, interior=kappa_minus, exterior=kappa_plus)
 
@@ -724,10 +712,6 @@ def get_grad_flux(
         actx = u_minus.array_context
         normal = actx.thaw(dcoll.normal(dd_bdry))
 
-#        # orthotropic materials
-#        if isinstance(kappa_minus, np.ndarray):
-#            kappa_minus = np.dot(normal, kappa_minus*normal)
-
         kappa_tpair = TracePair(
             dd_bdry, interior=kappa_minus, exterior=kappa_minus)
         u_tpair = TracePair(dd_bdry, interior=u_minus, exterior=u_minus)

From ee19d73855ec281b7e8ef70fae82fe0124e373ad Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 4 Mar 2024 10:29:40 -0600
Subject: [PATCH 05/15] order back to 1; remove CUDA flag

---
 test/test_multiphysics.py | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/test/test_multiphysics.py b/test/test_multiphysics.py
index b00fe4e37..b3e03b1c6 100644
--- a/test/test_multiphysics.py
+++ b/test/test_multiphysics.py
@@ -61,10 +61,6 @@
 logger = logging.getLogger(__name__)
 
 
-import os  # noqa
-os.environ["CUDA_VISIBLE_DEVICES"] = "-1"  # noqa
-
-
 @pytest.mark.parametrize("order", [1, 2, 3])
 def test_independent_volumes(actx_factory, order, visualize=False):
     """Check multi-volume machinery by setting up two independent volumes."""
@@ -457,7 +453,7 @@ def cv_from_temp(temp):
         or eoc_rec_wall.max_error() < 1e-11)
 
 
-@pytest.mark.parametrize("order", [2, 3])
+@pytest.mark.parametrize("order", [1, 3])
 @pytest.mark.parametrize("use_overintegration", [False, True])
 @pytest.mark.parametrize("orthotropic_kappa", [False, True])
 def test_thermally_coupled_fluid_wall_with_radiation(

From b08199880eefb7b11d3df0e32c596bddefeaeb43 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 5 Mar 2024 17:26:28 -0600
Subject: [PATCH 06/15] Dot the kappa_bc in the right place

---
 .../multiphysics/thermally_coupled_fluid_wall.py  | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
index 0636c152e..555483450 100644
--- a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
+++ b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
@@ -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 = self._t_plus.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(
@@ -276,8 +282,7 @@ 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, self.proj_kappa_plus(dcoll, dd_bdry))
 
     def temperature_plus(self, dcoll, dd_bdry):
         return project_from_base(dcoll, dd_bdry, self._t_plus)

From 7227285e649e941e08a77641fa2c9532fd7b4cbc Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 5 Mar 2024 17:32:07 -0600
Subject: [PATCH 07/15] pylint

---
 mirgecom/multiphysics/thermally_coupled_fluid_wall.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
index 555483450..4910e2754 100644
--- a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
+++ b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
@@ -236,7 +236,7 @@ def viscous_divergence_flux(
 
         if isinstance(state_bc.thermal_conductivity, np.ndarray):
             # orthotropic materials
-            actx = self._t_plus.array_context
+            actx = state_minus.array_context
             normal = actx.thaw(dcoll.normal(dd_bdry))
             kappa_bc = np.dot(normal, state_bc.thermal_conductivity*normal)
         else:

From 62b8bc009b6953925fd435969eb4b9a175ac6720 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 5 Mar 2024 20:16:27 -0600
Subject: [PATCH 08/15] remove zeros from test

---
 test/test_diffusion.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/test_diffusion.py b/test/test_diffusion.py
index f43bf32f6..244b3ab88 100644
--- a/test/test_diffusion.py
+++ b/test/test_diffusion.py
@@ -732,7 +732,7 @@ def test_orthotropic_diffusion(actx_factory):
     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 +753,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 +785,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

From b48506ec80c6f01559ab1932181606b22ebbf958 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 6 Mar 2024 10:38:33 -0600
Subject: [PATCH 09/15] Fix kappa_plus in kappa_bc; modify tests and now only
 asses the energy rhs shape

---
 .../thermally_coupled_fluid_wall.py           |  18 +--
 test/test_diffusion.py                        |   1 -
 test/test_multiphysics.py                     | 109 +++++++++++++++---
 3 files changed, 105 insertions(+), 23 deletions(-)

diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
index 4910e2754..ddb0184be 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::
@@ -82,6 +82,7 @@
 )
 import grudge.op as op
 
+from pytools.obj_array import make_obj_array
 from mirgecom.math import harmonic_mean
 from mirgecom.boundary import (
     MengaldoBoundaryCondition,
@@ -262,7 +263,7 @@ def __init__(
         self._t_plus = t_plus
         self._grad_t_plus = grad_t_plus
 
-    def proj_kappa_plus(self, dcoll, dd_bdry):
+    def get_normal_kappa_plus(self, dcoll, dd_bdry):
         # project kappa plus in case of overintegration
         if isinstance(self._kappa_plus, np.ndarray):
             # orthotropic materials
@@ -272,7 +273,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 get_normal_kappa_minus(self, dcoll, dd_bdry, kappa):
         # state minus is already in the quadrature domain
         if isinstance(kappa, np.ndarray):
             # orthotropic materials
@@ -282,7 +283,10 @@ def proj_kappa_minus(self, dcoll, dd_bdry, kappa):
         return kappa
 
     def kappa_bc(self, dcoll, dd_bdry, kappa_minus):
-        return harmonic_mean(kappa_minus, self.proj_kappa_plus(dcoll, dd_bdry))
+        kappa_plus = make_obj_array([
+            project_from_base(dcoll, dd_bdry, self._kappa_plus[i])
+            for i in range(self._kappa_plus.shape[0])])
+        return harmonic_mean(kappa_minus, kappa_plus)
 
     def temperature_plus(self, dcoll, dd_bdry):
         return project_from_base(dcoll, dd_bdry, self._t_plus)
@@ -290,9 +294,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 = 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.get_normal_kappa_plus(dcoll, dd_bdry)
+        kappa_minus = self.get_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,
diff --git a/test/test_diffusion.py b/test/test_diffusion.py
index 244b3ab88..08e4d7e10 100644
--- a/test/test_diffusion.py
+++ b/test/test_diffusion.py
@@ -731,7 +731,6 @@ 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, kappa1])
     u = 30.0*nodes[0] + 60.0*nodes[1]
 
diff --git a/test/test_multiphysics.py b/test/test_multiphysics.py
index b3e03b1c6..8647add97 100644
--- a/test/test_multiphysics.py
+++ b/test/test_multiphysics.py
@@ -150,10 +150,8 @@ def get_rhs(t, u):
 
 @pytest.mark.parametrize("order", [2, 3])
 @pytest.mark.parametrize("use_overintegration", [False, True])
-@pytest.mark.parametrize("orthotropic_kappa", [False, True])
 def test_thermally_coupled_fluid_wall(
-        actx_factory, order, use_overintegration, orthotropic_kappa,
-        visualize=False):
+        actx_factory, order, use_overintegration, visualize=False):
     """Check the thermally-coupled fluid/wall interface."""
     actx = actx_factory()
 
@@ -220,13 +218,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 = {
@@ -244,20 +236,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
@@ -552,10 +544,97 @@ 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("order", [1, 3])
+@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, order, use_overintegration, use_radiation, use_noslip,
+        visualize=False):
+    """Check the RHS shape for orthotropic kappa cases."""
+    actx = actx_factory()
+
+    dim = 2
+    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)
+
+
 if __name__ == "__main__":
     import sys
     if len(sys.argv) > 1:

From d89c33ff335f3f505b7b5e853c522046d42fd40c Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 6 Mar 2024 10:46:28 -0600
Subject: [PATCH 10/15] Assert wall_energy_rhs too.

---
 test/test_multiphysics.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/test_multiphysics.py b/test/test_multiphysics.py
index 8647add97..a00dd0eab 100644
--- a/test/test_multiphysics.py
+++ b/test/test_multiphysics.py
@@ -633,6 +633,7 @@ def test_orthotropic_flux(
 
     from meshmode.dof_array import DOFArray
     assert isinstance(fluid_rhs.energy, DOFArray)
+    assert isinstance(wall_energy_rhs, DOFArray)
 
 
 if __name__ == "__main__":

From 1395d90f963f20762d8550bc320f852f6325de2e Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 6 Mar 2024 13:07:15 -0600
Subject: [PATCH 11/15] Remove kappa + zeros

---
 mirgecom/diffusion.py | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py
index ec6abd34b..ab9759084 100644
--- a/mirgecom/diffusion.py
+++ b/mirgecom/diffusion.py
@@ -695,9 +695,6 @@ def grad_operator(
     dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag)
     dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL)
 
-    # Make sure it has an array context
-    kappa = kappa + dcoll.zeros(array_context=actx, dd=dd_vol)
-
     if kappa_tpairs is None:
         kappa_tpairs = interior_trace_pairs(
             dcoll, kappa, volume_dd=dd_vol, comm_tag=(_DiffusionKappaTag, comm_tag))
@@ -832,9 +829,6 @@ def diffusion_operator(
     dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag)
     dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL)
 
-    # Make sure it has an array context
-    kappa = kappa + dcoll.zeros(array_context=actx, dd=dd_vol)
-
     kappa_tpairs = interior_trace_pairs(
         dcoll, kappa, volume_dd=dd_vol, comm_tag=(_DiffusionKappaTag, comm_tag))
 

From 44c571ca924f6d3345900f9e3ed1cab86df086cd Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 6 Mar 2024 15:51:14 -0600
Subject: [PATCH 12/15] Remove kappa + zeros

---
 mirgecom/multiphysics/thermally_coupled_fluid_wall.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
index ddb0184be..fa77c7404 100644
--- a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
+++ b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
@@ -283,10 +283,8 @@ def get_normal_kappa_minus(self, dcoll, dd_bdry, kappa):
         return kappa
 
     def kappa_bc(self, dcoll, dd_bdry, kappa_minus):
-        kappa_plus = make_obj_array([
-            project_from_base(dcoll, dd_bdry, self._kappa_plus[i])
-            for i in range(self._kappa_plus.shape[0])])
-        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)

From 2cd9137d716936b0410aca5b2276acc13011718f Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Thu, 7 Mar 2024 10:18:07 -0600
Subject: [PATCH 13/15] flake8

---
 mirgecom/multiphysics/thermally_coupled_fluid_wall.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
index fa77c7404..ed237098d 100644
--- a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
+++ b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
@@ -82,7 +82,6 @@
 )
 import grudge.op as op
 
-from pytools.obj_array import make_obj_array
 from mirgecom.math import harmonic_mean
 from mirgecom.boundary import (
     MengaldoBoundaryCondition,

From 27c53617533dd08900f46a4f7958e486e41343c0 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Thu, 7 Mar 2024 12:44:05 -0600
Subject: [PATCH 14/15] Matt's review

---
 mirgecom/multiphysics/thermally_coupled_fluid_wall.py |  4 ++--
 test/test_multiphysics.py                             | 10 ++++------
 2 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
index ed237098d..3e76004d3 100644
--- a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
+++ b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
@@ -262,7 +262,7 @@ def __init__(
         self._t_plus = t_plus
         self._grad_t_plus = grad_t_plus
 
-    def get_normal_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
@@ -272,7 +272,7 @@ def get_normal_kappa_plus(self, dcoll, dd_bdry):
             return np.dot(normal, kappa_plus*normal)
         return project_from_base(dcoll, dd_bdry, self._kappa_plus)
 
-    def get_normal_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
diff --git a/test/test_multiphysics.py b/test/test_multiphysics.py
index a00dd0eab..48ea892cb 100644
--- a/test/test_multiphysics.py
+++ b/test/test_multiphysics.py
@@ -447,10 +447,8 @@ def cv_from_temp(temp):
 
 @pytest.mark.parametrize("order", [1, 3])
 @pytest.mark.parametrize("use_overintegration", [False, True])
-@pytest.mark.parametrize("orthotropic_kappa", [False, True])
 def test_thermally_coupled_fluid_wall_with_radiation(
-        actx_factory, order, use_overintegration, orthotropic_kappa,
-        visualize=False):
+        actx_factory, order, use_overintegration, visualize=False):
     """Check the thermally-coupled fluid/wall interface with radiation.
 
     Analytic solution prescribed as initial condition, then the RHS is assessed
@@ -509,7 +507,7 @@ def test_thermally_coupled_fluid_wall_with_radiation(
     # Made-up wall material
     wall_rho = 1.0
     wall_cp = 1000.0
-    wall_kappa = make_obj_array([1.0, 1.0]) if orthotropic_kappa else 1.0
+    wall_kappa = 1.0
     wall_emissivity = 1.0
 
     base_fluid_temp = 2.0
@@ -549,17 +547,17 @@ def test_thermally_coupled_fluid_wall_with_radiation(
     assert actx.to_numpy(op.norm(dcoll, solid_rhs, np.inf)) < 1e-4
 
 
-@pytest.mark.parametrize("order", [1, 3])
 @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, order, use_overintegration, use_radiation, use_noslip,
+        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)

From 058e3b1b8d644b8ef624745870f1ff53a71c60ff Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Thu, 7 Mar 2024 12:52:41 -0600
Subject: [PATCH 15/15] pylint

---
 mirgecom/multiphysics/thermally_coupled_fluid_wall.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
index 3e76004d3..392ce95ad 100644
--- a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
+++ b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py
@@ -291,9 +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 = self.get_normal_kappa_plus(dcoll, dd_bdry)
-        kappa_minus = self.get_normal_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,