Fixing heat flux corrective factor (Fusion-Power-Plant-Framework#3558)

* Fixing DN heat flux corrective factor

* Added no intersection region to SN

* Added notes for the perpendicular transport at the mid-plane in analyse().
Fixed test

bluemira/radiation_transport/advective_transport.py

@@ -19,8 +19,9 @@
 from bluemira.base.constants import EPS
 from bluemira.base.look_and_feel import bluemira_warn
 from bluemira.display.plotter import Zorder, plot_coordinates
-from bluemira.geometry.coordinates import Coordinates
+from bluemira.geometry.coordinates import Coordinates, coords_plane_intersect
 from bluemira.geometry.plane import BluemiraPlane
+from bluemira.geometry.tools import make_polygon
 from bluemira.radiation_transport.error import AdvectionTransportError
 from bluemira.radiation_transport.flux_surfaces_maker import _clip_flux_surfaces
 
@@ -126,8 +127,7 @@ def _check_params(self):
                 " in either the lower or upper directions."
             )
 
-    @staticmethod
-    def _process_first_wall(first_wall):
+    def _process_first_wall(self, first_wall):
         """
         Force working first wall geometry to be closed and counter-clockwise.
         """
@@ -142,7 +142,11 @@ def _process_first_wall(first_wall):
         if not first_wall.closed:
             bluemira_warn("First wall should be a closed geometry. Closing it.")
             first_wall.close()
-        return first_wall
+
+        int_intersection = coords_plane_intersect(first_wall, self._yz_plane)[0]
+        out_intersection = coords_plane_intersect(first_wall, self._yz_plane)[1]
+
+        return first_wall, int_intersection, out_intersection
 
     @staticmethod
     def _get_arrays(flux_surfaces):
@@ -217,6 +221,35 @@ def _clip_flux_surfaces(self, first_wall):
             ],
         )
 
+    def _no_wall_intersection_region(
+        self, x_up_inter, z_up_inter, x_down_inter, z_down_inter, *, lfs=True
+    ):
+        """
+        Get first wall mid-plane region with no flux line inetrsections.
+        """
+        up_end_i = self.first_wall.argmin(np.array([x_up_inter[-1], 0, z_up_inter[-1]]))
+        down_end_i = self.first_wall.argmin(
+            np.array([x_down_inter[-1], 0, z_down_inter[-1]])
+        )
+
+        reg_i = np.nonzero(
+            (self.first_wall.z < self.first_wall.z[up_end_i])
+            & (self.first_wall.z >= self.first_wall.z[down_end_i])
+            & (
+                (self.first_wall.x > self._o_point.x)
+                if lfs
+                else (self.first_wall.x < self._o_point.x)
+            )
+        )[0]
+
+        x_reg_inter = self.first_wall.x[reg_i]
+        z_reg_inter = self.first_wall.z[reg_i]
+
+        reg_wire = make_polygon(self.first_wall.T[reg_i])
+        wire_length = reg_wire.length
+
+        return x_reg_inter, z_reg_inter, wire_length
+
     def analyse(self, first_wall: Coordinates):
         """
         Perform the calculation to obtain charged particle heat fluxes on the
@@ -235,8 +268,23 @@ def analyse(self, first_wall: Coordinates):
             The z coordinates of the flux surface intersections
         heat_flux: np.array
             The perpendicular heat fluxes at the intersection points [MW/m^2]
+
+        Notes
+        -----
+        The heat flux model assumes pure parallel transport and fudges
+        the perpendicular transport via the power decay length, lambda.
+        This approach, while is widely used, leads to no power deposited
+        on the wall at the mid-plane.
+        _analyse_SN and _analyse_DN assume, in the area in proximity of
+        the mid-plane, where the outermost flux tube is open, that the
+        remaining power of from the exponential decay is deposited on
+        the wall perpendicularly.
         """
-        self.first_wall = self._process_first_wall(first_wall)
+        (
+            self.first_wall,
+            self.imp_int,
+            self.omp_int,
+        ) = self._process_first_wall(first_wall)
 
         if self.eq.is_double_null:
             x, z, hf = self._analyse_DN()
@@ -249,6 +297,7 @@ def analyse(self, first_wall: Coordinates):
     def _analyse_SN(self):
         """
         Calculation for the case of single nulls.
+
         """
         self._make_flux_surfaces_ob()
 
@@ -284,16 +333,42 @@ def _analyse_SN(self):
         heat_flux_lfs = self.params.f_lfs_lower_target * q_par_lfs * np.sin(alpha_lfs)
         heat_flux_hfs = self.params.f_hfs_lower_target * q_par_hfs * np.sin(alpha_hfs)
 
