From 88895de2f5824fec05b2a8c829f1007dfc28983b Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 6 Feb 2023 11:16:26 +0000 Subject: [PATCH] Switch to external FreeQDSK for *qdsk reading/writing --- docs/input_and_output.rst | 4 +- freegs/_aeqdsk.py | 329 -------------------------------------- freegs/_divgeo.py | 2 +- freegs/_fileutils.py | 126 --------------- freegs/_geqdsk.py | 306 ----------------------------------- freegs/geqdsk.py | 7 +- freegs/test_aeqdsk.py | 111 ------------- freegs/test_fileutils.py | 28 ---- freegs/test_geqdsk.py | 49 ------ pyproject.toml | 3 +- test_fileutils.py | 28 ---- 11 files changed, 9 insertions(+), 984 deletions(-) delete mode 100644 freegs/_aeqdsk.py delete mode 100644 freegs/_fileutils.py delete mode 100644 freegs/_geqdsk.py delete mode 100644 freegs/test_aeqdsk.py delete mode 100644 freegs/test_fileutils.py delete mode 100644 freegs/test_geqdsk.py delete mode 100644 test_fileutils.py diff --git a/docs/input_and_output.rst b/docs/input_and_output.rst index d609eb2..01eef4b 100644 --- a/docs/input_and_output.rst +++ b/docs/input_and_output.rst @@ -118,11 +118,11 @@ To convert a G-EQDSK file directly to DivGeo without any solves, and without needing to know the coil locations, use the low-level routines:: - from freegs import _geqdsk + from freegdsk import geqdsk from freegs import _divgeo with open("input.geqdsk") as f: - data = _geqdsk.read(f) + data = geqdsk.read(f) with open("output.equ", "w") as f: _divgeo.write(data, f) diff --git a/freegs/_aeqdsk.py b/freegs/_aeqdsk.py deleted file mode 100644 index cf36cde..0000000 --- a/freegs/_aeqdsk.py +++ /dev/null @@ -1,329 +0,0 @@ -""" - -fields - Lists the variables stored in the file, a default value, and a description - -""" - -from . import _fileutils as fu - -import warnings - -# List of file data variables, default values, and documentation -# This is used in both reader and writer -fields = [ - ( - "tsaisq", - 0.0, - "total chi2 from magnetic probes, flux loops, Rogowski and external coils", - ), - ("rcencm", 100.0, "major radius in cm for vacuum field BCENTR"), - ("bcentr", 1.0, "vacuum toroidal magnetic field in Tesla at RCENCM"), - ("pasmat", 1e6, "measured plasma toroidal current in Ampere"), - ("cpasma", 1e6, "fitted plasma toroidal current in Ampere-turn"), - ("rout", 100.0, "major radius of geometric center in cm"), - ("zout", 0.0, "Z of geometric center in cm"), - ("aout", 50.0, "plasma minor radius in cm"), - ("eout", 1.0, "Plasma boundary elongation"), - ("doutu", 1.0, "upper triangularity"), - ("doutl", 1.0, "lower triangularity"), - ("vout", 1000.0, "plasma volume in cm3"), - ("rcurrt", 100.0, "major radius in cm of current centroid"), - ("zcurrt", 0.0, "Z in cm at current centroid"), - ("qsta", 5.0, "equivalent safety factor q*"), - ("betat", 1.0, "toroidal beta in %"), - ( - "betap", - 1.0, - "poloidal beta with normalization average poloidal magnetic BPOLAV defined through Ampere's law", - ), - ( - "ali", - 0.0, - "li with normalization average poloidal magnetic defined through Ampere's law", - ), - ("oleft", 10.0, "plasma inner gap in cm"), - ("oright", 10.0, "plasma outer gap in cm"), - ("otop", 10.0, "plasma top gap in cm"), - ("obott", 10.0, "plasma bottom gap in cm"), - ("qpsib", 5.0, "q at 95% of poloidal flux"), - ("vertn", 1.0, "vacuum field (index? -- seems to be float) at current centroid"), - # fmt_1040 = r '^\s*' + 4 * r '([\s\-]\d+\.\d+[Ee][\+\-]\d\d)' - # read(neqdsk, 1040)(rco2v(k, jj), k = 1, mco2v) - (None, None, None), # New line - ( - "rco2v", - lambda data: [0.0] * data["mco2v"], - "1D array : path length in cm of vertical CO2 density chord", - ), - # read(neqdsk, 1040)(dco2v(jj, k), k = 1, mco2v) - (None, None, None), # New line - ( - "dco2v", - lambda data: [0.0] * data["mco2v"], - "line average electron density in cm3 from vertical CO2 chord", - ), - # read(neqdsk, 1040)(rco2r(k, jj), k = 1, mco2r) - (None, None, None), # New line - ( - "rco2r", - lambda data: [0.0] * data["mco2r"], - "path length in cm of radial CO2 density chord", - ), - # read(neqdsk, 1040)(dco2r(jj, k), k = 1, mco2r) - (None, None, None), # New line - ( - "dco2r", - lambda data: [0.0] * data["mco2r"], - "line average electron density in cm3 from radial CO2 chord", - ), - (None, None, None), # New line - ("shearb", 0.0, ""), - ( - "bpolav", - 1.0, - "average poloidal magnetic field in Tesla defined through Ampere's law", - ), - ("s1", 0.0, "Shafranov boundary line integrals"), - ("s2", 0.0, "Shafranov boundary line integrals"), - ("s3", 0.0, "Shafranov boundary line integrals"), - ("qout", 0.0, "q at plasma boundary"), - ("olefs", 0.0, ""), - ("orighs", 0.0, "outer gap of external second separatrix in cm"), - ("otops", 0.0, "top gap of external second separatrix in cm"), - ("sibdry", 1.0, ""), - ("areao", 100.0, "cross sectional area in cm2"), - ("wplasm", 0.0, ""), - ("terror", 0.0, "equilibrium convergence error"), - ("elongm", 0.0, "elongation at magnetic axis"), - ("qqmagx", 0.0, "axial safety factor q(0)"), - ("cdflux", 0.0, "computed diamagnetic flux in Volt-sec"), - ("alpha", 0.0, "Shafranov boundary line integral parameter"), - ("rttt", 0.0, "Shafranov boundary line integral parameter"), - ("psiref", 1.0, "reference poloidal flux in VS/rad"), - ( - "xndnt", - 0.0, - "vertical stability parameter, vacuum field index normalized to critical index value", - ), - ("rseps1", 1.0, "major radius of x point in cm"), - ("zseps1", -1.0, ""), - ("rseps2", 1.0, "major radius of x point in cm"), - ("zseps2", 1.0, ""), - ("sepexp", 0.0, "separatrix radial expansion in cm"), - ("obots", 0.0, "bottom gap of external second separatrix in cm"), - ("btaxp", 1.0, "toroidal magnetic field at magnetic axis in Tesla"), - ("btaxv", 1.0, "vacuum toroidal magnetic field at magnetic axis in Tesla"), - ("aaq1", 100.0, "minor radius of q=1 surface in cm, 100 if not found"), - ("aaq2", 100.0, "minor radius of q=2 surface in cm, 100 if not found"), - ("aaq3", 100.0, "minor radius of q=3 surface in cm, 100 if not found"), - ( - "seplim", - 0.0, - "> 0 for minimum gap in cm in divertor configurations, < 0 absolute value for minimum distance to external separatrix in limiter configurations", - ), - ("rmagx", 100.0, "major radius in cm at magnetic axis"), - ("zmagx", 0.0, ""), - ("simagx", 0.0, "Poloidal flux at the magnetic axis"), - ("taumhd", 0.0, "energy confinement time in ms"), - ("betapd", 0.0, "diamagnetic poloidal b"), - ("betatd", 0.0, "diamagnetic toroidal b in %"), - ("wplasmd", 0.0, "diamagnetic plasma stored energy in Joule"), - ("diamag", 0.0, "measured diamagnetic flux in Volt-sec"), - ("vloopt", 0.0, "measured loop voltage in volt"), - ("taudia", 0.0, "diamagnetic energy confinement time in ms"), - ( - "qmerci", - 0.0, - "Mercier stability criterion on axial q(0), q(0) > QMERCI for stability", - ), - ("tavem", 0.0, "average time in ms for magnetic and MSE data"), - # ishot > 91000 - # The next section is dependent on the EFIT version - # New version of EFIT on 05/24/97 writes aeqdsk that includes - # data values for parameters nsilop,magpri,nfcoil and nesum. - (None, True, None), # New line - ( - "nsilop", - lambda data: len(data.get("csilop", [])), - "Number of flux loop signals, len(csilop)", - ), - ( - "magpri", - lambda data: len(data.get("cmpr2", [])), - "Number of flux loop signals, len(cmpr2) (added to nsilop)", - ), - ( - "nfcoil", - lambda data: len(data.get("ccbrsp", [])), - "Number of calculated external coil currents, len(ccbrsp)", - ), - ( - "nesum", - lambda data: len(data.get("eccurt", [])), - "Number of measured E-coil currents", - ), - (None, None, None), # New line - ( - "csilop", - lambda data: [0.0] * data.get("nsilop", 0), - "computed flux loop signals in Weber", - ), - ("cmpr2", lambda data: [0.0] * data.get("magpri", 0), ""), - ( - "ccbrsp", - lambda data: [0.0] * data.get("nfcoil", 0), - "computed external coil currents in Ampere", - ), - ( - "eccurt", - lambda data: [0.0] * data.get("nesum", 0), - "measured E-coil current in Ampere", - ), - ("pbinj", 0.0, "neutral beam injection power in Watts"), - ("rvsin", 0.0, "major radius of vessel inner hit spot in cm"), - ("zvsin", 0.0, "Z of vessel inner hit spot in cm"), - ("rvsout", 0.0, "major radius of vessel outer hit spot in cm"), - ("zvsout", 0.0, "Z of vessel outer hit spot in cm"), - ("vsurfa", 0.0, "plasma surface loop voltage in volt, E EQDSK only"), - ("wpdot", 0.0, "time derivative of plasma stored energy in Watt, E EQDSK only"), - ("wbdot", 0.0, "time derivative of poloidal magnetic energy in Watt, E EQDSK only"), - ("slantu", 0.0, ""), - ("slantl", 0.0, ""), - ("zuperts", 0.0, ""), - ("chipre", 0.0, "total chi2 pressure"), - ("cjor95", 0.0, ""), - ("pp95", 0.0, "normalized P'(y) at 95% normalized poloidal flux"), - ("ssep", 0.0, ""), - ("yyy2", 0.0, "Shafranov Y2 current moment"), - ("xnnc", 0.0, ""), - ("cprof", 0.0, "current profile parametrization parameter"), - ("oring", 0.0, "not used"), - ( - "cjor0", - 0.0, - "normalized flux surface average current density at 99% of normalized poloidal flux", - ), - ("fexpan", 0.0, "flux expansion at x point"), - ("qqmin", 0.0, "minimum safety factor qmin"), - ("chigamt", 0.0, "total chi2 MSE"), - ("ssi01", 0.0, "magnetic shear at 1% of normalized poloidal flux"), - ("fexpvs", 0.0, "flux expansion at outer lower vessel hit spot"), - ( - "sepnose", - 0.0, - "radial distance in cm between x point and external field line at ZNOSE", - ), - ("ssi95", 0.0, "magnetic shear at 95% of normalized poloidal flux"), - ("rqqmin", 0.0, "normalized radius of qmin , square root of normalized volume"), - ("cjor99", 0.0, ""), - ( - "cj1ave", - 0.0, - "normalized average current density in plasma outer 5% normalized poloidal flux region", - ), - ("rmidin", 0.0, "inner major radius in m at Z=0.0"), - ("rmidout", 0.0, "outer major radius in m at Z=0.0"), -] - - -def write(data, fh): - """ - data [dict] - keys are given with documentation in the `fields` list. - Also includes - shot [int] - The shot number - time - in ms - - """ - # First line identification string - # Default to date > 1997 since that format includes nsilop etc. - fh.write("{0:11s}\n".format(data.get("header", " 26-OCT-98 09/07/98 "))) - - # Second line shot number - fh.write(" {:d} 1\n".format(data.get("shot", 0))) - - # Third line time - fh.write(" " + fu.f2s(data.get("time", 0.0)) + "\n") - - # Fourth line - # time(jj),jflag(jj),lflag,limloc(jj), mco2v,mco2r,qmflag - # jflag = 0 if error (? Seems to contradict example) - # lflag > 0 if error (? Seems to contradict example) - # limloc IN/OUT/TOP/BOT: limiter inside/outside/top/bot SNT/SNB: single null top/bottom DN: double null - # mco2v number of vertical CO2 density chords - # mco2r number of radial CO2 density chords - # qmflag axial q(0) flag, FIX if constrained and CLC for float - fh.write( - "*{:s} {:d} {:d} {:s} {:d} {:d} {:s}\n".format( - fu.f2s(data.get("time", 0.0)).strip(), - data.get("jflag", 1), - data.get("lflag", 0), - data.get("limloc", "DN"), - data.get("mco2v", 0), - data.get("mco2r", 0), - data.get("qmflag", "CLC"), - ) - ) - # Output data in lines of 4 values each - with fu.ChunkOutput(fh, chunksize=4) as output: - for key, default, description in fields: - if callable(default): - # Replace the default function with the value, which may depend on previously read data - default = default(data) - - if key is None: - output.newline() # Ensure on a new line - else: - output.write(data.get(key, default)) - - -def read(fh): - """ - Read an AEQDSK file, returning a dictionary of data - """ - # First line label. Date. - header = fh.readline() - - # Second line shot number - shot = int(fh.readline().split()[0]) - - # Third line time [ms] - time = float(fh.readline()) - - # Fourth line has (up to?) 9 entries - # time(jj),jflag(jj),lflag,limloc(jj), mco2v,mco2r,qmflag - words = fh.readline().split() - - # Dictionary to hold result - data = { - "header": header, - "shot": shot, - "time": time, - "jflag": int(words[1]), - "lflag": int(words[2]), - "limloc": words[3], # e.g. "SNB" - "mco2v": int(words[4]), - "mco2r": int(words[5]), - "qmflag": words[6], - } # e.g. "CLC" - - # Read each value from the file, and put into variables - values = fu.next_value(fh) - for key, default, doc in fields: - if key is None: - continue # skip - - if callable(default): - default = default(data) - - if isinstance(default, list): - # Read a list the same length as the default - data[key] = [next(values) for elt in default] - else: - value = next(values) - if isinstance(default, int) and not isinstance(value, int): - # Expecting an integer, but didn't get one - warnings.warn("Expecting an integer for '" + key + "' in aeqdsk file") - break - data[key] = value - - return data diff --git a/freegs/_divgeo.py b/freegs/_divgeo.py index 13c076f..6762c1b 100644 --- a/freegs/_divgeo.py +++ b/freegs/_divgeo.py @@ -21,7 +21,7 @@ """ import numpy as np -from ._fileutils import ChunkOutput, write_1d, write_2d +from freeqdsk._fileutils import ChunkOutput, write_1d, write_2d def write(data, fh, label=None): diff --git a/freegs/_fileutils.py b/freegs/_fileutils.py deleted file mode 100644 index 8410f85..0000000 --- a/freegs/_fileutils.py +++ /dev/null @@ -1,126 +0,0 @@ -""" -Utilities for writing and reading files compatible with Fortran - -""" - -import re - - -def f2s(f): - """ - Format a string containing a float - """ - s = "" - if f >= 0.0: - s += " " - return s + "%1.9E" % f - - -class ChunkOutput: - """ - This outputs values in lines, inserting - newlines when needed. - """ - - def __init__(self, filehandle, chunksize=5, extraspaces=0): - """ - filehandle output to write to - chunksize number of values on a line - extraspaces number of extra spaces between outputs - """ - self.fh = filehandle - self.counter = 0 - self.chunk = chunksize - self.extraspaces = extraspaces - - def write(self, value): - """ " - Write a value to the output, adding a newline if needed - - Distinguishes between: - - list : Iterates over the list and writes each element - - int : Converts using str - - float : Converts using f2s to Fortran-formatted string - - """ - if isinstance(value, list): - for elt in value: - self.write(elt) - return - - self.fh.write(" " * self.extraspaces) - - if isinstance(value, int): - self.fh.write(" " + str(value)) - else: - self.fh.write(f2s(value)) - - self.counter += 1 - if self.counter == self.chunk: - self.fh.write("\n") - self.counter = 0 - - def newline(self): - """ - Ensure that the file is at the start of a new line - """ - if self.counter != 0: - self.fh.write("\n") - self.counter = 0 - - def endblock(self): - """ - Make sure next block of data is on new line - """ - self.fh.write("\n") - self.counter = 0 - - def __enter__(self): - return self - - def __exit__(self, type, value, traceback): - """Ensure that the chunk finishes with a new line""" - if self.counter != 0: - self.counter = 0 - self.fh.write("\n") - - -def write_1d(val, out): - """ - Writes a 1D variable val to the file handle out - """ - for i in range(len(val)): - out.write(val[i]) - out.newline() - - -def write_2d(val, out): - """ - Writes a 2D array. Note that this transposes - the array, looping over the first index fastest - """ - nx, ny = val.shape - for y in range(ny): - for x in range(nx): - out.write(val[x, y]) - out.newline() - - -def next_value(fh): - """ - A generator which yields values from a file handle - - Checks if the value is a float or int, returning - the correct type depending on if '.' is in the string - - """ - pattern = re.compile(r"[ +\-]?\d+(?:\.\d+(?:[Ee][\+\-]\d\d)?)?") - - # Go through each line, extract values, then yield them one by one - for line in fh: - matches = pattern.findall(line) - for match in matches: - if "." in match: - yield float(match) - else: - yield int(match) diff --git a/freegs/_geqdsk.py b/freegs/_geqdsk.py deleted file mode 100644 index 1e91af1..0000000 --- a/freegs/_geqdsk.py +++ /dev/null @@ -1,306 +0,0 @@ -""" -Low level routines for reading and writing G-EQDSK files - -Copyright 2016 Ben Dudson, University of York. Email: benjamin.dudson@york.ac.uk - -This file is part of FreeGS. - -FreeGS is free software: you can redistribute it and/or modify -it under the terms of the GNU Lesser General Public License as published by -the Free Software Foundation, either version 3 of the License, or -(at your option) any later version. - -FreeGS is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU Lesser General Public License for more details. - -You should have received a copy of the GNU Lesser General Public License -along with FreeGS. If not, see . - -""" - -from datetime import date -from numpy import zeros, pi - -from ._fileutils import f2s, ChunkOutput, write_1d, write_2d, next_value - - -def write(data, fh, label=None, shot=None, time=None): - """ - Write a GEQDSK equilibrium file, given a dictionary of data - - data - dictionary - nx, ny Number of points in R (x), Z (y) - rdim, zdim Sizes of the R,Z dimensions - rcentr Reference value of R - bcentr Vacuum toroidal magnetic field at rcentr - rleft R at left (inner) boundary - zmid Z at middle of domain - rmagx, zmagx R,Z at magnetic axis (O-point) - simagx Poloidal flux psi at magnetic axis - sibdry Poloidal flux psi at plasma boundary - cpasma Plasma current [Amps] - - fpol 1D array of f(psi)=R*Bt [meter-Tesla] - pres 1D array of p(psi) [Pascals] - qpsi 1D array of q(psi) - - psi 2D array (nx,ny) of poloidal flux - - fh - file handle - - label - Text label to put in the file - """ - - nx = data["nx"] - ny = data["ny"] - - if not label: - label = "FREEGS" - if len(label) > 11: - label = label[0:12] - print("WARNING: label too long, it will be shortened to {}".format(label)) - - creation_date = date.today().strftime("%d/%m/%Y") - - if not shot: - shot = 0 - - if isinstance(shot, int): - shot = "# {:d}".format(shot) - - if not time: - time = 0 - - if isinstance(time, int): - time = " {:d}ms".format(time) - - # I have no idea what idum is, here it is set to 3 - idum = 3 - header = "{0:11s}{1:10s} {2:>8s}{3:16s}{4:4d}{5:4d}{6:4d}\n".format( - label, creation_date, shot, time, idum, nx, ny - ) - - # First line: Identification string, followed by resolution - fh.write(header) - - # Second line - fh.write( - f2s(data["rdim"]) - + f2s(data["zdim"]) - + f2s(data["rcentr"]) - + f2s(data["rleft"]) - + f2s(data["zmid"]) - + "\n" - ) - - # Third line - fh.write( - f2s(data["rmagx"]) - + f2s(data["zmagx"]) - + f2s(data["simagx"]) - + f2s(data["sibdry"]) - + f2s(data["bcentr"]) - + "\n" - ) - - # 4th line - fh.write( - f2s(data["cpasma"]) - + f2s(data["simagx"]) - + f2s(0.0) - + f2s(data["rmagx"]) - + f2s(0.0) - + "\n" - ) - - # 5th line - fh.write( - 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) - - write_1d(data["fpol"], co) - write_1d(data["pres"], co) - if "ffprime" in data: - write_1d(data["ffprime"], co) - else: - write_1d(workk, co) - if "pprime" in data: - write_1d(data["pprime"], co) - else: - write_1d(workk, co) - - write_2d(data["psi"], co) - write_1d(data["qpsi"], co) - - # Boundary / limiters - nbdry = 0 - nlim = 0 - if "rbdry" in data: - nbdry = len(data["rbdry"]) - if "rlim" in data: - nlim = len(data["rlim"]) - - co.newline() - fh.write("{0:5d}{1:5d}\n".format(nbdry, nlim)) - - if nbdry > 0: - for r, z in zip(data["rbdry"], data["zbdry"]): - co.write(r) - co.write(z) - co.newline() - - if nlim > 0: - for r, z in zip(data["rlim"], data["zlim"]): - co.write(r) - co.write(z) - co.newline() - - -def read(fh, cocos=1): - """ - Read a G-EQDSK formatted equilibrium file - - Format is specified here: - https://fusion.gat.com/theory/Efitgeqdsk - - cocos - COordinate COnventions. Not fully handled yet, - only whether psi is divided by 2pi or not. - if < 10 then psi is divided by 2pi, otherwise not. - - Returns - ------- - - A dictionary containing: - nx, ny Number of points in R (x), Z (y) - rdim, zdim Sizes of the R,Z dimensions - rcentr Reference value of R - bcentr Vacuum toroidal magnetic field at rcentr - rleft R at left (inner) boundary - zmid Z at middle of domain - rmagx, zmagx R,Z at magnetic axis (O-point) - simagx Poloidal flux psi at magnetic axis - sibdry Poloidal flux psi at plasma boundary - cpasma Plasma current [Amps] - - fpol 1D array of f(psi)=R*Bt [meter-Tesla] - pres 1D array of p(psi) [Pascals] - qpsi 1D array of q(psi) - - psi 2D array (nx,ny) of poloidal flux - - """ - - # Read the first line - header = fh.readline() - words = header.split() # Split on whitespace - if len(words) < 3: - raise ValueError("Expecting at least 3 numbers on first line") - - nx = int(words[-2]) - ny = int(words[-1]) - - print(" nx = {0}, ny = {1}".format(nx, ny)) - - # Dictionary to hold result - data = {"nx": nx, "ny": ny} - - # List of fields to read. None means discard value - fields = [ - "rdim", - "zdim", - "rcentr", - "rleft", - "zmid", - "rmagx", - "zmagx", - "simagx", - "sibdry", - "bcentr", - "cpasma", - "simagx", - None, - "rmagx", - None, - "zmagx", - None, - "sibdry", - None, - None, - ] - - values = next_value(fh) - - for f in fields: - val = next(values) - if f: - data[f] = val - - # Read arrays - - def read_1d(n): - """ - Read a 1D array of length n - """ - val = zeros(n) - for i in range(n): - val[i] = next(values) - return val - - def read_2d(n, m): - """ - Read a 2D (n,m) array in Fortran order - """ - val = zeros([n, m]) - for y in range(m): - for x in range(n): - val[x, y] = next(values) - return val - - data["fpol"] = read_1d(nx) - data["pres"] = read_1d(nx) - data["ffprime"] = read_1d(nx) - data["pprime"] = read_1d(nx) - - data["psi"] = read_2d(nx, ny) - - data["qpsi"] = read_1d(nx) - - # Ensure that psi is divided by 2pi - if cocos > 10: - for var in ["psi", "simagx", "sibdry"]: - data[var] /= 2 * pi - - nbdry = next(values) - nlim = next(values) - - # print(nbdry, nlim) - - if nbdry > 0: - # Read (R,Z) pairs - print(nbdry) - data["rbdry"] = zeros(nbdry) - data["zbdry"] = zeros(nbdry) - for i in range(nbdry): - data["rbdry"][i] = next(values) - data["zbdry"][i] = next(values) - - if nlim > 0: - # Read (R,Z) pairs - data["rlim"] = zeros(nlim) - data["zlim"] = zeros(nlim) - for i in range(nlim): - data["rlim"][i] = next(values) - data["zlim"][i] = next(values) - - return data diff --git a/freegs/geqdsk.py b/freegs/geqdsk.py index 142803d..c459eba 100755 --- a/freegs/geqdsk.py +++ b/freegs/geqdsk.py @@ -23,7 +23,6 @@ """ -from . import _geqdsk from . import critical from .equilibrium import Equilibrium from .machine import Wall @@ -32,6 +31,8 @@ from . import picard from .gradshafranov import mu0 + +from freeqdsk import geqdsk from scipy import interpolate from numpy import ( linspace, @@ -46,7 +47,7 @@ from scipy.integrate import romb -def write(eq, fh, label=None, oxpoints=None, fileformat=_geqdsk.write): +def write(eq, fh, label=None, oxpoints=None, fileformat=geqdsk.write): """ Write a GEQDSK equilibrium file, given a FreeGS Equilibrium object @@ -231,7 +232,7 @@ def read( show = True # Read the data as a Dictionary - data = _geqdsk.read(fh, cocos=cocos) + data = geqdsk.read(fh, cocos=cocos) # If data contains a limiter, set the machine wall if "rlim" in data: diff --git a/freegs/test_aeqdsk.py b/freegs/test_aeqdsk.py deleted file mode 100644 index e1982be..0000000 --- a/freegs/test_aeqdsk.py +++ /dev/null @@ -1,111 +0,0 @@ -import numpy -from io import StringIO - -from . import _aeqdsk - - -def test_writeread(): - """ - Test that data can be written then read back - """ - data = { - "shot": 66832, - "time": 2384.0, - "jflag": 1, - "lflag": 0, - "limloc": "SNB", - "mco2v": 3, - "mco2r": 1, - "qmflag": "CLC", - "tsaisq": 11.7513361, - "rcencm": 169.550003, - "bcentr": -2.06767821, - "pasmat": 1213135.38, - "cpasma": 1207042.13, - "rout": 168.491165, - "zout": -6.82398081, - "aout": 63.0725098, - "eout": 1.73637426, - "doutu": 0.160389453, - "doutl": 0.329588085, - "vout": 19912044.0, - "rcurrt": 170.800049, - "zcurrt": 7.52815676, - "qsta": 6.26168156, - "betat": 0.60095495, - "betap": 0.326897353, - "ali": 1.47733176, - "oleft": 3.73984718, - "oright": 4.84749842, - "otop": 32.4465942, - "obott": 20.2485809, - "qpsib": 4.39304399, - "vertn": -0.675418258, - "rco2v": [216.307495, 155.99646, 121.109322], - "dco2v": [27324300900000.0, 29309569900000.0, 29793563200000.0], - "rco2r": [125.545105], - "dco2r": [32812950400000.0], - "shearb": 4.02759838, - "bpolav": 0.282110274, - "s1": 2.155056, - "s2": 1.09512568, - "s3": 0.640428185, - "qout": 7.93196821, - "olefs": -50.0, - "orighs": -50.0, - "otops": -50.0, - "sibdry": -0.016445132, - "areao": 19292.2441, - "wplasm": 309183.625, - "terror": 0.000789557525, - "elongm": 1.18666041, - "qqmagx": 0.620565712, - "cdflux": -0.0285836495, - "alpha": 1.32450712, - "rttt": 155.478485, - "psiref": 0.0, - "xndnt": 0.0, - "rseps1": 147.703217, - "zseps1": -116.341461, - "rseps2": -999.0, - "zseps2": -999.0, - "sepexp": 5.94302845, - "obots": -50.0, - "btaxp": -2.18051338, - "btaxv": -2.03286076, - "aaq1": 30.0145092, - "aaq2": 46.8485107, - "aaq3": 54.2332726, - "seplim": 3.73984718, - "rmagx": 172.453949, - "zmagx": 9.105937, - "simagx": 0.380751818, - "taumhd": 64.3431244, - "betapd": 0.473303556, - "betatd": 0.870102286, - "wplasmd": 447656.406, - "diamag": -2.33697938e-05, - "vloopt": 0.110378414, - "taudia": 93.1602173, - "qmerci": 0.0, - "tavem": 10.0, - } - - output = StringIO() - - # Write to string - _aeqdsk.write(data, output) - - # Move to the beginning of the buffer - output.seek(0) - - # Read from string - data2 = _aeqdsk.read(output) - - # Check that data and data2 are the same - for key in data: - if isinstance(data[key], str): - assert data2[key] == data[key] - else: - # Number, or list of numbers - numpy.testing.assert_allclose(data2[key], data[key]) diff --git a/freegs/test_fileutils.py b/freegs/test_fileutils.py deleted file mode 100644 index f51a2a1..0000000 --- a/freegs/test_fileutils.py +++ /dev/null @@ -1,28 +0,0 @@ -from . import _fileutils - -try: - from io import StringIO -except: - # Python 2 - from StringIO import StringIO - - -def test_f2s(): - assert _fileutils.f2s(0.0) == " 0.000000000E+00" - assert _fileutils.f2s(1234) == " 1.234000000E+03" - assert _fileutils.f2s(-1.65281e12) == "-1.652810000E+12" - assert _fileutils.f2s(-1.65281e-2) == "-1.652810000E-02" - - -def test_ChunkOutput(): - output = StringIO() - co = _fileutils.ChunkOutput(output) - - for val in [1.0, -3.2, 6.2e5, 8.7654e-12, 42.0, -76]: - co.write(val) - - assert ( - output.getvalue() - == """ 1.000000000E+00-3.200000000E+00 6.200000000E+05 8.765400000E-12 4.200000000E+01 - -76""" - ) diff --git a/freegs/test_geqdsk.py b/freegs/test_geqdsk.py deleted file mode 100644 index 9010e8d..0000000 --- a/freegs/test_geqdsk.py +++ /dev/null @@ -1,49 +0,0 @@ -import numpy - -from io import StringIO - -from . import _geqdsk - - -def test_writeread(): - """ - Test that data can be written then read back - """ - nx = 65 - ny = 65 - - # Create a dummy dataset - data = { - "nx": nx, - "ny": ny, - "rdim": 2.0, - "zdim": 1.5, - "rcentr": 1.2, - "bcentr": 2.42, - "rleft": 0.5, - "zmid": 0.1, - "rmagx": 1.1, - "zmagx": 0.2, - "simagx": -2.3, - "sibdry": 0.21, - "cpasma": 1234521, - "fpol": numpy.random.rand(nx), - "pres": numpy.random.rand(nx), - "qpsi": numpy.random.rand(nx), - "psi": numpy.random.rand(nx, ny), - } - - output = StringIO() - - # Write to string - _geqdsk.write(data, output) - - # Move to the beginning of the buffer - output.seek(0) - - # Read from string - data2 = _geqdsk.read(output) - - # Check that data and data2 are the same - for key in data: - numpy.testing.assert_allclose(data2[key], data[key]) diff --git a/pyproject.toml b/pyproject.toml index 376e331..1f29448 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,7 +30,8 @@ dependencies = [ "matplotlib>=1.3", "h5py>=2.10.0", "Shapely>=1.7.1", - "importlib-metadata<4.3,>=1.1.0" + "importlib-metadata<4.3,>=1.1.0", + "freeqdsk>=0.1.0", ] dynamic = ["version"] diff --git a/test_fileutils.py b/test_fileutils.py deleted file mode 100644 index 9a35c18..0000000 --- a/test_fileutils.py +++ /dev/null @@ -1,28 +0,0 @@ -from freegs._fileutils import next_value -from numpy import allclose - - -def test_next_value(): - data = [ - '+9.500000000e-01+2.000000000e+00+1.000000000e+00+2.500000000e+01+0.000000000e+00\n', - ' 9.500000000E-01 2.000000000E+00 1.000000000E+00 2.500000000E+01 0.000000000E+00\n', - '+3.563359524e+02+2.337058846e-02-1.212203732e-02+1.953790839e-03-7.116987284e-03', - ' 3.563359524E+02 2.337058846e-02-1.212203732E-02 1.953790839E-03-7.116987284e-03', - '0 0 0\n', - ' 0 0\n', - ] - - expected = [ - 0.95, 2.0, 1.0, 25.0, 0.0, - 0.95, 2.0, 1.0, 25.0, 0.0, - 356.3359524, 0.02337058846, -0.01212203732, 0.001953790839, -0.007116987284, - 356.3359524, 0.02337058846, -0.01212203732, 0.001953790839, -0.007116987284, - 0, 0, 0, - 0, 0 - ] - - values = next_value(data) - - actual = [val for val in values] - - assert allclose(expected, actual)