From d235648ee1ce05a3721ff5b17c692d6a83ba1fd4 Mon Sep 17 00:00:00 2001 From: wearysebas Date: Tue, 3 Mar 2026 12:23:30 +0000 Subject: [PATCH 1/3] Added the working gridue to bout converter. --- .../gridue_to_bout_converter/README.md | 57 ++ .../gridue_to_bout.py | 946 ++++++++++++++++++ 2 files changed, 1003 insertions(+) create mode 100644 src/boutdata/gridue_to_bout_converter/README.md create mode 100644 src/boutdata/gridue_to_bout_converter/gridue_to_bout.py diff --git a/src/boutdata/gridue_to_bout_converter/README.md b/src/boutdata/gridue_to_bout_converter/README.md new file mode 100644 index 000000000..155ab5543 --- /dev/null +++ b/src/boutdata/gridue_to_bout_converter/README.md @@ -0,0 +1,57 @@ +This script was modified from B. Dudson's original gridue_to_bout.py script. It creates a BOUT grid file from a GRIDUE file. + +It can be ran from the terminal. + +To run it directly after the creation of a gridue file on INGRID the following code needs to be added to ingrid.py: + +```python + def ExportBOUTgrid(self, gridue_file, bout_grid_name: str = 'bout_from_in.grd.nc', plotting = True, verbose = True, ignore_checks = False): + """ + Export a BOUT grid file for the created grid. + + Parameters + ---------- + gridue_file : str, mandatory + Name of gridue file to convert to BOUT grid. + + bout_grid_name : str, optional + Name of BOUT grid file to save. + + plotting : bool, optional + If True, plot the gridue file before conversion. + + verbose : bool, optional + If True, print verbose output during conversion. + + ignore_checks : bool, optional + If True, ignore checks for gridue file format and structure. + + """ + bout_grid_name = gridue_file + "_" + bout_grid_name + Convert_grids(gridue_file,bout_grid_name,plotting, verbose, ignore_checks) + +``` + +And modifying the `ExportGridue` function to: + +```python + if type(self.CurrentTopology) in [SNL]: + if self.WriteGridueSNL(self.CurrentTopology.gridue_settings, fname): + print(f"# Successfully saved gridue file as '{fname}'") + self.ExportBOUTgrid(fname, 'bout_from_in.grd.nc', plotting = True, verbose = True, ignore_checks = False) + elif type(self.CurrentTopology) in [SF15, SF45, SF75, SF105, SF135, SF165, UDN, CDN]: + if self.WriteGridueDNL(self.CurrentTopology.gridue_settings, fname): + print(f"# Successfully saved gridue file as '{fname}'") + self.ExportBOUTgrid(fname, 'bout_from_in.grd.nc', plotting = True, verbose = True, ignore_checks = False) +``` + +*Note:* The created file will have +`dimensions: + x ; + y ; + z ; + x2 ; + y2 ; + t = UNLIMITED ; //` + +Where x has guard cells following INGRID's convention, and y doesn't have any. The converter takes the BOUT coordinates to be x -> y2 (and removes the guard cells) and y -> x2 (and adds the guard cells). So BOUT ends up using x2 as X and y2 as Y coordinates. \ No newline at end of file diff --git a/src/boutdata/gridue_to_bout_converter/gridue_to_bout.py b/src/boutdata/gridue_to_bout_converter/gridue_to_bout.py new file mode 100644 index 000000000..334959c33 --- /dev/null +++ b/src/boutdata/gridue_to_bout_converter/gridue_to_bout.py @@ -0,0 +1,946 @@ +#!/usr/bin/env python3 +# +# Convert UEDGE grid files (gridue) to BOUT++ grids +# +# Parts adapted from INGRID https://github.com/LLNL/INGRID +# Copyright (c) 2020, Lawrence Livermore National Security, LLC +# Parts adapted from Hypnotoad https://github.com/boutproject/hypnotoad/ +# Copyright 2019 J.T. Omotani +# +# +# Modified and finished by S. Ruiz, 2025 from the orinigal gridue_to_bout.py made by B. Dudson. +# +# + +import numpy as np +import matplotlib +import matplotlib.pyplot as plt + +from scipy import linalg + + +def _importBody(gridue_settings, f): + """ + Imports the body of a gridue file and returns a dictionary with the data. + + Parameters: + gridue_settings : dict + Dictionary with header information from the gridue file. + + f : file object + Opened file object to read the body of the gridue file. + + Returns: + gridue_settings : dict + A dictionary containing the gridue body data, with keys corresponding to + the items in BodyItems and values as numpy arrays. + """ + next(f) + BodyItems = ["rm", "zm", "psi", "br", "bz", "bpol", "bphi", "b"] + Str = {i: [] for i in BodyItems} + k = iter(Str.keys()) + Key = next(k) + + for line in f: + if line == "iogridue\n": + continue + if line == "\n": + try: + Key = next(k) + except: + continue + else: + Str[Key].append(line) + f.close() + nx = gridue_settings["nxm"] + 2 + ny = gridue_settings["nym"] + 2 + + for k, v in Str.items(): + L = ("".join(v).replace("\n", "").replace("D", "e")).split() + _l = iter(L) + vv = next(_l) + + data_ = np.zeros((nx, ny, 5)) + for n in range(5): + for j in range(ny): + for i in range(nx): + + data_[i][j][n] = float(vv) + + try: + vv = next(_l) + except: + continue + gridue_settings[k] = data_ + return gridue_settings + + +def _importSN(values, f): + """ + Import a single null file. Only used for Single null divertor topology. + """ + + HeaderItems = ["nxm", "nym", "ixpt1", "ixpt2", "iyseparatrix1"] + gridue_settings = dict(zip(HeaderItems, values)) + + # Add indices that are present in DN but not SN + gridue_settings.update( + { + "iyseparatrix2": gridue_settings["nym"] + 10, # Outside grid + "ix_cut1": gridue_settings["ixpt1"], + "ix_cut2": gridue_settings["nxm"] // 2, + "ix_inner": gridue_settings["nxm"] // 2, + "ix_cut3": gridue_settings["nxm"] // 2, + "ix_cut4": gridue_settings["ixpt2"], + } + ) + + return _importBody(gridue_settings, f) + + +def _importDN(values, f): + """ + Import a double null file. Used for most divertor topologies. + """ + + gridue_settings = dict(zip(["nxm", "nym"], values)) + + header_rows = [ + ["iyseparatrix1", "iyseparatrix2"], + ["ix_plate1", "ix_cut1", "_FILLER_", "ix_cut2", "ix_plate2"], + ["iyseparatrix3", "iyseparatrix4"], + ["ix_plate3", "ix_cut3", "_FILLER_", "ix_cut4", "ix_plate4"], + ] + + for row in header_rows: + values = [int(x) for x in next(f).split()] + if len(values) != len(row): + raise ValueError( + "Expected row with {} integers, found {}".format(len(row), len(values)) + ) + gridue_settings.update(zip(row, values)) + + gridue_settings.update({"ix_inner": gridue_settings["ix_plate2"]}) + + return _importBody(gridue_settings, f) + + +def importGridue(fname: str = "gridue") -> dict: + """ + Import UEDGE grid file as dictionary. + + Parameters + ---------- + fname : str, optional + Path/file name to gridue formatted file. + + Returns + ------- + A dict containing header and body information from the gridue file. + + """ + f = open(fname, mode="r") + values = [int(x) for x in next(f).split()] + + if len(values) == 5: + return _importSN(values, f) + if len(values) == 2: + return _importDN(values, f) + + raise ValueError("Unrecognised gridue format") + + +def plot(GridueParams: dict, edgecolor="black", ax: object = None, show=True): + """ + Plot UEDGE grid from 'dict' obtained from method 'ImportGridue' + + Parameters + ---------- + GridueParams : dict + Gridue header and body information as a dictionary. + + edgecolor : str, optional + Color of grid. + + ax : object, optional + Matplotlib axes to plot on. + + """ + r = GridueParams["rm"] + z = GridueParams["zm"] + Nx = len(r) + Ny = len(r[0]) + patches = [] + plt.figure(figsize=(6, 10)) + if ax is None: + ax = plt.gca() + idx = [np.array([1, 2, 4, 3, 1])] + for i in range(Nx): + for j in range(Ny): + p = matplotlib.patches.Polygon( + np.concatenate((r[i][j][idx], z[i][j][idx])).reshape(2, 5).T, + fill=False, + closed=True, + edgecolor=edgecolor, + ) + ax.add_patch(p) # draw the contours + ax.set_aspect("equal", adjustable="box") + ax.set_xlabel("R") + ax.set_ylabel("Z") + ax.set_ylim(z.min(), z.max()) + ax.set_xlim(r.min(), r.max()) + if show: + plt.show() + return ax + + +def calcHy( nx, ny, g: dict, dy: float): + """ + Calculate poloidal arc length metric from gridue dictionary + """ + rm = g["rm"] + zm = g["zm"] + + hy = np.zeros((nx, ny)) + for i in range(nx): + for j in range(ny): + r = rm[j, i, :] + z = zm[j, i, :] + # Find intersection with (1) -- (3) + R1 = 0.5 * (r[1] + r[3]) + Z1 = 0.5 * (z[1] + z[3]) + # Find intersection with (2) -- (4) + R2 = 0.5 * (r[2] + r[4]) + Z2 = 0.5 * (z[2] + z[4]) + + # Distance from one side to the other + dl = np.sqrt((R2 - R1) ** 2 + (Z2 - Z1) ** 2) + hy[i, j] = dl / dy + return hy + + +def calcGridAngle(g: dict): + Brxy = g["br"][:, :, 0].T + Bzxy = g["bz"][:, :, 0].T + Bpxy = g["bpol"][:, :, 0].T + + psi = g["psi"] + rm = g["rm"] + zm = g["zm"] + + delta_y = [Brxy / Bpxy, Bzxy / Bpxy] # Unit vector along e_y + + R_xlow = 0.5 * (rm[:, :, 1] + rm[:, :, 2]).T + Z_xlow = 0.5 * (zm[:, :, 1] + zm[:, :, 2]).T + R_xhigh = 0.5 * (rm[:, :, 3] + rm[:, :, 4]).T + Z_xhigh = 0.5 * (zm[:, :, 3] + zm[:, :, 4]).T + dR = R_xhigh - R_xlow + dZ = Z_xhigh - Z_xlow + dl = np.sqrt(dR**2 + dZ**2) + delta_x = [dR / dl, dZ / dl] # unit vector along e_x + + # Calculate angle between x and y unit vectors. + # sin(beta) = cos(pi/2 - beta) = e_x_hat.e_y_hat = delta_x.delta_y + sinBeta = delta_x[0] * delta_y[0] + delta_x[1] * delta_y[1] + + # Rotate delta_y by 90 degrees anticlockwise to get unit vector in psi direction + delta_psi = [-delta_y[1], delta_y[0]] + + # cosBeta = delta_x.delta_psi + cosBeta = delta_x[0] * delta_psi[0] + delta_x[1] * delta_psi[1] + + return sinBeta, cosBeta + + +def calcRZderivs(R, Z, var): + """ + Calculate derivatives of var w.r.t R and Z, using 5 points + + Inputs R, Z and var should be 3D [poloidal, radial, 5] + as stored in gridue files. + + Returns (dvar/dR, dvar/dZ) as 2D arrays [poloidal, radial] + at the location of point index 0 + """ + npol, nr, _ = R.shape + dR = np.zeros((npol, nr)) + dZ = np.zeros((npol, nr)) + for i in range(npol): + for j in range(nr): + A = np.zeros((5, 5)) + for p in range(5): + # Constraint on derivatives for this point + A[p, 0] = 1.0 # Constant + A[p, 1] = R[i, j, p] - R[i, j, 0] # dR + A[p, 2] = Z[i, j, p] - Z[i, j, 0] # dZ + A[p, 3] = (R[i, j, p] - R[i, j, 0]) ** 2 # dR^2 + A[p, 4] = (Z[i, j, p] - Z[i, j, 0]) ** 2 # dR^2 + # Values of the function being fitted + b = var[i, j, :].squeeze() + # Invert using LU for stability + lu, piv = linalg.lu_factor(A) + x = linalg.lu_solve((lu, piv), b) + dR[i, j] = x[1] + dZ[i, j] = x[2] + return dR, dZ # dvar/dR, dvar/dZ at cell center + + +def calcMetric(grd: dict, bpsign, verbose=False, ignore_checks=False): + """ + Calculate metric tensor given BOUT++ cell centers at each grid point. + + Parameters: + grd : dict + Dictionary containing grid data with keys: + "Rxy", "Zxy", "Brxy", "Bzxy", "Btxy", "Bpxy", "Bxy", "hy", + "cosBeta", "tanBeta", "curl_bOverB_Rhat", "curl_bOverB_Zhat", + "curl_bOverB_zetahat". + bpsign : int + Sign of the magnetic field (1 for normal, -1 for reversed). + verbose : bool, optional + If True, print detailed information about the metric tensor and Jacobian. + ignore_checks : bool, optional + If True, ignore checks on the Jacobian's relative error. + + Returns: + dict + A dictionary containing the metric tensor components, Jacobian, and other related quantities for each grid point. + --> Each component is a 2D array with shape [radial, poloidal]. + """ + + Rxy = grd["Rxy"] + Zxy = grd["Zxy"] + Brxy = grd["Brxy"] + Bzxy = grd["Bzxy"] + Btxy = grd["Btxy"] + Bpxy = grd["Bpxy"] + Bxy = grd["Bxy"] + hy = grd["hy"] + cosBeta = grd["cosBeta"] + tanBeta = grd["tanBeta"] + curl_bOverB_Rhat = grd["curl_bOverB_Rhat"] + curl_bOverB_Zhat = grd["curl_bOverB_Zhat"] + curl_bOverB_zetahat = grd["curl_bOverB_zetahat"] + + if verbose: + for var, name in [ + (cosBeta, "cos(beta)"), + (tanBeta, "tan(beta)"), + ]: + print( + "{} min {}, mean {}, max {}".format( + name, np.amin(var), np.mean(var), np.amax(var) + ) + ) + + dphidy = hy * Btxy / (Bpxy * Rxy) + + I = np.zeros(Rxy.shape) + + g11 = (Rxy * Bpxy) ** 2 + g22 = 1.0 / (hy * cosBeta) ** 2 + g33 = ( + 1.0 / Rxy**2 + + (Rxy * Bpxy * I) ** 2 + + (dphidy / (hy * cosBeta)) ** 2 + + 2.0 * Rxy * Bpxy * I * dphidy * tanBeta / hy + ) + g12 = Rxy * np.abs(Bpxy) * tanBeta / hy + g13 = -Rxy * Bpxy * dphidy * tanBeta / hy - I * (Rxy * Bpxy) ** 2 + g23 = -bpsign * dphidy / (hy * cosBeta) ** 2 - Rxy * np.abs(Bpxy) * I * tanBeta / hy + + J = hy / Bpxy + + g_11 = 1.0 / (Rxy * Bpxy * cosBeta) ** 2 + (I * Rxy) ** 2 + g_22 = hy**2 + (dphidy * Rxy) ** 2 + g_33 = Rxy**2 + g_12 = bpsign * I * dphidy * Rxy**2 - hy * tanBeta / (Rxy * np.abs(Bpxy)) + g_13 = I * Rxy**2 + g_23 = bpsign * dphidy * Rxy**2 + + Jcheck = ( + bpsign + * 1.0 + / np.sqrt( + g11 * g22 * g33 + + 2.0 * g12 * g13 * g23 + - g11 * g23**2 + - g22 * g13**2 + - g33 * g12**2 + ) + ) + + rel_error = (J - Jcheck) / J + + if verbose: + for var, name in [ + (hy, "hy"), + (J, "J"), + (Jcheck, "Jcheck"), + (J - Jcheck, "J - Jcheck"), + (rel_error, "(J - Jcheck)/J"), + ]: + print( + "{} min {}, mean {}, max {}".format( + name, np.amin(var), np.mean(var), np.amax(var) + ) + ) + + if np.max(np.abs(rel_error)) > 1e-6: + if ignore_checks: + print("WARNING: Relative error in Jacobian too large.") + #else: + # raise ValueError(f"Relative error in Jacobian too large: {np.max(np.abs(rel_error))}") + + # We want to output contravariant components of Curl(b/B) in the + # locally field-aligned coordinate system. + # The contravariant components of an arbitrary vector A are + # A^x = A.Grad(x) + # A^y = A.Grad(y) + # A^z = A.Grad(z) + + # Grad in cylindrical coordinates is + # Grad(f) = df/dR Rhat + 1/R df/dzeta zetahat + df/dZ Zhat + # https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates, + + # x = psi - psi_min + # dpsi/dR = -R*BZ + # dpsi/dZ = R*BR + # => Grad(x) = (dpsi/dR, 0, dpsi/dZ).(Rhat, zetahat, Zhat) + # => Grad(x) = (-R BZ, 0, R BR).(Rhat, zetahat, Zhat) + curl_bOverB_x = -Rxy * Bzxy * curl_bOverB_Rhat + Rxy * Brxy * curl_bOverB_Zhat + + # Grad(y) = (d_Z, 0, -d_R)/(hy*cosBeta) + # = (BR*cosBeta-BZ*sinBeta, 0, BZ*cosBeta+BR*sinBeta) + # /(Bp*hy*cosBeta) + # = (BR-BZ*tanBeta, 0, BZ+BR*tanBeta)/(Bp*hy) + curl_bOverB_y = ( + curl_bOverB_Rhat * (Brxy - Bzxy * tanBeta) + + curl_bOverB_Zhat * (Bzxy + Brxy * tanBeta) + ) / (Bpxy * hy) + + # Grad(z) = Grad(zeta) - Bt*hy/(Bp*R)*Grad(y) - I*Grad(x) + # Grad(z) = (0, 1/R, 0) - Bt*hy/(Bp*R)*Grad(y) - I*Grad(x) + curl_bOverB_z = ( + curl_bOverB_zetahat / Rxy + - Btxy * hy / (Bpxy * Rxy) * curl_bOverB_y + - I * curl_bOverB_x + ) + + bxcvx = Bxy / 2.0 * curl_bOverB_x + bxcvy = Bxy / 2.0 * curl_bOverB_y + bxcvz = Bxy / 2.0 * curl_bOverB_z + + if verbose: + for var, name in [ + (bxcvx, "bxcvx"), + (bxcvy, "bxcvy"), + (bxcvz, "bxcvz"), + ]: + print( + "{} min {}, mean {}, max {}".format( + name, np.amin(var), np.mean(var), np.amax(var) + ) + ) + + return { + "dphidy": dphidy, + # Metric tensor + "g11": g11, + "g22": g22, + "g33": g33, + "g12": g12, + "g13": g13, + "g23": g23, + # Inverse metric tensor + "g_11": g_11, + "g_22": g_22, + "g_33": g_33, + "g_12": g_12, + "g_13": g_13, + "g_23": g_23, + # Jacobian + "J": J, + # Integrated shear + "sinty": I, + # Curvature + "curl_bOverB_x": curl_bOverB_x, + "curl_bOverB_y": curl_bOverB_y, + "curl_bOverB_z": curl_bOverB_z, + "bxcvx": bxcvx, + "bxcvy": bxcvy, + "bxcvz": bxcvz, + } + + +def calcRZCurvature(g: dict): + """ + Calculate curvature in (R, Z, zeta) + Returns 2D arrays [radial, poloidal] + """ + # Variables in UEDGE format [poloidal, radial, 5] + BR = g["br"] + BZ = g["bz"] + Bzeta = g["bphi"] + B2 = g["b"] ** 2 + R = g["rm"] + Z = g["zm"] + + dBzetadR, dBzetadZ = calcRZderivs(R, Z, Bzeta) + dBRdR, dBRdZ = calcRZderivs(R, Z, BR) + dBZdR, dBZdZ = calcRZderivs(R, Z, BZ) + dB2dR, dB2dZ = calcRZderivs(R, Z, B2) + + # Select point at centre of cell + BR = BR[:, :, 0] + BZ = BZ[:, :, 0] + Bzeta = Bzeta[:, :, 0] + B2 = B2[:, :, 0] + R = R[:, :, 0] + Z = Z[:, :, 0] + + # In cylindrical coords + # curl(A) = (1/R*d(AZ)/dzeta - d(Azeta)/dZ) * Rhat + # + 1/R*(d(R Azeta)/dR - d(AR)/dzeta) * Zhat + # + (d(AR)/dZ - d(AZ)/dR) * zetahat + # Where AR, AZ and Azeta are the components on a basis of unit vectors, + # i.e. AR = A.Rhat; AZ = A.Zhat; Azeta = A.zetahat + # https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates, + # + # curl(b/B) = curl((BR/B2), (BZ/B2), (Bzeta/B2)) + # curl(b/B)_Rhat = 1/R d(BZ/B2)/dzeta - d(Bzeta/B2)/dZ + # = 1/(R*B2)*d(BZ)/dzeta - BZ/(R*B4)*d(B2)/dzeta + # - 1/B2*d(Bzeta)/dZ + Bzeta/B4*d(B2)/dZ + # = -1/B2*d(Bzeta)/dZ + Bzeta/B4*d(B2)/dZ + # curl(b/B)_Zhat = 1/R * (d(R Bzeta/B2)/dR - d(BR/B2)/dzeta) + # = Bzeta/(R*B2) + 1/B2*d(Bzeta)/dR - Bzeta/B4*d(B2)/dR + # - 1/(R*B2)*d(BR)/dzeta + BR/(R*B4)*d(B2)/dzeta + # = Bzeta/(R*B2) + 1/B2*d(Bzeta)/dR - Bzeta/B4*d(B2)/dR + # curl(b/B)_zetahat = d(BR/B2)/dZ - d(BZ/B2)/dR + # = 1/B2*d(BR)/dZ - BR/B4*d(B2)/dZ + # - 1/B2*d(BZ)/dR + BZ/B4*d(B2)/dR + # remembering d/dzeta=0 for axisymmetric equilibrium + + curl_bOverB_Rhat = -dBzetadZ / B2 + Bzeta / B2**2 * dB2dZ + curl_bOverB_Zhat = Bzeta / (R * B2) + dBzetadR / B2 - Bzeta / B2**2 * dB2dR + curl_bOverB_zetahat = ( + dBRdZ / B2 - BR / B2**2 * dB2dZ - dBZdR / B2 + BZ / B2**2 * dB2dR + ) + + # Return as [radial, poloidal] + return curl_bOverB_Rhat.T, curl_bOverB_Zhat.T, curl_bOverB_zetahat.T + +def main(): + import argparse + + parser = argparse.ArgumentParser( + description="""Converts UEDGE grid files (gridue) into BOUT++ grids. +Note that in most cases these grids are non-orthogonal.""" + ) + parser.add_argument("gridue_file", type=str) + parser.add_argument("-o", "--output", default="bout.grd.nc") + parser.add_argument("-p", "--plot", action="store_true", default=False) + parser.add_argument("-v", "--verbose", action="store_true", default=False) + parser.add_argument("-i", "--ignore-checks", action="store_true", default=False) + try: + import argcomplete + + argcomplete.autocomplete(parser) + except ImportError: + pass + + args = parser.parse_args() + gridue_file = args.gridue_file + output_filename = args.output + plotting = args.plot + verbose = args.verbose + ignore_checks = args.ignore_checks + + Convert_grids(gridue_file, output_filename, plotting, verbose, ignore_checks) + +def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False, verbose: bool = False, ignore_checks: bool = False): + """ + Convert UEDGE grid file to BOUT++ grid format. + + Parameters: + gridue_file : str + Path to the UEDGE grid file (gridue format). + output_filename : str + Path to save the converted BOUT++ grid file. + plotting : bool, optional + If True, plot the grid after conversion. + verbose : bool, optional + If True, print detailed information about the grid and conversion process. + ignore_checks : bool, optional + If True, ignore checks on the Jacobian's relative error. + Raises: + ------- + ValueError + If the magnetic field direction and poloidal magnetic field have opposite signs. + ValueError + If the relative error in the Jacobian is too large and ignore_checks is False. + ValueError + If the gridue file format is unrecognized. + ValueError + If the gridue file does not contain the expected header information. + ValueError + If the gridue file does not contain the expected body information. + ValueError + If the gridue file does not contain the expected single null or double null topology. + ValueError + If the gridue file does not contain the expected number of grid points. + ValueError + If the gridue file does not contain the expected number of poloidal or radial points. + ValueError + If the gridue file does not contain the expected number of boundary points. + ValueError + If the gridue file does not contain the expected number of inner or outer points. + ValueError + If the gridue file does not contain the expected number of inner or outer boundary points. + ValueError + If the gridue file does not contain the expected number of inner or outer plate points. + """ + + g = importGridue(gridue_file) + + if plotting: + plot(g) + + psi = g["psi"] + rm = g["rm"] + zm = g["zm"] + + Rxy = rm[:, :, 0].T + Zxy = zm[:, :, 0].T + Brxy = g["br"][:, :, 0].T + Bzxy = g["bz"][:, :, 0].T + Bpxy = g["bpol"][:, :, 0].T + Btxy = g["bphi"][:, :, 0].T + Bxy = g["b"][:, :, 0].T + psixy = psi[:, :, 0].T + nx, ny = Rxy.shape + + #idx = [np.array([1, 2, 4, 3, 1])] + + #pol = [] + #for i in range(nx): + # for j in range(ny): + # np.concatenate((rm[i][j][idx], zm[i][j][idx])).reshape(2, 5).T + + # Ordering + # (1) -- (3) + # | | + # | (0) | -> Radial, BOUT++ "x" + # | | + # (2) -- (4) + + # calculate change in psi across cell -> dx + + dx = np.zeros((nx, ny)) + for i in range(nx): + for j in range(ny): + if i > 1 and i < nx-2: + dx[i, j] = 0.5*(psi[j, i + 1, 0] - psi[j, i - 1, 0]) + else: + dx[i, j] = 0.5 * (psi[j, i, 3] + psi[j, i, 4] - psi[j, i, 1] - psi[j, i, 2]) + + # Note: UEDGE grids have narrow cells on the radial + # boundaries. BOUT++ applies boundary conditions half-way between + # cells, and usually expects two boundary cells. + + # Calculate direction of magnetic field Bp dot Grad(y) + Bp_dot_grady = Brxy[1, 1] * (Rxy[1, 1] - Rxy[1, 0]) + Bzxy[1, 1] * ( + Zxy[1, 1] - Zxy[1, 0] + ) + + if Bp_dot_grady * Bpxy[1, 1] < 0.0: + raise ValueError("Bp_dot_grady and Bpxy have opposite signs") + bpsign = np.sign(Bpxy[1, 1]) # Sign of the poloidal magnetic field + + # Choose angle dy across poloidal cell. This is somewhat arbitrary, + # but it is helpful if the y angle changes by 2pi for each poloidal transit of the core + + dy = 2 * np.pi / ny + + # Calculate hy, the arc length along the flux surface passing through + # the center of each cell. + hy = calcHy(nx,ny,g, dy) + + # Calculate angle between x and y coordinates. sinBeta = 0, cosBeta = 1 for an orthogonal mesh + sinBeta, cosBeta = calcGridAngle(g) + tanBeta = sinBeta / cosBeta + + # Calculate curvature + curl_bOverB_Rhat, curl_bOverB_Zhat, curl_bOverB_zetahat = calcRZCurvature(g) + + # Collect 2D variables for output + grd = { + "Rxy": Rxy, + "Zxy": Zxy, + "psixy": psixy, + "dx": dx, + "dy": np.full((nx, ny), dy), + # Magnetic field components + "Brxy": Brxy, + "Bzxy": Bzxy, + "Bpxy": Bpxy, + "Btxy": Btxy, + "Bxy": Bxy, + # Poloidal arc length + "hy": hy, + "hthe": hy, + # Curvature + "curl_bOverB_Rhat": curl_bOverB_Rhat, + "curl_bOverB_Zhat": curl_bOverB_Zhat, + "curl_bOverB_zetahat": curl_bOverB_zetahat, + # Grid angles + "cosBeta": cosBeta, + "tanBeta": tanBeta, + } + + # Remove Y (poloidal) boundary cells + for name in grd: + grd[name] = grd[name][:, 1:-1] + + Rxy = grd["Rxy"] + Zxy = grd["Zxy"] + nx, ny = Rxy.shape + + #Get Mesh Topology and remove guard cells accordingly + mesh_topology = getMeshTopology(g, nx, ny) + + if mesh_topology == "SF": + #SF case + ixseps1 = g["iyseparatrix1"] + 2 # Main X-point separatrix + ixseps2 = min(g["iyseparatrix3"] + 2, nx) # Secondary X-point separatrix + # Remove guard cells on either side of X-point. + ny_inner = g["ix_inner"] + for name in grd: + var = grd[name] + nx, ny = var.shape + newvar = np.zeros((nx, ny - 2)) + newvar[:, :ny_inner] = var[:, :ny_inner] + newvar[:, ny_inner:] = var[:, (ny_inner + 2) :] + grd[name] = newvar + g["ix_cut2"] = g["ix_cut2"] - 1 + g["ix_cut3"] = g["ix_cut3"] - 1 + g["ix_cut4"] = g["ix_cut4"] - 2 + else: + ixseps1 = g["iyseparatrix1"] + 2 # Lower X-point separatrix + ixseps2 = min(g["iyseparatrix2"] + 2, nx) # Upper X-point separatrix + # Double null -> Remove upper Y guard cells + ny_inner = g["ix_inner"] + for name in grd: + var = grd[name] + nx, ny = var.shape + newvar = np.zeros((nx, ny - 2)) + newvar[:, :ny_inner] = var[:, :ny_inner] + newvar[:, ny_inner:] = var[:, (ny_inner + 2) :] + grd[name] = newvar + g["ix_cut2"] = g["ix_cut2"] - 1 + g["ix_cut3"] = g["ix_cut3"] - 3 + g["ix_cut4"] = g["ix_cut4"] - 2 + + # Extrapolate X (radial) boundary cells + # Removing one cell, adding two on each X boundary + for name in grd: + var = grd[name] + nx, ny = var.shape + newvar = np.zeros((nx + 2, ny)) + newvar[2:-2, :] = var[1:-1, :] + if name in ["Rxy", "Zxy", "psixy"]: + # Linear extrapolation + newvar[1, :] = 2.0 * newvar[2, :] - newvar[3, :] + newvar[0, :] = 2.0 * newvar[1, :] - newvar[2, :] + newvar[-2, :] = 2.0 * newvar[-3, :] - newvar[-4, :] + newvar[-1, :] = 2.0 * newvar[-2, :] - newvar[-3, :] + else: + # Constant extrapolation + newvar[1, :] = newvar[2, :] + newvar[0, :] = newvar[2, :] + newvar[-2, :] = newvar[-3, :] + newvar[-1, :] = newvar[-3, :] + grd[name] = newvar + + # Re assign grid indices after removing guard cells + jyseps1_1 = g["ix_cut1"] - 1 + jyseps2_1 = g["ix_cut2"] + ny_inner = g["ix_inner"] + jyseps1_2 = g["ix_cut3"] + jyseps2_2 = g["ix_cut4"] - 1 + + # Calculate metric tensor + grd.update(calcMetric(grd, bpsign, verbose, ignore_checks)) + + dphidy = grd["dphidy"] + Rxy = grd["Rxy"] + Zxy = grd["Zxy"] + nx, ny = Rxy.shape + + # Calculate zShift and ShiftAngle + zShift = np.zeros((nx, ny)) + + # Inner core region + zShift[:, jyseps1_1 + 1] = 0.5 * dphidy[:, jyseps1_1 + 1] * dy + # Note: This goes to jyseps2_1 + 1, the first cell in the upper inner leg (including upper PF) + for jy in range(jyseps1_1 + 2, jyseps2_1 + 2): + zShift[:, jy] = ( + zShift[:, jy - 1] + 0.5 * (dphidy[:, jy] + dphidy[:, jy - 1]) * dy + ) + + # Outer core + # Note: ixseps2 points are connected from inner to outer core + # If single null, ixseps2 = nx, jyseps1_2 = jyseps2_1 + zShift[:, jyseps1_2 + 1] = ( + zShift[:, jyseps2_1] + + 0.5 * (dphidy[:, jyseps2_1] + dphidy[:, jyseps1_2 + 1]) * dy + ) + for jy in range(jyseps1_2 + 2, jyseps2_2 + 2): + zShift[:, jy] = ( + zShift[:, jy - 1] + 0.5 * (dphidy[:, jy] + dphidy[:, jy - 1]) * dy + ) + + if jyseps2_1 != jyseps1_2: + # Double null, with upper X-point + + # Upper inner leg. Note that upper PF region set from core at y index jyseps2_1 + 1 + for jy in range(jyseps2_1 + 2, ny_inner): + zShift[:, jy] = ( + zShift[:, jy - 1] + 0.5 * (dphidy[:, jy] + dphidy[:, jy - 1]) * dy + ) + + # Upper outer leg. Joins onto upper inner leg (jyseps2_1 + 1) for x < ixseps2 + zShift[:ixseps2, jyseps1_2] = ( + zShift[:ixseps2, jyseps2_1 + 1] + - 0.5 * (dphidy[:ixseps2, jyseps2_1 + 1] + dphidy[:ixseps2, jyseps1_2]) * dy + ) + # joins outer core/SOL region for x >= ixseps2 + zShift[ixseps2:, jyseps1_2] = ( + zShift[ixseps2:, jyseps1_2 + 1] + - 0.5 * (dphidy[ixseps2:, jyseps1_2 + 1] + dphidy[ixseps2:, jyseps1_2]) * dy + ) + # Iterate backwards along upper outer leg from X-point towards target + for jy in range(jyseps1_2 - 1, ny_inner - 1, -1): + zShift[:, jy] = ( + zShift[:, jy + 1] - 0.5 * (dphidy[:, jy] + dphidy[:, jy - 1]) * dy + ) + + # Lower outer leg + zShift[:ixseps1, jyseps2_2 + 1] = 0.5 * dphidy[:ixseps1, jyseps2_2 + 1] * dy + for jy in range(jyseps2_2 + 2, ny): + zShift[:, jy] = ( + zShift[:, jy - 1] + 0.5 * (dphidy[:, jy] + dphidy[:, jy - 1]) * dy + ) + # Lower inner leg. Going backwards in Y toward the plate + zShift[:ixseps1, jyseps1_1] = -0.5 * dphidy[:ixseps1, jyseps1_1] * dy + for jy in range(jyseps1_1 - 1, -1, -1): + zShift[:, jy] = ( + zShift[:, jy + 1] - 0.5 * (dphidy[:, jy] + dphidy[:, jy + 1]) * dy + ) + + ShiftAngle = np.zeros(nx) + ShiftAngle[:ixseps1] = np.sum( + dphidy[:ixseps1, (jyseps1_1 + 1) : (jyseps2_1 + 1)] * dy, axis=1 # Inner core + ) + np.sum( + dphidy[:ixseps1, (jyseps1_2 + 1) : (jyseps2_2 + 1)] * dy, axis=1 # Outer core + ) + + if verbose: + print( + "Safety factor: min {}, mean {}, max {}".format( + np.amin(ShiftAngle[:ixseps1]) / (2 * np.pi), + np.mean(ShiftAngle[:ixseps1]) / (2 * np.pi), + np.amax(ShiftAngle[:ixseps1]) / (2 * np.pi), + ) + ) + + if plotting: + plt.plot(Rxy, Zxy, "x") #needs to be transposed. + plt.plot(Rxy[ixseps1, :], Zxy[ixseps1, :], color="magenta", label="ixseps1") + if ixseps2 < nx: + plt.plot(Rxy[ixseps2, :], Zxy[ixseps2, :], color="r", label="ixseps2") + + plt.plot(Rxy[:, jyseps1_1], Zxy[:, jyseps1_1], color="k", label="jyseps1_1") + plt.plot(Rxy[:, jyseps1_2], Zxy[:, jyseps1_2], color="b", label="jyseps1_2") + plt.plot(Rxy[:, ny_inner], Zxy[:, ny_inner], color="c", label="ny_inner") + plt.plot(Rxy[:, jyseps2_1], Zxy[:, jyseps2_1], color="g", label="jyseps2_1") + plt.plot(Rxy[:, jyseps2_2], Zxy[:, jyseps2_2], color="r", label="jyseps2_2") + + ax = plt.gca() + ax.set_aspect("equal", adjustable="box") + ax.set_xlabel("R") + ax.set_ylabel("Z") + plt.legend() + plt.show() + + from boututils.datafile import DataFile + + if verbose: + print("Saving to " + output_filename) + + with DataFile(output_filename, create=True, format="NETCDF4") as f: + # Save unique ID for grid file + import uuid + + f.write_file_attribute("grid_id", str(uuid.uuid1())) + f.write_file_attribute("gridue", str(gridue_file)) + + f.write("nx", nx) + f.write("ny", ny) + f.write("ixseps1", ixseps1) + f.write("ixseps2", ixseps2) + f.write("jyseps1_1", jyseps1_1) + f.write("jyseps2_1", jyseps2_1) + f.write("ny_inner", ny_inner) + f.write("jyseps1_2", jyseps1_2) + f.write("jyseps2_2", jyseps2_2) + f.write("rm", rm) + f.write("zm", zm) + f.write("topology", mesh_topology) + + # 2D fields + for name in grd: + f.write(name, grd[name]) + + f.write("zShift", zShift) + f.write("ShiftAngle", ShiftAngle) + +def getMeshTopology(g, nx, ny): + """ + Get mesh topology from gridue data. + + Returns: + -------- + dict + A dictionary containing the mesh topology information, including: + - "single_null": bool indicating if the mesh is single null or double null. + - "ixseps1": int, index of the lower X-point separatrix. + - "ixseps2": int, index of the upper X-point separatrix (if double null). + - "jyseps1_1": int, index of the lower inner leg separatrix. + - "jyseps2_1": int, index of the upper inner leg separatrix. + - "ny_inner": int, number of inner poloidal points. + - "jyseps1_2": int, index of the lower outer leg separatrix. + - "jyseps2_2": int, index of the upper outer leg separatrix. + """ + + ixseps1 = g["iyseparatrix1"] + 2 # Lower X-point separatrix + ixseps2 = min(g["iyseparatrix2"] + 2, nx) # Upper X-point separatrix + jyseps1_1 = g["ix_cut1"] - 1 + jyseps2_1 = g["ix_cut2"] + ny_inner = g["ix_inner"] + jyseps1_2 = g["ix_cut3"] + jyseps2_2 = g["ix_cut4"] - 1 + + if (jyseps1_1 < 0 and jyseps2_2 >= ny - 1): + return "CFL" + elif (jyseps2_1 == jyseps1_2): + return "SN" + elif (jyseps1_2 <= ny_inner and ny_inner <= jyseps2_2): + return "SF" + elif (ixseps1 == ixseps2): + return "CDN" + else: + return "UDN"; + +if __name__ == "__main__": + main() From bca3491ade78467939df858adbf7a2abb88fda43 Mon Sep 17 00:00:00 2001 From: wearysebas Date: Mon, 9 Mar 2026 12:15:03 +0000 Subject: [PATCH 2/3] Fixed an issue with the branch cuts not being consistent between BOUT and INGRID for Single Null topology. --- .../gridue_to_bout_converter/gridue_to_bout.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/boutdata/gridue_to_bout_converter/gridue_to_bout.py b/src/boutdata/gridue_to_bout_converter/gridue_to_bout.py index 334959c33..6fabf6762 100644 --- a/src/boutdata/gridue_to_bout_converter/gridue_to_bout.py +++ b/src/boutdata/gridue_to_bout_converter/gridue_to_bout.py @@ -769,6 +769,13 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False jyseps1_2 = g["ix_cut3"] jyseps2_2 = g["ix_cut4"] - 1 + #jyseps2_1 should always be smaller that jyseps1_2. Only inconsistency found here is for SN. + if jyseps1_2 < jyseps2_1: + jyseps1_2 = jyseps2_1 + ny_inner = jyseps2_1 + #For Single Null (SN) cases, Ingrid sets this values to be different, but BOUT++ expects them to be the same. + print ("WARNING: Adjusting jyseps1_2 to be equal to jyseps2_1 for consistency with BOUT++ expectations. This is expected for Single Null cases.") + # Calculate metric tensor grd.update(calcMetric(grd, bpsign, verbose, ignore_checks)) @@ -930,6 +937,13 @@ def getMeshTopology(g, nx, ny): ny_inner = g["ix_inner"] jyseps1_2 = g["ix_cut3"] jyseps2_2 = g["ix_cut4"] - 1 + + #jyseps2_1 should always be smaller that jyseps1_2. Only inconsistency found here is for SN. + if jyseps1_2 < jyseps2_1: + jyseps1_2 = jyseps2_1 + ny_inner = jyseps2_1 + #For Single Null (SN) cases, Ingrid sets this values to be different, but BOUT++ expects them to be the same. + print ("WARNING: Adjusting jyseps1_2 to be equal to jyseps2_1 for consistency with BOUT++ expectations. This is expected for Single Null cases.") if (jyseps1_1 < 0 and jyseps2_2 >= ny - 1): return "CFL" From 7853c91db7b238a2b81a67695a812b5a73debf9a Mon Sep 17 00:00:00 2001 From: wearysebas Date: Mon, 9 Mar 2026 12:22:19 +0000 Subject: [PATCH 3/3] Apply ruff changes --- .../gridue_to_bout.py | 95 +++++++++++-------- 1 file changed, 57 insertions(+), 38 deletions(-) diff --git a/src/boutdata/gridue_to_bout_converter/gridue_to_bout.py b/src/boutdata/gridue_to_bout_converter/gridue_to_bout.py index 6fabf6762..737e9ad1c 100644 --- a/src/boutdata/gridue_to_bout_converter/gridue_to_bout.py +++ b/src/boutdata/gridue_to_bout_converter/gridue_to_bout.py @@ -12,10 +12,9 @@ # # -import numpy as np import matplotlib import matplotlib.pyplot as plt - +import numpy as np from scipy import linalg @@ -64,7 +63,6 @@ def _importBody(gridue_settings, f): for n in range(5): for j in range(ny): for i in range(nx): - data_[i][j][n] = float(vv) try: @@ -194,7 +192,7 @@ def plot(GridueParams: dict, edgecolor="black", ax: object = None, show=True): return ax -def calcHy( nx, ny, g: dict, dy: float): +def calcHy(nx, ny, g: dict, dy: float): """ Calculate poloidal arc length metric from gridue dictionary """ @@ -294,17 +292,17 @@ def calcMetric(grd: dict, bpsign, verbose=False, ignore_checks=False): Dictionary containing grid data with keys: "Rxy", "Zxy", "Brxy", "Bzxy", "Btxy", "Bpxy", "Bxy", "hy", "cosBeta", "tanBeta", "curl_bOverB_Rhat", "curl_bOverB_Zhat", - "curl_bOverB_zetahat". + "curl_bOverB_zetahat". bpsign : int Sign of the magnetic field (1 for normal, -1 for reversed). - verbose : bool, optional + verbose : bool, optional If True, print detailed information about the metric tensor and Jacobian. ignore_checks : bool, optional If True, ignore checks on the Jacobian's relative error. Returns: dict - A dictionary containing the metric tensor components, Jacobian, and other related quantities for each grid point. + A dictionary containing the metric tensor components, Jacobian, and other related quantities for each grid point. --> Each component is a 2D array with shape [radial, poloidal]. """ @@ -389,7 +387,7 @@ def calcMetric(grd: dict, bpsign, verbose=False, ignore_checks=False): if np.max(np.abs(rel_error)) > 1e-6: if ignore_checks: print("WARNING: Relative error in Jacobian too large.") - #else: + # else: # raise ValueError(f"Relative error in Jacobian too large: {np.max(np.abs(rel_error))}") # We want to output contravariant components of Curl(b/B) in the @@ -530,6 +528,7 @@ def calcRZCurvature(g: dict): # Return as [radial, poloidal] return curl_bOverB_Rhat.T, curl_bOverB_Zhat.T, curl_bOverB_zetahat.T + def main(): import argparse @@ -558,7 +557,14 @@ def main(): Convert_grids(gridue_file, output_filename, plotting, verbose, ignore_checks) -def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False, verbose: bool = False, ignore_checks: bool = False): + +def Convert_grids( + gridue_file: str, + output_filename: str, + plotting: bool = False, + verbose: bool = False, + ignore_checks: bool = False, +): """ Convert UEDGE grid file to BOUT++ grid format. @@ -620,10 +626,10 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False psixy = psi[:, :, 0].T nx, ny = Rxy.shape - #idx = [np.array([1, 2, 4, 3, 1])] + # idx = [np.array([1, 2, 4, 3, 1])] - #pol = [] - #for i in range(nx): + # pol = [] + # for i in range(nx): # for j in range(ny): # np.concatenate((rm[i][j][idx], zm[i][j][idx])).reshape(2, 5).T @@ -639,10 +645,12 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False dx = np.zeros((nx, ny)) for i in range(nx): for j in range(ny): - if i > 1 and i < nx-2: - dx[i, j] = 0.5*(psi[j, i + 1, 0] - psi[j, i - 1, 0]) + if i > 1 and i < nx - 2: + dx[i, j] = 0.5 * (psi[j, i + 1, 0] - psi[j, i - 1, 0]) else: - dx[i, j] = 0.5 * (psi[j, i, 3] + psi[j, i, 4] - psi[j, i, 1] - psi[j, i, 2]) + dx[i, j] = 0.5 * ( + psi[j, i, 3] + psi[j, i, 4] - psi[j, i, 1] - psi[j, i, 2] + ) # Note: UEDGE grids have narrow cells on the radial # boundaries. BOUT++ applies boundary conditions half-way between @@ -664,7 +672,7 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False # Calculate hy, the arc length along the flux surface passing through # the center of each cell. - hy = calcHy(nx,ny,g, dy) + hy = calcHy(nx, ny, g, dy) # Calculate angle between x and y coordinates. sinBeta = 0, cosBeta = 1 for an orthogonal mesh sinBeta, cosBeta = calcGridAngle(g) @@ -706,11 +714,11 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False Zxy = grd["Zxy"] nx, ny = Rxy.shape - #Get Mesh Topology and remove guard cells accordingly + # Get Mesh Topology and remove guard cells accordingly mesh_topology = getMeshTopology(g, nx, ny) if mesh_topology == "SF": - #SF case + # SF case ixseps1 = g["iyseparatrix1"] + 2 # Main X-point separatrix ixseps2 = min(g["iyseparatrix3"] + 2, nx) # Secondary X-point separatrix # Remove guard cells on either side of X-point. @@ -769,12 +777,14 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False jyseps1_2 = g["ix_cut3"] jyseps2_2 = g["ix_cut4"] - 1 - #jyseps2_1 should always be smaller that jyseps1_2. Only inconsistency found here is for SN. + # jyseps2_1 should always be smaller that jyseps1_2. Only inconsistency found here is for SN. if jyseps1_2 < jyseps2_1: jyseps1_2 = jyseps2_1 ny_inner = jyseps2_1 - #For Single Null (SN) cases, Ingrid sets this values to be different, but BOUT++ expects them to be the same. - print ("WARNING: Adjusting jyseps1_2 to be equal to jyseps2_1 for consistency with BOUT++ expectations. This is expected for Single Null cases.") + # For Single Null (SN) cases, Ingrid sets this values to be different, but BOUT++ expects them to be the same. + print( + "WARNING: Adjusting jyseps1_2 to be equal to jyseps2_1 for consistency with BOUT++ expectations. This is expected for Single Null cases." + ) # Calculate metric tensor grd.update(calcMetric(grd, bpsign, verbose, ignore_checks)) @@ -846,10 +856,15 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False ) ShiftAngle = np.zeros(nx) - ShiftAngle[:ixseps1] = np.sum( - dphidy[:ixseps1, (jyseps1_1 + 1) : (jyseps2_1 + 1)] * dy, axis=1 # Inner core - ) + np.sum( - dphidy[:ixseps1, (jyseps1_2 + 1) : (jyseps2_2 + 1)] * dy, axis=1 # Outer core + ShiftAngle[:ixseps1] = ( + np.sum( + dphidy[:ixseps1, (jyseps1_1 + 1) : (jyseps2_1 + 1)] * dy, + axis=1, # Inner core + ) + + np.sum( + dphidy[:ixseps1, (jyseps1_2 + 1) : (jyseps2_2 + 1)] * dy, + axis=1, # Outer core + ) ) if verbose: @@ -862,7 +877,7 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False ) if plotting: - plt.plot(Rxy, Zxy, "x") #needs to be transposed. + plt.plot(Rxy, Zxy, "x") # needs to be transposed. plt.plot(Rxy[ixseps1, :], Zxy[ixseps1, :], color="magenta", label="ixseps1") if ixseps2 < nx: plt.plot(Rxy[ixseps2, :], Zxy[ixseps2, :], color="r", label="ixseps2") @@ -888,7 +903,7 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False with DataFile(output_filename, create=True, format="NETCDF4") as f: # Save unique ID for grid file import uuid - + f.write_file_attribute("grid_id", str(uuid.uuid1())) f.write_file_attribute("gridue", str(gridue_file)) @@ -912,6 +927,7 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False f.write("zShift", zShift) f.write("ShiftAngle", ShiftAngle) + def getMeshTopology(g, nx, ny): """ Get mesh topology from gridue data. @@ -929,32 +945,35 @@ def getMeshTopology(g, nx, ny): - "jyseps1_2": int, index of the lower outer leg separatrix. - "jyseps2_2": int, index of the upper outer leg separatrix. """ - + ixseps1 = g["iyseparatrix1"] + 2 # Lower X-point separatrix - ixseps2 = min(g["iyseparatrix2"] + 2, nx) # Upper X-point separatrix + ixseps2 = min(g["iyseparatrix2"] + 2, nx) # Upper X-point separatrix jyseps1_1 = g["ix_cut1"] - 1 jyseps2_1 = g["ix_cut2"] ny_inner = g["ix_inner"] jyseps1_2 = g["ix_cut3"] jyseps2_2 = g["ix_cut4"] - 1 - #jyseps2_1 should always be smaller that jyseps1_2. Only inconsistency found here is for SN. + # jyseps2_1 should always be smaller that jyseps1_2. Only inconsistency found here is for SN. if jyseps1_2 < jyseps2_1: jyseps1_2 = jyseps2_1 ny_inner = jyseps2_1 - #For Single Null (SN) cases, Ingrid sets this values to be different, but BOUT++ expects them to be the same. - print ("WARNING: Adjusting jyseps1_2 to be equal to jyseps2_1 for consistency with BOUT++ expectations. This is expected for Single Null cases.") - - if (jyseps1_1 < 0 and jyseps2_2 >= ny - 1): + # For Single Null (SN) cases, Ingrid sets this values to be different, but BOUT++ expects them to be the same. + print( + "WARNING: Adjusting jyseps1_2 to be equal to jyseps2_1 for consistency with BOUT++ expectations. This is expected for Single Null cases." + ) + + if jyseps1_1 < 0 and jyseps2_2 >= ny - 1: return "CFL" - elif (jyseps2_1 == jyseps1_2): + elif jyseps2_1 == jyseps1_2: return "SN" - elif (jyseps1_2 <= ny_inner and ny_inner <= jyseps2_2): + elif jyseps1_2 <= ny_inner and ny_inner <= jyseps2_2: return "SF" - elif (ixseps1 == ixseps2): + elif ixseps1 == ixseps2: return "CDN" else: - return "UDN"; + return "UDN" + if __name__ == "__main__": main()