Feature/jeppe fixes

vita/modules/equilibrium/equdsk/equdsk.py

@@ -2,18 +2,13 @@
 # -*- coding: utf-8 -*-
 """
 Created on Tue Oct 22 09:53:55 2019
-
 Addapted from: ReadEQDSK
-
 Original AUTHOR: Leonardo Pigatto
 Maintainer: jmbols
-
 DESCRIPTION
 Python function to read EQDSK files
-
 CALLING SEQUENCE
 out = ReadEQDSK(in_filename)
-
 CHANGE HISTORY:
 -started on 29/09/2015 - porting from Matlab
 -16/10/2015 - generators introduced for reading input file
@@ -22,7 +17,6 @@
 -02/2018 - added  backwards compatibility with python 2
 -15/07/2020 - changed to comply with PEP8
 NOTES:
-
 """
 try:
     import builtins                  # <- python 3
@@ -36,192 +30,7 @@
 #from pyTokamak.formats.geqdsk import file_numbers,file_tokens
 
 # Other imports
-#from vita.modules.utils import intersection
-
-class eqdsk():
-	def __init__(self,comment,switch,nrbox,nzbox,rboxlength,zboxlength,R0EXP,rboxleft,zmid,Raxis,Zaxis,psiaxis,psiedge,B0EXP,Ip,T,p,TTprime,pprime,psi,q,nLCFS,nlimits,R,Z,R_limits,Z_limits,R_grid,Z_grid,psi_grid,rhopsi):
-		self.comment=comment
-		self.switch=switch
-		self.nrbox=nrbox
-		self.nzbox=nzbox
-		self.rboxlength=rboxlength
-		self.zboxlength=zboxlength
-		self.R0EXP=R0EXP
-		self.rboxleft=rboxleft
-		self.zmid = zmid
-		self.Raxis=Raxis
-		self.Zaxis=Zaxis
-		self.psiaxis=psiaxis
-		self.psiedge=psiedge
-		self.B0EXP=B0EXP
-		self.Ip=Ip
-		self.T=T
-		self.p=p
-		self.TTprime=TTprime
-		self.pprime=pprime
-		self.psi=psi
-		self.q=q
-		self.nLCFS=nLCFS
-		self.nlimits=nlimits
-		self.R=R
-		self.Z=Z
-		self.R_limits=R_limits
-		self.Z_limits=Z_limits
-		self.R_grid=R_grid
-		self.Z_grid=Z_grid
-		self.psi_grid=psi_grid
-		self.rhopsi=rhopsi
-
-
-#def file_tokens(fp):
-	#""" A generator to split a file into tokens
-	#"""
-	#toklist = []
-	#while True:
-		#line = fp.readline()
-		#if not line: break
-		#toklist = line.split()
-		#for tok in toklist:
-			#yield tok
-
-def file_numbers(fp):
-	"""Generator to get numbers from a text file"""
-	toklist = []
-	while True:
-		line = fp.readline()
-		if not line: break
-		# Match numbers in the line using regular expression
-		pattern = r'[+-]?\d*[\.]?\d+(?:[Ee][+-]?\d+)?'
-		toklist = re.findall(pattern, line)
-		for tok in toklist:
-			yield tok
-
-
-def ReadEQDSK(in_filename):
-
-	fin = open(in_filename,"r")
-
-	desc = fin.readline()
-	data = desc.split()
-	switch = int(data[-3])
-	nrbox = int(data[-2])
-	nzbox = int(data[-1])
-	comment = data[0:-3]
-
-	token = file_numbers(fin)
-
-	#first line
-	rboxlength = float(token.__next__())
-	zboxlength = float(token.__next__())
-	R0EXP = float(token.__next__())
-	rboxleft = float(token.__next__())
-	zmid   = float(token.__next__()) #(maxygrid+minygrid)/2
-
-	#second line
-	Raxis = float(token.__next__())
-	Zaxis = float(token.__next__())
-	psiaxis = float(token.__next__()) # psi_axis-psi_edge
-	psiedge = float(token.__next__()) # psi_edge-psi_edge (=0)
-	B0EXP = float(token.__next__()) # normalizing magnetic field in chease
-
-	#third line
-	#ip is first element, all others are already stored
-	Ip = float(token.__next__())
-
-	""" DEFINING USEFUL FUNCTIONS TO READ ARRAYS """
-	def consume(iterator,n):
-		#Advance iterator n steps ahead
-		next(islice(iterator,n,n),None)
-
-	def read_array(n,name="Unknown"):
-		data = np.zeros([n])
-		try:
-			for i in np.arange(n):
-				data[i] = float(token.__next__())
-		except:
-			raise IOError("Failed reading array '"+name+"' of size ", n)
-		return data
-
-	def read_2d(nr,nz,name="Unknown"):
-		data = np.zeros([nr, nz])
-		for i in np.arange(nr):
-			data[i,:] = read_array(nz, name+"["+str(i)+"]")
-		return data
-
-	#fourth line - nothing or already stored
-	#advance to next significant token
-	consume(token,9)
-
-	#T (or T - poloidal flux function)
-	T = read_array(nrbox,"T")
-
-	#p (pressure)
-	p = read_array(nrbox,"p")
-
-	#TT'
-	TTprime = read_array(nrbox,"TTprime")
-
-	#p'
-	pprime = read_array(nrbox,"pprime")
-
-	#psi
-	psi = read_2d(nrbox,nzbox,"psi")
-
-	#q safety factor
-	q = read_array(nrbox,"safety_factor")
-
-	#n of points for the lcfs and limiter boundary
-	nLCFS = int(token.__next__())
-	nlimits = int(token.__next__())
-
-	#rz lcfs coordinates
-	if nLCFS > 0:
-		R = np.zeros([nLCFS])
-		Z = np.zeros([nLCFS])
-		for ii in range(nLCFS):
-			R[ii] = float(token.__next__())
-			Z[ii] = float(token.__next__())
-	else:
-		R = [0]
-		Z = [0]
-
-	#rz limits
-	if nlimits > 0:
-		R_limits = np.zeros([nlimits])
-		Z_limits = np.zeros([nlimits])
-		for ii in range(nlimits):
-			R_limits[ii] = float(token.__next__())
-			Z_limits[ii] = float(token.__next__())
-	else:
-		R_limits = [0]
-		Z_limits = [0]
-
-	# construct r-z mesh
-	R_grid = np.zeros([nrbox,nzbox])
-	Z_grid = R_grid.copy()
-	for ii in range(nrbox):
-		R_grid[ii,:] = rboxleft + rboxlength*ii/float(nrbox-1)
-	for jj in range(nzbox):
-		Z_grid[:,jj] = (zmid-0.5*zboxlength) + zboxlength*jj/float(nzbox-1)
-
-	#psi grid for radial prifiles
-	psi_grid = np.linspace(psiaxis,psiedge,nrbox)
-
-	#corresponding s=sqrt(psi_bar)
-	rhopsi = np.sqrt(abs(psi_grid-psiaxis)/abs(psiedge-psiaxis))
-	fin.close()
-
-
-	##RZ grid for psi
-	#R_grid = np.linspace(rboxleft,rboxleft+rboxlength,nrbox)
-	#Z_grid = np.linspace(-zboxlength/2.,zboxlength/2.,nzbox)
-	##psi grid for radial prifiles
-	#psi_grid = np.linspace(psiaxis,psiedge,nrbox)
-	##corresponding rhopsi
-	#rhopsi = sqrt(abs(psi_grid-psiaxis)/abs(psiedge-psiaxis))
-
-	out = eqdsk(comment,switch,nrbox,nzbox,rboxlength,zboxlength,R0EXP,rboxleft,zmid,Raxis,Zaxis,psiaxis,psiedge,B0EXP,Ip,T,p,TTprime,pprime,psi,q,nLCFS,nlimits,R,Z,R_limits,Z_limits,R_grid,Z_grid,psi_grid,rhopsi)
-	return out
+from vita.modules.utils import intersection
 
 class Equdsk():
     '''
@@ -272,17 +81,14 @@ def __init__(self, in_filename):
     def read_eqdsk(self, in_filename):
         '''
         Function for reading the equdsk variables given a filename
