From 34a8af824c5d3e6e6115f8ea26c5743d236e31ed Mon Sep 17 00:00:00 2001 From: apearce Date: Mon, 3 Jun 2024 14:09:22 +0100 Subject: [PATCH 1/4] =?UTF-8?q?=F0=9F=94=A5=20remove=20bremsstrahlung=20an?= =?UTF-8?q?d=20line=20radiation=20calcuations?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/io/plot_proc.py | 6 +----- process/physics.py | 16 ---------------- process/physics_functions.py | 4 +--- process/stellarator.py | 2 -- source/fortran/physics_variables.f90 | 8 -------- tests/unit/test_impurity_radiation.py | 10 ---------- 6 files changed, 2 insertions(+), 44 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index b56a35cf87..665bc27a7a 100755 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -1042,13 +1042,10 @@ 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 # pradint = prad[1:] @ drho * 2.0e-5 @@ -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 diff --git a/process/physics.py b/process/physics.py index b60174fd27..bfd1862dee 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1183,8 +1183,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 @@ -3257,20 +3255,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)", diff --git a/process/physics_functions.py b/process/physics_functions.py index e0f337d377..c9a885e251 100644 --- a/process/physics_functions.py +++ b/process/physics_functions.py @@ -256,7 +256,7 @@ def radpwr(plasma_profile): pradpv = imp_rad.radtot + psyncpv return RadpwrData( - imp_rad.radb, imp_rad.radl, psyncpv, pcoreradpv, pedgeradpv, pradpv + psyncpv, pcoreradpv, pedgeradpv, pradpv ) @@ -462,8 +462,6 @@ class BoschHaleConstants: class RadpwrData: """DataClass which holds the output of the function radpwr""" - pbrempv: float - plinepv: float psyncpv: float pcoreradpv: float pedgeradpv: float diff --git a/process/stellarator.py b/process/stellarator.py index 8979be67ec..45050ca633 100644 --- a/process/stellarator.py +++ b/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 diff --git a/source/fortran/physics_variables.f90 b/source/fortran/physics_variables.f90 index 3d3e0150bc..ebb6789e5d 100644 --- a/source/fortran/physics_variables.f90 +++ b/source/fortran/physics_variables.f90 @@ -578,9 +578,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) @@ -638,9 +635,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) @@ -1010,7 +1004,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 @@ -1030,7 +1023,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 diff --git a/tests/unit/test_impurity_radiation.py b/tests/unit/test_impurity_radiation.py index ad5f6cac5e..cd400c38ec 100644 --- a/tests/unit/test_impurity_radiation.py +++ b/tests/unit/test_impurity_radiation.py @@ -10,16 +10,6 @@ def initialise_impurity_radiation(): impurity_radiation_module.init_impurity_radiation_module() 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 From 6060cafb851a51fa13de53cca59f158e3db071a9 Mon Sep 17 00:00:00 2001 From: apearce Date: Mon, 3 Jun 2024 14:18:04 +0100 Subject: [PATCH 2/4] =?UTF-8?q?=F0=9F=9A=A8=20flake?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/impurity_radiation.py | 46 ----------------------------------- process/io/plot_proc.py | 2 +- process/physics_functions.py | 4 +-- 3 files changed, 2 insertions(+), 50 deletions(-) diff --git a/process/impurity_radiation.py b/process/impurity_radiation.py index 604090c064..8003cceda6 100644 --- a/process/impurity_radiation.py +++ b/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. @@ -509,8 +486,6 @@ 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) @@ -545,33 +520,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): diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 665bc27a7a..4948593c73 100755 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -1045,7 +1045,7 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float: for l in range(imp_data.shape[0]): # noqa: E741 prad[k] = prad[k] + pimpden[l][k] * 2.0e-6 - + # benchmark prad again outfile so mod prad # pbremint = (rho[1:] * pbrem[1:]) @ drho # pradint = prad[1:] @ drho * 2.0e-5 diff --git a/process/physics_functions.py b/process/physics_functions.py index c9a885e251..403ff83acd 100644 --- a/process/physics_functions.py +++ b/process/physics_functions.py @@ -255,9 +255,7 @@ def radpwr(plasma_profile): # Total radiation power/volume. pradpv = imp_rad.radtot + psyncpv - return RadpwrData( - psyncpv, pcoreradpv, pedgeradpv, pradpv - ) + return RadpwrData(psyncpv, pcoreradpv, pedgeradpv, pradpv) def psync_albajar_fidone(): From 4272f8dccf666f936200b9e478e0b49e38fe8b5b Mon Sep 17 00:00:00 2001 From: apearce Date: Mon, 3 Jun 2024 14:28:37 +0100 Subject: [PATCH 3/4] =?UTF-8?q?=E2=9C=85=20fix=20tests?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- source/fortran/physics_variables.f90 | 2 +- tests/unit/test_impurity_radiation.py | 69 --------------------------- 2 files changed, 1 insertion(+), 70 deletions(-) diff --git a/source/fortran/physics_variables.f90 b/source/fortran/physics_variables.f90 index ebb6789e5d..381e914a5b 100644 --- a/source/fortran/physics_variables.f90 +++ b/source/fortran/physics_variables.f90 @@ -620,7 +620,7 @@ module physics_variables real(dp) :: pedgeradpv !! edge radiation power per volume (MW/m3) - real(dp) :: pfuscmw + real(dp) :: pfuscmw !! charged particle fusion power (MW) real(dp) :: phiint diff --git a/tests/unit/test_impurity_radiation.py b/tests/unit/test_impurity_radiation.py index cd400c38ec..397653e412 100644 --- a/tests/unit/test_impurity_radiation.py +++ b/tests/unit/test_impurity_radiation.py @@ -10,75 +10,6 @@ def initialise_impurity_radiation(): impurity_radiation_module.init_impurity_radiation_module() impurity_radiation.initialise_imprad() - """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): From 3383639870b96ece10a08efae36b08181e9c77ce Mon Sep 17 00:00:00 2001 From: apearce Date: Tue, 4 Jun 2024 14:42:49 +0100 Subject: [PATCH 4/4] =?UTF-8?q?=E2=9C=85=20fix=20integration=20tests?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/impurity_radiation.py | 22 +++------------------- source/fortran/physics_variables.f90 | 2 +- 2 files changed, 4 insertions(+), 20 deletions(-) diff --git a/process/impurity_radiation.py b/process/impurity_radiation.py index 8003cceda6..0c513a1e04 100644 --- a/process/impurity_radiation.py +++ b/process/impurity_radiation.py @@ -469,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): """ @@ -489,13 +488,9 @@ def __init__(self, plasma_profile): 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.""" @@ -544,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 ) diff --git a/source/fortran/physics_variables.f90 b/source/fortran/physics_variables.f90 index 381e914a5b..ebb6789e5d 100644 --- a/source/fortran/physics_variables.f90 +++ b/source/fortran/physics_variables.f90 @@ -620,7 +620,7 @@ module physics_variables real(dp) :: pedgeradpv !! edge radiation power per volume (MW/m3) - real(dp) :: pfuscmw + real(dp) :: pfuscmw !! charged particle fusion power (MW) real(dp) :: phiint