-        # Correct power (energy conservation)
-        q_omp_int = 2 * np.pi * np.sum(q_par_omp / (B_omp / Bp_omp) * self.dx_mp * x_omp)
-        f_correct_power = self.params.P_sep_particle / q_omp_int
+        # Find FW portion for perpendicular power
+        if self.first_wall.argmin(
+            np.array([x_hfs_inter[-1], 0, z_hfs_inter[-1]])
+        ) != self.first_wall.argmin(np.array([x_lfs_inter[-1], 0, z_lfs_inter[-1]])):
+            x_out_inter, z_out_inter, outb_length = self._no_wall_intersection_region(
+                x_hfs_inter, z_hfs_inter, x_lfs_inter, z_lfs_inter, lfs=True
+            )
+        else:
+            mid_i = self.first_wall.argmin(
+                np.array([x_lfs_inter[-1], 0, z_lfs_inter[-1]])
+            )
+            x_out_inter = self.first_wall.x[mid_i]
+            z_out_inter = self.first_wall.z[mid_i]
+            outb_length = 1
+
+        # Calculating missing power from parallel transport
+        q_omp_int = 2 * np.pi * np.sum(q_par_omp * Bp_omp / B_omp * self.dx_mp * x_omp)
+        miss_omp = self.params.P_sep_particle - q_omp_int
+        outb_surf = outb_length * 2 * np.pi * self.omp_int[0]
+
+        # Calculating mid-outboard and mid-inboard heat flux
+        heat_flux_x_outb = miss_omp / outb_surf
+        if outb_length != 1:
+            heat_flux_x_outb = [heat_flux_x_outb] * len(x_out_inter)
+
         return (
-            np.append(x_lfs_inter, x_hfs_inter),
-            np.append(z_lfs_inter, z_hfs_inter),
-            f_correct_power * np.append(heat_flux_lfs, heat_flux_hfs),
+            np.concatenate([np.atleast_1d(x_out_inter), x_lfs_inter, x_hfs_inter]),
+            np.concatenate([np.atleast_1d(z_out_inter), z_lfs_inter, z_hfs_inter]),
+            np.concatenate([
+                np.atleast_1d(heat_flux_x_outb),
+                heat_flux_lfs,
+                heat_flux_hfs,
+            ]),
         )
 
-    def _analyse_DN(self):
+    def _analyse_DN(self):  # noqa: PLR0914
         """
         Calculation for the case of double nulls.
         """
@@ -336,7 +411,8 @@ def _analyse_DN(self):
         Bt_imp = self.eq.Bt(x_imp)
         B_imp = np.hypot(Bp_imp, Bt_imp)
 
-        # Parallel power at the outboard and inboard midplane
+        # Parallel power set-up at the outboard and inboard midplane
+        # Note that the power is not split yet into hfs and lfs rates
         q_par_omp = self._q_par(x_omp, dx_omp, B_omp, Bp_omp)
         q_par_imp = self._q_par(x_imp, dx_imp, B_imp, Bp_imp, outboard=False)
 
@@ -347,13 +423,15 @@ def _analyse_DN(self):
         Bp_hfs_up = self.eq.Bp(x_hfs_up_inter, z_hfs_up_inter)
 
         # Calculate parallel power at the intersections
+        # Each q_par_* stores full P_sep_particle
         # Note that flux expansion terms cancel down to this
         q_par_lfs_down = q_par_omp * Bp_lfs_down / B_omp
         q_par_lfs_up = q_par_omp * Bp_lfs_up / B_omp
         q_par_hfs_down = q_par_imp * Bp_hfs_down / B_imp
         q_par_hfs_up = q_par_imp * Bp_hfs_up / B_imp
 
         # Calculate perpendicular heat fluxes
+        # Here P_sep_particle actually gets distributed over the four targets
         heat_flux_lfs_down = (
             self.params.f_lfs_lower_target * q_par_lfs_down * np.sin(alpha_lfs_down)
         )
@@ -367,44 +445,65 @@ def _analyse_DN(self):
             self.params.f_hfs_upper_target * q_par_hfs_up * np.sin(alpha_hfs_up)
         )
 
