From 27b919e73e4e3123e305ca6ebca7a5eb5bf8a572 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 5 Jul 2024 15:38:19 -0700 Subject: [PATCH 1/9] Save Jpar0 to output tokamak grids Reads ffprime, pprime from geqdsk file; passes through to regions, calculates Jpar0 and saves to the output grid --- hypnotoad/cases/tokamak.py | 64 ++++++++++++++++++++++++++++++++++--- hypnotoad/core/mesh.py | 24 +++++++++++++- hypnotoad/geqdsk/_geqdsk.py | 2 ++ 3 files changed, 84 insertions(+), 6 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index ea5c98b9..4d9d27a2 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -317,6 +317,8 @@ def __init__( psi1D, fpol1D, pressure=None, + fprime=None, + pprime=None, wall=None, psi_axis_gfile=None, psi_bdry_gfile=None, @@ -340,6 +342,8 @@ def __init__( -------- pressure[nf] = 1D array of pressure as a function of psi1D [Pa] + fprime[nf] = 1D array of df / dpsi + pprime[nf] = 1D array of dp / dpsi . If set then pressure must also be set. wall = [(R0,Z0), (R1, Z1), ...] A list of coordinate pairs, defining the vessel wall. @@ -362,6 +366,9 @@ def __init__( """ + if pprime is not None: + assert pressure is not None + self.user_options = self.user_options_factory.create(settings) if self.user_options.reverse_current: @@ -409,12 +416,22 @@ def __init__( # fpol constant in SOL fpol1D = np.concatenate([fpol1D, np.full(psiSOL.shape, fpol1D[-1])]) + if fprime is not None: + fprime = np.concatenate([fprime, np.full(psiSOL.shape, 0.0)]) + if pressure is not None: # Use an exponential decay for the pressure, based on # the value and gradient at the plasma edge p0 = pressure[-1] # p = p0 * exp( (psi - psi0) * dpdpsi / p0) - pressure = np.concatenate([pressure, p0 * np.exp(psiSOL * dpdpsi / p0)]) + pressure = np.concatenate( + [pressure, p0 * np.exp((psiSOL - psi1D[-1]) * dpdpsi / p0)] + ) + + if pprime is not None: + pprime = np.concatenate( + [pprime, dpdpsi * np.exp((psiSOL - psi1D[-1]) * dpdpsi / p0)] + ) self.magneticFunctionsFromGrid( R1D, Z1D, psi2D, self.user_options.psi_interpolation_method @@ -434,7 +451,12 @@ def __init__( # ext=3 specifies that boundary values are used outside range # Spline representing the derivative of f - self.fprime_spl = self.f_spl.derivative() + if fprime is not None: + self.fprime_spl = interpolate.InterpolatedUnivariateSpline( + psi1D * self.f_psi_sign, fprime, ext=3 + ) + else: + self.fprime_spl = self.f_spl.derivative() else: self.f_spl = lambda psi: 0.0 self.fprime_spl = lambda psi: 0.0 @@ -444,9 +466,17 @@ def __init__( self.p_spl = interpolate.InterpolatedUnivariateSpline( psi1D * self.f_psi_sign, pressure, ext=3 ) + + if pprime is not None: + self.pprime_spl = interpolate.InterpolatedUnivariateSpline( + psi1D * self.f_psi_sign, pprime, ext=3 + ) + else: + self.pprime_spl = self.p_spl.derivative() else: # If no pressure, then not output to grid file self.p_spl = None + self.pprime_spl = None # Find critical points (O- and X-points) R2D, Z2D = np.meshgrid(R1D, Z1D, indexing="ij") @@ -1621,9 +1651,14 @@ def createRegionObjects(self, all_regions, segments): eqreg.pressure = lambda psi: self.pressure( leg_psi + sign * abs(psi - leg_psi) ) + # Set pprime and fprime to zero so Jpar0 = 0 + eqreg.pprime = lambda psi: 0.0 + eqreg.fprime = lambda psi: 0.0 else: - # Core region, so use the core pressure + # Core region, so use the core profiles eqreg.pressure = self.pressure + eqreg.pprime = self.pprime + eqreg.fprime = self.fpolprime region_objects[name] = eqreg # The region objects need to be sorted, so that the @@ -1677,8 +1712,12 @@ def fpol(self, psi): @Equilibrium.handleMultiLocationArray def fpolprime(self, psi): - """psi-derivative of fpol""" - return self.fprime_spl(psi * self.f_psi_sign) + """psi-derivative of fpol + Note: Zero outside core.""" + fprime = self.fprime_spl(psi * self.f_psi_sign) + psinorm = (psi - self.psi_axis) / (self.psi_bdry - self.psi_axis) + fprime[psinorm > 1.0] = 0.0 + return fprime @Equilibrium.handleMultiLocationArray def pressure(self, psi): @@ -1687,6 +1726,17 @@ def pressure(self, psi): return None return self.p_spl(psi * self.f_psi_sign) + @Equilibrium.handleMultiLocationArray + def pprime(self, psi): + """psi-derivative of plasma pressure + Note: Zero outside the core""" + if self.pprime_spl is None: + return None + pprime = self.pprime_spl(psi * self.f_psi_sign) + psinorm = (psi - self.psi_axis) / (self.psi_bdry - self.psi_axis) + pprime[psinorm > 1.0] = 0.0 + return pprime + @property def Bt_axis(self): """Calculate toroidal field on axis""" @@ -1751,6 +1801,8 @@ def read_geqdsk( pressure = data["pres"] fpol = data["fpol"] + fprime = data["ffprime"] / fpol + pprime = data["pprime"] result = TokamakEquilibrium( R1D, @@ -1761,6 +1813,8 @@ def read_geqdsk( psi_bdry_gfile=psi_bdry_gfile, psi_axis_gfile=psi_axis_gfile, pressure=pressure, + fprime=fprime, + pprime=pprime, wall=wall, make_regions=make_regions, settings=settings, diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 8f1c1c38..dae2a148 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -902,6 +902,25 @@ def geometry1(self): self.Bxy = numpy.sqrt(self.Bpxy**2 + self.Btxy**2) + if hasattr( + self.meshParent.equilibrium.regions[self.equilibriumRegion.name], "pprime" + ) and hasattr( + self.meshParent.equilibrium.regions[self.equilibriumRegion.name], "fprime" + ): + # Calculate parallel current density from p' and f' + + pprime = self.meshParent.equilibrium.regions[ + self.equilibriumRegion.name + ].pprime(self.psixy) + fprime = self.meshParent.equilibrium.regions[ + self.equilibriumRegion.name + ].fprime(self.psixy) + + mu0 = 4e-7 * numpy.pi + self.Jpar0 = ( + self.Bxy * fprime / mu0 + self.Rxy * self.Btxy * pprime / self.Bxy + ) + def geometry2(self): """ Continuation of geometry1(), but needs neighbours to have calculated Bp so called @@ -3657,9 +3676,12 @@ def addFromRegionsXArray(name): addFromRegions("bxcvy") addFromRegions("bxcvz") - if hasattr(next(iter(self.equilibrium.regions.values())), "pressure"): + if hasattr(next(iter(self.regions.values())), "pressure"): addFromRegions("pressure") + if hasattr(next(iter(self.regions.values())), "Jpar0"): + addFromRegions("Jpar0") + # Penalty mask self.penalty_mask = BoutArray( numpy.zeros((self.nx, self.ny)), diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index ed643216..5d3695a2 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -195,6 +195,8 @@ def read(fh, cocos=1): fpol 1D array of f(psi)=R*Bt [meter-Tesla] pres 1D array of p(psi) [Pascals] + ffprime 1D array of ff'(psi) + pprime 1D array of p'(psi) qpsi 1D array of q(psi) psi 2D array (nx,ny) of poloidal flux From 88a8d3869b0e485b39603fd08a02c07736787732 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 31 Jul 2024 14:42:07 -0700 Subject: [PATCH 2/9] Fix sign of Jpar0 H.Seto noted that in the COCOS +1 convention that Hypnotoad assumes, Jpar has a minus sign. --- hypnotoad/core/mesh.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index dae2a148..5cb996fb 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -908,6 +908,8 @@ def geometry1(self): self.meshParent.equilibrium.regions[self.equilibriumRegion.name], "fprime" ): # Calculate parallel current density from p' and f' + # Note: COCOS +1 convention is assumed + # so mu0 * Jpar = -B f' - mu0 p' f / B pprime = self.meshParent.equilibrium.regions[ self.equilibriumRegion.name @@ -917,7 +919,7 @@ def geometry1(self): ].fprime(self.psixy) mu0 = 4e-7 * numpy.pi - self.Jpar0 = ( + self.Jpar0 = -( self.Bxy * fprime / mu0 + self.Rxy * self.Btxy * pprime / self.Bxy ) From 642e0d054a436941c12aa2e2f97d082828f9d539 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 31 Jul 2024 14:54:00 -0700 Subject: [PATCH 3/9] tokamak: Ignore missing psi_bdry in ffprime and pprime functions Used to detect when outside the core plasma. Test cases have no X-points so no psi_bdry. --- hypnotoad/cases/tokamak.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 4d9d27a2..a148455d 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -1715,8 +1715,11 @@ def fpolprime(self, psi): """psi-derivative of fpol Note: Zero outside core.""" fprime = self.fprime_spl(psi * self.f_psi_sign) - psinorm = (psi - self.psi_axis) / (self.psi_bdry - self.psi_axis) - fprime[psinorm > 1.0] = 0.0 + try: + psinorm = (psi - self.psi_axis) / (self.psi_bdry - self.psi_axis) + fprime[psinorm > 1.0] = 0.0 + except: + pass return fprime @Equilibrium.handleMultiLocationArray @@ -1733,8 +1736,11 @@ def pprime(self, psi): if self.pprime_spl is None: return None pprime = self.pprime_spl(psi * self.f_psi_sign) - psinorm = (psi - self.psi_axis) / (self.psi_bdry - self.psi_axis) - pprime[psinorm > 1.0] = 0.0 + try: + psinorm = (psi - self.psi_axis) / (self.psi_bdry - self.psi_axis) + pprime[psinorm > 1.0] = 0.0 + except: + pass return pprime @property From 2f024c37977f1df0cfafbd9a7afb7fced3b69723 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 31 Jul 2024 16:58:09 -0700 Subject: [PATCH 4/9] geqdsk writer: Calculate ffprime and pprime if not set Use Scipy.interpolate to differentiate w.r.t psi and save to the output geqdsk file. --- hypnotoad/cases/tokamak.py | 10 ++++------ hypnotoad/geqdsk/_geqdsk.py | 24 ++++++++++++++++++++++-- 2 files changed, 26 insertions(+), 8 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index a148455d..1bbcacdd 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -506,6 +506,8 @@ def __init__( if len(xpoints) == 0: warnings.warn("No X-points found in TokamakEquilibrium input") + self.psi_bdry = psi_bdry_gfile + self.psi_bdry_gfile = psi_bdry_gfile else: self.psi_bdry = xpoints[0][2] # Psi on primary X-point self.x_point = Point2D(xpoints[0][0], xpoints[0][1]) @@ -1715,11 +1717,9 @@ def fpolprime(self, psi): """psi-derivative of fpol Note: Zero outside core.""" fprime = self.fprime_spl(psi * self.f_psi_sign) - try: + if self.psi_bdry is not None: psinorm = (psi - self.psi_axis) / (self.psi_bdry - self.psi_axis) fprime[psinorm > 1.0] = 0.0 - except: - pass return fprime @Equilibrium.handleMultiLocationArray @@ -1736,11 +1736,9 @@ def pprime(self, psi): if self.pprime_spl is None: return None pprime = self.pprime_spl(psi * self.f_psi_sign) - try: + if self.psi_bdry is not None: psinorm = (psi - self.psi_axis) / (self.psi_bdry - self.psi_axis) pprime[psinorm > 1.0] = 0.0 - except: - pass return pprime @property diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index 5d3695a2..e67b31f1 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -22,6 +22,7 @@ from datetime import date from numpy import zeros, pi +import numpy as np from ._fileutils import f2s, ChunkOutput, write_1d, write_2d, next_value @@ -48,6 +49,11 @@ def write(data, fh, label=None, shot=None, time=None): psi 2D array (nx,ny) of poloidal flux + ffprime 1D array of f(psi) * f'(psi). If not present + then this is calculated from fpol + pprime 1D array of p'(psi). If not present then + this is calculated from pres + fh - file handle label - Text label to put in the file @@ -134,11 +140,25 @@ def write(data, fh, label=None, shot=None, time=None): if "ffprime" in data: write_1d(data["ffprime"], co) else: - write_1d(workk, co) + psi1D = np.linspace(data["simagx"], data["sibdry"], nx) + from scipy import interpolate + + fprime_spl = interpolate.InterpolatedUnivariateSpline( + psi1D, data["fpol"] + ).derivative() + ffprime = data["fpol"] * fprime_spl(psi1D) + write_1d(ffprime, co) + if "pprime" in data: write_1d(data["pprime"], co) else: - write_1d(workk, co) + psi1D = np.linspace(data["simagx"], data["sibdry"], nx) + from scipy import interpolate + + pprime_spl = interpolate.InterpolatedUnivariateSpline( + psi1D, data["pres"] + ).derivative() + write_1d(pprime_spl(psi1D), co) write_2d(data["psi"], co) write_1d(data["qpsi"], co) From f769c3bcc924b45777d2332131a01a33dfd552c5 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 31 Jul 2024 17:11:13 -0700 Subject: [PATCH 5/9] integrated test: Add Jpar0 to expected difference Jpar0 isn't in the reference grid files so integrated tests fail. --- integrated_tests/utils.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/integrated_tests/utils.py b/integrated_tests/utils.py index 7fa2ce95..56315611 100644 --- a/integrated_tests/utils.py +++ b/integrated_tests/utils.py @@ -22,6 +22,9 @@ "penalty_mask", "closed_wall_R", "closed_wall_Z", + "Jpar0", + "Jpar0_xlow", + "Jpar0_ylow", ] From d4653e65aa8535ca5302820c791d86a6dab1c45f Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 31 Jul 2024 17:22:10 -0700 Subject: [PATCH 6/9] flake8: Remove unused workk variable --- hypnotoad/geqdsk/_geqdsk.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index e67b31f1..d62e1d2d 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -126,11 +126,7 @@ def write(data, fh, label=None, shot=None, time=None): f2s(data["zmagx"]) + f2s(0.0) + f2s(data["sibdry"]) + f2s(0.0) + f2s(0.0) + "\n" ) - # SCENE has actual ff' and p' data so can use that # fill arrays - # Lukas Kripner (16/10/2018): uncommenting this, since you left there - # check for data existence bellow. This seems to as safer variant. - workk = zeros([nx]) # Write arrays co = ChunkOutput(fh) From e0b3d22f4ea71fbd198893e8f3421de7816cc62e Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 2 Aug 2024 16:53:53 -0700 Subject: [PATCH 7/9] tokamak: Handle case where psi is a scalar --- hypnotoad/cases/tokamak.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 1bbcacdd..48827659 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -1719,7 +1719,10 @@ def fpolprime(self, psi): fprime = self.fprime_spl(psi * self.f_psi_sign) if self.psi_bdry is not None: psinorm = (psi - self.psi_axis) / (self.psi_bdry - self.psi_axis) - fprime[psinorm > 1.0] = 0.0 + if np.isscalar(psi) and psinorm > 1.0: + fprime = 0.0 + else: + fprime[psinorm > 1.0] = 0.0 return fprime @Equilibrium.handleMultiLocationArray @@ -1738,7 +1741,10 @@ def pprime(self, psi): pprime = self.pprime_spl(psi * self.f_psi_sign) if self.psi_bdry is not None: psinorm = (psi - self.psi_axis) / (self.psi_bdry - self.psi_axis) - pprime[psinorm > 1.0] = 0.0 + if np.isscalar(psi) and psinorm > 1.0: + pprime = 0.0 + else: + pprime[psinorm > 1.0] = 0.0 return pprime @property From fe21db945a8c3812e79d9b5a564c80ecab707d8f Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 2 Aug 2024 16:59:09 -0700 Subject: [PATCH 8/9] Check for poloidal field sign before fitting spline InterpolatedUnivariateSpline requires monotonically increasing X axis, so reverse sign of psi in negative Bp case. --- hypnotoad/geqdsk/_geqdsk.py | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index d62e1d2d..39409ed9 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -139,9 +139,15 @@ def write(data, fh, label=None, shot=None, time=None): psi1D = np.linspace(data["simagx"], data["sibdry"], nx) from scipy import interpolate - fprime_spl = interpolate.InterpolatedUnivariateSpline( - psi1D, data["fpol"] - ).derivative() + sign = -1.0 if psi1D[1] < psi1D[0] else 1.0 + # spline fitting requires increasing X axis, so reverse if needed + + fprime_spl = ( + sign + * interpolate.InterpolatedUnivariateSpline( + sign * psi1D, data["fpol"] + ).derivative() + ) ffprime = data["fpol"] * fprime_spl(psi1D) write_1d(ffprime, co) @@ -151,9 +157,13 @@ def write(data, fh, label=None, shot=None, time=None): psi1D = np.linspace(data["simagx"], data["sibdry"], nx) from scipy import interpolate - pprime_spl = interpolate.InterpolatedUnivariateSpline( - psi1D, data["pres"] - ).derivative() + sign = -1.0 if psi1D[1] < psi1D[0] else 1.0 + pprime_spl = ( + sign + * interpolate.InterpolatedUnivariateSpline( + sign * psi1D, data["pres"] + ).derivative() + ) write_1d(pprime_spl(psi1D), co) write_2d(data["psi"], co) From 521ec7ed572ed0e499a72d28feab64929252bcae Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 2 Aug 2024 17:08:47 -0700 Subject: [PATCH 9/9] Fix spline signs --- hypnotoad/geqdsk/_geqdsk.py | 22 ++++++++-------------- 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index 39409ed9..b8653010 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -142,13 +142,10 @@ def write(data, fh, label=None, shot=None, time=None): sign = -1.0 if psi1D[1] < psi1D[0] else 1.0 # spline fitting requires increasing X axis, so reverse if needed - fprime_spl = ( - sign - * interpolate.InterpolatedUnivariateSpline( - sign * psi1D, data["fpol"] - ).derivative() - ) - ffprime = data["fpol"] * fprime_spl(psi1D) + fprime_spl = interpolate.InterpolatedUnivariateSpline( + sign * psi1D, sign * data["fpol"] + ).derivative() + ffprime = data["fpol"] * fprime_spl(sign * psi1D) write_1d(ffprime, co) if "pprime" in data: @@ -158,13 +155,10 @@ def write(data, fh, label=None, shot=None, time=None): from scipy import interpolate sign = -1.0 if psi1D[1] < psi1D[0] else 1.0 - pprime_spl = ( - sign - * interpolate.InterpolatedUnivariateSpline( - sign * psi1D, data["pres"] - ).derivative() - ) - write_1d(pprime_spl(psi1D), co) + pprime_spl = interpolate.InterpolatedUnivariateSpline( + sign * psi1D, sign * data["pres"] + ).derivative() + write_1d(pprime_spl(sign * psi1D), co) write_2d(data["psi"], co) write_1d(data["qpsi"], co)