Remove bremsstrahlung and line radiation calculations

process/impurity_radiation.py

@@ -326,29 +326,6 @@ def z2index(zimp):
     error_handling.report_error(33)
 
 
-def pbremden(imp_element_index, nprofile, tprofile):
-    """Calculates the Bremsstrahlung power per volume for a given set of profiles.
-
-    :param imp_element_index: Impurity profile index
-    :type imp_element_index: int
-    :param nprofile: density profile
-    :type nprofile: numpy.array
-    :param tprofile: temperature profile
-    :type tprofile: numpy.array
-    :return: pbremden - Bremsstrahlung power array
-    :rtype: numpy.array
-    """
-    pbremden = (
-        impurity_radiation_module.impurity_arr_frac[imp_element_index]
-        * nprofile
-        * nprofile
-        * numpy.square(zav_of_te(imp_element_index, tprofile))
-        * 5.355e-37
-        * numpy.sqrt(tprofile)
-    )
-    return pbremden
-
-
 def fradcore(rho, coreradius, coreradiationfraction):
     """Finds the fraction of radiation from the core that is subtracted in impurity radiation model.
 
@@ -492,10 +469,9 @@ def pimpden(imp_element_index, nprofile, tprofile):
 
 class ImpurityRadiation:
     """This class calculates the impurity radiation losses for given temperature and density profiles.
-    The considers the Bremsstrahlung (radb), line radiation (radl), total impurity radiation
-    from the core (radcore) and total impurity radiation (radtot) [MW/(m^3)].
-    The class is used to sum the impurity radiation loss from each impurity element to find the
-    total impurity radiation loss."""
+    The considers the  total impurity radiation from the core (radcore) and total impurity radiation
+    (radtot) [MW/(m^3)]. The class is used to sum the impurity radiation loss from each impurity
+    element to find the total impurity radiation loss."""
 
     def __init__(self, plasma_profile):
         """
@@ -509,18 +485,12 @@ def __init__(self, plasma_profile):
             0
         ]
 
-        self.pline_profile = numpy.zeros(self.plasma_profile.profile_size)
-        self.pbrem_profile = numpy.zeros(self.plasma_profile.profile_size)
         self.pimp_profile = numpy.zeros(self.plasma_profile.profile_size)
         self.radtot_profile = numpy.zeros(self.plasma_profile.profile_size)
         self.radcore_profile = numpy.zeros(self.plasma_profile.profile_size)
-        self.radb_profile = numpy.zeros(self.plasma_profile.profile_size)
-        self.radl_profile = numpy.zeros(self.plasma_profile.profile_size)
 
         self.radtot = 0.0
         self.radcore = 0.0
-        self.radb = 0.0
-        self.radl = 0.0
 
     def map_imprad_profile(self):
         """Map imprad_profile() over each impurity element index."""
@@ -545,33 +515,12 @@ def imprad_profile(self, imp_element_index):
         :type imp_element_index: Int
         """
 
-        pbrem = pbremden(
-            imp_element_index,
-            self.plasma_profile.neprofile.profile_y,
-            self.plasma_profile.teprofile.profile_y,
-        )
-        self.pbrem_profile = numpy.add(self.pbrem_profile, pbrem)
-
         pimp = pimpden(
             imp_element_index,
             self.plasma_profile.neprofile.profile_y,
             self.plasma_profile.teprofile.profile_y,
         )
 
-        pline = numpy.zeros(self.plasma_profile.profile_size)
-
-        # Using a 'mask' to apply logic to the array values in lieu of loops
-        pimp_greater_pbrem_mask = pimp >= pbrem
-        pline[pimp_greater_pbrem_mask] = (
-            pimp[pimp_greater_pbrem_mask] - pbrem[pimp_greater_pbrem_mask]
-        )
-        self.pline_profile = numpy.add(self.pline_profile, pline)
-
-        # There is a model inconsistency where pimp < pbrem, which should never be true.
-        # This fix is okay at high T but should be remedied.
-        pimp[numpy.invert(pimp_greater_pbrem_mask)] = pbrem[
-            numpy.invert(pimp_greater_pbrem_mask)
-        ]
         self.pimp_profile = numpy.add(self.pimp_profile, pimp)
 
     def calculate_radiation_loss_profiles(self):