-        # Correct power (energy conservation)
-        q_omp_int = 2 * np.pi * np.sum(q_par_omp * Bp_omp / B_omp * self.dx_mp * x_omp)
-        q_imp_int = 2 * np.pi * np.sum(q_par_imp * Bp_imp / B_imp * self.dx_mp * x_imp)
-
-        total_power = self.params.P_sep_particle
-        f_outboard = self.params.f_lfs_lower_target + self.params.f_lfs_upper_target
-        f_inboard = self.params.f_hfs_lower_target + self.params.f_hfs_upper_target
-        f_correct_lfs_down = (
-            total_power * self.params.f_lfs_lower_target / f_outboard
-        ) / q_omp_int
-        f_correct_lfs_up = (
-            total_power * self.params.f_lfs_upper_target / f_outboard
-        ) / q_omp_int
-        f_correct_hfs_down = (
-            total_power * self.params.f_hfs_lower_target / f_inboard
-        ) / q_imp_int
-        f_correct_hfs_up = (
-            total_power * self.params.f_hfs_upper_target / f_inboard
-        ) / q_imp_int
+        # Find FW portion for perpendicular power
+        x_out_inter, z_out_inter, outb_length = self._no_wall_intersection_region(
+            x_lfs_up_inter, z_lfs_up_inter, x_lfs_down_inter, z_lfs_down_inter, lfs=True
+        )
+        x_in_inter, z_in_inter, inb_length = self._no_wall_intersection_region(
+            x_hfs_up_inter, z_hfs_up_inter, x_hfs_down_inter, z_hfs_down_inter, lfs=False
+        )
+        # Calculating missing power from parallel transport
+        q_omp_int = (
+            2
+            * np.pi
+            * np.sum(q_par_omp * Bp_omp / B_omp * self.dx_mp * x_omp)
+            * (self.params.f_lfs_lower_target + self.params.f_lfs_upper_target)
+        )
+        q_imp_int = (
+            2
+            * np.pi
+            * np.sum(q_par_imp * Bp_imp / B_imp * self.dx_mp * x_imp)
+            * (self.params.f_hfs_lower_target + self.params.f_hfs_upper_target)
+        )
+        miss_omp = (
+            self.params.P_sep_particle
+            * (self.params.f_lfs_lower_target + self.params.f_lfs_upper_target)
+        ) - q_omp_int
+        miss_imp = (
+            self.params.P_sep_particle
+            * (self.params.f_hfs_lower_target + self.params.f_hfs_upper_target)
+        ) - q_imp_int
+        outb_surf = outb_length * 2 * np.pi * self.omp_int[0]
+        inb_surf = inb_length * 2 * np.pi * self.imp_int[0]
+
+        # Calculating mid-outboard and mid-inboard heat flux
+        heat_flux_x_outb = [miss_omp / outb_surf] * len(x_out_inter)
+        heat_flux_x_inb = [miss_imp / inb_surf] * len(x_in_inter)
 
         return (
             np.concatenate([
+                x_out_inter,
                 x_lfs_down_inter,
                 x_lfs_up_inter,
+                x_in_inter,
                 x_hfs_down_inter,
                 x_hfs_up_inter,
             ]),
             np.concatenate([
+                z_out_inter,
                 z_lfs_down_inter,
                 z_lfs_up_inter,
+                z_in_inter,
                 z_hfs_down_inter,
                 z_hfs_up_inter,
             ]),
             np.concatenate([
-                f_correct_lfs_down * heat_flux_lfs_down,
-                f_correct_lfs_up * heat_flux_lfs_up,
-                f_correct_hfs_down * heat_flux_hfs_down,
-                f_correct_hfs_up * heat_flux_hfs_up,
+                heat_flux_x_outb,
+                heat_flux_lfs_down,
+                heat_flux_lfs_up,
+                heat_flux_x_inb,
+                heat_flux_hfs_down,
+                heat_flux_hfs_up,
             ]),
         )
 

tests/radiation_transport/test_advective_transport.py

@@ -108,20 +108,22 @@ def test_single_null_warnings(self, caplog):
 
     def test_recursion(self):
         assert np.isclose(np.max(self.hf), 5.379, rtol=1e-2)
-        assert np.argmax(self.hf) == 0
+        assert np.argmax(self.hf) == 1
         assert np.isclose(np.sum(self.hf), 399, rtol=1e-2)
 
     def test_n_intersections(self):
         """
         Because it is a single null, we expect the same number of flux surfaces LFS and
         HFS.
+        The final number of intersection points may still be higher than the number of
+        flux surfaces as there are additional points in the mid-plane region
         """
         len_ob_lfs = len(self.solver.flux_surfaces_ob_down)
         len_ob_hfs = len(self.solver.flux_surfaces_ob_up)
 
         assert len_ob_hfs == len_ob_lfs
-        assert len_ob_lfs + len_ob_hfs == len(self.x)
-        assert len(self.solver.flux_surfaces) == len(self.x)
+        assert len_ob_lfs + len_ob_hfs <= len(self.x)
+        assert len(self.solver.flux_surfaces) <= len(self.x)
 
     def test_integrals(self):
         n_fs = len(self.solver.flux_surfaces)
@@ -235,8 +237,8 @@ def test_double_null_warnings(self, caplog):
         assert caplog.records[0].levelname == "WARNING"
 
     def test_recursion(self):
-        assert np.isclose(np.max(self.hf), 86.194, rtol=1e-2)
-        assert np.isclose(np.sum(self.hf), 830.6, rtol=1e-2)
+        assert np.isclose(np.max(self.hf), 130.992, rtol=1e-2)
+        assert np.isclose(np.sum(self.hf), 1254.656, rtol=1e-2)
 
     def test_analyse_DN(self, caplog):
         fw = deepcopy(self.solver.first_wall)