-
         Parameters
         ----------
         in_filename : string
             A string with the equdsk filename to load
-
         Raises
         ------
         IOError
             Raises error if the array is not loaded.
-
         Returns
         -------
         None
@@ -412,11 +218,9 @@ def read_2d(nr, nz, name="Unknown"):
     def get_midplane_lcfs(self, psi_p=1.0001):
         '''
         Function for getting the inner and outer radial position of the LCFS at the midplane
-
         input: self,  a reference to the object itself
                psi_p, the flux surface of the LCFS, standard is psi_p = 1.005
                (otherwise the field-line is located inside the LCFS)
-
         return: Rcross, a list with the outer and inner radial position of the mid-plane LCFS
         '''
 
@@ -456,4 +260,3 @@ def file_numbers(fp):
         toklist = re.findall(pattern, line)
         for tok in toklist:
             yield tok
-

vita/modules/projection/projection2D/interp_div.py

@@ -0,0 +1,46 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Dec 10 14:24:20 2020
+
+@author: jmbols
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+def get_div_coords(divertor_coords):
+    # Define the 1D polynomial that represents the divertor
+    for i in range(len(divertor_coords[0, :])-1):
+        divertor_x = [divertor_coords[0, i], divertor_coords[0, i+1]]
+        divertor_y = [divertor_coords[1, i], divertor_coords[1, i+1]]
+
+        divertor_func = np.polyfit(divertor_x, divertor_y, 1)
+        divertor_polyfit = np.poly1d(divertor_func)
+
+        # Define a 1D polynomial for a surface just above the divertor (to be used
+        # for evaluating the angle of incidence)
+        divertor_func_above = [divertor_func[0], divertor_func[1] + 0.001]
+        divertor_polyfit_above = np.poly1d(divertor_func_above)
+
+        # Define a 1D polynomial for a surface just below the divertor (to be used
+        # for evaluating the angle of incidence)
+        divertor_func_below = [divertor_func[0], divertor_func[1] - 0.001]
+        divertor_polyfit_below = np.poly1d(divertor_func_below)
+
+def get_z_from_r(r):
+    z = r
+    return z
+
+if __name__ == '__main__':
+    divertor_coords_x = [0.252, 0.332, 0.332, 0.368]
+    divertor_coords_y = [-0.426, -0.621, -0.710, -0.806]
+    divertor_coords = np.array([divertor_coords_x, divertor_coords_y])
+    #divertor_coords_x = [0.226, 0.314, 0.365, 0.406]
+    #divertor_coords_y = [-0.425, -0.621, -0.708, -0.811]
+
+    get_div_coords(divertor_coords)
+
+    plt.plot(divertor_coords_x, divertor_coords_y)
+
+    r = 5
+    function_z_from_r = get_z_from_r(r)

