New class to treat plasma data files

mb_scenario.py

@@ -5,10 +5,10 @@
 from dolfinx.fem.function import Constant
 from numpy.typing import NDArray
 
+from hisp import make_mb_model
+from hisp.plamsa_data_handling import PlasmaDataHandling
 from hisp.festim_models import make_mb_model
-from hisp.dina import get_particle_flux, heat
 from hisp.scenario import Scenario
-from hisp import CustomProblem
 
 # dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
 
@@ -26,6 +26,16 @@
 
 if __name__ == "__main__":
     my_scenario = Scenario.from_txt_file("scenario_test.txt", old_format=True)
+    plasma_data_handling = PlasmaDataHandling(
+        pulse_type_to_data = {
+            "FP": np.loadtxt("Binned_Flux_Data.dat", skiprows=1),
+            "ICWC": np.loadtxt("ICWC_data.dat", skiprows=1),
+            "GDC": np.loadtxt("GDC_data.dat", skiprows=1),
+        },
+        path_to_ROSP_data="ROSP_data",
+        path_to_RISP_data="RISP_data",
+        path_to_RISP_wall_data="RISP_Wall_data.dat",
+    )
 
     nb_mb = 64
     L = 6e-3  # m
@@ -49,7 +59,7 @@ def T_function(x: NDArray, t: Constant) -> float:
             T_bake = 483.15  # K
             flat_top_value = np.full_like(x[0], T_bake)
         else:
-            heat_flux = heat(pulse_type, nb_mb, t_rel)
+            heat_flux = plasma_data_handling.heat(pulse_type, nb_mb, t_rel)
             T_surface = 1.1e-4 * heat_flux + COOLANT_TEMP
             T_rear = 2.2e-5 * heat_flux + COOLANT_TEMP
             a = (T_rear - T_surface) / L
@@ -78,7 +88,7 @@ def deuterium_ion_flux(t: float) -> float:
         time_start_current_pulse = my_scenario.get_time_till_row(pulse_row)
         relative_time = t - time_start_current_pulse
 