@@ -590,23 +539,12 @@ def calculate_radiation_loss_profiles(self):
             )
         )
 
-        radb = self.pbrem_profile * self.rho
-        radl = self.pline_profile * self.rho
-
         self.radtot_profile = numpy.add(self.radtot_profile, radtot)
         self.radcore_profile = numpy.add(self.radcore_profile, radcore)
-        self.radb_profile = numpy.add(self.radb_profile, radb)
-        self.radl_profile = numpy.add(self.radl_profile, radl)
 
     def integrate_radiation_loss_profiles(self):
         """Integrate the radiation loss profiles using the Simpson rule.
         Store the total values for each aspect of impurity radiation loss."""
-        self.radb = 2.0e-6 * integrate.simpson(
-            self.radb_profile, x=self.rho, dx=self.rhodx
-        )
-        self.radl = 2.0e-6 * integrate.simpson(
-            self.radl_profile, x=self.rho, dx=self.rhodx
-        )
         self.radtot = 2.0e-6 * integrate.simpson(
             self.radtot_profile, x=self.rho, dx=self.rhodx
         )

process/io/plot_proc.py

@@ -1042,12 +1042,9 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float:
                         # Zav[i][k] = imp_data[i][j][2]
             # The impurity radiation
             pimpden[i][k] = imp_frac[i] * ne[k] * ne[k] * lz[i][k]
-            # The Bremsstrahlung
-            # pbremden[i][k] = imp_frac[i] * ne[k] * ne[k] * Zav[i][k] * Zav[i][k] * 5.355e-37 * np.sqrt(te[k])
 
         for l in range(imp_data.shape[0]):  # noqa: E741
             prad[k] = prad[k] + pimpden[l][k] * 2.0e-6
-            # pbrem[k] = pbrem[k] + pbremden[l][k] * 2.0e-6
 
     # benchmark prad again outfile so mod prad
     # pbremint = (rho[1:] * pbrem[1:]) @ drho
@@ -1058,7 +1055,6 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float:
     # print('pbrem = ',pbrem)
     # print(1.0e32*lz[12])
     # print('pradpv = ',pradint)
-    # print('pbrempv = ',pbremint)
     # print('pbremmw = ',pbremint*vol)
     # print('pradmw = ', pradint*vol, 'MW') # pimp = pline + pbrem
 

process/physics.py

@@ -1177,8 +1177,6 @@ def physics(self):
         # Calculate radiation power
 
         radpwrdata = physics_funcs.radpwr(self.plasma_profile)
-        physics_variables.pbrempv = radpwrdata.pbrempv
-        physics_variables.plinepv = radpwrdata.plinepv
         physics_variables.psyncpv = radpwrdata.psyncpv
         physics_variables.pcoreradpv = radpwrdata.pcoreradpv
         physics_variables.pedgeradpv = radpwrdata.pedgeradpv
@@ -3248,20 +3246,6 @@ def outplas(self):
         # po.ovarre(self.outfile,'Total power deposited in plasma (MW)','()',falpha*palpmw+pchargemw+pohmmw+pinjmw, 'OP ')
 
         po.osubhd(self.outfile, "Radiation Power (excluding SOL):")
