diff --git a/freegs4e/equilibrium.py b/freegs4e/equilibrium.py index ee5c2b5..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): """ @@ -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): """ @@ -1648,18 +1648,19 @@ 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, 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. - - R_geom = radial coords of geometric axis. - - geom_elon = geometric elongation. - - eff_elon = effective elongation. + - V = total plasma volume. + + This definition is also used in the Fiesta/Topeol equilibrium codes. Parameters ---------- @@ -1667,11 +1668,55 @@ def internalInductance1(self): Returns ------- float - Plasma internal inductance (li1). + 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 + + # 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 + + # first order integral + integral = np.sum(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), + + 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. + + This definition is also used in the Fiesta/Topeol equilibrium codes. + + Parameters + ---------- + + Returns + ------- + float + Plasma internal inductance. + """ # volume elements dV = 2.0 * np.pi * self.R * self.dR * self.dZ @@ -1687,21 +1732,79 @@ def internalInductance1(self): raise e # integrate B_pol squared over the volume - integral = 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 + + # 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^2 * Ip^2 * R_geom), + + where: + - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). + - dV = volume element. + - mu0 = magnetic permeability in vacuum. + - Ip = total plasma current. + - R_geom = radial coords of geometric axis. + + This definition is also used in the Fiesta/Topeol equilibrium codes. + + 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 squared over the volume + integral = np.sum((self.Bpol(self.R, self.Z) ** 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^2 * Ip^2 * R_mag), where: - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). @@ -1717,12 +1820,9 @@ def internalInductance2(self): Returns ------- float - Plasma internal inductance (li2). + 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 @@ -1737,7 +1837,7 @@ def internalInductance2(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 @@ -1745,12 +1845,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^2 * Ip^2 * R_geom)] * [(1 + geom_elon^2)/(2 * eff_elon)], where: - Bpol(R,Z) = poloidal magnetic field (at each (R,Z)). @@ -1758,6 +1857,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,12 +1866,9 @@ def internalInductance3(self): Returns ------- float - Plasma internal inductance (li3). + 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 @@ -1785,17 +1883,19 @@ 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 - * integral + (2 * integral) / ((mu0 * self.plasmaCurrent()) ** 2 * self.Rgeometric()) + ) * ( + (1 + self.geometricElongation() ** 2) + / (2.0 * self.effectiveElongation()) ) - def poloidalBeta(self): + def poloidalBeta1(self): """ - Calculates the poloidal beta from the following definition: + Calculates poloidal beta from the following definition: betap = (8 * pi / (mu0 * Ip^2)) * integral(p) dR dZ, @@ -1804,6 +1904,8 @@ def poloidalBeta(self): - Ip = total plasma current. - p = plasma pressure (at each (R,Z)). + This definition is also used in the Fiesta/Topeol equilibrium codes. + Parameters ---------- @@ -1813,12 +1915,12 @@ def poloidalBeta(self): Returns the poloidal beta value. """ - # 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,18 +1928,19 @@ def poloidalBeta(self): ) raise e - # calculate the poloidal beta by integrating pressure in 2D + pressure_integral = np.sum( + self.pressure(self.psiNRZ(self.R, self.Z)) * dV + ) + return ( ((8.0 * pi) / mu0) - * romb(romb(pressure)) - * self.dR - * self.dZ + * pressure_integral / (self.plasmaCurrent() ** 2) ) def poloidalBeta2(self): """ - Calculates an (alternative) poloidal beta from the following definition: + Calculates poloidal beta from the following definition: betap = 2 * mu0 * integral(p) dV / integral(Bpol^2) dV, @@ -1845,7 +1948,9 @@ def poloidalBeta2(self): - 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. + - dV = volume element + + This definition is also used in the Fiesta/Topeol equilibrium codes. Parameters ---------- @@ -1856,14 +1961,114 @@ def poloidalBeta2(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 + + # 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 + ) + field_integral_pol = np.sum((self.Bpol(self.R, self.Z) ** 2) * dV) + + return 2 * mu0 * pressure_integral / field_integral_pol + + def poloidalBeta3(self): + """ + Calculates poloidal beta from the following definition: + + 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 - # 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 + 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. + - R = radial grid points. + - V = total plasma volume. + + 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: @@ -1875,22 +2080,27 @@ def poloidalBeta2(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 + ) + R_avg = np.sum(self.R * dV) / self.plasmaVolume() - return 2 * mu0 * pressure_integral / field_integral_pol + 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 / integral(Btor^2) 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 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 ---------- @@ -1901,14 +2111,95 @@ def toroidalBeta(self): Returns the toroidal beta value. """ - # extract Btor squared - B_tor_2 = self.Btor(self.R, self.Z) ** 2 + # 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 + ---------- + + Returns + ------- + float + Returns the toroidal beta value. + """ # 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 = np.sum( + self.pressure(self.psiNRZ(self.R, self.Z)) * dV + ) + Btor2_int = np.sum((self.Btor(self.R, self.Z) ** 2) * dV) + + return 2 * mu0 * pressure_integral / Btor2_int + + def toroidalBeta3(self): + """ + Calculates a toroidal beta from the following definition: + + betat = 2 * mu0 * integral(p) dV / integral(Btor**2 + Bpol**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)). + - 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: @@ -1920,14 +2211,62 @@ def toroidalBeta(self): ) raise e - pressure_integral = romb(romb(pressure * dV)) + 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 / Btor_vac(R_geom,Z_geom)^2 * V, + + where: + - mu0 = magnetic permeability in vacuum. + - p = plasma pressure (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 + ---------- + + Returns + ------- + float + Returns the toroidal beta value. + """ + + # volume elements + dV = 2.0 * np.pi * self.R * self.dR * self.dZ - # correct for errors in Btor and core masking - np.nan_to_num(B_tor_2, copy=False) + # 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 - 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): """ @@ -1946,7 +2285,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( @@ -2254,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): """