From 4c9d789b6afec204dd6e4fd659f47a164138a36a Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 12 Jun 2025 09:40:07 +0100 Subject: [PATCH 1/3] patched triangularity issue --- freegs4e/equilibrium.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/freegs4e/equilibrium.py b/freegs4e/equilibrium.py index ee5c2b5..958a468 100644 --- a/freegs4e/equilibrium.py +++ b/freegs4e/equilibrium.py @@ -1438,7 +1438,7 @@ def triangularity_upper(self): R_geom = (self._sep_Rmax + self._sep_Rmin) / 2 R_minor = (self._sep_Rmax - self._sep_Rmin) / 2 - return (R_geom - self._sep_RZmin) / R_minor + return (R_geom - self._sep_RZmax) / R_minor def triangularity_lower(self): """ @@ -1459,7 +1459,7 @@ def triangularity_lower(self): R_geom = (self._sep_Rmax + self._sep_Rmin) / 2 R_minor = (self._sep_Rmax - self._sep_Rmin) / 2 - return (R_geom - self._sep_RZmax) / R_minor + return (R_geom - self._sep_RZmin) / R_minor def triangularity(self): """ @@ -1480,7 +1480,7 @@ def triangularity(self): R_geom = (self._sep_Rmax + self._sep_Rmin) / 2 R_minor = (self._sep_Rmax - self._sep_Rmin) / 2 - return (R_geom - (self._sep_RZmax + self._sep_RZmin) / 2) / R_minor + return (R_geom - 0.5 * (self._sep_RZmax + self._sep_RZmin)) / R_minor def squareness(self): """ From 1dd8ad6993bae3a1697f09e2f4a26fb4765f53f0 Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 12 Jun 2025 19:06:01 +0100 Subject: [PATCH 2/3] patching betap, betat, and li params --- freegs4e/equilibrium.py | 241 +++++++++++++++++++++++++++++++++------- 1 file changed, 200 insertions(+), 41 deletions(-) diff --git a/freegs4e/equilibrium.py b/freegs4e/equilibrium.py index 958a468..163a480 100644 --- a/freegs4e/equilibrium.py +++ b/freegs4e/equilibrium.py @@ -1648,9 +1648,124 @@ def closest_wall_point(self): def internalInductance1(self): """ - Calculates li1 plasma internal inductance according to: + Calculates plasma internal inductance according to: - li_1 = [(2 * integral(Bpol^2) dV) / (mu0 * Ip^2 * R_geom)] * [(1 + geom_elon^2)/(2 * eff_elon)], + li = ( (integral() dl * integral(Bpol) dV) / (mu0 * Ip * V) )^2, + + This definition is also used in the Fiesta equilibrium code. + + where: + - L = length of the separatrix (LCFS). + - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). + - dV = volume element. + - mu0 = magnetic permeability in vacuum. + - Ip = total plasma current. + - V = total plasma volume. + + + Parameters + ---------- + + Returns + ------- + float + Plasma internal inductance. + """ + + # volume elements + dV = 2.0 * np.pi * self.R * self.dR * self.dZ + + # mask with the core plasma + try: + dV *= self._profiles.limiter_core_mask + except AttributeError as e: + print(e) + warnings.warn( + "The core mask is not in place. You need to solve for the equilibrium first!" + ) + raise e + + # integrate B_pol over the plasma volume + integral = romb(romb(self.Bpol(self.R, self.Z) * dV)) + + return ( + (self.separatrix_length() * integral) + / (mu0 * self.plasmaCurrent() * self.plasmaVolume()) + ) ** 2 + + def internalInductance2(self): + """ + Calculates plasma internal inductance according to: + + li = ( integral(Bpol**2) dV ) / ( ( integral(Bpol) dl / L )**2 * V), + + This definition is also used in the Fiesta equilibrium code. + + where: + - L = length of the separatrix (LCFS). + - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). + - dV = volume element. + - dl = line element (along LCFS). + - mu0 = magnetic permeability in vacuum. + - Ip = total plasma current. + - V = total plasma volume. + + + Parameters + ---------- + + Returns + ------- + float + Plasma internal inductance. + """ + + # extract Bpol + B_pol = self.Bpol(self.R, self.Z) + + # volume elements + dV = 2.0 * np.pi * self.R * self.dR * self.dZ + + # mask with the core plasma + try: + dV *= self._profiles.limiter_core_mask + except AttributeError as e: + print(e) + warnings.warn( + "The core mask is not in place. You need to solve for the equilibrium first!" + ) + raise e + + # integrate B_pol squared over the volume + integral_Bpol_2 = romb(romb((B_pol**2) * dV)) + + # integrate B_pol over the LCFS + + # find arc length + separatrix = self.separatrix() + dl = np.sqrt( + np.diff(separatrix[:, 0]) ** 2 + np.diff(separatrix[:, 1]) ** 2 + ) # element-wise arc length + l = np.concatenate(([0], np.cumsum(dl))) # cumulative arc length + + # evaluate Bp on the LCFS + Bp = self.Bpol(separatrix[:, 0], separatrix[:, 1]) + + # integrate + integral_Bpol_lcfs = np.trapz(Bp, l) + + return integral_Bpol_2 / ( + (integral_Bpol_lcfs / self.separatrix_length()) ** 2 + * self.plasmaVolume() + ) + + def internalInductance3(self): + """ + Calculates plasma internal inductance according to: + + li = (2 * integral(Bpol^2) dV) / (mu0 * Ip^2 * R_geom), + + This definition is also used in the Fiesta equilibrium code. where: - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). @@ -1658,8 +1773,6 @@ def internalInductance1(self): - mu0 = magnetic permeability in vacuum. - Ip = total plasma current. - R_geom = radial coords of geometric axis. - - geom_elon = geometric elongation. - - eff_elon = effective elongation. Parameters ---------- @@ -1667,7 +1780,7 @@ def internalInductance1(self): Returns ------- float - Plasma internal inductance (li1). + Plasma internal inductance. """ # extract Bpol squared @@ -1690,18 +1803,16 @@ def internalInductance1(self): integral = romb(romb(B_pol_2 * dV)) return ( - (2 * integral) + 2 + * integral / ((mu0 * self.plasmaCurrent()) ** 2 * self.Rgeometric()) - ) * ( - (1 + self.geometricElongation() ** 2) - / (2.0 * self.effectiveElongation()) ) - def internalInductance2(self): + def internalInductance4(self): """ - Calculates li2 plasma internal inductance according to: + Calculates plasma internal inductance according to: - li_2 = (2 * integral(Bpol^2) dV) / (mu0 * Ip^2 * R_mag), + li = (2 * integral(Bpol^2) dV) / (mu0 * Ip^2 * R_mag), where: - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). @@ -1717,7 +1828,7 @@ def internalInductance2(self): Returns ------- float - Plasma internal inductance (li2). + Plasma internal inductance. """ # extract Bpol squared @@ -1745,12 +1856,11 @@ def internalInductance2(self): / ((mu0 * self.plasmaCurrent()) ** 2 * self.Rmagnetic()) ) - def internalInductance3(self): + def internalInductance5(self): """ - Calculates li3 plasma internal inductance according to: - - li_3 = (2 * integral(Bpol^2) dV) / (mu0 * Ip^2 * R_geom), + Calculates plasma internal inductance according to: + li = [(2 * integral(Bpol^2) dV) / (mu0 * Ip^2 * R_geom)] * [(1 + geom_elon^2)/(2 * eff_elon)], where: - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). @@ -1758,6 +1868,8 @@ def internalInductance3(self): - mu0 = magnetic permeability in vacuum. - Ip = total plasma current. - R_geom = radial coords of geometric axis. + - geom_elon = geometric elongation. + - eff_elon = effective elongation. Parameters ---------- @@ -1765,7 +1877,7 @@ def internalInductance3(self): Returns ------- float - Plasma internal inductance (li3). + Plasma internal inductance. """ # extract Bpol squared @@ -1788,14 +1900,61 @@ def internalInductance3(self): integral = romb(romb(B_pol_2 * dV)) return ( - 2 - * integral + (2 * integral) / ((mu0 * self.plasmaCurrent()) ** 2 * self.Rgeometric()) + ) * ( + (1 + self.geometricElongation() ** 2) + / (2.0 * self.effectiveElongation()) ) def poloidalBeta(self): """ - Calculates the poloidal beta from the following definition: + Calculates poloidal beta from the following definition: + + betap = 2 * mu0 * integral(p) dV / integral(Bpol^2) dV, + + where: + - mu0 = magnetic permeability in vacuum. + - p = plasma pressure (at each (R,Z)). + - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). + - dV = volume element + + Parameters + ---------- + + Returns + ------- + float + Returns the poloidal beta value. + """ + + # extract Bpol squared + B_pol_2 = self.Bpol(self.R, self.Z) ** 2 + + # volume elements + dV = 2.0 * np.pi * self.R * self.dR * self.dZ + + # plasma pressure + pressure = self.pressure(self.psiNRZ(self.R, self.Z)) + + # mask with the core plasma + try: + dV *= self._profiles.limiter_core_mask + except AttributeError as e: + print(e) + warnings.warn( + "The core mask is not in place. You need to solve for the equilibrium first!" + ) + raise e + + pressure_integral = romb(romb(pressure * dV)) + field_integral_pol = romb(romb(B_pol_2 * dV)) + + return 2 * mu0 * pressure_integral / field_integral_pol + + def poloidalBeta2(self): + """ + Calculates poloidal beta from the following definition: betap = (8 * pi / (mu0 * Ip^2)) * integral(p) dR dZ, @@ -1816,9 +1975,12 @@ def poloidalBeta(self): # plasma pressure pressure = self.pressure(self.psiNRZ(self.R, self.Z)) + # volume elements + dV = self.dR * self.dZ + # mask with the core plasma try: - pressure *= self._profiles.limiter_core_mask + dV *= self._profiles.limiter_core_mask except AttributeError as e: print(e) warnings.warn( @@ -1826,25 +1988,24 @@ def poloidalBeta(self): ) raise e - # calculate the poloidal beta by integrating pressure in 2D + pressure_integral = romb(romb(pressure * dV)) + return ( ((8.0 * pi) / mu0) - * romb(romb(pressure)) - * self.dR - * self.dZ + * pressure_integral / (self.plasmaCurrent() ** 2) ) - def poloidalBeta2(self): + def toroidalBeta(self): """ - Calculates an (alternative) poloidal beta from the following definition: + Calculates a toroidal beta from the following definition: - betap = 2 * mu0 * integral(p) dV / integral(Bpol^2) dV, + betat = 2 * mu0 * integral(p) dV / Btor(R_geom,Z_geom)^2 integral() dV, where: - mu0 = magnetic permeability in vacuum. - p = plasma pressure (at each (R,Z)). - - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). + - Btor(R,Z) = toroidal magnetic field (at an (R,Z)). - dV = volume element. Parameters @@ -1853,12 +2014,9 @@ def poloidalBeta2(self): Returns ------- float - Returns the poloidal beta value. + Returns the toroidal beta value. """ - # extract Bpol squared - B_pol_2 = self.Bpol(self.R, self.Z) ** 2 - # volume elements dV = 2.0 * np.pi * self.R * self.dR * self.dZ @@ -1876,11 +2034,15 @@ def poloidalBeta2(self): raise e pressure_integral = romb(romb(pressure * dV)) - field_integral_pol = romb(romb(B_pol_2 * dV)) - return 2 * mu0 * pressure_integral / field_integral_pol + # extract Btor squared + B_tor_2 = self.Btor(self.Rgeometric(), self.Zgeometric()) ** 2 - def toroidalBeta(self): + field_integral_tor = romb(romb(B_tor_2 * dV)) + + return 2 * mu0 * pressure_integral / field_integral_tor + + def toroidalBeta2(self): """ Calculates a toroidal beta from the following definition: @@ -1922,9 +2084,6 @@ def toroidalBeta(self): pressure_integral = romb(romb(pressure * dV)) - # correct for errors in Btor and core masking - np.nan_to_num(B_tor_2, copy=False) - field_integral_tor = romb(romb(B_tor_2 * dV)) return 2 * mu0 * pressure_integral / field_integral_tor @@ -1946,7 +2105,7 @@ def normalised_total_Beta(self): """ return 1.0 / ( - (1.0 / self.poloidalBeta2()) + (1.0 / self.toroidalBeta()) + (1.0 / self.poloidalBeta()) + (1.0 / self.toroidalBeta()) ) def strikepoints( From fd88168e5e68df0a41bab9700e22eb22eeac2374 Mon Sep 17 00:00:00 2001 From: Kamran Pentland Date: Mon, 16 Jun 2025 13:40:23 +0100 Subject: [PATCH 3/3] finalising definitions of betap, betat, and li --- freegs4e/equilibrium.py | 330 +++++++++++++++++++++++++++++++--------- 1 file changed, 255 insertions(+), 75 deletions(-) diff --git a/freegs4e/equilibrium.py b/freegs4e/equilibrium.py index 163a480..1a3fbf7 100644 --- a/freegs4e/equilibrium.py +++ b/freegs4e/equilibrium.py @@ -27,7 +27,7 @@ import shapely as sh from numpy import array, exp, linspace, meshgrid, pi from scipy import interpolate -from scipy.integrate import cumulative_trapezoid, romb +from scipy.integrate import cumulative_trapezoid from scipy.spatial.distance import pdist, squareform from . import critical, machine, multigrid, polygons # multigrid solver @@ -370,7 +370,7 @@ def plasmaVolume(self): raise e # integrate volume in 2D - return romb(romb(dV)) + return np.sum(dV) def plasmaBr(self, R, Z): """ @@ -1652,8 +1652,6 @@ def internalInductance1(self): li = ( (integral() dl * integral(Bpol) dV) / (mu0 * Ip * V) )^2, - This definition is also used in the Fiesta equilibrium code. - where: - L = length of the separatrix (LCFS). - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). @@ -1662,6 +1660,7 @@ def internalInductance1(self): - Ip = total plasma current. - V = total plasma volume. + This definition is also used in the Fiesta/Topeol equilibrium codes. Parameters ---------- @@ -1685,8 +1684,8 @@ def internalInductance1(self): ) raise e - # integrate B_pol over the plasma volume - integral = romb(romb(self.Bpol(self.R, self.Z) * dV)) + # first order integral + integral = np.sum(self.Bpol(self.R, self.Z) * dV) return ( (self.separatrix_length() * integral) @@ -1699,8 +1698,6 @@ def internalInductance2(self): li = ( integral(Bpol**2) dV ) / ( ( integral(Bpol) dl / L )**2 * V), - This definition is also used in the Fiesta equilibrium code. - where: - L = length of the separatrix (LCFS). - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). @@ -1710,6 +1707,7 @@ def internalInductance2(self): - Ip = total plasma current. - V = total plasma volume. + This definition is also used in the Fiesta/Topeol equilibrium codes. Parameters ---------- @@ -1720,9 +1718,6 @@ def internalInductance2(self): Plasma internal inductance. """ - # extract Bpol - B_pol = self.Bpol(self.R, self.Z) - # volume elements dV = 2.0 * np.pi * self.R * self.dR * self.dZ @@ -1737,7 +1732,7 @@ def internalInductance2(self): raise e # integrate B_pol squared over the volume - integral_Bpol_2 = romb(romb((B_pol**2) * dV)) + integral_Bpol_2 = np.sum((self.Bpol(self.R, self.Z) ** 2) * dV) # integrate B_pol over the LCFS @@ -1763,9 +1758,7 @@ def internalInductance3(self): """ Calculates plasma internal inductance according to: - li = (2 * integral(Bpol^2) dV) / (mu0 * Ip^2 * R_geom), - - This definition is also used in the Fiesta equilibrium code. + li = (2 * integral(Bpol^2) dV) / (mu0^2 * Ip^2 * R_geom), where: - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). @@ -1774,6 +1767,8 @@ def internalInductance3(self): - Ip = total plasma current. - R_geom = radial coords of geometric axis. + This definition is also used in the Fiesta/Topeol equilibrium codes. + Parameters ---------- @@ -1783,9 +1778,6 @@ def internalInductance3(self): Plasma internal inductance. """ - # extract Bpol squared - B_pol_2 = self.Bpol(self.R, self.Z) ** 2 - # volume elements dV = 2.0 * np.pi * self.R * self.dR * self.dZ @@ -1800,7 +1792,7 @@ def internalInductance3(self): raise e # integrate B_pol squared over the volume - integral = romb(romb(B_pol_2 * dV)) + integral = np.sum((self.Bpol(self.R, self.Z) ** 2) * dV) return ( 2 @@ -1812,7 +1804,7 @@ def internalInductance4(self): """ Calculates plasma internal inductance according to: - li = (2 * integral(Bpol^2) dV) / (mu0 * Ip^2 * R_mag), + li = (2 * integral(Bpol^2) dV) / (mu0^2 * Ip^2 * R_mag), where: - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). @@ -1831,9 +1823,6 @@ def internalInductance4(self): Plasma internal inductance. """ - # extract Bpol squared - B_pol_2 = self.Bpol(self.R, self.Z) ** 2 - # volume elements dV = 2.0 * np.pi * self.R * self.dR * self.dZ @@ -1848,7 +1837,7 @@ def internalInductance4(self): raise e # integrate B_pol squared over the volume - integral = romb(romb(B_pol_2 * dV)) + integral = np.sum((self.Bpol(self.R, self.Z) ** 2) * dV) return ( 2 @@ -1860,7 +1849,7 @@ def internalInductance5(self): """ Calculates plasma internal inductance according to: - li = [(2 * integral(Bpol^2) dV) / (mu0 * Ip^2 * R_geom)] * [(1 + geom_elon^2)/(2 * eff_elon)], + li = [(2 * integral(Bpol^2) dV) / (mu0^2 * Ip^2 * R_geom)] * [(1 + geom_elon^2)/(2 * eff_elon)], where: - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). @@ -1880,9 +1869,6 @@ def internalInductance5(self): Plasma internal inductance. """ - # extract Bpol squared - B_pol_2 = self.Bpol(self.R, self.Z) ** 2 - # volume elements dV = 2.0 * np.pi * self.R * self.dR * self.dZ @@ -1897,7 +1883,7 @@ def internalInductance5(self): raise e # integrate B_pol squared over the volume - integral = romb(romb(B_pol_2 * dV)) + integral = np.sum((self.Bpol(self.R, self.Z) ** 2) * dV) return ( (2 * integral) @@ -1907,7 +1893,52 @@ def internalInductance5(self): / (2.0 * self.effectiveElongation()) ) - def poloidalBeta(self): + def poloidalBeta1(self): + """ + Calculates poloidal beta from the following definition: + + betap = (8 * pi / (mu0 * Ip^2)) * integral(p) dR dZ, + + where: + - mu0 = magnetic permeability in vacuum. + - Ip = total plasma current. + - p = plasma pressure (at each (R,Z)). + + This definition is also used in the Fiesta/Topeol equilibrium codes. + + Parameters + ---------- + + Returns + ------- + float + Returns the poloidal beta value. + """ + + # volume elements + dV = self.dR * self.dZ + + # mask with the core plasma + try: + dV *= self._profiles.limiter_core_mask + except AttributeError as e: + print(e) + warnings.warn( + "The core mask is not in place. You need to solve for the equilibrium first!" + ) + raise e + + pressure_integral = np.sum( + self.pressure(self.psiNRZ(self.R, self.Z)) * dV + ) + + return ( + ((8.0 * pi) / mu0) + * pressure_integral + / (self.plasmaCurrent() ** 2) + ) + + def poloidalBeta2(self): """ Calculates poloidal beta from the following definition: @@ -1919,6 +1950,8 @@ def poloidalBeta(self): - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). - dV = volume element + This definition is also used in the Fiesta/Topeol equilibrium codes. + Parameters ---------- @@ -1928,15 +1961,9 @@ def poloidalBeta(self): Returns the poloidal beta value. """ - # extract Bpol squared - B_pol_2 = self.Bpol(self.R, self.Z) ** 2 - # volume elements dV = 2.0 * np.pi * self.R * self.dR * self.dZ - # plasma pressure - pressure = self.pressure(self.psiNRZ(self.R, self.Z)) - # mask with the core plasma try: dV *= self._profiles.limiter_core_mask @@ -1947,21 +1974,89 @@ def poloidalBeta(self): ) raise e - pressure_integral = romb(romb(pressure * dV)) - field_integral_pol = romb(romb(B_pol_2 * dV)) + pressure_integral = np.sum( + self.pressure(self.psiNRZ(self.R, self.Z)) * dV + ) + field_integral_pol = np.sum((self.Bpol(self.R, self.Z) ** 2) * dV) return 2 * mu0 * pressure_integral / field_integral_pol - def poloidalBeta2(self): + def poloidalBeta3(self): """ Calculates poloidal beta from the following definition: - betap = (8 * pi / (mu0 * Ip^2)) * integral(p) dR dZ, + betap = ( integral(p) dV / V) / (integral(Bpol) dl / L)**2, where: + - p = pressure field. + - dV = volume element. + - V = total plasma volume. + - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). + - dl = line element (along LCFS). + - L = length of the separatrix (LCFS). + + This definition is also used in the Fiesta/Topeol equilibrium codes. + + Parameters + ---------- + + Returns + ------- + float + Returns the poloidal beta value. + """ + + # volume elements + dV = 2.0 * np.pi * self.R * self.dR * self.dZ + + # mask with the core plasma + try: + dV *= self._profiles.limiter_core_mask + except AttributeError as e: + print(e) + warnings.warn( + "The core mask is not in place. You need to solve for the equilibrium first!" + ) + raise e + + # pressure integral + pressure_integral = np.sum( + self.pressure(self.psiNRZ(self.R, self.Z)) * dV + ) + + # integrate B_pol over the LCFS + # find arc length + separatrix = self.separatrix() + dl = np.sqrt( + np.diff(separatrix[:, 0]) ** 2 + np.diff(separatrix[:, 1]) ** 2 + ) # element-wise arc length + l = np.concatenate(([0], np.cumsum(dl))) # cumulative arc length + + # evaluate Bp on the LCFS + Bp = self.Bpol(separatrix[:, 0], separatrix[:, 1]) + + # integrate + integral_Bpol_lcfs = np.trapz(Bp, l) + + return (pressure_integral / self.plasmaVolume()) / ( + integral_Bpol_lcfs / self.separatrix_length() + ) ** 2 + + def poloidalBeta4(self): + """ + Calculates poloidal beta from the following definition: + + betap = 4 * integral(p) dV / (mu0 * Ip^2 * integral(R) dV / V), + + where: + - p = plasma pressure (at each (R,Z)). + - dV = volume element - mu0 = magnetic permeability in vacuum. - Ip = total plasma current. - - p = plasma pressure (at each (R,Z)). + - R = radial grid points. + - V = total plasma volume. + + This definition is also used in the Fiesta/Topeol equilibrium codes. Parameters ---------- @@ -1972,11 +2067,8 @@ def poloidalBeta2(self): Returns the poloidal beta value. """ - # plasma pressure - pressure = self.pressure(self.psiNRZ(self.R, self.Z)) - # volume elements - dV = self.dR * self.dZ + dV = 2.0 * np.pi * self.R * self.dR * self.dZ # mask with the core plasma try: @@ -1988,25 +2080,70 @@ def poloidalBeta2(self): ) raise e - pressure_integral = romb(romb(pressure * dV)) + pressure_integral = np.sum( + self.pressure(self.psiNRZ(self.R, self.Z)) * dV + ) + R_avg = np.sum(self.R * dV) / self.plasmaVolume() - return ( - ((8.0 * pi) / mu0) - * pressure_integral - / (self.plasmaCurrent() ** 2) + return (4 * pressure_integral) / ( + mu0 * self.plasmaCurrent() ** 2 * R_avg ) - def toroidalBeta(self): + def toroidalBeta1(self): """ Calculates a toroidal beta from the following definition: - betat = 2 * mu0 * integral(p) dV / Btor(R_geom,Z_geom)^2 integral() dV, + betat = 2 * mu0 * integral(p) dR dZ / integral(Btor**2) dR dZ, where: - mu0 = magnetic permeability in vacuum. - p = plasma pressure (at each (R,Z)). - Btor(R,Z) = toroidal magnetic field (at an (R,Z)). - - dV = volume element. + + This definition is also used in the Fiesta/Topeol equilibrium codes. + + Parameters + ---------- + + Returns + ------- + float + Returns the toroidal beta value. + """ + + # plasma pressure and toroidal fields + pressure = self.pressure(self.psiNRZ(self.R, self.Z)) + Btor2 = self.Btor(self.R, self.Z) ** 2 + + # mask with the core plasma + try: + pressure *= self._profiles.limiter_core_mask + Btor2 *= self._profiles.limiter_core_mask + except AttributeError as e: + print(e) + warnings.warn( + "The core mask is not in place. You need to solve for the equilibrium first!" + ) + raise e + + pressure_integral = np.sum(pressure * self.dR * self.dZ) + Btor2_integral = np.sum(Btor2 * self.dR * self.dZ) + + return 2 * mu0 * pressure_integral / Btor2_integral + + def toroidalBeta2(self): + """ + Calculates a toroidal beta from the following definition: + + betat = 2 * mu0 * integral(p) dV / integral(Btor**2) dV, + + where: + - mu0 = magnetic permeability in vacuum. + - p = plasma pressure (at each (R,Z)). + - dV = volume element + - Btor(R,Z) = toroidal magnetic field (at an (R,Z)). + + This definition is also used in the Fiesta/Topeol equilibrium codes. Parameters ---------- @@ -2020,9 +2157,6 @@ def toroidalBeta(self): # volume elements dV = 2.0 * np.pi * self.R * self.dR * self.dZ - # plasma pressure - pressure = self.pressure(self.psiNRZ(self.R, self.Z)) - # mask with the core plasma try: dV *= self._profiles.limiter_core_mask @@ -2033,26 +2167,75 @@ def toroidalBeta(self): ) raise e - pressure_integral = romb(romb(pressure * dV)) + pressure_integral = np.sum( + self.pressure(self.psiNRZ(self.R, self.Z)) * dV + ) + Btor2_int = np.sum((self.Btor(self.R, self.Z) ** 2) * dV) - # extract Btor squared - B_tor_2 = self.Btor(self.Rgeometric(), self.Zgeometric()) ** 2 + return 2 * mu0 * pressure_integral / Btor2_int - field_integral_tor = romb(romb(B_tor_2 * dV)) + def toroidalBeta3(self): + """ + Calculates a toroidal beta from the following definition: - return 2 * mu0 * pressure_integral / field_integral_tor + betat = 2 * mu0 * integral(p) dV / integral(Btor**2 + Bpol**2) dV, - def toroidalBeta2(self): + where: + - mu0 = magnetic permeability in vacuum. + - p = plasma pressure (at each (R,Z)). + - dV = volume element + - Btor(R,Z) = toroidal magnetic field (at an (R,Z)). + - Bpol(R,Z) = poloidal magnetic field (at an (R,Z)). + + This definition is also used in the Fiesta/Topeol equilibrium codes. + + Parameters + ---------- + + Returns + ------- + float + Returns the toroidal beta value. + """ + + # volume elements + dV = 2.0 * np.pi * self.R * self.dR * self.dZ + + # mask with the core plasma + try: + dV *= self._profiles.limiter_core_mask + except AttributeError as e: + print(e) + warnings.warn( + "The core mask is not in place. You need to solve for the equilibrium first!" + ) + raise e + + pressure_integral = np.sum( + self.pressure(self.psiNRZ(self.R, self.Z)) * dV + ) + B_int = np.sum( + (self.Btor(self.R, self.Z) ** 2 + self.Bpol(self.R, self.Z) ** 2) + * dV + ) + + return 2 * mu0 * pressure_integral / B_int + + def toroidalBeta4(self): """ Calculates a toroidal beta from the following definition: - betat = 2 * mu0 * integral(p) dV / integral(Btor^2) dV, + betat = 2 * mu0 * integral(p) dV / Btor_vac(R_geom,Z_geom)^2 * V, where: - mu0 = magnetic permeability in vacuum. - p = plasma pressure (at each (R,Z)). - - Btor(R,Z) = toroidal magnetic field (at each (R,Z)). - dV = volume element. + - Btor_vac(R,Z) = vacuum toroidal magnetic field (at an (R,Z)). + - (Rgeom,Zgeom) = geoemtric axis. + - V = total plasma volume. + + This definition is also used in the Fiesta/Topeol equilibrium codes. Parameters ---------- @@ -2063,15 +2246,9 @@ def toroidalBeta2(self): Returns the toroidal beta value. """ - # extract Btor squared - B_tor_2 = self.Btor(self.R, self.Z) ** 2 - # volume elements dV = 2.0 * np.pi * self.R * self.dR * self.dZ - # plasma pressure - pressure = self.pressure(self.psiNRZ(self.R, self.Z)) - # mask with the core plasma try: dV *= self._profiles.limiter_core_mask @@ -2082,11 +2259,14 @@ def toroidalBeta2(self): ) raise e - pressure_integral = romb(romb(pressure * dV)) - - field_integral_tor = romb(romb(B_tor_2 * dV)) + pressure_integral = np.sum( + self.pressure(self.psiNRZ(self.R, self.Z)) * dV + ) + Btor_vac = self.fvac() / self.Rgeometric() - return 2 * mu0 * pressure_integral / field_integral_tor + return ( + 2 * mu0 * pressure_integral / (Btor_vac**2 * self.plasmaVolume()) + ) def normalised_total_Beta(self): """ @@ -2413,7 +2593,7 @@ def solve(self, profiles, Jtor=None, psi=None, psi_bndry=None): self._updatePlasmaPsi(plasma_psi) # update plasma current - self._current = romb(romb(Jtor)) * self.dR * self.dZ + self._current = np.sum(Jtor) * self.dR * self.dZ def _updatePlasmaPsi(self, plasma_psi): """