1212#
1313#
1414
15- import numpy as np
1615import matplotlib
1716import matplotlib .pyplot as plt
18-
17+ import numpy as np
1918from scipy import linalg
2019
2120
@@ -64,7 +63,6 @@ def _importBody(gridue_settings, f):
6463 for n in range (5 ):
6564 for j in range (ny ):
6665 for i in range (nx ):
67-
6866 data_ [i ][j ][n ] = float (vv )
6967
7068 try :
@@ -194,7 +192,7 @@ def plot(GridueParams: dict, edgecolor="black", ax: object = None, show=True):
194192 return ax
195193
196194
197- def calcHy ( nx , ny , g : dict , dy : float ):
195+ def calcHy (nx , ny , g : dict , dy : float ):
198196 """
199197 Calculate poloidal arc length metric from gridue dictionary
200198 """
@@ -294,17 +292,17 @@ def calcMetric(grd: dict, bpsign, verbose=False, ignore_checks=False):
294292 Dictionary containing grid data with keys:
295293 "Rxy", "Zxy", "Brxy", "Bzxy", "Btxy", "Bpxy", "Bxy", "hy",
296294 "cosBeta", "tanBeta", "curl_bOverB_Rhat", "curl_bOverB_Zhat",
297- "curl_bOverB_zetahat".
295+ "curl_bOverB_zetahat".
298296 bpsign : int
299297 Sign of the magnetic field (1 for normal, -1 for reversed).
300- verbose : bool, optional
298+ verbose : bool, optional
301299 If True, print detailed information about the metric tensor and Jacobian.
302300 ignore_checks : bool, optional
303301 If True, ignore checks on the Jacobian's relative error.
304302
305303 Returns:
306304 dict
307- A dictionary containing the metric tensor components, Jacobian, and other related quantities for each grid point.
305+ A dictionary containing the metric tensor components, Jacobian, and other related quantities for each grid point.
308306 --> Each component is a 2D array with shape [radial, poloidal].
309307 """
310308
@@ -389,7 +387,7 @@ def calcMetric(grd: dict, bpsign, verbose=False, ignore_checks=False):
389387 if np .max (np .abs (rel_error )) > 1e-6 :
390388 if ignore_checks :
391389 print ("WARNING: Relative error in Jacobian too large." )
392- #else:
390+ # else:
393391 # raise ValueError(f"Relative error in Jacobian too large: {np.max(np.abs(rel_error))}")
394392
395393 # We want to output contravariant components of Curl(b/B) in the
@@ -530,6 +528,7 @@ def calcRZCurvature(g: dict):
530528 # Return as [radial, poloidal]
531529 return curl_bOverB_Rhat .T , curl_bOverB_Zhat .T , curl_bOverB_zetahat .T
532530
531+
533532def main ():
534533 import argparse
535534
@@ -558,7 +557,14 @@ def main():
558557
559558 Convert_grids (gridue_file , output_filename , plotting , verbose , ignore_checks )
560559
561- def Convert_grids (gridue_file : str , output_filename : str , plotting : bool = False , verbose : bool = False , ignore_checks : bool = False ):
560+
561+ def Convert_grids (
562+ gridue_file : str ,
563+ output_filename : str ,
564+ plotting : bool = False ,
565+ verbose : bool = False ,
566+ ignore_checks : bool = False ,
567+ ):
562568 """
563569 Convert UEDGE grid file to BOUT++ grid format.
564570
@@ -620,10 +626,10 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False
620626 psixy = psi [:, :, 0 ].T
621627 nx , ny = Rxy .shape
622628
623- #idx = [np.array([1, 2, 4, 3, 1])]
629+ # idx = [np.array([1, 2, 4, 3, 1])]
624630
625- #pol = []
626- #for i in range(nx):
631+ # pol = []
632+ # for i in range(nx):
627633 # for j in range(ny):
628634 # np.concatenate((rm[i][j][idx], zm[i][j][idx])).reshape(2, 5).T
629635
@@ -639,10 +645,12 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False
639645 dx = np .zeros ((nx , ny ))
640646 for i in range (nx ):
641647 for j in range (ny ):
642- if i > 1 and i < nx - 2 :
643- dx [i , j ] = 0.5 * (psi [j , i + 1 , 0 ] - psi [j , i - 1 , 0 ])
648+ if i > 1 and i < nx - 2 :
649+ dx [i , j ] = 0.5 * (psi [j , i + 1 , 0 ] - psi [j , i - 1 , 0 ])
644650 else :
645- dx [i , j ] = 0.5 * (psi [j , i , 3 ] + psi [j , i , 4 ] - psi [j , i , 1 ] - psi [j , i , 2 ])
651+ dx [i , j ] = 0.5 * (
652+ psi [j , i , 3 ] + psi [j , i , 4 ] - psi [j , i , 1 ] - psi [j , i , 2 ]
653+ )
646654
647655 # Note: UEDGE grids have narrow cells on the radial
648656 # boundaries. BOUT++ applies boundary conditions half-way between
@@ -664,7 +672,7 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False
664672
665673 # Calculate hy, the arc length along the flux surface passing through
666674 # the center of each cell.
667- hy = calcHy (nx ,ny ,g , dy )
675+ hy = calcHy (nx , ny , g , dy )
668676
669677 # Calculate angle between x and y coordinates. sinBeta = 0, cosBeta = 1 for an orthogonal mesh
670678 sinBeta , cosBeta = calcGridAngle (g )
@@ -706,11 +714,11 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False
706714 Zxy = grd ["Zxy" ]
707715 nx , ny = Rxy .shape
708716
709- #Get Mesh Topology and remove guard cells accordingly
717+ # Get Mesh Topology and remove guard cells accordingly
710718 mesh_topology = getMeshTopology (g , nx , ny )
711719
712720 if mesh_topology == "SF" :
713- #SF case
721+ # SF case
714722 ixseps1 = g ["iyseparatrix1" ] + 2 # Main X-point separatrix
715723 ixseps2 = min (g ["iyseparatrix3" ] + 2 , nx ) # Secondary X-point separatrix
716724 # 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
769777 jyseps1_2 = g ["ix_cut3" ]
770778 jyseps2_2 = g ["ix_cut4" ] - 1
771779
772- # jyseps2_1 should always be smaller that jyseps1_2. Only inconsistency found here is for SN.
780+ # jyseps2_1 should always be smaller that jyseps1_2. Only inconsistency found here is for SN.
773781 if jyseps1_2 < jyseps2_1 :
774782 jyseps1_2 = jyseps2_1
775783 ny_inner = jyseps2_1
776- #For Single Null (SN) cases, Ingrid sets this values to be different, but BOUT++ expects them to be the same.
777- print ("WARNING: Adjusting jyseps1_2 to be equal to jyseps2_1 for consistency with BOUT++ expectations. This is expected for Single Null cases." )
784+ # For Single Null (SN) cases, Ingrid sets this values to be different, but BOUT++ expects them to be the same.
785+ print (
786+ "WARNING: Adjusting jyseps1_2 to be equal to jyseps2_1 for consistency with BOUT++ expectations. This is expected for Single Null cases."
787+ )
778788
779789 # Calculate metric tensor
780790 grd .update (calcMetric (grd , bpsign , verbose , ignore_checks ))
@@ -846,10 +856,15 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False
846856 )
847857
848858 ShiftAngle = np .zeros (nx )
849- ShiftAngle [:ixseps1 ] = np .sum (
850- dphidy [:ixseps1 , (jyseps1_1 + 1 ) : (jyseps2_1 + 1 )] * dy , axis = 1 # Inner core
851- ) + np .sum (
852- dphidy [:ixseps1 , (jyseps1_2 + 1 ) : (jyseps2_2 + 1 )] * dy , axis = 1 # Outer core
859+ ShiftAngle [:ixseps1 ] = (
860+ np .sum (
861+ dphidy [:ixseps1 , (jyseps1_1 + 1 ) : (jyseps2_1 + 1 )] * dy ,
862+ axis = 1 , # Inner core
863+ )
864+ + np .sum (
865+ dphidy [:ixseps1 , (jyseps1_2 + 1 ) : (jyseps2_2 + 1 )] * dy ,
866+ axis = 1 , # Outer core
867+ )
853868 )
854869
855870 if verbose :
@@ -862,7 +877,7 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False
862877 )
863878
864879 if plotting :
865- plt .plot (Rxy , Zxy , "x" ) # needs to be transposed.
880+ plt .plot (Rxy , Zxy , "x" ) # needs to be transposed.
866881 plt .plot (Rxy [ixseps1 , :], Zxy [ixseps1 , :], color = "magenta" , label = "ixseps1" )
867882 if ixseps2 < nx :
868883 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
888903 with DataFile (output_filename , create = True , format = "NETCDF4" ) as f :
889904 # Save unique ID for grid file
890905 import uuid
891-
906+
892907 f .write_file_attribute ("grid_id" , str (uuid .uuid1 ()))
893908 f .write_file_attribute ("gridue" , str (gridue_file ))
894909
@@ -912,6 +927,7 @@ def Convert_grids(gridue_file: str, output_filename: str, plotting: bool = False
912927 f .write ("zShift" , zShift )
913928 f .write ("ShiftAngle" , ShiftAngle )
914929
930+
915931def getMeshTopology (g , nx , ny ):
916932 """
917933 Get mesh topology from gridue data.
@@ -929,32 +945,35 @@ def getMeshTopology(g, nx, ny):
929945 - "jyseps1_2": int, index of the lower outer leg separatrix.
930946 - "jyseps2_2": int, index of the upper outer leg separatrix.
931947 """
932-
948+
933949 ixseps1 = g ["iyseparatrix1" ] + 2 # Lower X-point separatrix
934- ixseps2 = min (g ["iyseparatrix2" ] + 2 , nx ) # Upper X-point separatrix
950+ ixseps2 = min (g ["iyseparatrix2" ] + 2 , nx ) # Upper X-point separatrix
935951 jyseps1_1 = g ["ix_cut1" ] - 1
936952 jyseps2_1 = g ["ix_cut2" ]
937953 ny_inner = g ["ix_inner" ]
938954 jyseps1_2 = g ["ix_cut3" ]
939955 jyseps2_2 = g ["ix_cut4" ] - 1
940956
941- #jyseps2_1 should always be smaller that jyseps1_2. Only inconsistency found here is for SN.
957+ # jyseps2_1 should always be smaller that jyseps1_2. Only inconsistency found here is for SN.
942958 if jyseps1_2 < jyseps2_1 :
943959 jyseps1_2 = jyseps2_1
944960 ny_inner = jyseps2_1
945- #For Single Null (SN) cases, Ingrid sets this values to be different, but BOUT++ expects them to be the same.
946- print ("WARNING: Adjusting jyseps1_2 to be equal to jyseps2_1 for consistency with BOUT++ expectations. This is expected for Single Null cases." )
947-
948- if (jyseps1_1 < 0 and jyseps2_2 >= ny - 1 ):
961+ # For Single Null (SN) cases, Ingrid sets this values to be different, but BOUT++ expects them to be the same.
962+ print (
963+ "WARNING: Adjusting jyseps1_2 to be equal to jyseps2_1 for consistency with BOUT++ expectations. This is expected for Single Null cases."
964+ )
965+
966+ if jyseps1_1 < 0 and jyseps2_2 >= ny - 1 :
949967 return "CFL"
950- elif ( jyseps2_1 == jyseps1_2 ) :
968+ elif jyseps2_1 == jyseps1_2 :
951969 return "SN"
952- elif ( jyseps1_2 <= ny_inner and ny_inner <= jyseps2_2 ) :
970+ elif jyseps1_2 <= ny_inner and ny_inner <= jyseps2_2 :
953971 return "SF"
954- elif ( ixseps1 == ixseps2 ) :
972+ elif ixseps1 == ixseps2 :
955973 return "CDN"
956974 else :
957- return "UDN" ;
975+ return "UDN"
976+
958977
959978if __name__ == "__main__" :
960979 main ()
0 commit comments