From 8b37602149d8ba465e1e953ccdf225be35706c37 Mon Sep 17 00:00:00 2001 From: mattngc Date: Tue, 15 Dec 2020 16:34:37 +0000 Subject: [PATCH] Manual merge in of Jeppe's fixes that don't involve the interface surface. Will handle these separately at the same time as refactoring API interface. --- vita/modules/equilibrium/equdsk/equdsk.py | 199 +----------------- .../projection/projection2D/interp_div.py | 46 ++++ .../projection2D/psi_map_projection.py | 4 +- vita/modules/sol_heat_flux/eich/eich.py | 4 +- .../sol_heat_flux/mid_plane_heat_flux.py | 5 +- 5 files changed, 55 insertions(+), 203 deletions(-) create mode 100644 vita/modules/projection/projection2D/interp_div.py diff --git a/vita/modules/equilibrium/equdsk/equdsk.py b/vita/modules/equilibrium/equdsk/equdsk.py index c40a100..85fe05a 100644 --- a/vita/modules/equilibrium/equdsk/equdsk.py +++ b/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 - diff --git a/vita/modules/projection/projection2D/interp_div.py b/vita/modules/projection/projection2D/interp_div.py new file mode 100644 index 0000000..18e3560 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/vita/modules/projection/projection2D/psi_map_projection.py b/vita/modules/projection/projection2D/psi_map_projection.py index 8df5908..66b2202 100644 --- a/vita/modules/projection/projection2D/psi_map_projection.py +++ b/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): diff --git a/vita/modules/sol_heat_flux/eich/eich.py b/vita/modules/sol_heat_flux/eich/eich.py index 42e1d2c..c89e822 100644 --- a/vita/modules/sol_heat_flux/eich/eich.py +++ b/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"): diff --git a/vita/modules/sol_heat_flux/mid_plane_heat_flux.py b/vita/modules/sol_heat_flux/mid_plane_heat_flux.py index 45916c1..fa564ac 100644 --- a/vita/modules/sol_heat_flux/mid_plane_heat_flux.py +++ b/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 :