-        ion_flux = get_particle_flux(
+        ion_flux = plasma_data_handling.get_particle_flux(
             pulse_type=pulse_type, nb_mb=nb_mb, t_rel=relative_time, ion=True
         )
         tritium_fraction = PULSE_TYPE_TO_TRITIUM_FRACTION[pulse_type]
@@ -101,7 +111,7 @@ def tritium_ion_flux(t: float) -> float:
         time_start_current_pulse = my_scenario.get_time_till_row(pulse_row)
         relative_time = t - time_start_current_pulse
 
-        ion_flux = get_particle_flux(
+        ion_flux = plasma_data_handling.get_particle_flux(
             pulse_type=pulse_type, nb_mb=nb_mb, t_rel=relative_time, ion=True
         )
 
@@ -124,7 +134,7 @@ def deuterium_atom_flux(t: float) -> float:
         total_time_pulse = my_scenario.get_pulse_duration(pulse_row)
         time_start_current_pulse = my_scenario.get_time_till_row(pulse_row)
         relative_time = t - time_start_current_pulse
-        atom_flux = get_particle_flux(
+        atom_flux = plasma_data_handling.get_particle_flux(
             pulse_type=pulse_type, nb_mb=nb_mb, t_rel=relative_time, ion=False
         )
 
@@ -147,7 +157,7 @@ def tritium_atom_flux(t: float) -> float:
         time_start_current_pulse = my_scenario.get_time_till_row(pulse_row)
         relative_time = t - time_start_current_pulse
 
-        atom_flux = get_particle_flux(
+        atom_flux = plasma_data_handling.get_particle_flux(
             pulse_type=pulse_type, nb_mb=nb_mb, t_rel=relative_time, ion=False
         )
         tritium_fraction = PULSE_TYPE_TO_TRITIUM_FRACTION[pulse_type]

src/hisp/plamsa_data_handling/__init__.py

@@ -0,0 +1 @@
+from .main import PlasmaDataHandling

src/hisp/plamsa_data_handling/main.py

@@ -0,0 +1,129 @@
+import numpy as np
+from numpy.typing import NDArray
+
+from typing import Dict
+
+
+class PlasmaDataHandling:
+    def __init__(
+        self,
+        pulse_type_to_data: Dict[str, NDArray],
+        path_to_RISP_data: str,
+        path_to_ROSP_data: str,
+        path_to_RISP_wall_data: str,
+    ):
+        self.pulse_type_to_data = pulse_type_to_data
+        self.path_to_RISP_data = path_to_RISP_data
+        self.path_to_ROSP_data = path_to_ROSP_data
+        self.path_to_RISP_wall_data = path_to_RISP_wall_data
+
+    def get_particle_flux(
+        self, pulse_type: str, nb_mb: int, t_rel: float, ion=True
+    ) -> float:
+        """Returns the particle flux for a given pulse type
+
+        Args:
+            pulse_type: pulse type (eg. FP, ICWC, RISP, GDC, BAKE)
+            nb_mb: monoblock number
+            t_rel: t_rel as an integer (in seconds).
+                t_rel = t - t_pulse_start where t_pulse_start is the start of the pulse in seconds
+            ion (bool, optional): _description_. Defaults to True.
+
+        Returns:
+            float: particle flux in part/m2/s
+        """
+        if ion:
+            FP_index = 2
+            other_index = 0
+        if not ion:
+            FP_index = 3
+            other_index = 1
+
+        if pulse_type == "FP":
+            flux = self.pulse_type_to_data[pulse_type][:, FP_index][nb_mb - 1]
+        elif pulse_type == "RISP":
+            assert isinstance(
+                t_rel, float
+            ), f"t_rel should be a float, not {type(t_rel)}"
+            flux = self.RISP_data(nb_mb=nb_mb, t_rel=t_rel)[other_index]
+        elif pulse_type == "BAKE":
+            flux = 0.0
+        else:
+            flux = self.pulse_type_to_data[pulse_type][:, other_index][nb_mb - 1]
+
+        return flux
+
+    def RISP_data(self, nb_mb: int, t_rel: float | int) -> NDArray:
+        """Returns the correct RISP data file for indicated monoblock
+
+        Args:
+            nb_mb: mb number
+            t_rel: t_rel as an integer(in seconds).
+                t_rel = t - t_pulse_start where t_pulse_start is the start of the pulse in seconds
+
+        Returns:
+            data: data from correct file as a numpy array
+        """
+        inner_swept_bins = list(range(46, 65))
+        outer_swept_bins = list(range(19, 34))
+
+        if nb_mb in inner_swept_bins:
+            folder = self.path_to_RISP_data
+            div = True
+            offset_mb = 46
+        elif nb_mb in outer_swept_bins:
+            folder = self.path_to_ROSP_data
+            div = True
+            offset_mb = 19
+        else:
+            div = False
+            offset_mb = 0
+
+        t_rel = int(t_rel)
+        # NOTE: what is the point of this test since it takes nb_mb as an argument?
+        if div:
+            if 1 <= t_rel <= 9:
+                data = np.loadtxt(f"{folder}/time0.dat", skiprows=1)
+            elif 10 <= t_rel <= 98:
+                data = np.loadtxt(f"{folder}/time10.dat", skiprows=1)
+            elif 100 <= t_rel <= 260:
+                data = np.loadtxt(f"{folder}/time{t_rel}.dat", skiprows=1)
+            elif 261 <= t_rel <= 269:
+                data = np.loadtxt(f"{folder}/time270.dat", skiprows=1)
+            else:  # NOTE: so if time is too large a MB transforms into a FW element???
+                data = np.loadtxt(self.path_to_RISP_wall_data, skiprows=1)
+        else:
+            data = np.loadtxt(self.path_to_RISP_wall_data, skiprows=1)
+
+        return data[nb_mb - offset_mb, :]
+
+    def heat(self, pulse_type: str, nb_mb: int, t_rel: float) -> float:
+        """Returns the surface heat flux (W/m2) for a given pulse type
+
+        Args:
+            pulse_type: pulse type (eg. FP, ICWC, RISP, GDC, BAKE)
+            nb_mb: monoblock number
+            t_rel: t_rel as an integer (in seconds).
+                t_rel = t - t_pulse_start where t_pulse_start is the start of the pulse in seconds
+
+        Raises:
+            ValueError: if the pulse type is unknown
+
+        Returns:
+            the surface heat flux in W/m2
+        """
+        if pulse_type == "RISP":
+            data = self.RISP_data(nb_mb, t_rel=t_rel)
+        elif pulse_type in self.pulse_type_to_data.keys():
+            data = self.pulse_type_to_data[pulse_type]
+        else:
+            raise ValueError(f"Invalid pulse type {pulse_type}")
+
+        if pulse_type == "FP":
+            heat_val = data[:, -2][nb_mb - 1]
+        elif pulse_type == "RISP":
+            heat_val = data[-1]
+        else:
+            heat_val = data[:, -1][nb_mb - 1]
+
+        return heat_val