From 1de69e38d7ba26425456d36e39bc2a5a2d57e5eb Mon Sep 17 00:00:00 2001 From: ym1906 Date: Mon, 9 Sep 2024 16:58:23 +0100 Subject: [PATCH 1/8] Upated bootstrap sauter to use existing plasma profile instead of calc one --- process/physics.py | 43 +++++++++++-------------------------------- 1 file changed, 11 insertions(+), 32 deletions(-) diff --git a/process/physics.py b/process/physics.py index e616abff60..76f97d0e6a 100644 --- a/process/physics.py +++ b/process/physics.py @@ -4552,8 +4552,8 @@ def bootstrap_fraction_nevins( aibs = 2.5 * betae0 * rmajor * bt * q95 * ainteg return 1.0e6 * aibs / plascur - @staticmethod - def bootstrap_fraction_sauter(): + + def bootstrap_fraction_sauter(self): """Bootstrap current fraction from Sauter et al scaling author: P J Knight, CCFE, Culham Science Centre None @@ -4566,37 +4566,19 @@ def bootstrap_fraction_sauter(): O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) Physics of Plasmas 9 (2002) 5140 """ - NR = 200 - - roa = np.arange(1, NR + 1, step=1) / NR - + NR = self.plasma_profile.profile_size + + roa=self.plasma_profile.neprofile.profile_x rho = np.sqrt(physics_variables.xarea / np.pi) * roa + sqeps = np.sqrt(roa * (physics_variables.rminor / physics_variables.rmajor)) - - ne = 1e-19 * np.vectorize( - lambda r: profiles_module.nprofile( - r, - physics_variables.rhopedn, - physics_variables.ne0, - physics_variables.neped, - physics_variables.nesep, - physics_variables.alphan, - ) - )(roa) + + ne=self.plasma_profile.neprofile.profile_y*1e-19 ni = (physics_variables.dnitot / physics_variables.dene) * ne - tempe = np.vectorize( - lambda r: profiles_module.tprofile( - r, - physics_variables.rhopedt, - physics_variables.te0, - physics_variables.teped, - physics_variables.tesep, - physics_variables.alphat, - physics_variables.tbeta, - ) - )(roa) + + tempe=self.plasma_profile.teprofile.profile_y # tempi = (physics_variables.ti / physics_variables.te) * tempe - + zef = np.full_like(tempi, physics_variables.zeff) # Flat Zeff profile assumed # mu = 1/safety factor @@ -4620,13 +4602,11 @@ def bootstrap_fraction_sauter(): # Looping from 2 because dcsa etc should return 0 @ j == 1 nr_rng = np.arange(2, NR) nr_rng_1 = nr_rng - 1 - drho = rho[nr_rng] - rho[nr_rng_1] da = 2 * np.pi * rho[nr_rng_1] * drho # area of annulus dlogte_drho = (np.log(tempe[nr_rng]) - np.log(tempe[nr_rng_1])) / drho dlogti_drho = (np.log(tempi[nr_rng]) - np.log(tempi[nr_rng_1])) / drho dlogne_drho = (np.log(ne[nr_rng]) - np.log(ne[nr_rng_1])) / drho - jboot = ( 0.5 * ( @@ -4689,7 +4669,6 @@ def bootstrap_fraction_sauter(): * mu[nr_rng_1] ) ) # A/m2 - return np.sum(da * jboot, axis=0) / physics_variables.plascur @staticmethod From 1f7b0e638d4c740cf083d61fa3a32442752aad72 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Mon, 9 Sep 2024 17:04:05 +0100 Subject: [PATCH 2/8] Remove unused profiles_module import --- process/physics.py | 1 - 1 file changed, 1 deletion(-) diff --git a/process/physics.py b/process/physics.py index 76f97d0e6a..4c3939bee3 100644 --- a/process/physics.py +++ b/process/physics.py @@ -24,7 +24,6 @@ numerics, stellarator_variables, process_output as po, - profiles_module, ) From 311aac311e86744f74281be604b18efa64290dee Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 10 Sep 2024 08:11:40 +0100 Subject: [PATCH 3/8] black fixes --- process/physics.py | 64 ++++++++++++++++------------------------------ 1 file changed, 22 insertions(+), 42 deletions(-) diff --git a/process/physics.py b/process/physics.py index 4c3939bee3..97265d2818 100644 --- a/process/physics.py +++ b/process/physics.py @@ -258,9 +258,7 @@ def vscalc( aeps = (1.0 + 1.81 * np.sqrt(eps) + 2.05 * eps) * np.log(8.0 / eps) - ( 2.0 + 9.25 * np.sqrt(eps) - 1.21 * eps ) - beps = ( - 0.73 * np.sqrt(eps) * (1.0 + 2.0 * eps**4 - 6.0 * eps**5 + 3.7 * eps**6) - ) + beps = 0.73 * np.sqrt(eps) * (1.0 + 2.0 * eps**4 - 6.0 * eps**5 + 3.7 * eps**6) rlpext = rmajor * rmu0 * aeps * (1.0 - eps) / (1.0 - eps + beps * kappa) rlp = rlpext + rlpint @@ -932,9 +930,7 @@ def physics(self): ) betat = ( - physics_variables.beta - * physics_variables.btot**2 - / physics_variables.bt**2 + physics_variables.beta * physics_variables.btot**2 / physics_variables.bt**2 ) current_drive_variables.bscf_nevins = ( current_drive_variables.cboot @@ -1628,18 +1624,14 @@ def culdlm( # This applies to the density at the plasma edge, so must be scaled # to give the density limit applying to the average plasma density. - dlimit[0] = ( - 1.54e20 * qperp**0.43 * bt**0.31 / (q95 * rmajor) ** 0.45 - ) / prn1 + dlimit[0] = (1.54e20 * qperp**0.43 * bt**0.31 / (q95 * rmajor) ** 0.45) / prn1 # Borrass density limit model for ITER (I) # This applies to the density at the plasma edge, so must be scaled # to give the density limit applying to the average plasma density. # Borrass et al, ITER-TN-PH-9-6 (1989) - dlimit[1] = ( - 1.8e20 * qperp**0.53 * bt**0.31 / (q95 * rmajor) ** 0.22 - ) / prn1 + dlimit[1] = (1.8e20 * qperp**0.53 * bt**0.31 / (q95 * rmajor) ** 0.22) / prn1 # Borrass density limit model for ITER (II) # This applies to the density at the plasma edge, so must be scaled @@ -1647,9 +1639,7 @@ def culdlm( # This formula is (almost) identical to that in the original routine # denlim (now deleted). - dlimit[2] = ( - 0.5e20 * qperp**0.57 * bt**0.31 / (q95 * rmajor) ** 0.09 - ) / prn1 + dlimit[2] = (0.5e20 * qperp**0.57 * bt**0.31 / (q95 * rmajor) ** 0.09) / prn1 # JET edge radiation density limit model # This applies to the density at the plasma edge, so must be scaled @@ -2118,10 +2108,7 @@ def culcur( 0.5 * (1.17 - 0.65 * eps) / ((1.0 - eps * eps) ** 2) - * ( - 1.0 - + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3) - ) + * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) ) elif icurr in [5, 6]: # Todd empirical scalings fq = ( @@ -2910,9 +2897,9 @@ def outplas(self): for imp in range(impurity_radiation_module.nimp): # MDK Update fimp, as this will make the ITV output work correctly. - impurity_radiation_module.fimp[ - imp - ] = impurity_radiation_module.impurity_arr_frac[imp] + impurity_radiation_module.fimp[imp] = ( + impurity_radiation_module.impurity_arr_frac[imp] + ) str1 = ( f2py_compatible_to_string( impurity_radiation_module.impurity_arr_label[imp] @@ -4551,7 +4538,6 @@ def bootstrap_fraction_nevins( aibs = 2.5 * betae0 * rmajor * bt * q95 * ainteg return 1.0e6 * aibs / plascur - def bootstrap_fraction_sauter(self): """Bootstrap current fraction from Sauter et al scaling author: P J Knight, CCFE, Culham Science Centre @@ -4566,25 +4552,24 @@ def bootstrap_fraction_sauter(self): Physics of Plasmas 9 (2002) 5140 """ NR = self.plasma_profile.profile_size - - roa=self.plasma_profile.neprofile.profile_x + + roa = self.plasma_profile.neprofile.profile_x rho = np.sqrt(physics_variables.xarea / np.pi) * roa - + sqeps = np.sqrt(roa * (physics_variables.rminor / physics_variables.rmajor)) - - ne=self.plasma_profile.neprofile.profile_y*1e-19 + + ne = self.plasma_profile.neprofile.profile_y * 1e-19 ni = (physics_variables.dnitot / physics_variables.dene) * ne - - tempe=self.plasma_profile.teprofile.profile_y # + + tempe = self.plasma_profile.teprofile.profile_y tempi = (physics_variables.ti / physics_variables.te) * tempe - - zef = np.full_like(tempi, physics_variables.zeff) # Flat Zeff profile assumed + # Flat Zeff profile assumed + zef = np.full_like(tempi, physics_variables.zeff) # mu = 1/safety factor # Parabolic q profile assumed mu = 1 / ( - physics_variables.q0 - + (physics_variables.q - physics_variables.q0) * roa**2 + physics_variables.q0 + (physics_variables.q - physics_variables.q0) * roa**2 ) amain = np.full_like(mu, physics_variables.afuel) zmain = np.full_like(mu, 1.0 + physics_variables.fhe3) @@ -4668,6 +4653,7 @@ def bootstrap_fraction_sauter(self): * mu[nr_rng_1] ) ) # A/m2 + return np.sum(da * jboot, axis=0) / physics_variables.plascur @staticmethod @@ -5079,11 +5065,7 @@ def pcond( * np.sqrt(kappa95) * denfac / powerht**0.4e0 - * ( - zeff**2 - * pcur**4 - / (rmajor * rminor * qstar**3 * kappa95**1.5e0) - ) + * (zeff**2 * pcur**4 / (rmajor * rminor * qstar**3 * kappa95**1.5e0)) ** 0.08e0 ) @@ -5538,9 +5520,7 @@ def pcond( # Table 4. (Issue #311) # Note that aspect ratio and M (afuel) do not appear, and B (bt) only # appears in the "saturation factor" h. - h = dnla19**0.448e0 / ( - 1.0e0 + np.exp(-9.403e0 * (bt / dnla19) ** 1.365e0) - ) + h = dnla19**0.448e0 / (1.0e0 + np.exp(-9.403e0 * (bt / dnla19) ** 1.365e0)) tauee = ( hfact * 0.0367e0 From c3d183ef054590be9435a180857a68e8fc1365fe Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 10 Sep 2024 09:38:02 +0100 Subject: [PATCH 4/8] reverted to static method, added plasma_profile requirement --- process/physics.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/process/physics.py b/process/physics.py index 97265d2818..2dce524978 100644 --- a/process/physics.py +++ b/process/physics.py @@ -982,7 +982,8 @@ def physics(self): current_drive_variables.pscf_scene = ps_fraction_scene(physics_variables.beta) current_drive_variables.bscf_sauter = ( - current_drive_variables.cboot * self.bootstrap_fraction_sauter() + current_drive_variables.cboot + * self.bootstrap_fraction_sauter(self.plasma_profile) ) current_drive_variables.bscf_sakai = ( @@ -4538,7 +4539,8 @@ def bootstrap_fraction_nevins( aibs = 2.5 * betae0 * rmajor * bt * q95 * ainteg return 1.0e6 * aibs / plascur - def bootstrap_fraction_sauter(self): + @staticmethod + def bootstrap_fraction_sauter(plasma_profile): """Bootstrap current fraction from Sauter et al scaling author: P J Knight, CCFE, Culham Science Centre None @@ -4551,17 +4553,17 @@ def bootstrap_fraction_sauter(self): O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) Physics of Plasmas 9 (2002) 5140 """ - NR = self.plasma_profile.profile_size + NR = plasma_profile.profile_size - roa = self.plasma_profile.neprofile.profile_x + roa = plasma_profile.neprofile.profile_x rho = np.sqrt(physics_variables.xarea / np.pi) * roa sqeps = np.sqrt(roa * (physics_variables.rminor / physics_variables.rmajor)) - ne = self.plasma_profile.neprofile.profile_y * 1e-19 + ne = plasma_profile.neprofile.profile_y * 1e-19 ni = (physics_variables.dnitot / physics_variables.dene) * ne - tempe = self.plasma_profile.teprofile.profile_y + tempe = plasma_profile.teprofile.profile_y tempi = (physics_variables.ti / physics_variables.te) * tempe # Flat Zeff profile assumed zef = np.full_like(tempi, physics_variables.zeff) @@ -4653,7 +4655,6 @@ def bootstrap_fraction_sauter(self): * mu[nr_rng_1] ) ) # A/m2 - return np.sum(da * jboot, axis=0) / physics_variables.plascur @staticmethod From e4a9b1a9e58dfe62ef36fe47059ba6831f4e86c5 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 10 Sep 2024 09:38:30 +0100 Subject: [PATCH 5/8] updated bootstrap to use new version of function --- tests/unit/test_physics.py | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 32bbac8ce6..31aa7d06bb 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -350,33 +350,33 @@ class BootstrapFractionSauterParam(NamedTuple): "bootstrapfractionsauterparam", ( BootstrapFractionSauterParam( - dnitot=6.6125550702454276e19, + dnitot=7.1297522422781575e+19, rminor=2.6666666666666665, tesep=0.10000000000000001, - ti=12, + ti=12.570861186498382, triang=0.5, q0=1, afuel=2.5, - zeff=2.0909945616489103, - rhopedn=0.94000000000000006, - bt=5.7000000000000002, - plascur=18398455.678867526, + zeff=2.5211399464385624, + rhopedn=0.9400000000000001, + bt=5.326133750416047, + plascur=16528278.760008096, xarea=38.39822223637151, fhe3=0, teped=5.5, - dene=7.5e19, - te=12, + dene=8.016748468651018e+19, + te=12.570861186498382, rmajor=8, q=3.5, - nesep=4.1177885154594193e19, - te0=24.402321098330372, - neped=7.000240476281013e19, + nesep=3.6992211545476006e+19, + te0=25.986118047669795, + neped=6.2886759627309195e+19, tbeta=2, - ne0=8.515060981068918e19, + ne0=1.054474759840606e+20, alphan=1, - rhopedt=0.94000000000000006, + rhopedt=0.9400000000000001, alphat=1.45, - expected_bfs=0.27635918746616817, + expected_bfs=0.4110838247346975, ), ), ) @@ -462,8 +462,8 @@ def test_bootstrap_fraction_sauter(bootstrapfractionsauterparam, monkeypatch, ph monkeypatch.setattr( physics_variables, "alphat", bootstrapfractionsauterparam.alphat ) - - bfs = physics.bootstrap_fraction_sauter() + physics.plasma_profile.run() + bfs = physics.bootstrap_fraction_sauter(physics.plasma_profile) assert bfs == pytest.approx(bootstrapfractionsauterparam.expected_bfs) From b0c68110aa776113b9c842c55df6d5d2c55f54f0 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 10 Sep 2024 10:03:29 +0100 Subject: [PATCH 6/8] Black fix --- tests/unit/test_physics.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 31aa7d06bb..dd95d87dae 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -350,7 +350,7 @@ class BootstrapFractionSauterParam(NamedTuple): "bootstrapfractionsauterparam", ( BootstrapFractionSauterParam( - dnitot=7.1297522422781575e+19, + dnitot=7.1297522422781575e19, rminor=2.6666666666666665, tesep=0.10000000000000001, ti=12.570861186498382, @@ -364,15 +364,15 @@ class BootstrapFractionSauterParam(NamedTuple): xarea=38.39822223637151, fhe3=0, teped=5.5, - dene=8.016748468651018e+19, + dene=8.016748468651018e19, te=12.570861186498382, rmajor=8, q=3.5, - nesep=3.6992211545476006e+19, + nesep=3.6992211545476006e19, te0=25.986118047669795, - neped=6.2886759627309195e+19, + neped=6.2886759627309195e19, tbeta=2, - ne0=1.054474759840606e+20, + ne0=1.054474759840606e20, alphan=1, rhopedt=0.9400000000000001, alphat=1.45, From 56460af55727f5ec99b4018896b915ffb151d0cc Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 10 Sep 2024 10:20:22 +0100 Subject: [PATCH 7/8] trying to fix the black fail --- process/physics.py | 1 + 1 file changed, 1 insertion(+) diff --git a/process/physics.py b/process/physics.py index 2dce524978..98b1995e3e 100644 --- a/process/physics.py +++ b/process/physics.py @@ -4655,6 +4655,7 @@ def bootstrap_fraction_sauter(plasma_profile): * mu[nr_rng_1] ) ) # A/m2 + return np.sum(da * jboot, axis=0) / physics_variables.plascur @staticmethod From f7891f771419c2f87368fa17c722fc59f7535323 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 10 Sep 2024 13:01:02 +0100 Subject: [PATCH 8/8] black fixed --- process/physics.py | 44 ++++++++++++++++++++++++++++++++------------ 1 file changed, 32 insertions(+), 12 deletions(-) diff --git a/process/physics.py b/process/physics.py index 98b1995e3e..0aa4345617 100644 --- a/process/physics.py +++ b/process/physics.py @@ -258,7 +258,9 @@ def vscalc( aeps = (1.0 + 1.81 * np.sqrt(eps) + 2.05 * eps) * np.log(8.0 / eps) - ( 2.0 + 9.25 * np.sqrt(eps) - 1.21 * eps ) - beps = 0.73 * np.sqrt(eps) * (1.0 + 2.0 * eps**4 - 6.0 * eps**5 + 3.7 * eps**6) + beps = ( + 0.73 * np.sqrt(eps) * (1.0 + 2.0 * eps**4 - 6.0 * eps**5 + 3.7 * eps**6) + ) rlpext = rmajor * rmu0 * aeps * (1.0 - eps) / (1.0 - eps + beps * kappa) rlp = rlpext + rlpint @@ -930,7 +932,9 @@ def physics(self): ) betat = ( - physics_variables.beta * physics_variables.btot**2 / physics_variables.bt**2 + physics_variables.beta + * physics_variables.btot**2 + / physics_variables.bt**2 ) current_drive_variables.bscf_nevins = ( current_drive_variables.cboot @@ -1625,14 +1629,18 @@ def culdlm( # This applies to the density at the plasma edge, so must be scaled # to give the density limit applying to the average plasma density. - dlimit[0] = (1.54e20 * qperp**0.43 * bt**0.31 / (q95 * rmajor) ** 0.45) / prn1 + dlimit[0] = ( + 1.54e20 * qperp**0.43 * bt**0.31 / (q95 * rmajor) ** 0.45 + ) / prn1 # Borrass density limit model for ITER (I) # This applies to the density at the plasma edge, so must be scaled # to give the density limit applying to the average plasma density. # Borrass et al, ITER-TN-PH-9-6 (1989) - dlimit[1] = (1.8e20 * qperp**0.53 * bt**0.31 / (q95 * rmajor) ** 0.22) / prn1 + dlimit[1] = ( + 1.8e20 * qperp**0.53 * bt**0.31 / (q95 * rmajor) ** 0.22 + ) / prn1 # Borrass density limit model for ITER (II) # This applies to the density at the plasma edge, so must be scaled @@ -1640,7 +1648,9 @@ def culdlm( # This formula is (almost) identical to that in the original routine # denlim (now deleted). - dlimit[2] = (0.5e20 * qperp**0.57 * bt**0.31 / (q95 * rmajor) ** 0.09) / prn1 + dlimit[2] = ( + 0.5e20 * qperp**0.57 * bt**0.31 / (q95 * rmajor) ** 0.09 + ) / prn1 # JET edge radiation density limit model # This applies to the density at the plasma edge, so must be scaled @@ -2109,7 +2119,10 @@ def culcur( 0.5 * (1.17 - 0.65 * eps) / ((1.0 - eps * eps) ** 2) - * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) + * ( + 1.0 + + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3) + ) ) elif icurr in [5, 6]: # Todd empirical scalings fq = ( @@ -2898,9 +2911,9 @@ def outplas(self): for imp in range(impurity_radiation_module.nimp): # MDK Update fimp, as this will make the ITV output work correctly. - impurity_radiation_module.fimp[imp] = ( - impurity_radiation_module.impurity_arr_frac[imp] - ) + impurity_radiation_module.fimp[ + imp + ] = impurity_radiation_module.impurity_arr_frac[imp] str1 = ( f2py_compatible_to_string( impurity_radiation_module.impurity_arr_label[imp] @@ -4571,7 +4584,8 @@ def bootstrap_fraction_sauter(plasma_profile): # mu = 1/safety factor # Parabolic q profile assumed mu = 1 / ( - physics_variables.q0 + (physics_variables.q - physics_variables.q0) * roa**2 + physics_variables.q0 + + (physics_variables.q - physics_variables.q0) * roa**2 ) amain = np.full_like(mu, physics_variables.afuel) zmain = np.full_like(mu, 1.0 + physics_variables.fhe3) @@ -5067,7 +5081,11 @@ def pcond( * np.sqrt(kappa95) * denfac / powerht**0.4e0 - * (zeff**2 * pcur**4 / (rmajor * rminor * qstar**3 * kappa95**1.5e0)) + * ( + zeff**2 + * pcur**4 + / (rmajor * rminor * qstar**3 * kappa95**1.5e0) + ) ** 0.08e0 ) @@ -5522,7 +5540,9 @@ def pcond( # Table 4. (Issue #311) # Note that aspect ratio and M (afuel) do not appear, and B (bt) only # appears in the "saturation factor" h. - h = dnla19**0.448e0 / (1.0e0 + np.exp(-9.403e0 * (bt / dnla19) ** 1.365e0)) + h = dnla19**0.448e0 / ( + 1.0e0 + np.exp(-9.403e0 * (bt / dnla19) ** 1.365e0) + ) tauee = ( hfact * 0.0367e0