From 6f6e5463e9d11f3fb8a2328dfeaed137126b65a2 Mon Sep 17 00:00:00 2001 From: h_seto Date: Thu, 1 Aug 2024 11:32:55 -0700 Subject: [PATCH 01/60] fix fprime for the negative Bp case --- hypnotoad/cases/tokamak.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index ea5c98b9..cef489fb 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -434,7 +434,17 @@ 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: + # fprime is derivative with respect to psi rather than + # increasing radial label psi*f_psi_sign + fpol_f_psi_sign_spl = interpolate.InterpolatedUnivariateSpline( + psi1D * self.f_psi_sign, fpol1D * self.f_psi_sign, ext=3 + ) + self.fprime_spl = fpol_f_psi_sign_spl.derivative() else: self.f_spl = lambda psi: 0.0 self.fprime_spl = lambda psi: 0.0 From 731a330f5804a17d87fa04730cbc6712f670a897 Mon Sep 17 00:00:00 2001 From: h_seto Date: Thu, 1 Aug 2024 12:13:00 -0700 Subject: [PATCH 02/60] define increasing radial label xcoord --- hypnotoad/cases/tokamak.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index cef489fb..ba93df4c 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -421,28 +421,30 @@ def __init__( ) self.f_psi_sign = 1.0 - if len(fpol1D) > 0: - # Spline for interpolation of f = R*Bt - # Note: psi1D must be increasing - if psi1D[-1] < psi1D[0]: - self.f_psi_sign = -1.0 + # Note: radial label must be increasing + if psi1D[-1] < psi1D[0]: + self.f_psi_sign = -1.0 + xcoord = self.f_psi_sign* psi1D + if len(fpol1D) > 0: + # Spline for interpolation of f = R*Bt + self.f_spl = interpolate.InterpolatedUnivariateSpline( - psi1D * self.f_psi_sign, fpol1D, ext=3 + xcoord, fpol1D, ext=3 ) # ext=3 specifies that boundary values are used outside range # Spline representing the derivative of f if fprime is not None: self.fprime_spl = interpolate.InterpolatedUnivariateSpline( - psi1D * self.f_psi_sign, fprime, ext=3 + xcoord, fprime, ext=3 ) else: # fprime is derivative with respect to psi rather than # increasing radial label psi*f_psi_sign fpol_f_psi_sign_spl = interpolate.InterpolatedUnivariateSpline( - psi1D * self.f_psi_sign, fpol1D * self.f_psi_sign, ext=3 + xcoord, fpol1D * self.f_psi_sign, ext=3 ) self.fprime_spl = fpol_f_psi_sign_spl.derivative() else: @@ -452,7 +454,7 @@ def __init__( # Optional pressure profile if pressure is not None: self.p_spl = interpolate.InterpolatedUnivariateSpline( - psi1D * self.f_psi_sign, pressure, ext=3 + xcoord, pressure, ext=3 ) else: # If no pressure, then not output to grid file From 107c5dbb562a755432312f53695615891d5c154b Mon Sep 17 00:00:00 2001 From: h_seto Date: Fri, 2 Aug 2024 09:51:18 -0700 Subject: [PATCH 03/60] pressure extrapolation fix --- hypnotoad/cases/tokamak.py | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index ba93df4c..7cf69c86 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -403,19 +403,24 @@ def __init__( ): # if psi_outer is not beyond the last point of psi1D, no need to extend # Exclude first point since duplicates last point in core - psiSOL = np.linspace(psi1D[-1], psi_outer, 50)[1:] + psi_lcfs = psi1D[-1] + psiSOL = np.linspace(psi_lcfs, psi_outer, 50)[1:] psi1D = np.concatenate([psi1D, psiSOL]) - + # fpol constant in SOL fpol1D = np.concatenate([fpol1D, np.full(psiSOL.shape, fpol1D[-1])]) - 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)]) - + if pressure is not None: + # Use an exponential decay for the pressure, based on + # the value and gradient at the plasma edge + p_lcfs = pressure[-1] + # p_SOL = p_lcfs * exp( (psi - psi_lcfs) * dpdpsi / p_lcfs) + p_SOL = p_lcfs * np.exp((psiSOL-psi_lcfs)* dpdpsi / p_lcfs) + pressure = np.concatenate([pressure,p_SOL]) + + if pprime is not None: + pprime = np.concatenate([pprime, dpdpsi*p_SOL/p_lcfs]) + self.magneticFunctionsFromGrid( R1D, Z1D, psi2D, self.user_options.psi_interpolation_method ) @@ -442,7 +447,7 @@ def __init__( ) else: # fprime is derivative with respect to psi rather than - # increasing radial label psi*f_psi_sign + # increasing radial label xcoord fpol_f_psi_sign_spl = interpolate.InterpolatedUnivariateSpline( xcoord, fpol1D * self.f_psi_sign, ext=3 ) From 0c6a08491c051f1c6da8d2b0bcff9b0d191c0199 Mon Sep 17 00:00:00 2001 From: h_seto Date: Fri, 2 Aug 2024 10:03:13 -0700 Subject: [PATCH 04/60] pressure and pprime has been added in __init__ --- hypnotoad/cases/tokamak.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 7cf69c86..1b60d05e 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,7 +342,9 @@ 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. The wall is closed, so the last point connects to the first. @@ -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: From aecafcba51e9eba6234269fce0842fd1cda69e09 Mon Sep 17 00:00:00 2001 From: h_seto Date: Fri, 2 Aug 2024 10:29:50 -0700 Subject: [PATCH 05/60] slight modification on f_psi_sign --- hypnotoad/cases/tokamak.py | 29 ++++++++++++++++++++++++----- 1 file changed, 24 insertions(+), 5 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 1b60d05e..e2140585 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -416,7 +416,10 @@ 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 @@ -432,11 +435,13 @@ def __init__( R1D, Z1D, psi2D, self.user_options.psi_interpolation_method ) - self.f_psi_sign = 1.0 # Note: radial label must be increasing if psi1D[-1] < psi1D[0]: self.f_psi_sign = -1.0 + else: + self.f_psi_sign = 1.0 + xcoord = self.f_psi_sign* psi1D if len(fpol1D) > 0: @@ -453,8 +458,8 @@ def __init__( xcoord, fprime, ext=3 ) else: - # fprime is derivative with respect to psi rather than - # increasing radial label xcoord + # fprime is derivative of fpol with respect to psi + # rather than increasing radial label xcoord fpol_f_psi_sign_spl = interpolate.InterpolatedUnivariateSpline( xcoord, fpol1D * self.f_psi_sign, ext=3 ) @@ -468,10 +473,24 @@ def __init__( self.p_spl = interpolate.InterpolatedUnivariateSpline( xcoord, pressure, ext=3 ) + if pprime is not None: + self.pprime_spl = interpolate.InterpolatedUnivariateSpline( + xcoord, pprime, ext=3 + ) + else: + # pprime is derivative of pressure with respect to psi + # rather than increasing radial label xcoord + + pressure_f_psi_sign_spl = interpolate.InterpolatedUnivariateSpline( + xcoord, pressure * self.f_psi_sign, ext=3 + ) + self.pprime_spl = pressure_f_psi_sign_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") opoints, xpoints = critical.find_critical( From 16c5697f7612b228c043ce617de8b2ab11ce221b Mon Sep 17 00:00:00 2001 From: h_seto Date: Fri, 2 Aug 2024 10:47:38 -0700 Subject: [PATCH 06/60] fprime and pprime have been added in read_geqdsk --- hypnotoad/cases/tokamak.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index e2140585..912a96e0 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -1794,6 +1794,8 @@ def read_geqdsk( pressure = data["pres"] fpol = data["fpol"] + fprime = data["ffprime"] / fpol + pprime = data["pprime"] result = TokamakEquilibrium( R1D, @@ -1804,6 +1806,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, From 2816d33f1903732430c21b097a0c813f309d828b Mon Sep 17 00:00:00 2001 From: h_seto Date: Fri, 2 Aug 2024 15:01:58 -0700 Subject: [PATCH 07/60] extrapolation rule for fprime/pprime should be ext=1 for consistency with f_pol and pressure with ext=3 --- hypnotoad/cases/tokamak.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 912a96e0..8480e323 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -455,7 +455,7 @@ def __init__( # Spline representing the derivative of f if fprime is not None: self.fprime_spl = interpolate.InterpolatedUnivariateSpline( - xcoord, fprime, ext=3 + xcoord, fprime, ext=1 ) else: # fprime is derivative of fpol with respect to psi @@ -475,7 +475,7 @@ def __init__( ) if pprime is not None: self.pprime_spl = interpolate.InterpolatedUnivariateSpline( - xcoord, pprime, ext=3 + xcoord, pprime, ext=1 ) else: # pprime is derivative of pressure with respect to psi @@ -1664,9 +1664,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 From 6f38a0957e9ce882aab7fc6f5dc3e63378aa21dc Mon Sep 17 00:00:00 2001 From: h_seto Date: Fri, 2 Aug 2024 15:51:02 -0700 Subject: [PATCH 08/60] function pprime has been added --- hypnotoad/cases/tokamak.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 8480e323..7eee6f65 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -1735,6 +1735,13 @@ 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""" + if self.pprime_spl is None: + return None + return self.pprime_spl(psi * self.f_psi_sign) + @property def Bt_axis(self): """Calculate toroidal field on axis""" From 4c1fb53604028ca88da46c8c97ae2091478e12a4 Mon Sep 17 00:00:00 2001 From: h_seto Date: Fri, 2 Aug 2024 16:33:51 -0700 Subject: [PATCH 09/60] add pprime and ffprime in doc of read --- hypnotoad/geqdsk/_geqdsk.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index ed643216..c96d9863 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 44c2778483f26c85d4a8f04a842dbed0e88c1450 Mon Sep 17 00:00:00 2001 From: h_seto Date: Fri, 2 Aug 2024 17:18:09 -0700 Subject: [PATCH 10/60] fix interpolation for negative Bp case --- hypnotoad/geqdsk/_geqdsk.py | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index c96d9863..ed2dac92 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -22,6 +22,8 @@ from datetime import date from numpy import zeros, pi +import numpy as np +import scipy as sp from ._fileutils import f2s, ChunkOutput, write_1d, write_2d, next_value @@ -131,14 +133,36 @@ def write(data, fh, label=None, shot=None, time=None): write_1d(data["fpol"], co) write_1d(data["pres"], co) + if "ffprime" in data: write_1d(data["ffprime"], co) else: - write_1d(workk, co) + psi_axis = data["simagx"] + psi_bdry = data["sibdry"] + fpol = data["fpol"] + + sign_dpsi = np.sign(psi_bdry-psi_axis) + xcrd = np.linspace(psi_axis,psi_bdry,nx)*sign_dpsi + fprime_spl = sp.interpolate.InterpolatedUnivariateSpline( + xrd, fpol*sign_dpsi).derivative() + ffprime = fpol*fprime_spl(xrcd) + + write_1d(ffprime,co) + + if "pprime" in data: write_1d(data["pprime"], co) else: - write_1d(workk, co) + psi_axis = data["simagx"] + psi_bdry = data["sibdry"] + + sign_dpsi = np.sign(psi_bdry-psi_axis) + xcrd = np.linspace(psi_axis,psi_bdry,nx)*sign_dpsi + + pprime_spl = sp.interpolate.InterpolatedUnivariateSpline( + xrd, data["pres"]*sign_dpsi).derivative() + write_1d(pprime(xcrd),co) + write_2d(data["psi"], co) write_1d(data["qpsi"], co) From f3a364c7394f3912793207ebb36db371a5deef39 Mon Sep 17 00:00:00 2001 From: h_seto Date: Mon, 12 Aug 2024 14:57:40 -0700 Subject: [PATCH 11/60] cocos=+1 convention fix in find_critical --- hypnotoad/utils/critical.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/hypnotoad/utils/critical.py b/hypnotoad/utils/critical.py index a2e0a1c9..4a86541f 100644 --- a/hypnotoad/utils/critical.py +++ b/hypnotoad/utils/critical.py @@ -110,9 +110,8 @@ def find_critical(R, Z, psi, atol, maxits, discard_xpoints=False): count = 0 while True: - Br = -f(R1, Z1, dy=1, grid=False) / R1 - Bz = f(R1, Z1, dx=1, grid=False) / R1 - + Br = f(R1, Z1, dy=1, grid=False) / R1 + Bz = -f(R1, Z1, dx=1, grid=False) / R1 if Br**2 + Bz**2 < atol: # Found a minimum. Classify as either # O-point or X-point @@ -181,10 +180,10 @@ def find_critical(R, Z, psi, atol, maxits, discard_xpoints=False): # J = ( dBr/dR, dBr/dZ ) # ( dBz/dR, dBz/dZ ) - J[0, 0] = -Br / R1 - f(R1, Z1, dy=1, dx=1)[0][0] / R1 - J[0, 1] = -f(R1, Z1, dy=2)[0][0] / R1 - J[1, 0] = -Bz / R1 + f(R1, Z1, dx=2) / R1 - J[1, 1] = f(R1, Z1, dx=1, dy=1)[0][0] / R1 + J[0, 0] = -Br / R1 + f(R1, Z1, dy=1, dx=1)[0][0] / R1 + J[0, 1] = f(R1, Z1, dy=2)[0][0] / R1 + J[1, 0] = -Bz / R1 - f(R1, Z1, dx=2) / R1 + J[1, 1] = -f(R1, Z1, dx=1, dy=1)[0][0] / R1 d = dot(inv(J), [Br, Bz]) From 9c86ab1fff424ae62566d71209c369e8d7f5023f Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 13 Aug 2024 10:05:24 -0700 Subject: [PATCH 12/60] fix xrd in _geqdsk.py --- hypnotoad/geqdsk/_geqdsk.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index ed2dac92..daa37217 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -144,8 +144,8 @@ def write(data, fh, label=None, shot=None, time=None): sign_dpsi = np.sign(psi_bdry-psi_axis) xcrd = np.linspace(psi_axis,psi_bdry,nx)*sign_dpsi fprime_spl = sp.interpolate.InterpolatedUnivariateSpline( - xrd, fpol*sign_dpsi).derivative() - ffprime = fpol*fprime_spl(xrcd) + xcrd, fpol*sign_dpsi).derivative() + ffprime = fpol*fprime_spl(xcrd) write_1d(ffprime,co) From 58489ae40c66eb19eda9f059914919861fe4eb54 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 13 Aug 2024 10:09:13 -0700 Subject: [PATCH 13/60] fix xrd in _geqdsk.py --- hypnotoad/geqdsk/_geqdsk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index daa37217..648d92f2 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -160,7 +160,7 @@ def write(data, fh, label=None, shot=None, time=None): xcrd = np.linspace(psi_axis,psi_bdry,nx)*sign_dpsi pprime_spl = sp.interpolate.InterpolatedUnivariateSpline( - xrd, data["pres"]*sign_dpsi).derivative() + xcrd, data["pres"]*sign_dpsi).derivative() write_1d(pprime(xcrd),co) From 59dabcab8b92e85e8d61b058b823843a8ca0210e Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 13 Aug 2024 10:22:54 -0700 Subject: [PATCH 14/60] fix pprime in _geqdsk.py --- hypnotoad/geqdsk/_geqdsk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index 648d92f2..c6d1e1e5 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -161,7 +161,7 @@ def write(data, fh, label=None, shot=None, time=None): pprime_spl = sp.interpolate.InterpolatedUnivariateSpline( xcrd, data["pres"]*sign_dpsi).derivative() - write_1d(pprime(xcrd),co) + write_1d(pprime_spl(xcrd),co) write_2d(data["psi"], co) From 4b26b8033b84c3938334fa034b5b23af1724a1f4 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 13 Aug 2024 10:45:35 -0700 Subject: [PATCH 15/60] flake8 fix in tokamak.py --- hypnotoad/cases/tokamak.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 7eee6f65..e744fce7 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -344,7 +344,7 @@ 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. The wall is closed, so the last point connects to the first. @@ -368,7 +368,7 @@ 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: @@ -429,24 +429,23 @@ def __init__( pressure = np.concatenate([pressure,p_SOL]) if pprime is not None: - pprime = np.concatenate([pprime, dpdpsi*p_SOL/p_lcfs]) + pprime = np.concatenate([pprime, dpdpsi * p_SOL / p_lcfs]) self.magneticFunctionsFromGrid( R1D, Z1D, psi2D, self.user_options.psi_interpolation_method ) - # Note: radial label must be increasing if psi1D[-1] < psi1D[0]: self.f_psi_sign = -1.0 else: self.f_psi_sign = 1.0 - - xcoord = self.f_psi_sign* psi1D + + xcoord = self.f_psi_sign * psi1D if len(fpol1D) > 0: + # Spline for interpolation of f = R*Bt - self.f_spl = interpolate.InterpolatedUnivariateSpline( xcoord, fpol1D, ext=3 ) @@ -485,12 +484,12 @@ def __init__( xcoord, pressure * self.f_psi_sign, ext=3 ) self.pprime_spl = pressure_f_psi_sign_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") opoints, xpoints = critical.find_critical( @@ -1741,7 +1740,7 @@ def pprime(self, psi): if self.pprime_spl is None: return None return self.pprime_spl(psi * self.f_psi_sign) - + @property def Bt_axis(self): """Calculate toroidal field on axis""" From 93bdecf5a95d8ee389d56b9f6f22d2470fa12e85 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 13 Aug 2024 10:47:14 -0700 Subject: [PATCH 16/60] flake8 fix in _geqdsk.py --- hypnotoad/geqdsk/_geqdsk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index c6d1e1e5..0f354caf 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -126,7 +126,7 @@ def write(data, fh, label=None, shot=None, time=None): # 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]) + #workk = zeros([nx]) # Write arrays co = ChunkOutput(fh) From 0bf561ff69659925a09c12cf477b408533fa8c74 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 13 Aug 2024 10:53:04 -0700 Subject: [PATCH 17/60] flake8 fix in _geqdsk.py --- hypnotoad/geqdsk/_geqdsk.py | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index 0f354caf..6cb4c338 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -140,16 +140,14 @@ def write(data, fh, label=None, shot=None, time=None): psi_axis = data["simagx"] psi_bdry = data["sibdry"] fpol = data["fpol"] - + sign_dpsi = np.sign(psi_bdry-psi_axis) - xcrd = np.linspace(psi_axis,psi_bdry,nx)*sign_dpsi + xcrd = np.linspace(psi_axis, psi_bdry, nx) * sign_dpsi fprime_spl = sp.interpolate.InterpolatedUnivariateSpline( - xcrd, fpol*sign_dpsi).derivative() - ffprime = fpol*fprime_spl(xcrd) - + xcrd, fpol * sign_dpsi).derivative() + ffprime = fpol * fprime_spl(xcrd) write_1d(ffprime,co) - - + if "pprime" in data: write_1d(data["pprime"], co) else: @@ -157,12 +155,11 @@ def write(data, fh, label=None, shot=None, time=None): psi_bdry = data["sibdry"] sign_dpsi = np.sign(psi_bdry-psi_axis) - xcrd = np.linspace(psi_axis,psi_bdry,nx)*sign_dpsi + xcrd = np.linspace(psi_axis, psi_bdry, nx) * sign_dpsi pprime_spl = sp.interpolate.InterpolatedUnivariateSpline( - xcrd, data["pres"]*sign_dpsi).derivative() - write_1d(pprime_spl(xcrd),co) - + xcrd, data["pres"] * sign_dpsi).derivative() + write_1d(pprime_spl(xcrd), co) write_2d(data["psi"], co) write_1d(data["qpsi"], co) From efbd96188b84a72b85ee73494eb25af44611a5f2 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 13 Aug 2024 11:00:06 -0700 Subject: [PATCH 18/60] flake8 fix in _geqdsk.py --- hypnotoad/geqdsk/_geqdsk.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index 6cb4c338..20dce969 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -126,7 +126,7 @@ def write(data, fh, label=None, shot=None, time=None): # 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]) + # workk = zeros([nx]) # Write arrays co = ChunkOutput(fh) @@ -141,7 +141,7 @@ def write(data, fh, label=None, shot=None, time=None): psi_bdry = data["sibdry"] fpol = data["fpol"] - sign_dpsi = np.sign(psi_bdry-psi_axis) + sign_dpsi = np.sign(psi_bdry - psi_axis) xcrd = np.linspace(psi_axis, psi_bdry, nx) * sign_dpsi fprime_spl = sp.interpolate.InterpolatedUnivariateSpline( xcrd, fpol * sign_dpsi).derivative() @@ -154,8 +154,8 @@ def write(data, fh, label=None, shot=None, time=None): psi_axis = data["simagx"] psi_bdry = data["sibdry"] - sign_dpsi = np.sign(psi_bdry-psi_axis) - xcrd = np.linspace(psi_axis, psi_bdry, nx) * sign_dpsi + sign_dpsi = np.sign(psi_bdry - psi_axis) + xcrd = np.linspace(psi_axis, psi_bdry, nx) * sign_dpsi pprime_spl = sp.interpolate.InterpolatedUnivariateSpline( xcrd, data["pres"] * sign_dpsi).derivative() @@ -217,7 +217,7 @@ 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) + pprime 1D array of p'(psi) qpsi 1D array of q(psi) psi 2D array (nx,ny) of poloidal flux From be66d650dd407c05e9bedc688f23de0cc60dd9d9 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 13 Aug 2024 11:02:31 -0700 Subject: [PATCH 19/60] flake8 fix in tokamak.py --- hypnotoad/cases/tokamak.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index e744fce7..75aff168 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -413,24 +413,24 @@ def __init__( psi_lcfs = psi1D[-1] psiSOL = np.linspace(psi_lcfs, psi_outer, 50)[1:] psi1D = np.concatenate([psi1D, psiSOL]) - + # 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 p_lcfs = pressure[-1] # p_SOL = p_lcfs * exp( (psi - psi_lcfs) * dpdpsi / p_lcfs) - p_SOL = p_lcfs * np.exp((psiSOL-psi_lcfs)* dpdpsi / p_lcfs) - pressure = np.concatenate([pressure,p_SOL]) + p_SOL = p_lcfs * np.exp((psiSOL-psi_lcfs) * dpdpsi / p_lcfs) + pressure = np.concatenate([pressure, p_SOL]) if pprime is not None: pprime = np.concatenate([pprime, dpdpsi * p_SOL / p_lcfs]) - + self.magneticFunctionsFromGrid( R1D, Z1D, psi2D, self.user_options.psi_interpolation_method ) From 0d3ae2626ea39973fca42a84554dbb5f3fd9154a Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 13 Aug 2024 11:05:03 -0700 Subject: [PATCH 20/60] flake8 fix in tokamak.py --- hypnotoad/cases/tokamak.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 75aff168..631717b5 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -416,7 +416,7 @@ 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)]) @@ -425,7 +425,7 @@ def __init__( # the value and gradient at the plasma edge p_lcfs = pressure[-1] # p_SOL = p_lcfs * exp( (psi - psi_lcfs) * dpdpsi / p_lcfs) - p_SOL = p_lcfs * np.exp((psiSOL-psi_lcfs) * dpdpsi / p_lcfs) + p_SOL = p_lcfs * np.exp((psiSOL - psi_lcfs) * dpdpsi / p_lcfs) pressure = np.concatenate([pressure, p_SOL]) if pprime is not None: From e806fb4f7363a933e222e26838a5d244b09a4fb1 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 13 Aug 2024 11:15:14 -0700 Subject: [PATCH 21/60] black fix in tokamak.py --- hypnotoad/cases/tokamak.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 631717b5..1c30a776 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -446,9 +446,7 @@ def __init__( if len(fpol1D) > 0: # Spline for interpolation of f = R*Bt - self.f_spl = interpolate.InterpolatedUnivariateSpline( - xcoord, fpol1D, ext=3 - ) + self.f_spl = interpolate.InterpolatedUnivariateSpline(xcoord, fpol1D, ext=3) # ext=3 specifies that boundary values are used outside range # Spline representing the derivative of f From 8f5b7c6d3b1b9e7d7f47e69380f610d913e149a8 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 13 Aug 2024 11:16:42 -0700 Subject: [PATCH 22/60] black fixs in _geqdsk.py --- hypnotoad/geqdsk/_geqdsk.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/hypnotoad/geqdsk/_geqdsk.py b/hypnotoad/geqdsk/_geqdsk.py index 20dce969..753483df 100644 --- a/hypnotoad/geqdsk/_geqdsk.py +++ b/hypnotoad/geqdsk/_geqdsk.py @@ -144,9 +144,10 @@ def write(data, fh, label=None, shot=None, time=None): sign_dpsi = np.sign(psi_bdry - psi_axis) xcrd = np.linspace(psi_axis, psi_bdry, nx) * sign_dpsi fprime_spl = sp.interpolate.InterpolatedUnivariateSpline( - xcrd, fpol * sign_dpsi).derivative() + xcrd, fpol * sign_dpsi + ).derivative() ffprime = fpol * fprime_spl(xcrd) - write_1d(ffprime,co) + write_1d(ffprime, co) if "pprime" in data: write_1d(data["pprime"], co) @@ -158,7 +159,8 @@ def write(data, fh, label=None, shot=None, time=None): xcrd = np.linspace(psi_axis, psi_bdry, nx) * sign_dpsi pprime_spl = sp.interpolate.InterpolatedUnivariateSpline( - xcrd, data["pres"] * sign_dpsi).derivative() + xcrd, data["pres"] * sign_dpsi + ).derivative() write_1d(pprime_spl(xcrd), co) write_2d(data["psi"], co) From bf27fa46c0eb3cdc57b4eff279d51ad0d984cf21 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 13 Aug 2024 17:14:59 -0700 Subject: [PATCH 23/60] cocos class has been added --- hypnotoad/cocos/cocos.py | 109 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 hypnotoad/cocos/cocos.py diff --git a/hypnotoad/cocos/cocos.py b/hypnotoad/cocos/cocos.py new file mode 100644 index 00000000..0ee057f6 --- /dev/null +++ b/hypnotoad/cocos/cocos.py @@ -0,0 +1,109 @@ +# +# COCOS convention class +# +# Coordinate conventions based on +# COCOS paper: +# O. Sauter and S. Yu. Medvevdev, "Tokamak Coordinate Conventions: COCOS" +# Comput. Phys. Commun. 184 (2013) 293 +# +# and cocos_module.f90 in CHEASE code: https://gitlab.epfl.ch/spc/chease +# +# written by Haruki SETO seto.haruki@qst.go.jp +# + +import numpy as np + + +class Cocos: + """ + cocos_index: +1/+11, +2/+12, +3/+13, +4/+14, +5/+15, +6/+16, +7/+17, +8/+18 + exp_Bp: 0 for poloidal flux devided by 2pi, and +1 for full poloidal flux + sigma_Bp: sign of dpsi/drho (+1 or -1) for Ip > 0 and Bt > 0 + sigma_cylind: +1 for (R, phi, Z), and -1 for (R, Z, phi) + sigma_poloid: +1 for (rho, theta, phi), and -1 for (rho, phi, theta) + sigma_qsafe: sign of qsafe (+1/-1) for Ip > 0 and Bt > 0 + sigma_pprime: sign of pprime (+1/-1) for Ip > 0 and Bt > 0 + theta_clockwize: True for clockwize and False for cnt-clockwize + """ + + def __init__(self, cocos_index): + + cocos_index_sign = np.sign(cocos_index) + + if cocos_index_sign < 0: + raise ValueError("negative inputs are not supported yet:", cocos_index) + + cocos_index_2nd_digit = cocos_index // 10 + + if cocos_index_2nd_digit == 0: + self.exp_Bp = 0.0 + + elif cocos_index_2nd_digit == 1: + self.exp_Bp = 1.0 + + else: + raise ValueError("Invalid COCOS input:", cocos_index) + + cocos_index_1st_digit = cocos_index % 10 + + if cocos_index_1st_digit == 1: + self.sigma_Bp = +1.0 + self.sigma_cylind = +1.0 + self.sigma_poloid = +1.0 + self.sigma_qsafe = +1.0 + self.sigma_pprime = -1.0 + + elif cocos_index_1st_digit == 2: + self.sigma_Bp = +1.0 + self.sigma_cylind = -1.0 + self.sigma_poloid = +1.0 + self.sigma_qsafe = +1.0 + self.sigma_pprime = -1.0 + + elif cocos_index_1st_digit == 3: + self.sigma_Bp = -1.0 + self.sigma_cylind = +1.0 + self.sigma_poloid = -1.0 + self.sigma_qsafe = -1.0 + self.sigma_pprime = +1.0 + + elif cocos_index_1st_digit == 4: + self.sigma_Bp = -1.0 + self.sigma_cylind = -1.0 + self.sigma_poloid = -1.0 + self.sigma_qsafe = -1.0 + self.sigma_pprime = +1.0 + + elif cocos_index_1st_digit == 5: + self.sigma_Bp = +1.0 + self.sigma_cylind = +1.0 + self.sigma_poloid = -1.0 + self.sigma_qsafe = -1.0 + self.sigma_pprime = -1.0 + + elif cocos_index_1st_digit == 6: + self.sigma_Bp = +1.0 + self.sigma_cylind = -1.0 + self.sigma_poloid = -1.0 + self.sigma_qsafe = -1.0 + self.sigma_pprime = -1.0 + + elif cocos_index_1st_digit == 7: + self.sigma_Bp = -1.0 + self.sigma_cylind = +1.0 + self.sigma_poloid = +1.0 + self.sigma_qsafe = +1.0 + self.sigma_pprime = +1.0 + + elif cocos_index_1st_digit == 8: + self.sigma_Bp = -1.0 + self.sigma_cylind = -1.0 + self.sigma_poloid = +1.0 + self.sigma_qsafe = +1.0 + self.sigma_pprime = +1.0 + + else: + raise ValueError("Invalid COCOS input:", cocos_index) + + self.cocos_index = cocos_index + self.theta_clockwize = self.sigma_cylind * self.sigma_poloid > 0 From ac8593d3a93814970864c638c54f9747a998e19c Mon Sep 17 00:00:00 2001 From: HarukiST Date: Thu, 15 Aug 2024 14:33:18 -0700 Subject: [PATCH 24/60] flake8 and black fixes --- hypnotoad/cocos/cocos.py | 218 ++++++++++++++++++++++++++++++--------- 1 file changed, 168 insertions(+), 50 deletions(-) diff --git a/hypnotoad/cocos/cocos.py b/hypnotoad/cocos/cocos.py index 0ee057f6..3a1b1fa7 100644 --- a/hypnotoad/cocos/cocos.py +++ b/hypnotoad/cocos/cocos.py @@ -17,20 +17,21 @@ class Cocos: """ cocos_index: +1/+11, +2/+12, +3/+13, +4/+14, +5/+15, +6/+16, +7/+17, +8/+18 - exp_Bp: 0 for poloidal flux devided by 2pi, and +1 for full poloidal flux - sigma_Bp: sign of dpsi/drho (+1 or -1) for Ip > 0 and Bt > 0 - sigma_cylind: +1 for (R, phi, Z), and -1 for (R, Z, phi) - sigma_poloid: +1 for (rho, theta, phi), and -1 for (rho, phi, theta) - sigma_qsafe: sign of qsafe (+1/-1) for Ip > 0 and Bt > 0 - sigma_pprime: sign of pprime (+1/-1) for Ip > 0 and Bt > 0 - theta_clockwize: True for clockwize and False for cnt-clockwize + + exp_Bp: 0 for poloidal flux devided by 2pi, and +1 for full poloidal flux + sign_Bp: sign of dpsi/drho (+1 or -1) for Ip > 0 and Bt > 0 + sign_cylind: +1 for (R, phi, Z), and -1 for (R, Z, phi) + sign_poloid: +1 for (rho, theta, phi), and -1 for (rho, phi, theta) + sign_qsafe: sign of qsafe (+1/-1) for Ip > 0 and Bt > 0 + sign_pprime: sign of pprime (+1/-1) for Ip > 0 and Bt > 0 + theta_clockwize: True for clockwize and False for cnt-clockwize """ def __init__(self, cocos_index): - cocos_index_sign = np.sign(cocos_index) + sign_cocos_index = np.sign(cocos_index) - if cocos_index_sign < 0: + if sign_cocos_index < 0: raise ValueError("negative inputs are not supported yet:", cocos_index) cocos_index_2nd_digit = cocos_index // 10 @@ -47,63 +48,180 @@ def __init__(self, cocos_index): cocos_index_1st_digit = cocos_index % 10 if cocos_index_1st_digit == 1: - self.sigma_Bp = +1.0 - self.sigma_cylind = +1.0 - self.sigma_poloid = +1.0 - self.sigma_qsafe = +1.0 - self.sigma_pprime = -1.0 + self.sign_Bp = +1.0 + self.sign_cylind = +1.0 + self.sign_poloid = +1.0 + self.sign_qsafe = +1.0 + self.sign_pprime = -1.0 elif cocos_index_1st_digit == 2: - self.sigma_Bp = +1.0 - self.sigma_cylind = -1.0 - self.sigma_poloid = +1.0 - self.sigma_qsafe = +1.0 - self.sigma_pprime = -1.0 + self.sign_Bp = +1.0 + self.sign_cylind = -1.0 + self.sign_poloid = +1.0 + self.sign_qsafe = +1.0 + self.sign_pprime = -1.0 elif cocos_index_1st_digit == 3: - self.sigma_Bp = -1.0 - self.sigma_cylind = +1.0 - self.sigma_poloid = -1.0 - self.sigma_qsafe = -1.0 - self.sigma_pprime = +1.0 + self.sign_Bp = -1.0 + self.sign_cylind = +1.0 + self.sign_poloid = -1.0 + self.sign_qsafe = -1.0 + self.sign_pprime = +1.0 elif cocos_index_1st_digit == 4: - self.sigma_Bp = -1.0 - self.sigma_cylind = -1.0 - self.sigma_poloid = -1.0 - self.sigma_qsafe = -1.0 - self.sigma_pprime = +1.0 + self.sign_Bp = -1.0 + self.sign_cylind = -1.0 + self.sign_poloid = -1.0 + self.sign_qsafe = -1.0 + self.sign_pprime = +1.0 elif cocos_index_1st_digit == 5: - self.sigma_Bp = +1.0 - self.sigma_cylind = +1.0 - self.sigma_poloid = -1.0 - self.sigma_qsafe = -1.0 - self.sigma_pprime = -1.0 + self.sign_Bp = +1.0 + self.sign_cylind = +1.0 + self.sign_poloid = -1.0 + self.sign_qsafe = -1.0 + self.sign_pprime = -1.0 elif cocos_index_1st_digit == 6: - self.sigma_Bp = +1.0 - self.sigma_cylind = -1.0 - self.sigma_poloid = -1.0 - self.sigma_qsafe = -1.0 - self.sigma_pprime = -1.0 + self.sign_Bp = +1.0 + self.sign_cylind = -1.0 + self.sign_poloid = -1.0 + self.sign_qsafe = -1.0 + self.sign_pprime = -1.0 elif cocos_index_1st_digit == 7: - self.sigma_Bp = -1.0 - self.sigma_cylind = +1.0 - self.sigma_poloid = +1.0 - self.sigma_qsafe = +1.0 - self.sigma_pprime = +1.0 + self.sign_Bp = -1.0 + self.sign_cylind = +1.0 + self.sign_poloid = +1.0 + self.sign_qsafe = +1.0 + self.sign_pprime = +1.0 elif cocos_index_1st_digit == 8: - self.sigma_Bp = -1.0 - self.sigma_cylind = -1.0 - self.sigma_poloid = +1.0 - self.sigma_qsafe = +1.0 - self.sigma_pprime = +1.0 + self.sign_Bp = -1.0 + self.sign_cylind = -1.0 + self.sign_poloid = +1.0 + self.sign_qsafe = +1.0 + self.sign_pprime = +1.0 else: raise ValueError("Invalid COCOS input:", cocos_index) + self.theta_clockwize = self.sign_cylind * self.sign_poloid > 0 self.cocos_index = cocos_index - self.theta_clockwize = self.sigma_cylind * self.sigma_poloid > 0 + + +class CocosConversion(Cocos): + + def __init__( + self, cocos_index, geqdsk_dict, positive_definite_qsafe=False, show_signs=False + ): + + super().__init__(cocos_index) + + self.check_consistency( + cocos_index, geqdsk_dict, positive_definite_qsafe, show_signs + ) + + def check_consistency( + self, cocos_index, geqdsk_dict, positive_definite_qsafe, show_signs + ): + """ + check consistency between input GEQDSK and COCOS convension based on + Section 5 in COCOS paper. + + geqdsk_dict - dictionary of GEQDSK read by _geqdsk.py + + positive_definite_qsafe - always use |q| rather than q (e.g. JET) + """ + + xpos = geqdsk_dict["nx"] // 2 # radial index for checking sign of 1D quantities + + sign_Ip_geqdsk = np.sign(geqdsk_dict["cpasma"]) + sign_B0_geqdsk = np.sign(geqdsk_dict["fpol"][xpos]) + sign_dpsi_geqdsk = np.sign(geqdsk_dict["sibdry"] - geqdsk_dict["simagx"]) + sign_pprime_geqdsk = np.sign(geqdsk_dict["pprime"][xpos]) + sign_qsafe_geqdsk = np.sign(geqdsk_dict["qpsi"][xpos]) + + if show_signs: + print("sign_Ip_geqdsk %d" % sign_Ip_geqdsk) + print("sign_B0_geqdsk %d" % sign_B0_geqdsk) + print("sign_dpsi_geqdsk %d" % sign_dpsi_geqdsk) + print("sign_pprime_geqdsk %d" % sign_pprime_geqdsk) + print("sign_qsafe_geqdsk %d" % sign_qsafe_geqdsk) + + if sign_dpsi_geqdsk == np.sign(sign_Ip_geqdsk * self.sign_Bp): + self.sign_dpsi = sign_dpsi_geqdsk + + else: + raise ValueError("geqdsk is not consistent with COCOS (dpsi)") + + if sign_pprime_geqdsk == np.sign(-1 * sign_Ip_geqdsk * self.sign_Bp): + self.sign_pprime = sign_pprime_geqdsk + + else: + raise ValueError("geqdsk is not consistent with COCOS (pprime)") + + sign_qsafe_cocos = np.sign(sign_Ip_geqdsk * sign_B0_geqdsk * self.sign_poloid) + + if sign_qsafe_geqdsk == sign_qsafe_cocos: + self.sign_qsafe = sign_qsafe_cocos + + elif positive_definite_qsafe and sign_qsafe_cocos < 0: + print("geqdsk has positive definite qsafe but qsafe should be negative") + self.sign_qsafe = sign_qsafe_cocos + geqdsk_dict["qpsi"] *= -1 + + else: + raise ValueError("geqdsk is not consistent with COCOS input (qsafe)") + + self.sign_B0 = sign_B0_geqdsk + self.sign_Ip = sign_Ip_geqdsk + self.geqdsk_dict = geqdsk_dict + + def convert_geqdsk_to_cocos_out( + self, cocos_index_out=1, sign_B0_out=None, sign_Ip_out=None + ): + """ + convert geqdsk file to cocos_index_out convension. + + Based on Appendix C in COCOS paper but, l_d = 1[m], l_b = 1[T], exp_mu = 1 + are assumed for both input and output. + """ + + twopi = np.pi * 2.0 + + cocos_out = Cocos(cocos_index_out) + + sign_cylind_eff = np.sign(cocos_out.sign_cylind * self.sign_cylind) + sign_poloid_eff = np.sign(cocos_out.sign_poloid * self.sign_poloid) + sign_Bp_eff = np.sign(cocos_out.sign_Bp * self.sign_Bp) + exp_Bp_eff = cocos_out.exp_Bp - self.exp_Bp + + if sign_B0_out is not None: + sign_B0_eff = np.sign(sign_B0_out * self.sign_B0) + else: + sign_B0_eff = sign_cylind_eff + + if sign_Ip_out is not None: + sign_Ip_eff = np.sign(sign_Ip_out * self.sign_Ip) + else: + sign_Ip_eff = sign_cylind_eff + + # convert self.geqdsk_dict to COCOS_out convension rule + self.geqdsk_dict["cpasma"] *= sign_Ip_eff + self.geqdsk_dict["bcentr"] *= sign_B0_eff + + factor_psi = sign_Ip_eff * sign_Bp_eff * np.power(twopi, exp_Bp_eff) + self.geqdsk_dict["simagx"] *= factor_psi + self.geqdsk_dict["sibdry"] *= factor_psi + self.geqdsk_dict["psi"] *= factor_psi + + self.geqdsk_dict["fpol"] *= sign_B0_eff + + factor_dpsi = sign_Ip_eff * sign_Bp_eff / np.power(twopi, exp_Bp_eff) + self.geqdsk_dict["pprime"] *= factor_dpsi + self.geqdsk_dict["ffprime"] *= factor_dpsi + + self.geqdsk_dict["qpsi"] *= sign_Ip_eff * sign_B0_eff * sign_poloid_eff + + return self.geqdsk_dict From c6bf6dd96513dedeca0118991df44a1db4038324 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 18 Sep 2024 10:35:02 +0900 Subject: [PATCH 25/60] add signed Bp function in module equilibrium --- hypnotoad/core/equilibrium.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/hypnotoad/core/equilibrium.py b/hypnotoad/core/equilibrium.py index d4b75f04..3fad1c59 100644 --- a/hypnotoad/core/equilibrium.py +++ b/hypnotoad/core/equilibrium.py @@ -4040,6 +4040,13 @@ def Bzeta(self, R, Z): """ return self.fpol(self.psi(R, Z)) / R + def Bpol(self, R, Z): + """ + signed poloidal magnetic field as a function of R and Z + added by H. Seto + """ + return self.f_psi_sign * numpy.sqrt(self.Bp_R(R, Z) **2 + self.Bp_Z(R, Z) ** 2) + def B2(self, R, Z): """ B^2 as a function of R and Z From 8458d981abd21430839772fa253b45c2ae8779bf Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 18 Sep 2024 10:35:37 +0900 Subject: [PATCH 26/60] use signed Bp in zShift calculation --- hypnotoad/core/mesh.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 8f1c1c38..f9dc0d1b 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -1,4 +1,4 @@ -# Copyright 2019 J.T. Omotani +# CopyrigAht 2019 J.T. Omotani # # Contact John Omotani john.omotani@ukaea.uk # @@ -1558,10 +1558,12 @@ def integrand_func(R, Z): ) / R ) - Bp = numpy.sqrt( - self.meshParent.equilibrium.Bp_R(R, Z) ** 2 - + self.meshParent.equilibrium.Bp_Z(R, Z) ** 2 - ) + # Bp = numpy.sqrt( + # self.meshParent.equilibrium.Bp_R(R, Z) ** 2 + # + self.meshParent.equilibrium.Bp_Z(R, Z) ** 2 + # ) + Bp = self.meshParent.equilibrium.Bpol(R, Z) + return Bt / (R * Bp) integrand = integrand_func( From de512ca3c39e23df7dc07efa67e534569c2cb0f1 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 18 Sep 2024 10:46:34 +0900 Subject: [PATCH 27/60] black reformatting forsigned Bp --- hypnotoad/core/equilibrium.py | 4 ++-- hypnotoad/core/mesh.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/hypnotoad/core/equilibrium.py b/hypnotoad/core/equilibrium.py index 3fad1c59..d4ed8153 100644 --- a/hypnotoad/core/equilibrium.py +++ b/hypnotoad/core/equilibrium.py @@ -4045,8 +4045,8 @@ def Bpol(self, R, Z): signed poloidal magnetic field as a function of R and Z added by H. Seto """ - return self.f_psi_sign * numpy.sqrt(self.Bp_R(R, Z) **2 + self.Bp_Z(R, Z) ** 2) - + return self.f_psi_sign * numpy.sqrt(self.Bp_R(R, Z) ** 2 + self.Bp_Z(R, Z) ** 2) + def B2(self, R, Z): """ B^2 as a function of R and Z diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index f9dc0d1b..79687a30 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -1563,7 +1563,7 @@ def integrand_func(R, Z): # + self.meshParent.equilibrium.Bp_Z(R, Z) ** 2 # ) Bp = self.meshParent.equilibrium.Bpol(R, Z) - + return Bt / (R * Bp) integrand = integrand_func( From 7d849d312c90837c224cf6adfdd9970934ec9645 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 18 Sep 2024 14:10:11 +0900 Subject: [PATCH 28/60] f_psi_sign has been added for signed Bp --- hypnotoad/cases/circular.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/hypnotoad/cases/circular.py b/hypnotoad/cases/circular.py index 29c8e24f..48ad47f3 100644 --- a/hypnotoad/cases/circular.py +++ b/hypnotoad/cases/circular.py @@ -124,6 +124,8 @@ def __init__(self, settings=None, nonorthogonal_settings=None): # Core-only geometry: add connection so domain is periodic in y self.makeConnection("circular", 0, "circular", 0) + self.f_psi_sign = np.sign(self.psi_r(1.0) - self.psi_r(0.0)) + def makeRegion(self): """ Create the EquilibriumRegion to build the grid around From f7e788249abc9f069ff5c5af4981624eede08224 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 18 Sep 2024 14:50:26 +0900 Subject: [PATCH 29/60] f_psi_sign has been added in torpex.py for signed Bp --- hypnotoad/cases/torpex.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/hypnotoad/cases/torpex.py b/hypnotoad/cases/torpex.py index d03f8a7b..14370955 100644 --- a/hypnotoad/cases/torpex.py +++ b/hypnotoad/cases/torpex.py @@ -206,6 +206,9 @@ def __init__(self, equilibOptions, meshOptions): "direction of plasma current should be clockwise to be " "consistent with sign of grad(psi)" ) + + self.f_psi_sign = numpy.sign(psi_bndry - psi_axis) + # index of a point close to the magnetic axis i_axis = numpy.searchsorted(R, R_axis) j_axis = numpy.searchsorted(Z, Z_axis) From 7653245d779f6a64b50d9f2889f2472b7bc0c515 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 18 Sep 2024 15:07:59 +0900 Subject: [PATCH 30/60] f_psi_sign has been added for matfile in torpex.py for signed Bp --- hypnotoad/cases/torpex.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/hypnotoad/cases/torpex.py b/hypnotoad/cases/torpex.py index 14370955..9efaea2c 100644 --- a/hypnotoad/cases/torpex.py +++ b/hypnotoad/cases/torpex.py @@ -297,6 +297,9 @@ def __init__(self, equilibOptions, meshOptions): f"{eqfile['Z'][0, 0][ZindMid, RindMid]}" ) self.Bt_axis = Bt[ZindMid, RindMid] + + self.f_psi_sign = numpy.sign(0.0 - psi[ZindMid, RindMid]) # H.seto + else: raise ValueError("Failed to initialise psi function from inputs") From a69274d9ea707ec411f3a47aebef88137b9f13a2 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 18 Sep 2024 15:35:50 +0900 Subject: [PATCH 31/60] f_psi_sign has been added for coil in torpex.py for signed Bp --- hypnotoad/cases/torpex.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/hypnotoad/cases/torpex.py b/hypnotoad/cases/torpex.py index 9efaea2c..0a5dee96 100644 --- a/hypnotoad/cases/torpex.py +++ b/hypnotoad/cases/torpex.py @@ -162,6 +162,9 @@ def __init__(self, equilibOptions, meshOptions): self.Rmax = float("inf") self.Zmin = -float("inf") self.Zmax = float("inf") + + self.f_psi_sign = 1.0 # H.seto (dummy) + elif "gfile" in equilibOptions: # load a g-file with open(equilibOptions["gfile"], "rt") as fh: @@ -207,7 +210,7 @@ def __init__(self, equilibOptions, meshOptions): "consistent with sign of grad(psi)" ) - self.f_psi_sign = numpy.sign(psi_bndry - psi_axis) + self.f_psi_sign = numpy.sign(psi_bndry - psi_axis) # H.seto (dummy) # index of a point close to the magnetic axis i_axis = numpy.searchsorted(R, R_axis) @@ -298,7 +301,7 @@ def __init__(self, equilibOptions, meshOptions): ) self.Bt_axis = Bt[ZindMid, RindMid] - self.f_psi_sign = numpy.sign(0.0 - psi[ZindMid, RindMid]) # H.seto + self.f_psi_sign = numpy.sign(0.0 - psi[ZindMid, RindMid]) # H.seto (dummy) else: raise ValueError("Failed to initialise psi function from inputs") From bf4eebb2311f68fd6708866322536579db1da2fd Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 18 Sep 2024 16:53:54 +0900 Subject: [PATCH 32/60] signed Bp with self.bpsign rather than self.f_psi_sign in geometry1 for COCOS=+1 --- hypnotoad/core/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 79687a30..31dfde63 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -853,7 +853,7 @@ def geometry1(self): self.Brxy = self.meshParent.equilibrium.Bp_R(self.Rxy, self.Zxy) self.Bzxy = self.meshParent.equilibrium.Bp_Z(self.Rxy, self.Zxy) - self.Bpxy = numpy.sqrt(self.Brxy**2 + self.Bzxy**2) + self.Bpxy = self.bpsign * numpy.sqrt(self.Brxy**2 + self.Bzxy**2) self.calcPoloidalDistance() From dabf4d0befc49bee29c68f06224b10f203da99af Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 18 Sep 2024 17:12:27 +0900 Subject: [PATCH 33/60] flipping sign of Bp in if statement with Bp_dot_grady has been commented out --- hypnotoad/core/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 31dfde63..fb3a3c3e 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -882,7 +882,7 @@ def geometry1(self): print( "Poloidal field is in opposite direction to Grad(theta) -> Bp negative" ) - self.Bpxy = -self.Bpxy + # self.Bpxy = -self.Bpxy if self.bpsign > 0.0: raise ValueError( "Sign of Bp should be negative? (note this check will raise an " From bc2e8935eba98e51e68bf40b39fbe7ecefbafa96 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 18 Sep 2024 17:26:08 +0900 Subject: [PATCH 34/60] note on capBpYlowXpoint has been added --- hypnotoad/core/mesh.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index fb3a3c3e..bf727411 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -912,6 +912,9 @@ def geometry2(self): # (By default do _not_ do this) # Get rid of minimum in Bpxy.ylow field, because it can cause large spikes in # some metric coefficients, which may cause numerical problems in simulations + # + # note: capBpYlowXpoint assumes a positive definie Bp and doesnt work + # with a negative definit Bp (H.Seto) self.capBpYlowXpoint() self.hy = self.calcHy() From 575226c2cb0af8fb5423f4038c314b9593a4864e Mon Sep 17 00:00:00 2001 From: HarukiST Date: Fri, 20 Sep 2024 09:51:28 +0900 Subject: [PATCH 35/60] use positive definit Jacobian J=h/|Bp| --- hypnotoad/core/mesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index bf727411..74617935 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -1029,7 +1029,7 @@ def calcMetric(self): self.g13 = -self.I * self.g11 self.g23 = -self.dphidy / self.hy**2 - self.J = self.hy / self.Bpxy + self.J = self.hy / numpy.abs(self.Bpxy) self.g_11 = 1.0 / self.g11 + (self.I * self.Rxy) ** 2 self.g_22 = self.hy**2 + (self.Rxy * self.dphidy) ** 2 @@ -1062,7 +1062,7 @@ def calcMetric(self): - self.Rxy * numpy.abs(self.Bpxy) * self.I * self.tanBeta / self.hy ) - self.J = self.hy / self.Bpxy + self.J = self.hy / numpy.abs(self.Bpxy) self.g_11 = ( 1.0 / (self.Rxy * self.Bpxy * self.cosBeta) ** 2 From a25caacb3f265c86c31c2e6ac4c833ef3ec5208b Mon Sep 17 00:00:00 2001 From: HarukiST Date: Fri, 20 Sep 2024 09:55:13 +0900 Subject: [PATCH 36/60] modify Jcheck for Jacobian J=h/|Bp| --- hypnotoad/core/mesh.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 74617935..600e5390 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -1078,16 +1078,12 @@ def calcMetric(self): self.g_23 = self.bpsign * self.dphidy * self.Rxy**2 # check Jacobian is OK - Jcheck = ( - self.bpsign - * 1.0 - / numpy.sqrt( - self.g11 * self.g22 * self.g33 - + 2.0 * self.g12 * self.g13 * self.g23 - - self.g11 * self.g23**2 - - self.g22 * self.g13**2 - - self.g33 * self.g12**2 - ) + Jcheck = 1.0 / numpy.sqrt( + self.g11 * self.g22 * self.g33 + + 2.0 * self.g12 * self.g13 * self.g23 + - self.g11 * self.g23**2 + - self.g22 * self.g13**2 + - self.g33 * self.g12**2 ) # ignore grid points at X-points as J should diverge there (as Bp->0) if Jcheck._corners_array is not None: From 8453924c79ceac83b9900b88cea0ba637de2ce78 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Fri, 20 Sep 2024 10:40:57 +0900 Subject: [PATCH 37/60] add Jpar0 based on output-jpar0 branch --- hypnotoad/core/mesh.py | 43 ++++++++++++++++++++++++++++++++++++++---- 1 file changed, 39 insertions(+), 4 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 600e5390..6d96f6db 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -902,6 +902,27 @@ 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' + # Note: COCOS +1 convention is assumed + # so mu0 * Jpar = -B f' - mu0 p' f / B + + 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 @@ -1338,10 +1359,21 @@ def curl_bOverB_zetahat(R, Z): self.bxcvy = self.Bxy / 2.0 * self.curl_bOverB_y self.bxcvz = self.Bxy / 2.0 * self.curl_bOverB_z elif self.user_options.curvature_type == "bxkappa": - raise ValueError("bxkappa form of curvature not implemented yet") - self.bxcvx = float("nan") - self.bxcvy = float("nan") - self.bxcvz = float("nan") + + # bxkappa for cocos1-extension by H. Seto (QST) + + bxcvu = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDX("#Bxy") + bxcvv = -self.DDX("#Btxy*#Rxy/#Bxy") / self.J + bxcvw = self.DDX( + "#Bxy*#hy/#Bpxy" + ) / self.J - self.Btxy * self.Rxy / self.J / self.Bxy * self.DDX( + "#Btxy*#hy/#Bpxy/#Rxy" + ) + + self.bxcvx = bxcvu + self.bxcvy = bxcvv + self.bxcvz = bxcvw - self.I * bxcvu + else: raise ValueError( "Unrecognized option '" @@ -3661,6 +3693,9 @@ def addFromRegionsXArray(name): if hasattr(next(iter(self.equilibrium.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)), From 7705e2cac194721ed24aa2ada061a2809408852a Mon Sep 17 00:00:00 2001 From: HarukiST Date: Thu, 26 Sep 2024 16:35:11 +0900 Subject: [PATCH 38/60] add dummy sinty for backward compatibility --- hypnotoad/core/mesh.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 6d96f6db..d0f75fe9 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -1034,6 +1034,7 @@ def calcMetric(self): # derivatives. This means different (radial) regions can use different # locations for where zShift=0. self.I = MultiLocationArray(self.nx, self.ny).zero() + self.sinty = self.I # Here ShiftTorsion = d2phidxdy # Haven't checked this is exactly the quantity needed by BOUT++... @@ -1361,7 +1362,7 @@ def curl_bOverB_zetahat(R, Z): elif self.user_options.curvature_type == "bxkappa": # bxkappa for cocos1-extension by H. Seto (QST) - + bxcvu = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDX("#Bxy") bxcvv = -self.DDX("#Btxy*#Rxy/#Bxy") / self.J bxcvw = self.DDX( @@ -3664,8 +3665,9 @@ def addFromRegionsXArray(name): # IntShiftTorsion should never be used. It is only for some 'BOUT-06 style # differencing'. IntShiftTorsion is not written by Hypnotoad1, so don't write # here. /JTO 19/5/2019 - if not self.user_options.shiftedmetric: - addFromRegions("sinty") + # if not self.user_options.shiftedmetric: + # addFromRegions("sinty") + addFromRegions("sinty") addFromRegions("g11") addFromRegions("g22") addFromRegions("g33") From ed2c33ca4fe4a5a1718c39d641a3180154a01eb4 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Thu, 26 Sep 2024 16:35:32 +0900 Subject: [PATCH 39/60] add dummy sinty for backward compatibility (black) --- hypnotoad/core/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index d0f75fe9..1914bd09 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -1362,7 +1362,7 @@ def curl_bOverB_zetahat(R, Z): elif self.user_options.curvature_type == "bxkappa": # bxkappa for cocos1-extension by H. Seto (QST) - + bxcvu = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDX("#Bxy") bxcvv = -self.DDX("#Btxy*#Rxy/#Bxy") / self.J bxcvw = self.DDX( From f5b8f87421ffa5d763e7219874b9783047791418 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Thu, 26 Sep 2024 16:37:13 +0900 Subject: [PATCH 40/60] fix bxcvx in bxkappae --- hypnotoad/core/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 1914bd09..31689a73 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -1363,7 +1363,7 @@ def curl_bOverB_zetahat(R, Z): # bxkappa for cocos1-extension by H. Seto (QST) - bxcvu = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDX("#Bxy") + bxcvu = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDY("#Bxy") bxcvv = -self.DDX("#Btxy*#Rxy/#Bxy") / self.J bxcvw = self.DDX( "#Bxy*#hy/#Bpxy" From de25952e8f607f02fea7f2009cbc4aabaee41ae1 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 1 Oct 2024 14:28:27 +0900 Subject: [PATCH 41/60] fix curvature terms in bxkappae --- hypnotoad/core/mesh.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 31689a73..1e69ddb0 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -1363,14 +1363,18 @@ def curl_bOverB_zetahat(R, Z): # bxkappa for cocos1-extension by H. Seto (QST) - bxcvu = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDY("#Bxy") - bxcvv = -self.DDX("#Btxy*#Rxy/#Bxy") / self.J - bxcvw = self.DDX( + curlxb0u = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDY("#Bxy") + curlxb0w = self.DDX( "#Bxy*#hy/#Bpxy" ) / self.J - self.Btxy * self.Rxy / self.J / self.Bxy * self.DDX( "#Btxy*#hy/#Bpxy/#Rxy" ) + coefv = -self.Btxy * self.Bpxy * self.Rxy / (self.hy * self.Bxy**2) + bxcvu = curlxb0u + bxcvv = curlxb0w * coefv + bxcvw = curlxb0w + self.bxcvx = bxcvu self.bxcvy = bxcvv self.bxcvz = bxcvw - self.I * bxcvu From aeaddb931548fd79d0aec11ea8d2a4e4b63e06ad Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 1 Oct 2024 14:28:41 +0900 Subject: [PATCH 42/60] fix curvature terms in bxkappae --- hypnotoad/core/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 1e69ddb0..f6d79150 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -1372,7 +1372,7 @@ def curl_bOverB_zetahat(R, Z): coefv = -self.Btxy * self.Bpxy * self.Rxy / (self.hy * self.Bxy**2) bxcvu = curlxb0u - bxcvv = curlxb0w * coefv + bxcvv = curlxb0w * coefv bxcvw = curlxb0w self.bxcvx = bxcvu From 2617c5ddf8a2b362e1ee8ebe8686622e9f38580c Mon Sep 17 00:00:00 2001 From: HarukiST Date: Fri, 4 Oct 2024 11:42:15 +0900 Subject: [PATCH 43/60] add in bxkappa2 and bxkappa3 --- hypnotoad/core/mesh.py | 51 ++++++++++++++++++++++++++++++++++++++---- 1 file changed, 47 insertions(+), 4 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index f6d79150..98f612bf 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -82,7 +82,13 @@ class MeshRegion: "curl(b/B)", doc="Expression used to calculate curvature operator 'bxcv'", value_type=str, - allowed=["curl(b/B)", "curl(b/B) with x-y derivatives", "bxkappa"], + allowed=[ + "curl(b/B)", + "curl(b/B) with x-y derivatives", + "bxkappa", + "bxkappa2", + "bxkappa3", + ], ), curvature_smoothing=WithMeta( None, @@ -1084,7 +1090,7 @@ def calcMetric(self): - self.Rxy * numpy.abs(self.Bpxy) * self.I * self.tanBeta / self.hy ) - self.J = self.hy / numpy.abs(self.Bpxy) + self.J = self.bpsign * self.hy / self.Bpxy self.g_11 = ( 1.0 / (self.Rxy * self.Bpxy * self.cosBeta) ** 2 @@ -1359,10 +1365,10 @@ def curl_bOverB_zetahat(R, Z): self.bxcvx = self.Bxy / 2.0 * self.curl_bOverB_x self.bxcvy = self.Bxy / 2.0 * self.curl_bOverB_y self.bxcvz = self.Bxy / 2.0 * self.curl_bOverB_z - elif self.user_options.curvature_type == "bxkappa": - # bxkappa for cocos1-extension by H. Seto (QST) + elif self.user_options.curvature_type == "bxkappa": + # b0 x kappa terms for cocos1-extension by H. Seto (QST) curlxb0u = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDY("#Bxy") curlxb0w = self.DDX( "#Bxy*#hy/#Bpxy" @@ -1379,6 +1385,43 @@ def curl_bOverB_zetahat(R, Z): self.bxcvy = bxcvv self.bxcvz = bxcvw - self.I * bxcvu + elif self.user_options.curvature_type == "bxkappa2": + + curlxb0u = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDY("#Bxy") + curlxb0w = self.bpsign * ( + self.bpsign * self.Bpxy**2 / self.Bxy / self.J * self.DDX("#hy/#Bpxy") + + self.DDX("#Bxy") + - self.Btxy * self.Rxy / self.Bxy * self.DDX("#Btxy/#Rxy") + ) + + coefv = -self.Btxy * self.Bpxy * self.Rxy / (self.hy * self.Bxy**2) + bxcvu = curlxb0u + bxcvv = curlxb0w * coefv + bxcvw = curlxb0w + + self.bxcvx = bxcvu + self.bxcvy = bxcvv + self.bxcvz = bxcvw - self.I * bxcvu + + elif self.user_options.curvature_type == "bxkappa3": + + pprime = self.meshParent.equilibrium.regions[ + self.equilibriumRegion.name + ].pprime(self.psixy) + + mu0 = 4e-7 * numpy.pi + curlxb0u = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDY("#Bxy") + curlxb0w = -self.bpsign * (mu0 / self.Bxy * pprime + self.DDX("#Bxy")) + + coefv = -self.Btxy * self.Bpxy * self.Rxy / (self.hy * self.Bxy**2) + bxcvu = curlxb0u + bxcvv = curlxb0w * coefv + bxcvw = curlxb0w + + self.bxcvx = bxcvu + self.bxcvy = bxcvv + self.bxcvz = bxcvw - self.I * bxcvu + else: raise ValueError( "Unrecognized option '" From 36fb632621882ca40a495d1c21e05d51f9ff2d75 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Fri, 4 Oct 2024 14:37:26 +0900 Subject: [PATCH 44/60] re-define BOUT-3 style curvature with pre-cancellation as bxkappa --- hypnotoad/core/mesh.py | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 98f612bf..200f3a25 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -1367,13 +1367,13 @@ def curl_bOverB_zetahat(R, Z): self.bxcvz = self.Bxy / 2.0 * self.curl_bOverB_z elif self.user_options.curvature_type == "bxkappa": - - # b0 x kappa terms for cocos1-extension by H. Seto (QST) + # b0 x kappa terms with pre-cancellation for cocos1-extension + # by H. Seto (QST) curlxb0u = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDY("#Bxy") - curlxb0w = self.DDX( - "#Bxy*#hy/#Bpxy" - ) / self.J - self.Btxy * self.Rxy / self.J / self.Bxy * self.DDX( - "#Btxy*#hy/#Bpxy/#Rxy" + curlxb0w = self.bpsign * ( + self.bpsign * self.Bpxy**2 / self.Bxy / self.J * self.DDX("#hy/#Bpxy") + + self.DDX("#Bxy") + - self.Btxy * self.Rxy / self.Bxy * self.DDX("#Btxy/#Rxy") ) coefv = -self.Btxy * self.Bpxy * self.Rxy / (self.hy * self.Bxy**2) @@ -1385,13 +1385,16 @@ def curl_bOverB_zetahat(R, Z): self.bxcvy = bxcvv self.bxcvz = bxcvw - self.I * bxcvu + elif self.user_options.curvature_type == "bxkappa2": + # b0 x kappa terms for cocos1-extension by H. Seto (QST) + curlxb0u = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDY("#Bxy") - curlxb0w = self.bpsign * ( - self.bpsign * self.Bpxy**2 / self.Bxy / self.J * self.DDX("#hy/#Bpxy") - + self.DDX("#Bxy") - - self.Btxy * self.Rxy / self.Bxy * self.DDX("#Btxy/#Rxy") + curlxb0w = self.DDX( + "#Bxy*#hy/#Bpxy" + ) / self.J - self.Btxy * self.Rxy / self.J / self.Bxy * self.DDX( + "#Btxy*#hy/#Bpxy/#Rxy" ) coefv = -self.Btxy * self.Bpxy * self.Rxy / (self.hy * self.Bxy**2) @@ -1405,6 +1408,9 @@ def curl_bOverB_zetahat(R, Z): elif self.user_options.curvature_type == "bxkappa3": + # b0 x kappa terms with pprime for cocos1-extension + # by H. Seto (QST) + pprime = self.meshParent.equilibrium.regions[ self.equilibriumRegion.name ].pprime(self.psixy) From b5331ed9fd8710678c3bc15cf10e4fa2f86640f4 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 8 Oct 2024 17:59:20 +0900 Subject: [PATCH 45/60] dx should be always positive definie --- hypnotoad/core/mesh.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 200f3a25..5490e1e0 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -839,9 +839,6 @@ def geometry1(self): """ self.psixy = self.meshParent.equilibrium.psi(self.Rxy, self.Zxy) - self.dx = MultiLocationArray(self.nx, self.ny) - self.dx.centre = (self.psi_vals[2::2] - self.psi_vals[:-2:2])[:, numpy.newaxis] - self.dx.ylow = (self.psi_vals[2::2] - self.psi_vals[:-2:2])[:, numpy.newaxis] if self.psi_vals[0] > self.psi_vals[-1]: # x-coordinate is -psixy so x always increases radially across grid @@ -851,6 +848,11 @@ def geometry1(self): self.bpsign = 1.0 self.xcoord = self.psixy + # dx must be positive definite + self.dx = MultiLocationArray(self.nx, self.ny) + self.dx.centre = self.bpsign*(self.psi_vals[2::2] - self.psi_vals[:-2:2])[:, numpy.newaxis] + self.dx.ylow = self.bpsign*(self.psi_vals[2::2] - self.psi_vals[:-2:2])[:, numpy.newaxis] + self.dy = MultiLocationArray(self.nx, self.ny) self.dy.centre = self.meshParent.dy_scalar self.dy.ylow = self.meshParent.dy_scalar From ffb98cd30bde78bac40402c6b7da4c9169fdd67a Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 13 Nov 2024 08:28:17 +0900 Subject: [PATCH 46/60] comment fix --- hypnotoad/core/equilibrium.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hypnotoad/core/equilibrium.py b/hypnotoad/core/equilibrium.py index d4ed8153..7ccba0d7 100644 --- a/hypnotoad/core/equilibrium.py +++ b/hypnotoad/core/equilibrium.py @@ -3950,7 +3950,7 @@ def magneticFunctionsFromGrid(self, R, Z, psiRZ, option): @Equilibrium.handleMultiLocationArray def psi(self, R, Z): - "Return the poloidal flux at the given (R,Z) location" + """Return the poloidal flux at the given (R,Z) location""" return self.psi_func(R, Z, grid=False) # The __get__(self) call converts the function to a 'bound From 37f12e3ed3119eb03d2a2b0ad6dd1cda1708b226 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 13 Nov 2024 08:31:13 +0900 Subject: [PATCH 47/60] add netcdf3_64bit for BOUTv3.1 --- hypnotoad/core/mesh.py | 144 +++++++++++++++++++++++++---------------- 1 file changed, 88 insertions(+), 56 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 5490e1e0..d070a9a4 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -73,6 +73,14 @@ class MeshRegion: user_options_factory = OptionsFactory( EquilibriumRegion.user_options_factory, + output_format=WithMeta( + "NETCDF4", + doc="format of grid file", + value_type=str, + allowed=[ + "NETCDF4", + "NETCDF3_64BIT"], + ), shiftedmetric=WithMeta( True, doc="Is grid generated for paralleltransform=ShiftedMetric?", @@ -871,7 +879,17 @@ def geometry1(self): self.pressure = self.meshParent.equilibrium.regions[ self.equilibriumRegion.name ].pressure(self.psixy) - + """ + if numpy.min(self.pressure) < 0.0: + print("minimum pressure is negative, add vacuum pressure") + vacuum_pressure = -numpy.min(self.pressure) + self.pressure += vacuum_pressure + + if numpy.min(self.pressure) < numpy.max(self.pressure)*1.e-3: + print("minimum pressure is too small, add vacuum pressure") + vacuum_pressure = numpy.max(self.pressure)*1.e-3 + self.pressure += vacuum_pressure + """ # determine direction - dot Bp with Grad(y) vector # evaluate in 'sol' at outer radial boundary Bp_dot_grady = self.Brxy.centre[-1, self.ny // 2] * ( @@ -947,7 +965,7 @@ def geometry2(self): self.capBpYlowXpoint() self.hy = self.calcHy() - + if not self.user_options.orthogonal: # Calculate beta (angle between e_x and Grad(x), also the angle between e_y # and Grad(y)), used for non-orthogonal grid @@ -1048,8 +1066,10 @@ def calcMetric(self): # Haven't checked this is exactly the quantity needed by BOUT++... # ShiftTorsion is only used in Curl operator - Curl is rarely used. self.ShiftTorsion = self.DDX("#dphidy") - + self.hthe = self.hy/1.0 + if self.user_options.orthogonal: + self.g11 = (self.Rxy * self.Bpxy) ** 2 self.g22 = 1.0 / self.hy**2 self.g33 = ( @@ -3701,14 +3721,15 @@ def addFromRegionsXArray(name): addFromRegions("psixy") addFromRegions("dx") addFromRegions("dy") - addFromRegions("poloidal_distance") - addFromRegionsXArray("total_poloidal_distance") + #addFromRegions("poloidal_distance") + #addFromRegionsXArray("total_poloidal_distance") addFromRegions("Brxy") addFromRegions("Bzxy") addFromRegions("Bpxy") addFromRegions("Btxy") addFromRegions("Bxy") addFromRegions("hy") + addFromRegions("hthe") # if not self.user_options.orthogonal: # addFromRegions('beta') # addFromRegions('eta') @@ -3723,19 +3744,20 @@ def addFromRegionsXArray(name): # if not self.user_options.shiftedmetric: # addFromRegions("sinty") addFromRegions("sinty") - addFromRegions("g11") - addFromRegions("g22") - addFromRegions("g33") - addFromRegions("g12") - addFromRegions("g13") - addFromRegions("g23") - addFromRegions("J") - addFromRegions("g_11") - addFromRegions("g_22") - addFromRegions("g_33") - addFromRegions("g_12") - addFromRegions("g_13") - addFromRegions("g_23") + if self.user_options.output_format == "NETCDF4": + addFromRegions("g11") + addFromRegions("g22") + addFromRegions("g33") + addFromRegions("g12") + addFromRegions("g13") + addFromRegions("g23") + addFromRegions("J") + addFromRegions("g_11") + addFromRegions("g_22") + addFromRegions("g_33") + addFromRegions("g_12") + addFromRegions("g_13") + addFromRegions("g_23") if self.user_options.curvature_type in ( "curl(b/B) with x-y derivatives", "curl(b/B)", @@ -3746,13 +3768,15 @@ def addFromRegionsXArray(name): addFromRegions("bxcvx") addFromRegions("bxcvy") addFromRegions("bxcvz") - + """ if hasattr(next(iter(self.equilibrium.regions.values())), "pressure"): addFromRegions("pressure") if hasattr(next(iter(self.regions.values())), "Jpar0"): addFromRegions("Jpar0") - + """ + addFromRegions("pressure") + addFromRegions("Jpar0") # Penalty mask self.penalty_mask = BoutArray( numpy.zeros((self.nx, self.ny)), @@ -3798,8 +3822,11 @@ def writeArrayXDirection(self, name, array, f): def writeGridfile(self, filename): from boututils.datafile import DataFile + from netCDF4 import stringtochar - with DataFile(filename, create=True, format="NETCDF4") as f: + + # with DataFile(filename, create=True, format="NETCDF4") as f: + with DataFile(filename, create=True, format=self.user_options.output_format) as f: # Save unique ID for grid file import uuid @@ -3808,7 +3835,10 @@ def writeGridfile(self, filename): # ny for BOUT++ excludes boundary guard cells f.write("ny", self.ny_noguards) f.write("y_boundary_guards", self.user_options.y_boundary_guards) - f.write("curvature_type", self.user_options.curvature_type) + + if self.user_options.output_format == "NETCDF4": + f.write("curvature_type", self.user_options.curvature_type) + f.write("Bt_axis", self.equilibrium.Bt_axis) if hasattr(self.equilibrium, "psi_axis"): @@ -3835,9 +3865,9 @@ def writeGridfile(self, filename): for name in ["Rxy", "Zxy"]: self.writeCorners(name, self.__dict__[name], f) - if self.user_options.orthogonal: - # Also write hy as "hthe" for backward compatibility - self.writeArray("hthe", self.hy, f) + #if self.user_options.orthogonal: + # # Also write hy as "hthe" for backward compatibility + # self.writeArray("hthe", self.hy, f) # write the 1d fields for name in self.arrayXDirection_to_output: @@ -4004,42 +4034,44 @@ def writeGridfile(self, filename): + "\nMesh options:\n" + self.user_options.as_table() ) - f.write("hypnotoad_inputs", inputs_string) - options_dict = dict(self.equilibrium.user_options) - options_dict.update(self.equilibrium.nonorthogonal_options) - options_dict.update(self.user_options) - f.write("hypnotoad_inputs_yaml", yaml.dump(options_dict)) - f.write_file_attribute("hypnotoad_version", self.version) - if self.git_hash is not None: - f.write_file_attribute("hypnotoad_git_hash", self.git_hash) - f.write_file_attribute( - "hypnotoad_git_diff", - self.git_diff if self.git_diff is not None else "", - ) + if self.user_options.output_format == "NETCDF4": + f.write("hypnotoad_inputs", inputs_string) + options_dict = dict(self.equilibrium.user_options) + options_dict.update(self.equilibrium.nonorthogonal_options) + options_dict.update(self.user_options) + f.write("hypnotoad_inputs_yaml", yaml.dump(options_dict)) + + f.write_file_attribute("hypnotoad_version", self.version) + if self.git_hash is not None: + f.write_file_attribute("hypnotoad_git_hash", self.git_hash) + f.write_file_attribute( + "hypnotoad_git_diff", + self.git_diff if self.git_diff is not None else "", + ) - if hasattr(self.equilibrium, "geqdsk_filename"): - # If grid was created from a geqdsk file, save the file name - f.write_file_attribute( + if hasattr(self.equilibrium, "geqdsk_filename"): + # If grid was created from a geqdsk file, save the file name + f.write_file_attribute( "hypnotoad_geqdsk_filename", self.equilibrium.geqdsk_filename - ) + ) - if hasattr(self.equilibrium, "geqdsk_input"): - # If grid was created from a geqdsk file, save the file contents - # - # Write as string variable and not attribute because the string will be - # long and attributes are printed by 'ncdump -h' or by ncdump when - # looking at a different variable, which would be inconvenient. It is - # not likely that we need to load the geqdsk file contents in BOUT++, so - # no reason to save as an attribute. - f.write( - "hypnotoad_input_geqdsk_file_contents", - self.equilibrium.geqdsk_input, - ) + if hasattr(self.equilibrium, "geqdsk_input"): + # If grid was created from a geqdsk file, save the file contents + # + # Write as string variable and not attribute because the string will be + # long and attributes are printed by 'ncdump -h' or by ncdump when + # looking at a different variable, which would be inconvenient. It is + # not likely that we need to load the geqdsk file contents in BOUT++, so + # no reason to save as an attribute. + f.write( + "hypnotoad_input_geqdsk_file_contents", + self.equilibrium.geqdsk_input, + ) - # save Python and module versions to enable reproducibility - f.write("Python_version", sys.version) - f.write("module_versions", module_versions_formatted()) + # save Python and module versions to enable reproducibility + f.write("Python_version", sys.version) + f.write("module_versions", module_versions_formatted()) def plot2D(self, f, title=None): from matplotlib import pyplot From 0dcac982cc2410e586758aa03d740d77e7e50284 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 13 Nov 2024 08:33:52 +0900 Subject: [PATCH 48/60] black fix --- hypnotoad/cases/tokamak.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 1c30a776..769f4e3d 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -427,7 +427,6 @@ def __init__( # p_SOL = p_lcfs * exp( (psi - psi_lcfs) * dpdpsi / p_lcfs) p_SOL = p_lcfs * np.exp((psiSOL - psi_lcfs) * dpdpsi / p_lcfs) pressure = np.concatenate([pressure, p_SOL]) - if pprime is not None: pprime = np.concatenate([pprime, dpdpsi * p_SOL / p_lcfs]) @@ -470,6 +469,7 @@ def __init__( self.p_spl = interpolate.InterpolatedUnivariateSpline( xcoord, pressure, ext=3 ) + if pprime is not None: self.pprime_spl = interpolate.InterpolatedUnivariateSpline( xcoord, pprime, ext=1 From d7af0736f9e2b2d7735f87715ac9e14833e63587 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Wed, 13 Nov 2024 08:36:13 +0900 Subject: [PATCH 49/60] black fix --- hypnotoad/core/mesh.py | 41 +++++++++++++++++++++-------------------- 1 file changed, 21 insertions(+), 20 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index d070a9a4..b630d260 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -77,9 +77,7 @@ class MeshRegion: "NETCDF4", doc="format of grid file", value_type=str, - allowed=[ - "NETCDF4", - "NETCDF3_64BIT"], + allowed=["NETCDF4", "NETCDF3_64BIT"], ), shiftedmetric=WithMeta( True, @@ -847,7 +845,6 @@ def geometry1(self): """ self.psixy = self.meshParent.equilibrium.psi(self.Rxy, self.Zxy) - if self.psi_vals[0] > self.psi_vals[-1]: # x-coordinate is -psixy so x always increases radially across grid self.bpsign = -1.0 @@ -858,8 +855,12 @@ def geometry1(self): # dx must be positive definite self.dx = MultiLocationArray(self.nx, self.ny) - self.dx.centre = self.bpsign*(self.psi_vals[2::2] - self.psi_vals[:-2:2])[:, numpy.newaxis] - self.dx.ylow = self.bpsign*(self.psi_vals[2::2] - self.psi_vals[:-2:2])[:, numpy.newaxis] + self.dx.centre = ( + self.bpsign * (self.psi_vals[2::2] - self.psi_vals[:-2:2])[:, numpy.newaxis] + ) + self.dx.ylow = ( + self.bpsign * (self.psi_vals[2::2] - self.psi_vals[:-2:2])[:, numpy.newaxis] + ) self.dy = MultiLocationArray(self.nx, self.ny) self.dy.centre = self.meshParent.dy_scalar @@ -889,7 +890,7 @@ def geometry1(self): print("minimum pressure is too small, add vacuum pressure") vacuum_pressure = numpy.max(self.pressure)*1.e-3 self.pressure += vacuum_pressure - """ + """ # determine direction - dot Bp with Grad(y) vector # evaluate in 'sol' at outer radial boundary Bp_dot_grady = self.Brxy.centre[-1, self.ny // 2] * ( @@ -965,7 +966,7 @@ def geometry2(self): self.capBpYlowXpoint() self.hy = self.calcHy() - + if not self.user_options.orthogonal: # Calculate beta (angle between e_x and Grad(x), also the angle between e_y # and Grad(y)), used for non-orthogonal grid @@ -1066,8 +1067,8 @@ def calcMetric(self): # Haven't checked this is exactly the quantity needed by BOUT++... # ShiftTorsion is only used in Curl operator - Curl is rarely used. self.ShiftTorsion = self.DDX("#dphidy") - self.hthe = self.hy/1.0 - + self.hthe = self.hy / 1.0 + if self.user_options.orthogonal: self.g11 = (self.Rxy * self.Bpxy) ** 2 @@ -1407,11 +1408,10 @@ def curl_bOverB_zetahat(R, Z): self.bxcvy = bxcvv self.bxcvz = bxcvw - self.I * bxcvu - elif self.user_options.curvature_type == "bxkappa2": # b0 x kappa terms for cocos1-extension by H. Seto (QST) - + curlxb0u = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDY("#Bxy") curlxb0w = self.DDX( "#Bxy*#hy/#Bpxy" @@ -1432,7 +1432,7 @@ def curl_bOverB_zetahat(R, Z): # b0 x kappa terms with pprime for cocos1-extension # by H. Seto (QST) - + pprime = self.meshParent.equilibrium.regions[ self.equilibriumRegion.name ].pprime(self.psixy) @@ -3721,8 +3721,8 @@ def addFromRegionsXArray(name): addFromRegions("psixy") addFromRegions("dx") addFromRegions("dy") - #addFromRegions("poloidal_distance") - #addFromRegionsXArray("total_poloidal_distance") + # addFromRegions("poloidal_distance") + # addFromRegionsXArray("total_poloidal_distance") addFromRegions("Brxy") addFromRegions("Bzxy") addFromRegions("Bpxy") @@ -3824,9 +3824,10 @@ def writeGridfile(self, filename): from boututils.datafile import DataFile from netCDF4 import stringtochar - # with DataFile(filename, create=True, format="NETCDF4") as f: - with DataFile(filename, create=True, format=self.user_options.output_format) as f: + with DataFile( + filename, create=True, format=self.user_options.output_format + ) as f: # Save unique ID for grid file import uuid @@ -3838,7 +3839,7 @@ def writeGridfile(self, filename): if self.user_options.output_format == "NETCDF4": f.write("curvature_type", self.user_options.curvature_type) - + f.write("Bt_axis", self.equilibrium.Bt_axis) if hasattr(self.equilibrium, "psi_axis"): @@ -3865,7 +3866,7 @@ def writeGridfile(self, filename): for name in ["Rxy", "Zxy"]: self.writeCorners(name, self.__dict__[name], f) - #if self.user_options.orthogonal: + # if self.user_options.orthogonal: # # Also write hy as "hthe" for backward compatibility # self.writeArray("hthe", self.hy, f) @@ -4053,7 +4054,7 @@ def writeGridfile(self, filename): if hasattr(self.equilibrium, "geqdsk_filename"): # If grid was created from a geqdsk file, save the file name f.write_file_attribute( - "hypnotoad_geqdsk_filename", self.equilibrium.geqdsk_filename + "hypnotoad_geqdsk_filename", self.equilibrium.geqdsk_filename ) if hasattr(self.equilibrium, "geqdsk_input"): From 598fd8ed9e38f24a5d9116811732944c4e310443 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Fri, 15 Nov 2024 21:54:39 +0900 Subject: [PATCH 50/60] delete useless duplication of hthe and hy --- hypnotoad/core/mesh.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index b630d260..0c2588f7 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -1067,7 +1067,6 @@ def calcMetric(self): # Haven't checked this is exactly the quantity needed by BOUT++... # ShiftTorsion is only used in Curl operator - Curl is rarely used. self.ShiftTorsion = self.DDX("#dphidy") - self.hthe = self.hy / 1.0 if self.user_options.orthogonal: @@ -3729,7 +3728,7 @@ def addFromRegionsXArray(name): addFromRegions("Btxy") addFromRegions("Bxy") addFromRegions("hy") - addFromRegions("hthe") + # if not self.user_options.orthogonal: # addFromRegions('beta') # addFromRegions('eta') @@ -3866,9 +3865,9 @@ def writeGridfile(self, filename): for name in ["Rxy", "Zxy"]: self.writeCorners(name, self.__dict__[name], f) - # if self.user_options.orthogonal: - # # Also write hy as "hthe" for backward compatibility - # self.writeArray("hthe", self.hy, f) + if self.user_options.orthogonal: + # Also write hy as "hthe" for backward compatibility + self.writeArray("hthe", self.hy, f) # write the 1d fields for name in self.arrayXDirection_to_output: From 4a5dde043b0ea77714cc663fdeddb925e97ea492 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Thu, 16 Jan 2025 10:34:38 +0900 Subject: [PATCH 51/60] set dx at separatrix without pf regions --- hypnotoad/cases/tokamak.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 769f4e3d..dbc22244 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -856,8 +856,9 @@ def describeSingleNull(self): dpsidi_pf = (self.psi_sep[0] - self.psi_pf_upper) / self.user_options.nx_pf # Get the smallest absolute grid spacing for the separatrix - dpsidi_sep = min([dpsidi_sol, dpsidi_core, dpsidi_pf], key=abs) - + #dpsidi_sep = min([dpsidi_sol, dpsidi_core, dpsidi_pf], key=abs) + dpsidi_sep = min([dpsidi_sol, dpsidi_core], key=abs) # set dx at separatrix + # decrease (assuming the factor is <1) the spacing around the separatrix by the # factor psi_spacing_separatrix_multiplier if self.user_options.psi_spacing_separatrix_multiplier is not None: From 4ebd3346d6c4b491ff13ee84398a99bcd72440c4 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Mon, 19 May 2025 10:35:15 +0900 Subject: [PATCH 52/60] set minimum grid width around sepatatrix --- hypnotoad/cases/tokamak.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index dbc22244..912916cd 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -856,9 +856,14 @@ def describeSingleNull(self): dpsidi_pf = (self.psi_sep[0] - self.psi_pf_upper) / self.user_options.nx_pf # Get the smallest absolute grid spacing for the separatrix + # dpsidi_pf_minimum is introduced to avoid very narrow radial grid width + # around separatrix in core and SOL region and 0.3 is an emprical factor. + #dpsidi_sep = min([dpsidi_sol, dpsidi_core, dpsidi_pf], key=abs) dpsidi_sep = min([dpsidi_sol, dpsidi_core], key=abs) # set dx at separatrix - + #dpsidi_pf_minimum = max([0.3*dpsidi_sep, 0.3*dpsidi_pf], key=abs) + dpsidi_pf_minimum = max([0.3*dpsidi_sep, dpsidi_pf], key=abs) + # decrease (assuming the factor is <1) the spacing around the separatrix by the # factor psi_spacing_separatrix_multiplier if self.user_options.psi_spacing_separatrix_multiplier is not None: @@ -900,6 +905,7 @@ def describeSingleNull(self): "psi_start": self.psi_pf_lower, "psi_end": psi_sep, "grad_end": dpsidi_sep, + "grad_start": dpsidi_pf_minimum, } leg_regions = { @@ -950,6 +956,7 @@ def describeSingleNull(self): "psi_start": self.psi_pf_upper, "psi_end": psi_sep, "grad_end": dpsidi_sep, + "grad_start": dpsidi_pf_minimum, } leg_regions = { From 218c5014ae453a88919e611244ce6c5f98aa7605 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Mon, 19 May 2025 18:06:12 +0900 Subject: [PATCH 53/60] clean pressure and Jpar --- hypnotoad/core/mesh.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index 0c2588f7..e0c0c0a3 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -3767,15 +3767,10 @@ def addFromRegionsXArray(name): addFromRegions("bxcvx") addFromRegions("bxcvy") addFromRegions("bxcvz") - """ - if hasattr(next(iter(self.equilibrium.regions.values())), "pressure"): - addFromRegions("pressure") - if hasattr(next(iter(self.regions.values())), "Jpar0"): - addFromRegions("Jpar0") - """ addFromRegions("pressure") addFromRegions("Jpar0") + # Penalty mask self.penalty_mask = BoutArray( numpy.zeros((self.nx, self.ny)), From 80b01e57ecb05d7bc09838af1bf61e1f97063c82 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 20 May 2025 13:19:57 +0900 Subject: [PATCH 54/60] add dummy data for kinetic profiles --- hypnotoad/core/mesh.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index e0c0c0a3..fb41ed49 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -950,6 +950,18 @@ def geometry1(self): self.Bxy * fprime / mu0 + self.Rxy * self.Btxy * pprime / self.Bxy ) + + # set dummy data for kinetic quantities if kinetic_data is True. H. Seto (QST) + # they will be added in the post-script with peqdsk file + # since peqdsk file doesn't have uniform file format. + # + # ni0: bulk ion number density + # resistivity: Spitzer resistivity + # + if self.user_options.kinetic_data: + self.ni0 = MultiLocationArray(self.nx, self.ny).zero() + self.resistivity = MultiLocationArray(self.nx, self.ny).zero() + def geometry2(self): """ Continuation of geometry1(), but needs neighbours to have calculated Bp so called @@ -3770,7 +3782,11 @@ def addFromRegionsXArray(name): addFromRegions("pressure") addFromRegions("Jpar0") - + + if self.user_options.kinetic_data: + addFromRegions("ni0") + addFromRegions("resistivity") + # Penalty mask self.penalty_mask = BoutArray( numpy.zeros((self.nx, self.ny)), From 756be55393ba64b74024df4c17ecf7de839703c2 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 20 May 2025 13:48:19 +0900 Subject: [PATCH 55/60] add user option of kinetic profiles --- hypnotoad/cases/tokamak.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 912916cd..ed0efdd8 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -307,6 +307,16 @@ class TokamakEquilibrium(Equilibrium): value_type=int, check_all=is_positive, ), + # + # set dummy data of kinetic profile for post-process + # + kinetic_profiles=WithMeta( + False, + doc=( + "Set dummy data (zero-fill) of kinetic profiles for post-process" + ), + value_type=bool, + ), ) def __init__( From 30ee20a7897e7d955d1b3b8a3b0feafb7e74dd66 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Tue, 20 May 2025 13:49:09 +0900 Subject: [PATCH 56/60] name of option of kinetic profiles has been changed --- hypnotoad/core/mesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index fb41ed49..c9020b65 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -958,7 +958,7 @@ def geometry1(self): # ni0: bulk ion number density # resistivity: Spitzer resistivity # - if self.user_options.kinetic_data: + if self.user_options.kinetic_profiles: self.ni0 = MultiLocationArray(self.nx, self.ny).zero() self.resistivity = MultiLocationArray(self.nx, self.ny).zero() @@ -3783,7 +3783,7 @@ def addFromRegionsXArray(name): addFromRegions("pressure") addFromRegions("Jpar0") - if self.user_options.kinetic_data: + if self.user_options.kinetic_profiles: addFromRegions("ni0") addFromRegions("resistivity") From edbe03e4c337c361ab42501f9c99e5618a84b2ed Mon Sep 17 00:00:00 2001 From: HarukiST Date: Fri, 23 May 2025 15:16:41 +0900 Subject: [PATCH 57/60] delete options for kinetic profiles --- hypnotoad/core/mesh.py | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index c9020b65..d46e2279 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -950,17 +950,6 @@ def geometry1(self): self.Bxy * fprime / mu0 + self.Rxy * self.Btxy * pprime / self.Bxy ) - - # set dummy data for kinetic quantities if kinetic_data is True. H. Seto (QST) - # they will be added in the post-script with peqdsk file - # since peqdsk file doesn't have uniform file format. - # - # ni0: bulk ion number density - # resistivity: Spitzer resistivity - # - if self.user_options.kinetic_profiles: - self.ni0 = MultiLocationArray(self.nx, self.ny).zero() - self.resistivity = MultiLocationArray(self.nx, self.ny).zero() def geometry2(self): """ @@ -3783,10 +3772,6 @@ def addFromRegionsXArray(name): addFromRegions("pressure") addFromRegions("Jpar0") - if self.user_options.kinetic_profiles: - addFromRegions("ni0") - addFromRegions("resistivity") - # Penalty mask self.penalty_mask = BoutArray( numpy.zeros((self.nx, self.ny)), From 85741836a6fb660f1a71826d2ceac9536a60ad08 Mon Sep 17 00:00:00 2001 From: HarukiST Date: Thu, 4 Sep 2025 15:43:36 +0900 Subject: [PATCH 58/60] delete options for kinetic profiles --- hypnotoad/cases/tokamak.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index ed0efdd8..912916cd 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -307,16 +307,6 @@ class TokamakEquilibrium(Equilibrium): value_type=int, check_all=is_positive, ), - # - # set dummy data of kinetic profile for post-process - # - kinetic_profiles=WithMeta( - False, - doc=( - "Set dummy data (zero-fill) of kinetic profiles for post-process" - ), - value_type=bool, - ), ) def __init__( From 06a8fcada214de90ce9ff7b161eecc7a8e2f4a9f Mon Sep 17 00:00:00 2001 From: HarukiST Date: Thu, 4 Sep 2025 15:51:01 +0900 Subject: [PATCH 59/60] formatted with black --- hypnotoad/core/mesh.py | 1 - 1 file changed, 1 deletion(-) diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py index d46e2279..8cbf9bbb 100644 --- a/hypnotoad/core/mesh.py +++ b/hypnotoad/core/mesh.py @@ -950,7 +950,6 @@ def geometry1(self): 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 From 2707fb6c1302463848c637753e6259e9c556fa6d Mon Sep 17 00:00:00 2001 From: HarukiST Date: Thu, 4 Sep 2025 15:55:39 +0900 Subject: [PATCH 60/60] formatted with black --- hypnotoad/cases/tokamak.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 912916cd..8f17b5d8 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -858,11 +858,11 @@ def describeSingleNull(self): # Get the smallest absolute grid spacing for the separatrix # dpsidi_pf_minimum is introduced to avoid very narrow radial grid width # around separatrix in core and SOL region and 0.3 is an emprical factor. - - #dpsidi_sep = min([dpsidi_sol, dpsidi_core, dpsidi_pf], key=abs) - dpsidi_sep = min([dpsidi_sol, dpsidi_core], key=abs) # set dx at separatrix - #dpsidi_pf_minimum = max([0.3*dpsidi_sep, 0.3*dpsidi_pf], key=abs) - dpsidi_pf_minimum = max([0.3*dpsidi_sep, dpsidi_pf], key=abs) + + # dpsidi_sep = min([dpsidi_sol, dpsidi_core, dpsidi_pf], key=abs) + dpsidi_sep = min([dpsidi_sol, dpsidi_core], key=abs) # set dx at separatrix + # dpsidi_pf_minimum = max([0.3*dpsidi_sep, 0.3*dpsidi_pf], key=abs) + dpsidi_pf_minimum = max([0.3 * dpsidi_sep, dpsidi_pf], key=abs) # decrease (assuming the factor is <1) the spacing around the separatrix by the # factor psi_spacing_separatrix_multiplier