Add radiation tools tests (#3394)

* ✅ Add few tests

* ✅ Add test for make_wall_detectors()

* ✅ Add test for FirstWallRadiationSolver

* 🎨 Allow more customisation in FirstwallSolver

* test_FirstWallRadiationSolver: test the total power

* increase tolerance to test the total power

---------

Co-authored-by: Athoy Nilima <athoy.nilima@ukaea.uk>
Co-authored-by: james <james.cook1@ukaea.uk>

bluemira/radiation_transport/radiation_tools.py

@@ -1095,7 +1095,15 @@ def __init__(self, source_func: Callable, firstwall_shape: BluemiraWire):
         self.rad_source = source_func
         self.fw_shape = firstwall_shape
 
-    def solve(self, *, plot=True):
+    def solve(
+        self,
+        max_wall_len: float = 1e-3,
+        x_width: float = 0.01,
+        n_samples: int = 500,
+        *,
+        plot: bool = True,
+        verbose: bool = False,
+    ):
         """Solve first wall radiation problem"""
         shift = translate(0, 0, np.min(self.fw_shape.z))
         height = np.max(self.fw_shape.z) - np.min(self.fw_shape.z)
@@ -1113,12 +1121,10 @@ def solve(self, *, plot=True):
             parent=world,
             material=emitter,
         )
-        max_wall_len = 10.0e-2
-        X_WIDTH = 0.01
         wall_detectors = make_wall_detectors(
-            self.fw_shape.x, self.fw_shape.z, max_wall_len, X_WIDTH
+            self.fw_shape.x, self.fw_shape.z, max_wall_len, x_width, debug=verbose
         )
-        wall_loads = detect_radiation(wall_detectors, 500, world)
+        wall_loads = detect_radiation(wall_detectors, n_samples, world, verbose=verbose)
 
         if plot:
             plot_radiation_loads(

tests/radiation_transport/test_radiation_profile.py

@@ -10,6 +10,7 @@
 import pytest
 
 from bluemira.base import constants
+from bluemira.base.constants import raw_uc
 from bluemira.base.file import get_bluemira_path
 from bluemira.codes.process import api
 from bluemira.equilibria.equilibrium import Equilibrium
@@ -20,8 +21,15 @@
 )
 from bluemira.radiation_transport.radiation_profile import RadiationSource
 from bluemira.radiation_transport.radiation_tools import (
+    DetectedRadiation,
+    FirstWallRadiationSolver,
     electron_density_and_temperature_sol_decay,
+    grid_interpolator,
+    interpolated_field_values,
     ion_front_distance,
+    linear_interpolator,
+    make_wall_detectors,
+    pfr_filter,
     radiative_loss_function_values,
     target_temperature,
     upstream_temperature,
@@ -92,6 +100,7 @@ def setup_class(cls):
             sol_impurities=cls.config["f_imp_sol"],
         )
         source.analyse(firstwall_geom=fw_shape)
+        source.rad_map(fw_shape)
 
         cls.profiles = profiles
         cls.source = source
@@ -296,3 +305,53 @@ def test_radiative_loss_function_values(self):
         lvals = np.array([imp_data_l_ref[1], imp_data_l_ref[3]])
         l1 = radiative_loss_function_values(tvals, t_ref, l_ref)
         np.testing.assert_allclose(l1, lvals, rtol=2e-1)
+
+    def test_pfr_filter(self):
+        x_point_z = self.source.sol_rad.points["x_point"]["z_low"]
+        pfr_x_down, pfr_z_down = pfr_filter(self.source.sol_rad.separatrix, x_point_z)
+        assert pfr_x_down.shape == (59,)
+        assert pfr_z_down.shape == (59,)
+
+        assert np.all(pfr_z_down < x_point_z - 0.01)
+
+    def test_make_wall_detectors(self):
+        max_wall_len = 10.0e-2
+        X_WIDTH = 0.01
+        wall_detectors = make_wall_detectors(
+            self.fw_shape.x, self.fw_shape.z, max_wall_len, X_WIDTH
+        )
+        assert all(detector[2] <= max_wall_len for detector in wall_detectors)
+        assert all(detector[1] == X_WIDTH for detector in wall_detectors)
+        assert len(wall_detectors) == 532
+
+    def test_FirstWallRadiationSolver(self):
+        # Coversion required for CHERAB
+        f_sol = linear_interpolator(
+            self.source.sol_rad.x_tot,
+            self.source.sol_rad.z_tot,
+            raw_uc(self.source.sol_rad.rad_tot, "MW", "W"),
+        )
+
+        # SOL radiation grid
+        x_sol = np.linspace(min(self.fw_shape.x), max(self.fw_shape.x), 4)
+        z_sol = np.linspace(min(self.fw_shape.z), max(self.fw_shape.z), 4)
+
+        rad_sol_grid = interpolated_field_values(x_sol, z_sol, f_sol)
+        func = grid_interpolator(x_sol, z_sol, rad_sol_grid)
+        solver = FirstWallRadiationSolver(
+            source_func=func, firstwall_shape=self.fw_shape
+        )
+        assert solver.fw_shape == self.fw_shape
+
+        wall_loads = solver.solve(50, 1, 10, plot=False)
+
+        # check return types
+        assert isinstance(wall_loads, DetectedRadiation)
+        assert isinstance(wall_loads.total_power, float)
+
+        # check if the arrays are of same length
+        assert len(wall_loads.detector_numbers) == len(wall_loads.detected_power)
+
+        # the solver gives slightly different powers in each run
+        # so just asserting the order of total power
+        assert np.isclose(wall_loads.total_power, 3.14e8, rtol=0.05)

tests/radiation_transport/test_radiation_tools.py

@@ -11,7 +11,10 @@
     calculate_line_radiation_loss,
     calculate_z_species,
     exponential_decay,
+    filtering_in_or_out,
     gaussian_decay,
+    interpolated_field_values,
+    linear_interpolator,
 )
 
 
@@ -46,3 +49,34 @@ def test_calculate_line_radiation_loss():
     frac = 0.01
     rad = calculate_line_radiation_loss(ne, p_loss, frac)
     assert rad == pytest.approx(0.796, abs=1e-3)
+
+
+def test_interpolated_field_values():
+    x = np.array([0, 1, 1, 0])
+    z = np.array([0, 0, 1, 1])
+    field = np.array([0, 1, 2, 3])
+
+    interpolator = linear_interpolator(x, z, field)
+    x_new = np.array([0.5, 1.5])
+    z_new = np.array([0.5, 1.5])
+
+    field_grid = interpolated_field_values(x_new, z_new, interpolator)
+
+    assert field_grid.shape == (2, 2)
+    assert np.isclose(field_grid[0, 0], interpolator(0.5, 0.5))
+    assert np.isclose(field_grid[1, 1], interpolator(1.5, 1.5))
+
+
+def test_filtering_in_or_out():
+    test_domain_x = [1.0, 2.0, 3.0]
+    test_domain_z = [4.0, 5.0, 6.0]
+
+    include_func = filtering_in_or_out(test_domain_x, test_domain_z, include_points=True)
+    assert include_func([2.0, 5.0])
+    assert not include_func([0.0, 0.0])
+
+    exclude_func = filtering_in_or_out(
+        test_domain_x, test_domain_z, include_points=False
+    )
+    assert not exclude_func([2.0, 5.0])
+    assert exclude_func([0.0, 0.0])