vita/modules/projection/projection2D/psi_map_projection.py

@@ -58,7 +58,7 @@ def _shooting_algorithm(x_lims, function, function_z_from_r,
     '''
     x_lower = x_lims[0]
     x_upper = x_lims[1]
-    for _ in range(n_max):
+    for i in range(n_max):
         x_mid = (x_lower + x_upper)/2
         psi = function(x_mid, function_z_from_r(x_mid))
         if abs(psi - psi_target) > tol:
@@ -74,6 +74,8 @@ def _shooting_algorithm(x_lims, function, function_z_from_r,
                     x_lower = x_mid
         else:
             break
+    if i==49:
+        print("Warning: map_psi_omp_to_divertor._shooting_algorithm\n Convergence not reached")
     return x_mid
 
 def _flux_expansion(b_pol, points_omp, points_div):

vita/modules/sol_heat_flux/eich/eich.py

@@ -19,7 +19,7 @@ class Eich(HeatLoad):
         calculate_heat_flux_density: Function for calculating the heat-flux given
         a set of input parameters.
     '''
-    def __init__(self, lambda_q=1.5e-3, S=0., r0_lfs=0., r0_hfs=0.):
+    def __init__(self, lambda_q=1.5e-3, S=0., r0_lfs=0., r0_hfs=0., q_0 = 1.):
         '''
         Inputs for the class are
 
@@ -42,7 +42,7 @@ def __init__(self, lambda_q=1.5e-3, S=0., r0_lfs=0., r0_hfs=0.):
         HeatLoad.__init__(self, r0_lfs, r0_hfs)
         self.lambda_q = lambda_q
         self.S = S
-        self.q_0 = 1.
+        self.q_0 = q_0
         self.model_type = "Eich"
 
     def calculate_heat_flux_density(self, where="lfs-mp"):

vita/modules/sol_heat_flux/mid_plane_heat_flux.py

@@ -157,8 +157,9 @@ def plot_heat_power_density(self):
 
         '''
         plt.plot(self._s, self._q)
-        plt.xlabel('$s$')
-        plt.ylabel('$q(s)$')
+        plt.xlabel('$s$ - Distance (m)')
+        plt.ylabel('$q(s)$ - Normalised power')
+        plt.title('Outer mid-plane heat profile')
 
         imageFile = getOption('imageFile')
         if imageFile :