-        po.ovarre(
-            self.outfile,
-            "Bremsstrahlung radiation power (MW)",
-            "(pbrempv*vol)",
-            physics_variables.pbrempv * physics_variables.vol,
-            "OP ",
-        )
-        po.ovarre(
-            self.outfile,
-            "Line radiation power (MW)",
-            "(plinepv*vol)",
-            physics_variables.plinepv * physics_variables.vol,
-            "OP ",
-        )
         po.ovarre(
             self.outfile,
             "Synchrotron radiation power (MW)",

process/physics_functions.py

@@ -255,9 +255,7 @@ def radpwr(plasma_profile):
     #  Total radiation power/volume.
     pradpv = imp_rad.radtot + psyncpv
 
-    return RadpwrData(
-        imp_rad.radb, imp_rad.radl, psyncpv, pcoreradpv, pedgeradpv, pradpv
-    )
+    return RadpwrData(psyncpv, pcoreradpv, pedgeradpv, pradpv)
 
 
 def psync_albajar_fidone():
@@ -461,8 +459,6 @@ class BoschHaleConstants:
 class RadpwrData:
     """DataClass which holds the output of the function radpwr"""
 
-    pbrempv: float
-    plinepv: float
     psyncpv: float
     pcoreradpv: float
     pedgeradpv: float

process/stellarator.py

@@ -4074,8 +4074,6 @@ def stphys(self, output):
 
         #  Calculate radiation power
         radpwr_data = physics_funcs.radpwr(self.plasma_profile)
-        physics_variables.pbrempv = radpwr_data.pbrempv
-        physics_variables.plinepv = radpwr_data.plinepv
         physics_variables.psyncpv = radpwr_data.psyncpv
         physics_variables.pcoreradpv = radpwr_data.pcoreradpv
         physics_variables.pedgeradpv = radpwr_data.pedgeradpv

source/fortran/physics_variables.f90

@@ -577,9 +577,6 @@ module physics_variables
   real(dp) :: palpnb
   !! alpha power from hot neutral beam ions (MW)
 
-  real(dp) :: pbrempv
-  !! bremsstrahlung power per volume (MW/m3)
-
   real(dp) :: pchargemw
   !! non-alpha charged particle fusion power (MW)
 
@@ -637,9 +634,6 @@ module physics_variables
   real(dp) :: plascur
   !! plasma current (A)
 
-  real(dp) :: plinepv
-  !! line radiation power per volume (MW/m3)
-
   real(dp) :: pneutmw
   !! neutron fusion power (MW)
 
@@ -1009,7 +1003,6 @@ subroutine init_physics_variables
     palpipv = 0.0D0
     palpmw = 0.0D0
     palpnb = 0.0D0
-    pbrempv = 0.0D0
     pchargemw = 0.0D0
     pchargepv = 0.0D0
     pcoef = 0.0D0
@@ -1029,7 +1022,6 @@ subroutine init_physics_variables
     photon_wall = 0.0D0
     piepv = 0.0D0
     plascur = 0.0D0
-    plinepv = 0.0D0
     pneutmw = 0.0D0
     pneutpv = 0.0D0
     pohmmw = 0.0D0

tests/unit/test_impurity_radiation.py

@@ -12,85 +12,6 @@ def initialise_impurity_radiation():
     impurity_radiation.initialise_imprad()
 
 
-class PbremdenParam(NamedTuple):
-    imp_element_index: int = 0
-    ne: np.array = np.array
-    te: np.array = np.array
-    expected_pbremden: np.array = np.array
-
-
-def test_pbremden():
-    """Tests `pbremden` function.
-
-    :param imp_element_index: impurity element
-    :type imp_element_index: float
-
-    :param ne:  electron density (/m3).
-    :type ne: np.array
-
-    :param te: electron temperature (keV)
-    :type te: np.array
-
-    :param expected_pbremden: Bremsstrahlung radiation density (W/m3)
-    :type expected_pbremden: np.array
-    """
-    pbrem_parameters = PbremdenParam(
-        imp_element_index=0,
-        ne=np.array(
-            [
-                9.42593370e19,
-                9.37237672e19,
-                9.21170577e19,
-                8.94392086e19,
-                8.56902197e19,
-                8.08700913e19,
-                7.49788231e19,
-                6.80164153e19,
-                5.99828678e19,
-                3.28986749e19,
-            ]
-        ),
-        te=np.array(
-            [
-                27.73451868,
-                27.25167194,
-                25.82164396,
-                23.50149071,
-                20.39190536,
-                16.64794796,
-                12.50116941,
-                8.31182764,
-                4.74643357,
-                0.1,
-            ]
-        ),
-        expected_pbremden=np.array(
-            [
-                25056.39306004834,
-                24555.88124974145,
-                23090.40418076286,
-                20766.488978621197,
-                17756.23528821658,
-                14289.456043405115,
-                10644.175682502506,
-                7142.2528057621585,
-                4197.5669332439475,
-                183.28051080066396,
-            ]
-        ),
-    )
-    assert (
-        pytest.approx(
-            impurity_radiation.pbremden(
-                pbrem_parameters.imp_element_index,
-                pbrem_parameters.ne,
-                pbrem_parameters.te,
-            )
-        )
-        == pbrem_parameters.expected_pbremden
-    )
-
-
 class PimpdenParam(NamedTuple):
     imp_element_index: int = 0
     ne: np.array = np.array