python git

Scripts_Matthieu/AnalysisFunctions.py

+# Functions to extract the information from the fortran output and compare with analytical results
+
+import Constants
+import math
+import Utilities
+
+def extractValues(positions, drifts, efield, edescr, bdescr, ExB):
+    pos = []
+    enum = []
+    dnum = []
+    enumperp = []
+    dnumperp = []
+    ea = []
+    ba = []
+    da = []
+    for i in range(len(drifts)):
+        # Convert to float values to allow computations and interpretation
+        dn = [float(x) for x in drifts[i].split()]
+        dpn = [float(x) for x in positions[i].split()]
+        en = [float(x) for x in efield[i].split()]
+
+        # filter our boundary cells => These have zero field and drift and are not useful for interpretation!
+        if dn != [0,0,0] and not math.isnan(dn[0]):
+            dnum.append(dn)
+            pos.append(dpn)
+            enum.append(en)
+            # Compute analytical efield and drifts at these positions
+            efa = edescr(dpn[0], dpn[2], dpn[1], Constants.Rmajor)
+            bfa = bdescr(dpn[0], dpn[2], dpn[1], Constants.Rmajor, Constants.Rminor)
+            ea.append(efa)
+            ba.append(ba)
+            da.append(ExB(efa, bfa).tolist())
+	    dnumperp.append(transformCyltoRadPerp(dnum[i],pos[i]))
+	    enumperp.append(transformCyltoRadPerp(enum[i],pos[i]))
+    return [pos, enum, dnum, enumperp,dnumperp, ea , ba, da]
+
+def extractVectorInfo(numeric,analytical):
+    res_mag = []
+    an_mag = []
+    er_mag = []
+    er_vec = []
+    er_rel_vec = []
+    er_rel_mag = []
+    for i in range(len(numeric)):
+        res_mag.append(math.sqrt(sum([x**2 for x in numeric[i]])))
+        an_mag.append(math.sqrt(sum([x**2 for x in analytical[i]])))
+        er_mag.append(math.sqrt((res_mag[i]-an_mag[i])**2))
+        er_vec.append(Utilities.listdif(numeric[i],analytical[i]))
+        er_rel_vec.append([math.sqrt(x**2)/an_mag[i] for x in er_vec[i]])
+        er_rel_mag.append(er_mag[i]/an_mag[i])
+
+    return [res_mag,an_mag,er_mag,er_vec,er_rel_mag,er_rel_vec]
+
+def slice(positions,info,direction, coordinate):
+    # Get a slice of the full range of information along a certain direction
+    # positions are the coordinates corresponding to the info values
+    # direction: between 0 and 2: along which direction to take a slice!
+    # coordinate specifies which surface to keep as slice
+    pos_slice = []
+    info_slice = []
+    index1 = direction-1 %3
+    index2 = direction+1 %3
+    for i in range(len(positions)):
+        if positions[i][direction] == coordinate:
+            pos_slice.append([positions[i][index1],positions[i][index2]])
+            info_slice.append(info[i])
+
+    return [pos_slice, info_slice]
+
+def getAllSlices(positions,allInfo,direction,coordinate):
+    pos_slice = []
+    result = []
+    for i in range(len(allInfo)):
+        [pos_slice, info_slice] = slice(positions,allInfo[i],direction,coordinate)
+        result.append(info_slice)
+    result.append(pos_slice)
+    return result
+
+def getAllInfo(positions,drifts,efield,edescr,bdescr,ExB):
+    [pos, enum, dnum,enumperp,dnumperp, ea, ba, da] = extractValues(positions, drifts, efield, edescr, bdescr, ExB)
+
+    [dres_mag, dan_mag, der_mag, der_vec, der_rel_mag, der_rel_vec] =  extractVectorInfo(dnum,da)
+    [eres_mag, ean_mag, eer_mag, eer_vec, eer_rel_mag, eer_rel_vec] = extractVectorInfo(enum, ea)
+
+    torsurf = pos[0][1]
+
+    torSliceDict = {}
+    keys = ['drifts', 'efield','drifts_perp','efield_perp' 'ef_mag', 'ef_an_mag', 'drift_mag', 'drift_an_mag', 'ef_vec_er', 'ef_mag_er',
+            'drift_vec_er', 'drift_er_mag',
+            'efield_er_rel_mag', 'efield_er_rel_vec', 'drift_er_rel_mag', 'drift_er_rel_vec', 'pos']
+    #torslices = getAllSlices(pos,[dres_mag, dan_mag, der_mag, der_vec, der_rel_mag, der_rel_vec,eres_mag, ean_mag, eer_mag, eer_vec, eer_rel_mag, eer_rel_vec],1,torsurf)
+
+    values = [dnum,enum,dnumperp,enumperp,eres_mag,ean_mag,dres_mag,dan_mag,eer_vec,eer_mag,der_vec,der_mag,eer_rel_mag,eer_rel_vec,der_rel_mag,der_rel_vec,pos]
+    for i in range(len(values)):
+        torSliceDict[keys[i]] = values[i]
+
+    return torSliceDict
+
+
+def readText(filepath):
+    file = open(filepath,'r')
+    result = []
+    for line in file:
+        result.append(line)
+    return result
+    
+    
+def transformCyltoRadPerp(v,pos):
+    result = v
+    phi = math.atan2(pos[2],pos[0])
+    result[0] = v[0]*math.cos(phi)+v[2]*math.sin(phi)
+    result[2] = v[2]*math.cos(phi)-v[0]*math.sin(phi)
+    return result
+
+
+
+
+

Scripts_Matthieu/AnalyticalPrescriptions.py

+# This file contains the analytical functions describing potential, efield and magnetic field for different cases
+import math
+# Magnetic Fields
+def b_field_description1(r,z,theta, Rmajor, Rminor):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
+    bval = 0.01/R*100  # R in cm
+    return [0,bval,0]
+def bfieldstrength1(r,z,theta,Rmajor,Rminor):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
+    bval = 0.01/R*100  # R in cm
+    return bval
+
+def b_field_description2(r,z,theta, Rmajor, Rminor):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
+    bval = 0.01/R*100  # R in cm
+    return [0,bval,0]
+
+# Potentials
+def potential_description1(r,z,theta,Rmajor,Rminor):
+    r = r/100
+    Rmajor = Rmajor/100
+    z = z/100
+    D = math.sqrt((r-Rmajor)**2+z**2)
+    potval = math.exp(-D)*500
+    return potval
+
+def potential_description2(r,z,theta,Rmajor,Rminor):
+    R = (r-Rmajor)/100+0.001
+    potval = math.exp(-R)*100
+    return potval
+
+# Electric Fields
+def e_descr1(r,z,theta, Rmajor):
+    r = r/100
+    Rmajor = Rmajor/100
+    z = z/100
+    D = math.sqrt((r - Rmajor) ** 2 + z ** 2)
+    Er = -500*math.exp(-D)*1/D*(r-Rmajor)
+    Ephi = 0
+    Ez = -500*math.exp(-D)*1/D*z
+    return [Er,Ephi, Ez]
+
+def e_descr2(r,z,theta, Rmajor):
+    D = (r - Rmajor)/100+0.001
+    Er = -100*math.exp(-D)
+    Ephi = 0
+    Ez = 0
+    return [Er,Ephi, Ez]
\ No newline at end of file

Scripts_Matthieu/AnalyzeNumericalResults.py

+import math
+import driftFunctions
+import AnalyticalPrescriptions
+import Constants
+import Utilities
+#import pandas as pd
+import matplotlib.pyplot as plt
+import PlottingFunctions
+import AnalysisFunctions
+
+# Constants
+drifts_file = open('/u/mbj/Drifts/data-for-drift-computations/Description1/Case50/OutputGG/ExB_RESULTS','r')
+driftpos_file = open('/u/mbj/Drifts/data-for-drift-computations/Description1/Case50/OutputGG/DRIFT_POS','r')
+efield_file = open('/u/mbj/Drifts/data-for-drift-computations/Description1/Case50/OutputGG/EFIELD','r')
+drifts = []
+drift_pos = []
+efield = []
+
+# Read numeric drift values
+for line in drifts_file:
+    drifts.append(line)
+
+# Read drift positions
+for line in driftpos_file:
+    drift_pos.append(line)
+
+# Electric field
+for line in efield_file:
+    efield.append(line)
+
+# Convert values from fortran to float, filter boundary cells and compute analytical profiles for cell center locations
+[drift_pos_num, efield_num, drifts_num, efield_an, bfield_an, drifts_an] = AnalysisFunctions.extractValues(drift_pos,drifts, efield,AnalyticalPrescriptions.e_descr1,AnalyticalPrescriptions.b_field_description1,driftFunctions.ExBdrift)
+
+# Compute vectors and magnitude + errors for efield and drift
+[drift_res_mag,drift_an_mag,drift_er_mag,drift_er_vec,drift_er_rel_mag,drift_er_rel_vec] = AnalysisFunctions.extractVectorInfo(drifts_num,drifts_an)
+[efield_res_mag,efield_an_mag,efield_er_mag,efield_er_vec,efield_er_rel_mag,efield_er_rel_vec] = AnalysisFunctions.extractVectorInfo(efield_num,efield_an)
+
+# Take slices in different directions to allow for extra analysis
+torsurf = drift_pos_num[0][1]
+radsurf = drift_pos_num[150][0] # Based on inspection
+
+# Create a dictionary containing the different slices
+torSliceDict = {}
+keys = ['drifts','efield','drifts_perp','efield_perp','ef_mag','ef_an_mag','drift_mag','drift_an_mag','ef_vec_er','ef_mag_er','drift_vec_er','drift_er_mag',
+ 'efield_er_rel_mag','efield_er_rel_vec','drift_er_rel_mag','drift_er_rel_vec','pos']
+torslices = AnalysisFunctions.getAllSlices(drift_pos_num,[drifts_num,efield_num,efield_res_mag,efield_an_mag,drift_res_mag,
+                                            drift_an_mag,efield_er_vec,efield_er_mag,drift_er_vec,drift_er_mag,efield_er_rel_mag,
+                                            efield_er_rel_vec,drift_er_rel_mag,drift_er_rel_vec],1,torsurf)
+for i in range(len(torslices)):
+    torSliceDict[keys[i]] = torslices[i]
+
+print(torSliceDict.get('efield_er_rel_mag'))
+print(sum(torSliceDict.get('efield_er_rel_mag')))
+print(len(torSliceDict.get('efield_er_rel_mag')))
+# Plot results
+
+x = [r[0] for r in torSliceDict.get('pos')]
+y = [r[1] for r in torSliceDict.get('pos')]
+er = [x[0] for x in torSliceDict.get('efield')]
+ez = [x[1] for x in torSliceDict.get('efield')]
+
+
+plot1 = plt.figure(1)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+plt.title('Magnitude of Electic Field')
+sc1 = plt.scatter(x, y, s=20, c=torSliceDict.get('ef_mag'), cmap='plasma')
+plt.colorbar(sc1)
+
+
+plot2 = plt.figure(2)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+plt.title('Theta component of Electric Field')
+sc2 = plt.scatter(x, y, s=20, c = [x[1] for x in torSliceDict.get('efield')], cmap='plasma')
+plt.colorbar(sc2)
+
+plot3 = plt.figure(3)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+plt.title('Relative magnitude error of Electric Field')
+sc3 = plt.scatter(x, y, s=20, c=torSliceDict.get('efield_er_rel_mag'), cmap='plasma')
+plt.colorbar(sc3)
+
+plot4 = plt.figure(4)
+plt.title('Potential')
+PlottingFunctions.plot2DContour(Constants.ns,Constants.ns,Constants.Rmajor-Constants.Rminor,Constants.Rmajor+Constants.Rminor,-Constants.Rminor,Constants.Rminor,Constants.Rmajor,Constants.Rminor,AnalyticalPrescriptions.potential_description1,'X [cm]', 'Y [cm]','Potential [V]',False)
+
+plot5 = plt.figure(5)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+plt.title('Magnitude of ExB Drift (Numerical)')
+sc5 = plt.scatter(x, y, s=20, c=torSliceDict.get('drift_mag'), cmap='plasma')
+plt.colorbar(sc5)
+
+plot6 = plt.figure(6)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+plt.title('Magnitude of ExB Drift (Analytical)')
+sc6 = plt.scatter(x, y, s=20, c=torSliceDict.get('drift_an_mag'), cmap='plasma')
+plt.colorbar(sc6)
+
+plot7 = plt.figure(7)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+plt.title('Relative magnitude error of ExB Drift')
+sc7 = plt.scatter(x, y, s=20, c=torSliceDict.get('drift_er_rel_mag'), cmap='plasma')
+plt.colorbar(sc7)
+
+plot8 = plt.figure(8)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+plt.title('Magnitude of Efield (Rad)')
+sc8 = plt.scatter(x, y, s=20, c= [x[0] for x in torSliceDict.get('efield_perp')], cmap='plasma')
+plt.colorbar(sc8)
+
+plot9 = plt.figure(9)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+plt.title('Magnitude of Efield (perp)')
+sc9 = plt.scatter(x, y, s=20, c= [x[2] for x in torSliceDict.get('efield_perp')], cmap='plasma')
+plt.colorbar(sc9)
+
+
+plt.show()
+

Scripts_Matthieu/AnalyzeResults.py

+# Read the results from the EMC3 computation of drifts and compare to analytical
+import math
+from driftFunctions import *
+drifts_file = open('C:/Users/matth/github/EMC3/TestCase/DRIFTS','r')
+driftpos_file = open('C:/Users/matth/github/EMC3/TestCase/DRIFT_POS','r')
+efield_file = open('C:/Users/matth/github/EMC3/TestCase/EFIELD','r')
+drifts = []
+drift_pos = []
+efield = []
+# Read numeric drift values
+for line in drifts_file:
+    drifts.append(line)
+
+# Read drift positions
+for line in driftpos_file:
+    drift_pos.append(line)
+
+# Electric field
+for line in efield_file:
+    efield.append(line)
+
+print('length',len(drifts))
+print(drift_pos)
+first_drift = drifts[250]
+second_drift = drifts[150]
+first_pos = drift_pos[250]
+second_pos = drift_pos[150]
+first_el = efield[250]
+second_el = efield[150]
+
+fd = [float(x) for x in first_drift.split()]
+fp = [float(x) for x in first_pos.split()]
+fe = [float(x) for x in first_el.split()]
+sd = [float(x) for x in second_drift.split()]
+sp = [float(x) for x in second_pos.split()]
+se = [float(x) for x in second_el.split()]
+
+print('fp', fp)
+print('sp', sp)
+print('fd',fd)
+print('sd', sd)
+print('fe', fe)
+print('se', se)
+
+# Compute drift obtained analytically for the positions corresponding to the numerical drifts
+# Analytical prescription for the potential should correspond with the description used to evaluate potential at
+# plasma cell centers => match fortran and python!!!
+def b_field_description(r,z,theta, Rmajor, Rminor):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
+    bval = 0.01/R*100  # R in cm
+    return [0,bval,0]
+def b_field_description(r,z,theta, Rmajor, Rminor):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
+    bval = 0.01/R*100  # R in cm
+    return [0,bval,0]
+# Create the potential data
+# At the moment we use the same grid for magnetic and plasma cells
+def potential_description(r,z,theta,Rmajor,Rminor):
+    r = r/100
+    Rmajor = Rmajor/100
+    z = z/100
+    D = math.sqrt((r-Rmajor)**2+z**2)+0.001
+    potval = math.exp(-D)*500
+    return [potval]
+def potential_description(r,z,theta,Rmajor,Rminor):
+    R = (r-Rmajor)/100+0.001
+    potval = math.exp(-R)*100
+    return potval
+
+def e_descr(r,z,theta, Rmajor):
+    r = r/100
+    Rmajor = Rmajor/100
+    z = z/100
+    D = math.sqrt((r - Rmajor) ** 2 + z ** 2)
+    Er = -500*math.exp(-D)*1/D*(r-Rmajor)
+    Ephi = 0
+    Ez = -500*math.exp(-D)*1/D*z
+    return [Er,Ephi, Ez]
+
+def e_descr(r,z,theta, Rmajor):
+    D = (r - Rmajor)/100+0.001
+    Er = -100*math.exp(-D)
+    Ephi = 0
+    Ez = 0
+    return [Er,Ephi, Ez]
+fda = ExBdrift(e_descr(fp[0],fp[2],fp[1],500),b_field_description(fp[0],fp[2],fp[1],500,100))
+print('fda',fda)
+
+sda = ExBdrift(e_descr(sp[0],sp[2],sp[1],500),b_field_description(sp[0],sp[2],sp[1],500,100))
+print('sda',sda)
+
+
+# Quick checks
+pos = [  553.255879500000,0.104719755119660,8.43490275000000 ]
+E = e_descr(pos[0],pos[2],pos[1],500)
+B = b_field_description(pos[0],pos[2],pos[1],500,100)
+D = ExBdrift(E,B)
+potcheck = potential_description(562.81309, 0.1047197551199660, 20.40921, 500,100)
+print(potcheck)
+print("ALL", E, B, D)
\ No newline at end of file

Scripts_Matthieu/Constants.py

+# This file contains constant values, to be used in other files
+Rmajor = 500
+Rminor = 100
+ns = 100
\ No newline at end of file

Scripts_Matthieu/CreateMesh.py

+# This script calls python functions allowing to create a mesh for some simple computations
+
+# Import the required functions from "MeshFunctions"
+from MeshFunctions import *
+
+[tor,rad,hor] = create_toroidal_vertices(6,5,10,6,500,100)
+writeGEOMETRY3DDATA(tor,rad,hor)
+
+imposeAnalytic(tor,rad,hor,b_field_description,"BFIELDSTRENGTH")
+

Scripts_Matthieu/ErrorConvergence.py

+# Compute an error measure for different mesh resolutions and compare with expected decrease
+import AnalysisFunctions
+import driftFunctions
+import AnalyticalPrescriptions
+import matplotlib.pyplot as plt
+
+LSQfiles = []
+LSQfiles.append('TestCase2010050')
+LSQfiles.append('TestCase105020')
+LSQfiles.append('TestCase5500500')
+LSQfiles.append('TestCase10100100')
+
+GGfiles = []
+GGfiles.append('TestCase2010050')
+GGfiles.append('TestCase105020')
+GGfiles.append('TestCase5500500')
+GGfiles.append('TestCase10100100')
+
+errorsLSQ = []
+nbCellsLSQ = []
+errorsGG = []
+nbCellsGG = []
+path = '/u/mbj/Drifts/data-for-drift-computations/'
+
+for i in range(len(LSQfiles)):
+    torSliceDict = AnalysisFunctions.getAllInfo(AnalysisFunctions.readText(path + LSQfiles[i] + '/Output/DRIFT_POS'),
+                                                AnalysisFunctions.readText(path+LSQfiles[i]+'/Output/ExB_RESULTS'),
+                                                AnalysisFunctions.readText(path + LSQfiles[i] + '/Output/EFIELD'),
+                                                AnalyticalPrescriptions.e_descr2, AnalyticalPrescriptions.b_field_description2,
+                                                driftFunctions.ExBdrift)
+    relerrors = torSliceDict.get('efield_er_rel_mag')
+    errorsLSQ.append(sum(relerrors)/len(relerrors))
+    nbCellsLSQ.append(len(relerrors))
+    print(i)
+    print(sum(relerrors))
+    print(len(relerrors))
+    print(sum(relerrors)/len(relerrors))
+
+
+for i in range(len(GGfiles)):
+    torSliceDict = AnalysisFunctions.getAllInfo(AnalysisFunctions.readText(path + GGfiles[i] + '/OutputLSQ/DRIFT_POS'),
+                                                AnalysisFunctions.readText(path+GGfiles[i]+'/OutputLSQ/ExB_RESULTS'),
+                                                AnalysisFunctions.readText(path + GGfiles[i] + '/OutputLSQ/EFIELD'),
+                                                AnalyticalPrescriptions.e_descr2, AnalyticalPrescriptions.b_field_description2,
+                                                driftFunctions.ExBdrift)
+    relerrors = torSliceDict.get('efield_er_rel_mag')
+    errorsGG.append(sum(relerrors)/len(relerrors))
+    nbCellsGG.append(len(relerrors))
+    print(i)
+    print(sum(relerrors))
+    print(len(relerrors))
+    print(sum(relerrors)/len(relerrors))
+
+# For both methods, plot the error as a function of the number of cells used
+plot1 = plt.figure(1)
+sc1 = plt.scatter(nbCellsGG, errorsGG)
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+plot2 = plt.figure(2)
+sc2 = plt.scatter(nbCellsLSQ,errorsLSQ)
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+
+
+
+plt.show()

Scripts_Matthieu/MeshFunctions.py

+import numpy as np
+import math
+import os
+
+def generate_vertices(nb_tor,nb_rad,nb_pol,toroidal,):
+
+    arr = np.array([1, 2, 3, 4, 5])
+
+    toroidal
+
+
+
+    # Write the required geometry information to a text file
+    gd = open("GEOMETRY_3D_DATA","w+")
+    gd.write("teststring")
+
+    return
+
+
+def create_toroidal_vertices(nb_tor,nb_rad,nb_pol,toroidal_range, R, Rmin):
+    # nb_tor, nb_rad, nb_pol are the number of surfaces in each direction
+    # R (major radius) and a (minor radius) give the dimension of the torus
+    # Create array containing angles of toroidal surfaces
+    toroidal_surfaces = np.linspace(0, toroidal_range, nb_tor)  # Include first and last surface
+    # Initialize arrays for the radial and horizontal positions
+    # Radial position R = R0 (major radius) + a (minor radius)*cos(phi)
+    radial_positions = np.empty(shape=(nb_tor,nb_rad*nb_pol),dtype='float')
+    horizontal_positions = np.empty(shape=(nb_tor,nb_rad*nb_pol),dtype='float')
+    r = np.linspace(Rmin/2,Rmin,nb_rad)
+    for i in range(nb_tor):
+        for j in range(nb_pol):
+            for k in range(nb_rad):
+                radial_positions[i,j*nb_rad+k] = R+r[k]*math.cos(j*2*math.pi/nb_pol)
+                horizontal_positions[i,j*nb_rad+k] = r[k]*math.sin(j*2*math.pi/nb_pol)
+    return [toroidal_surfaces,radial_positions,horizontal_positions]
+
+
+def writeGEOMETRY3DDATA(tor_surf,rad_pos,hor_pos, ntheta, nr, nz, path):
+    # Open file to write the geometry information
+    filepath = path+"\\GEOMETRY_3D_DATA"
+    gd = open(filepath, "w+")
+    gd.write(str(nr)+"   ")
+    gd.write(str(nz)+"   ")
+    gd.write(str(ntheta)+"   ")
+    gd.write("\n")
+    for i in range(len(tor_surf)):
+        gd.write(str(tor_surf[i]))
+        gd.write('\n')
+        for j in range(math.ceil(len(rad_pos[i])/6)):
+            for k in range(min(len(rad_pos[i,j*6:]),6)):
+                fs = "{:.6f}".format(rad_pos[i,j*6+k])
+                gd.write(str(fs+'  '))
+            gd.write('\n')
+        for j in range(math.ceil(len(hor_pos[i])/6)):
+            for k in range(min(len(hor_pos[i,j*6:]),6)):
+                fs = "{:.6f}".format(hor_pos[i,j*6+k])
+                gd.write(str(fs+'  '))
+            gd.write('\n')
+
+    return
+
+
+def imposeAnalytic(tor_surf,rad_pos,hor_pos,analytic_description,Rmajor, Rminor, path):
+    bfs = open(path, "w+")
+    for i in range(len(tor_surf)):
+        for j in range(math.ceil(len(rad_pos[i])/10)):
+            for k in range(min(len(rad_pos[i,j*10:]),10)):
+                value = analytic_description(rad_pos[i,j*10+k],hor_pos[i,j*10+k],tor_surf[i],Rmajor,Rminor)
+                fs = "{:.6f}".format(value)
+                bfs.write(str(fs+'  '))
+            bfs.write('\n')
+    return
+
+def b_field_description(r,z,theta):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-500)**2 + z**2)+0.01+0*theta
+    bval = 0.01/R*100  # R in cm
+    return bval
+
+def Bfieldvector(r,z,theta):
+    R = math.sqrt((r - 500) ** 2 + z ** 2) + 0.01
+    Br = 0
+    Bz = math.cos(theta/180*1/10*math.pi)*0.01
+    Btheta = 0.01/R*100
+    return [Br,Bz,Btheta]
+
+def CreateGeneralInfo(filepath, nz,r,z,theta):
+    # nz is the number of zones considered
+    # r,z,theta are the number of surfaces in all directions
+    # They should be lists with length nz
+    gi = open(filepath,"w+")
+    gi.write("* input file for EMC3 code")
+    gi.write('\n')
+    gi.write("*** 1. fein grid mesh representing magnetic coordinates")
+    gi.write('\n')
+    for i in range(nz):
+        gi.write(str(i+1))
+        gi.write("\n")
+        gi.write(str(r[i])+'   ')
+        gi.write(str(z[i])+'   ')
+        gi.write(str(theta[i])+'   ')
+        gi.write("\n")
+        gi.write("*** 2. Non default, Non Transparant, Plate surface")
+        gi.write("\n")
+        gi.write("*** 2.1 Non default")
+        gi.write("\n")
+        gi.write("* radial")
+        gi.write("\n")
+        gi.write("0")
+        gi.write("\n")
+        gi.write("* poloidal")
+        gi.write("\n")
+        gi.write("2")
+        gi.write("\n")
+        gi.write("0  "+str(i)+"  1")
+        gi.write("\n")
+        gi.write("0  "+str(r[i]-2)+"  0  "+str(theta[i]-2))
+        gi.write("\n")
+        gi.write(str(z[i]-1)+"  "+str(i)+"  1")
+        gi.write("\n")
+        gi.write("0  "+str(r[i]-2)+"  0  "+str(theta[i]-2))
+        gi.write("\n")
+        gi.write("* toroidal")
+        gi.write("\n")
+        gi.write("2")
+        gi.write("\n")
+        gi.write("0  "+str(i)+"  3")
+        gi.write("\n")
+        gi.write("0  "+str(r[i]-2)+"  0  "+str(z[i]-2))
+        gi.write("\n")
+        gi.write(str(theta[i]-1)+"  "+str(i)+"  3")
+        gi.write("\n")
+        gi.write("0  "+str(r[i]-2)+"  0  "+str(z[i]-2))
+        gi.write("\n")
+        gi.write("*** 2.2 Non transparant")
+        gi.write("\n")
+        gi.write("* radial")
+        gi.write("\n")
+        gi.write("2")
+        gi.write("\n")
+        gi.write("3  "+str(i)+"  1")
+        gi.write("\n")
+        gi.write("0  "+str(z[i]-2)+"  0  "+str(theta[i]-2))
+        gi.write("\n")
+        gi.write(str(r[i]-2)+"  "+str(i)+"  -1")
+        gi.write("\n")
+        gi.write("0  "+str(z[i]-2)+"  0  "+str(theta[i]-2))
+        gi.write("\n")
+        gi.write("* poloidal")
+        gi.write("\n")
+        gi.write("0")
+        gi.write("\n")
+        gi.write("* toroidal")
+        gi.write("\n")
+        gi.write("0")
+        gi.write("\n")
+        gi.write("*** 2.3 Plate surface")
+        gi.write("\n")
+        gi.write("* radial")
+        gi.write("\n")
+        gi.write("0")
+        gi.write("\n")
+        gi.write("* poloidal")
+        gi.write("\n")
+        gi.write("0")
+        gi.write("\n")
+        gi.write("* toroidal")
+        gi.write("\n")
+        gi.write("2")
+        gi.write("\n")
+        gi.write(str(theta[i]-2)+"  "+str(i)+ "  1  ")
+        gi.write("\n")
+        gi.write(str(r[i]-3)+"  "+str(r[i]-2)+"  0  "+str(z[i]-2) )
+        gi.write("\n")
+        gi.write(str(theta[i]-2)+"  "+str(i)+ "  -1  ")
+        gi.write("\n")
+        gi.write(str(r[i]-3)+"  "+str(r[i]-2)+"  0  "+str(z[i]-2))
+        gi.write("\n")
+    gi.write("*** 3. Physical cell")
+    gi.write("\n")
+    gi.write("10  1") # Geometric cell = Physical cell
+    gi.write("\n")
+    gi.write("* check cell ?")
+    gi.write("\n")
+    gi.write("F")
+    return None
+
+
+
+
+
+
+

Scripts_Matthieu/NumericalAnalysisDescription1.py

+import math
+import driftFunctions
+import AnalyticalPrescriptions
+import Constants
+import Utilities
+#import pandas as pd
+import matplotlib.pyplot as plt
+import PlottingFunctions
+import AnalysisFunctions
+
+# Constants
+path = '/u/mbj/Drifts/data-for-drift-computations/Description1/'
+
+LSQfiles = []
+LSQfiles.append('Case50')
+LSQfiles.append('Case100')
+LSQfiles.append('Case200')
+LSQfiles.append('Case400')
+LSQfiles.append('Case1000')
+
+
+
+errorsLSQ = []
+nbCellsLSQ = []
+errorsGG = []
+nbCellsGG = []
+
+for i in range(len(LSQfiles)):
+    torSliceDict = AnalysisFunctions.getAllInfo(AnalysisFunctions.readText(path + LSQfiles[i] + '/OutputGG/DRIFT_POS'),
+                                                AnalysisFunctions.readText(path+LSQfiles[i]+'/OutputGG/ExB_RESULTS'),
+                                                AnalysisFunctions.readText(path + LSQfiles[i] + '/OutputGG/EFIELD'),
+                                                AnalyticalPrescriptions.e_descr1, AnalyticalPrescriptions.b_field_description1,
+                                                driftFunctions.ExBdrift)
+    relerrors = torSliceDict.get('efield_er_rel_mag')
+    y = len(relerrors)
+    errorsLSQ.append(sum([x/y for x in relerrors]))
+    nbCellsLSQ.append(len(relerrors))
+    print(i)
+    print(sum(relerrors))
+    print(len(relerrors))
+    print(sum(relerrors)/len(relerrors))
+
+print(errorsLSQ)
+print(nbCellsLSQ)
+# For both methods, plot the error as a function of the number of cells used
+#plot1 = plt.figure(1)
+#sc1 = plt.scatter(nbCellsGG, errorsGG)
+#ax = plt.gca()
+#ax.set_yscale('log')
+#ax.set_xscale('log')
+plot2 = plt.figure(2)
+sc2 = plt.scatter(nbCellsLSQ,errorsLSQ)
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+
+
+
+plt.show()

Scripts_Matthieu/PlotCase1.py

+from mpl_toolkits import mplot3d
+import numpy as nump
+import matplotlib
+import matplotlib.pyplot as plt
+from MeshFunctions import *
+import math
+
+from driftFunctions import *
+
+
+# General Info about the case to plot
+nt =1; nr = 10; np = 10
+Rmajor = 500
+rminor = 100
+Trange = 36
+# Create the vertices
+[tor,rad,hor] = create_toroidal_vertices(nt,nr,np,Trange,Rmajor,rminor)
+[xlines_pol,xlines_rad,xlines_tor,ylines_pol,ylines_rad,ylines_tor,zlines_pol,zlines_rad,zlines_tor] = plotMesh(tor,rad,hor,nt,nr,np,None,show=True)
+
+
+def b_field_description(r,z,theta, Rmajor, Rminor):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
+    R = 0.01*R
+    bval = math.exp(-R) # R in cm
+    return bval
+x,z,Bf = plot2DContour(100,100,Rmajor-rminor,Rmajor+rminor,-rminor,rminor,Rmajor,rminor,b_field_description,'','','',False)
+
+[xdata,ydata,zdata] = plotGridPoints(tor,rad,hor,show = False)
+Bdata = []
+for i in range(len(xdata)):
+    Bdata.append(b_field_description(xdata[i],zdata[i],0,Rmajor,rminor))
+
+# Evaluate the electric field at all points:
+def potential(r,z,theta,Rmaor,rminor):
+    R = math.sqrt((r - Rmajor) ** 2 + z ** 2)
+    return Rmajor/R
+def e_field_description(r,z,theta, Rmajor, Rminor):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-Rmajor)**2 + z**2)
+    # Potential as 1/r => Efield: -1/r^2
+    er = -(r-Rmajor)/(((r-Rmajor)**2+z**2))**(3/2)*Rmajor
+    etheta = 0
+    ez = -z/Rmajor/(((r-Rmajor)**2+z**2)/Rmajor)**(3/2)*Rmajor
+    return [er, etheta, ez]
+
+def gradB(r,z,theta,Rmajor,rminor):
+    gbr = -(r-Rmajor)/(Rmajor*math.sqrt(z**2+(r-Rmajor)**2))*math.exp(-math.sqrt(z**2+(r-Rmajor)**2)/Rmajor)
+    gbz = -z/(Rmajor*math.sqrt(z**2+(r-Rmajor)**2))*math.exp(-math.sqrt(z**2+(r-Rmajor)**2)/Rmajor)
+    return [gbr, 0, gbz]
+
+def bvector(r,z,theta, Rmajor, Rminor):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
+    R = R/Rmajor
+    btheta = math.exp(-R) # R in cm
+    return [0,btheta,0]
+
+# Get two points, where the drifts (exb or gradb) will be computed and the vector plotted
+p1 = [xlines_tor[math.floor(np/2)*nr-1][math.floor(nt/3)],ylines_tor[math.floor(np/2)*nr-1][math.floor(nt/3)],zlines_tor[math.floor(np/2)*nr-1][math.floor(nt/3)]]
+p2 = [xlines_tor[math.floor(np/4)*nr-1][math.floor(2*nt/3)],ylines_tor[math.floor(np/4)*nr-1][math.floor(2*nt/3)],zlines_tor[math.floor(np/4)*nr-1][math.floor(2*nt/3)]]
+r1 = math.sqrt(p1[0]**2+p1[1]**2)
+theta1 = math.atan2(p1[1],p1[0])
+z1 = p1[2]
+r2 = math.sqrt(p2[0]**2+p2[1]**2)
+theta2 = math.atan2(p2[1],p2[0])
+z2 = p2[2]
+
+E1 = e_field_description(r1,z1,theta1,Rmajor,rminor)
+B1 = bvector(r1,z1,theta1,Rmajor,rminor)
+B2 = bvector(r2,z2,theta2,Rmajor,rminor)
+gradB2 = gradB(r2,z2,theta2,Rmajor,rminor)
+vexb = ExBdrift(E1,B1)
+vgradb = gradBdrift(B2,gradB2)
+# Plot only the first outer surface and 2 toroidal lines
+fig = plt.figure()
+ax = plt.axes(projection='3d')
+scale = 100
+E1 = cyltocart(E1,p1)
+B1 = cyltocart(B1,p1)
+B2 = cyltocart(B2,p2)
+gradB2 = cyltocart(gradB2,p2)
+E1 = 10*scale*nump.array(E1)
+B1 = scale*nump.array(B1)
+B2 = scale*nump.array(B2)
+gradB2 = 100*scale*nump.array(gradB2)
+vexb = 10*scale*vexb
+vgradb = scale*vgradb
+
+ax.quiver(p1[0],p1[1],p1[2],E1[0], E1[1],E1[2],color = 'yellow',label='Electric Field Vector')
+ax.quiver(p1[0],p1[1],p1[2],B1[0], B1[1],B1[2],color = 'blue',label = 'Magnetic Field Vector')
+ax.quiver(p1[0],p1[1],p1[2],vexb[0],vexb[1],vexb[2],color = 'red',label = 'E x B Drift Velocity')
+ax.quiver(p2[0],p2[1],p2[2],gradB2[0], gradB2[1],gradB2[2],color = 'black',label = 'Magnetic Field Strength Gradient')
+ax.quiver(p2[0],p2[1],p2[2],B2[0], B2[1],B2[2],color = 'blue')
+ax.quiver(p2[0],p2[1],p2[2],vgradb[0],vgradb[1],vgradb[2],color = 'purple',label = 'Grad B Velocity')
+ax.plot(xlines_pol[nr-1], ylines_pol[nr-1], zlines_pol[nr-1],color='green')
+ax.plot(xlines_pol[(nt)*nr-1], ylines_pol[(nt)*nr-1], zlines_pol[(nt)*nr-1],color='green')
+ax.plot(xlines_tor[math.floor(np/4)*nr-1], ylines_tor[math.floor(np/4)*nr-1], zlines_tor[math.floor(np/4)*nr-1],color='green')
+ax.plot(xlines_tor[math.floor(np/2)*nr-1], ylines_tor[math.floor(np/2)*nr-1], zlines_tor[math.floor(np/2)*nr-1],color='green')
+#ax.contourf(z,x,Bf,100,zdir='z')
+ax.set_xlabel('X [cm]')
+ax.set_ylabel('Y [cm]')
+ax.set_zlabel('Z [cm]')
+ax.set_title('Test Case 1')
+ax.legend()
+plt.show()
+
+# Colorplot of magnetic field and potential
+# Plot 2d contour of Bfield
+ns = 100
+plot2DContour(ns,ns,Rmajor-rminor,Rmajor+rminor,-rminor,rminor,Rmajor,rminor,b_field_description,'X [cm]', 'Y [cm]','Magnetic field strength [T]',True)
+plot2DContour(ns,ns,Rmajor-rminor,Rmajor+rminor,-rminor,rminor,Rmajor,rminor,potential,'X [cm]', 'Y [cm]','Potential [V]',True)
\ No newline at end of file

Scripts_Matthieu/PlotCase3.py

+from mpl_toolkits import mplot3d
+import numpy as nump
+import matplotlib
+import matplotlib.pyplot as plt
+from MeshFunctions import *
+import math
+from PlottingFunctions import *
+from TracingFieldLines import *
+from plotSingleFluxTube import *
+from driftFunctions import *
+# General Info about the case to plot
+nt =10; nr = 5; np = 20
+Rmajor = 500
+rminor = 100
+Trange = 36
+# Create the vertices
+m2d = create_2D_mesh(nr,np,Rmajor,rminor,'toroidal')
+
+# Create the analytical description of the B field
+def bvectorcase3(r,theta,z, Rmajor, Rminor):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-Rmajor)**2 + z**2)
+    btheta = math.exp(-R*0.01)  # R in cm
+    br = -z/R*1.5
+    bz = (r-Rmajor)/R*1.5
+    return [br,btheta,bz]
+# Create the 3D mesh using the field lines
+precision = 100
+m3d = extend_mesh_3D(m2d,Trange,bvectorcase3,nt,precision, Rmajor, rminor)
+rad = m3d[0]
+hor = m3d[1]
+tor = m3d[2]
+field_lines = m3d[3]
+
+# Get the field lines in x,y,z
+lines = []
+for i in range(len(field_lines)):
+    field = []
+    x = []
+    y = []
+    z = []
+    for j in range(len(field_lines[i])):
+        x.append(field_lines[i][j][0] * math.cos(field_lines[i][j][2] * math.pi / 180))
+        y.append(field_lines[i][j][0] * math.sin(field_lines[i][j][2] * math.pi / 180))
+        z.append(field_lines[i][j][1])
+    field = [x,y,z]
+    lines.append(field)
+
+
+# Plot lines connecting the vertices in all different directions, this will be used to plot the mesh
+[xlines_pol,xlines_rad,xlines_tor,ylines_pol,ylines_rad,ylines_tor,zlines_pol,zlines_rad,zlines_tor] = plotMesh(tor,rad,hor,nt,nr,np,lines,show=True)
+
+# plot a single flux tube
+plotfluxTube(tor,rad,hor,nt,nr,np,None,True)
+
+#def b_field_description(r,z,theta, Rmajor, Rminor):
+#    # Radial B-field (linear and max at central surfaces of the torus
+#    R = math.sqrt((r-Rmajor)**2 + z**2)
+#    R = 0.01*R
+#    bval = math.exp(-R) # R in cm
+#    return bval
+#x,z,Bf = plot2DContour(100,100,Rmajor-rminor,Rmajor+rminor,-rminor,rminor,Rmajor,rminor,b_field_description,'','','',False)
+
+#[xdata,ydata,zdata] = plotGridPoints(tor,rad,hor,show = False)
+
+# Evaluate the electric field at all points:
+def e_field_description(r,theta,z, Rmajor, Rminor):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-Rmajor)**2 + z**2)
+    # Potential as 1/r => Efield: -1/r^2
+    er = -(r-Rmajor)/Rmajor/(((r-Rmajor)**2+z**2)/Rmajor)**(3/2)
+    etheta = 0
+    ez = -z/Rmajor/(((r-Rmajor)**2+z**2)/Rmajor)**(3/2)
+    return [er, etheta, ez]
+
+def gradB(r,theta,z,Rmajor,rminor):
+    gbr = -(r-Rmajor)/(Rmajor*math.sqrt(z**2+(r-Rmajor)**2))*math.exp(-math.sqrt(z**2+(r-Rmajor)**2)/Rmajor)
+    gbz = -z/(Rmajor*math.sqrt(z**2+(r-Rmajor)**2))*math.exp(-math.sqrt(z**2+(r-Rmajor)**2)/Rmajor)
+    return [gbr, 0, gbz]
+
+
+
+# Plot only the first outer surface and 2 toroidal lines
+#plt.figure()
+#ax = plt.figure().add_subplot(projection='3d')
+#ax.plot(xlines_pol[nr-1], ylines_pol[nr-1], zlines_pol[nr-1],color='green')
+#ax.plot(xlines_pol[(nt)*nr-1], ylines_pol[(nt)*nr-1], zlines_pol[(nt)*nr-1],color='green')
+#ax.plot(xlines_tor[math.floor(np/4)*nr-1], ylines_tor[math.floor(np/4)*nr-1], zlines_tor[math.floor(np/4)*nr-1],color='green')
+#ax.plot(xlines_tor[math.floor(np/2)*nr-1], ylines_tor[math.floor(np/2)*nr-1], zlines_tor[math.floor(np/2)*nr-1],color='green')
+#ax.contourf(z,x,Bf,100,zdir='z')
+#ax.set_xlabel('X [cm]')
+#ax.set_ylabel('Y [cm]')
+#ax.set_zlabel('Z [cm]')
+#ax.set_title('Test Case 3')
+#plt.show()
+
+# Get two points, where the drifts (exb or gradb) will be computed and the vector plotted
+p1 = [xlines_tor[math.floor(np/2)*nr-1][math.floor(nt/3)],ylines_tor[math.floor(np/2)*nr-1][math.floor(nt/3)],zlines_tor[math.floor(np/2)*nr-1][math.floor(nt/3)]]
+p2 = [xlines_tor[math.floor(np/4)*nr-1][math.floor(2*nt/3)],ylines_tor[math.floor(np/4)*nr-1][math.floor(2*nt/3)],zlines_tor[math.floor(np/4)*nr-1][math.floor(2*nt/3)]]
+p3= [xlines_tor[math.floor(np/2)*nr-1][math.floor(nt/3)+1],ylines_tor[math.floor(np/2)*nr-1][math.floor(nt/3)+1],zlines_tor[math.floor(np/2)*nr-1][math.floor(nt/3)+1]]
+p4 = [xlines_tor[math.floor(np/4)*nr-1][math.floor(2*nt/3)+1],ylines_tor[math.floor(np/4)*nr-1][math.floor(2*nt/3)+1],zlines_tor[math.floor(np/4)*nr-1][math.floor(2*nt/3)+1]]
+r1 = math.sqrt(p1[0]**2+p1[1]**2)
+theta1 = math.atan2(p1[1],p1[0])
+z1 = p1[2]
+p1cyl = [r1,theta1,z1]
+r2 = math.sqrt(p2[0]**2+p2[1]**2)
+theta2 = math.atan2(p2[1],p2[0])
+z2 = p2[2]
+p2cyl = [r2,theta2,z2]
+
+E1 = e_field_description(r1,theta1,z1,Rmajor,rminor)
+B1 = bvectorcase3(r1,theta1,z1,Rmajor,rminor)
+B2 = bvectorcase3(r2,theta2,z2,Rmajor,rminor)
+gradB2 = gradB(r2,theta2,z2,Rmajor,rminor)
+vexb = ExBdrift(E1,B1)
+vgradb = gradBdrift(B2,gradB2)
+# Plot only the first outer surface and 2 toroidal lines
+fig = plt.figure()
+ax = plt.axes(projection='3d')
+scale = 25
+E1 = cyltocart(E1,p1cyl)
+print('E1', E1)
+B1 = cylvec(B1,p1cyl)
+B1 = nump.array(p3)-nump.array(p1)
+print('B1', B1)
+B2 = cylvec(B2,p2cyl)
+B2 = nump.array(p4)-nump.array(p2)
+print('B2', B2)
+gradB2 = cyltocart(gradB2,p2cyl)
+E1 = 500*scale*nump.array(E1)
+B1 = 2*nump.array(B1)
+B2 = 2*nump.array(B2)
+gradB2 = 1000*scale*nump.array(gradB2)
+vexb = 1000*scale*vexb
+vgradb = 10*scale*vgradb
+
+ax.quiver(p1[0],p1[1],p1[2],E1[0], E1[1],E1[2],color = 'yellow',label='Electric Field Vector')
+ax.quiver(p1[0],p1[1],p1[2],B1[0], B1[1],B1[2],color = 'blue',label = 'Magnetic Field Vector')
+ax.quiver(p1[0],p1[1],p1[2],vexb[0],vexb[1],vexb[2],color = 'red',label = 'E x B Drift Velocity')
+ax.quiver(p2[0],p2[1],p2[2],gradB2[0], gradB2[1],gradB2[2],color = 'black',label = 'Magnetic Field Strength Gradient')
+ax.quiver(p2[0],p2[1],p2[2],B2[0], B2[1],B2[2],color = 'blue')
+ax.quiver(p2[0],p2[1],p2[2],vgradb[0],vgradb[1],vgradb[2],color = 'purple',label = 'Grad B Velocity')
+ax.plot(xlines_pol[nr-1], ylines_pol[nr-1], zlines_pol[nr-1],color='green')
+ax.plot(xlines_pol[(nt)*nr-1], ylines_pol[(nt)*nr-1], zlines_pol[(nt)*nr-1],color='green')
+ax.plot(xlines_tor[math.floor(np/4)*nr-1], ylines_tor[math.floor(np/4)*nr-1], zlines_tor[math.floor(np/4)*nr-1],color='green')
+ax.plot(xlines_tor[math.floor(np/2)*nr-1], ylines_tor[math.floor(np/2)*nr-1], zlines_tor[math.floor(np/2)*nr-1],color='green')
+#ax.contourf(z,x,Bf,100,zdir='z')
+ax.set_xlabel('X [cm]')
+ax.set_ylabel('Y [cm]')
+ax.set_zlabel('Z [cm]')
+ax.set_title('Test Case 3')
+ax.legend()
+plt.show()
\ No newline at end of file

Scripts_Matthieu/PlotMesh.py

+# Plot the vertices generated for the creation of a mesh in EMC3
+from mpl_toolkits import mplot3d
+import numpy as nump
+import matplotlib
+import matplotlib.pyplot as plt
+from MeshFunctions import *
+import math
+from PlottingFunctions import *
+
+# General Info about the case to plot
+nt =1; nr = 10; np = 100
+Rmajor = 500
+rminor = 100
+Trange = 36
+# Create the vertices
+[tor,rad,hor] = create_toroidal_vertices(nt,nr,np,Trange,Rmajor,rminor)
+# PLot the grid points
+plt.figure(1)
+[xdata,ydata,zdata] = plotGridPoints(tor,rad,hor,show = True)
+
+# Plot lines connecting the vertices in all different directions, this will be used to plot the mesh
+plt.figure(2)
+[xlines_pol,xlines_rad,xlines_tor,ylines_pol,ylines_rad,ylines_tor,zlines_pol,zlines_rad,zlines_tor] = plotMesh(tor,rad,hor,np,nr,nt,show=True)
+#plt.show()
+# Make a 2D plot
+for i in range(len(xlines_tor)):
+    plt.plot(xlines_tor[i],zlines_tor[i],color='green')
+for i in range(len(xlines_pol)):
+    plt.plot(xlines_pol[i],zlines_pol[i],color='green')
+for i in range(len(xlines_rad)):
+    plt.plot(xlines_rad[i],zlines_rad[i],color='green')
+plt.xlabel("X [cm]")
+plt.ylabel("Z [cm]")
+plt.title(" Toroidal section of the mesh")
+plt.show()
+# Plot the magnetic and potential values
+
+def b_field_description(r,z,theta, Rmajor, Rminor):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
+    bval = math.exp(-R/Rmajor) # R in cm
+    return bval
+Bfield = [0]*nr*nt*np
+rdata = [0]*nr*nt*np
+for i in range(nr*nt*np):
+    Bfield[i] = b_field_description(math.sqrt(xdata[i]**2+ydata[i]**2),zdata[i], math.atan2(ydata[i],xdata[i]),Rmajor,rminor)
+    rdata[i] = math.sqrt(xdata[i]**2+ydata[i]**2)
+
+#plt.figure(2)
+#ax = plt.axes(projection='3d')
+#ax.scatter3D(xdata, ydata, zdata, Bfield, cmap='plasma')
+#matplotlib.colorbar.ColorbarBase(ax=ax, values=sorted(Bfield),
+ #                                orientation="horizontal")
+#plt.show()
+print(Bfield)
+plt.figure(3)
+# Plot 2d contour of Bfield
+ns = 100
+plot2DContour(ns,ns,Rmajor-rminor,Rmajor+rminor,-rminor,rminor,Rmajor,rminor,b_field_description,'X [cm]', 'Y [cm]','Magnetic field strength',True)
\ No newline at end of file

Scripts_Matthieu/PlotMeshFollowingFieldLines.py

+import numpy as nump
+#import scipy as sp
+import matplotlib.pyplot as plt
+from MeshFunctions import *
+from TracingFieldLines import *
+from PlottingFunctions import *
+from plotSingleFluxTube import *
+import math
+# Create the 2D mesh to start with
+nt =6; nr = 2; np = 30
+Rmajor = 500
+rminor = 100
+Trange = 36
+m2d = create_2D_mesh(nr,np,Rmajor,rminor,'toroidal')
+# Create the analytical description of the B field
+def B(r,z,theta,Rmajor,rminor):
+    R = math.sqrt((r - Rmajor) ** 2 + z ** 2) + 0.01
+    R = 0.01*R
+    # Find poloidal direction => helical field
+    alpha = math.atan2(z,r-Rmajor)
+    Bphimag = 3
+    Br = Bphimag*-math.sin(alpha)
+    Bz = Bphimag*math.cos(alpha)
+    Btheta = 1 #math.exp(-R)
+    return [Br,Bz,Btheta]
+# Create the 3D mesh using the field lines
+precision = 100
+m3d = extend_mesh_3D(m2d,Trange,B,nt,precision, Rmajor, rminor)
+
+rad = m3d[0]
+hor = m3d[1]
+tor = m3d[2]
+field_lines = m3d[3]
+
+
+
+
+# Get the field lines in x,y,z
+lines = []
+for i in range(len(field_lines)):
+    field = []
+    x = []
+    y = []
+    z = []
+    for j in range(len(field_lines[i])):
+        x.append(field_lines[i][j][0] * math.cos(field_lines[i][j][2] * math.pi / 180))
+        y.append(field_lines[i][j][0] * math.sin(field_lines[i][j][2] * math.pi / 180))
+        z.append(field_lines[i][j][1])
+    field = [x,y,z]
+    lines.append(field)
+
+
+# Plot lines connecting the vertices in all different directions, this will be used to plot the mesh
+plt.figure(2)
+plotMesh(tor,rad,hor,nt,nr,np,lines,show=True)
+
+# Plot the field lines
+plt.figure()
+ax = plt.axes(projection='3d')
+for i in range(len(field_lines)):
+    ax.plot(lines[i][0],lines[i][1],lines[i][2])
+plt.show()
+
+
+
+# plot a single flux tube
+plotfluxTube(tor,rad,hor,nt,nr,np,None,True)

Scripts_Matthieu/PlottingFunctions.py

+# Contains functions for plotting
+import math
+import matplotlib
+import matplotlib.pyplot as plt
+import numpy as nump
+# Plot the grid points
+def plotGridPoints(tor,rad,hor,show = True):
+    xdata = []
+    ydata = []
+    zdata = []
+    cdata = []
+    for i in range(len(tor)):
+        for j in range(len(rad[i])):
+            xdata.append(rad[i, j] * math.cos(tor[i] * math.pi / 180))
+            ydata.append(rad[i, j] * math.sin(tor[i] * math.pi / 180))
+            zdata.append(hor[i, j])
+            cdata.append(1)
+
+    # fig = plt.figure()
+    ax = plt.axes(projection='3d')
+    ax.scatter3D(xdata, ydata, zdata, cdata, cmap='plasma')
+    if show:
+        plt.show()
+    return [xdata,ydata,zdata]
+
+def plotMesh(tor,rad,hor,nt,nr,np,field_lines,show=True):
+    plt.figure()
+    ax = plt.axes(projection='3d')
+    # Lines following toroidal direction
+
+    xlines_tor = [None] * np * nr
+    ylines_tor = [None] * np * nr
+    zlines_tor = [None] * np * nr
+    for i in range(len(rad[0])):
+        xtemp = []
+        ytemp = []
+        ztemp = []
+        for k in range(0,nt):
+            xtemp.append(rad[k, i] * math.cos(tor[k] * math.pi / 180))
+            ytemp.append(rad[k, i] * math.sin(tor[k] * math.pi / 180))
+            ztemp.append(hor[k, i])
+        xlines_tor[i] = xtemp
+        ylines_tor[i] = ytemp
+        zlines_tor[i] = ztemp
+
+    if field_lines == None:
+        for i in range(len(rad[0])):
+            ax.plot(xlines_tor[i], ylines_tor[i], zlines_tor[i], color='green')
+    else:
+        for i in range(len(field_lines)):
+            ax.plot(field_lines[i][0], field_lines[i][1], field_lines[i][2],color='green')
+
+    # Poloidal lines
+
+    xlines_pol = [None] * nt * nr
+    ylines_pol = [None] * nt * nr
+    zlines_pol = [None] * nt * nr
+    for i in range(0,nt):
+        for j in range(nr):
+            xtemp = []
+            ytemp = []
+            ztemp = []
+            for k in range(np):
+                xtemp.append(rad[i, k * nr + j] * math.cos(tor[i] * math.pi / 180))
+                ytemp.append(rad[i, k * nr + j] * math.sin(tor[i] * math.pi / 180))
+                ztemp.append(hor[i, k * nr + j])
+            xtemp.append(xtemp[0])
+            ytemp.append(ytemp[0])
+            ztemp.append(ztemp[0])
+            xlines_pol[i * nr + j] = xtemp
+            ylines_pol[i *nr + j] = ytemp
+            zlines_pol[i* nr + j] = ztemp
+            ax.plot(xlines_pol[i * nr + j], ylines_pol[i  * nr + j], zlines_pol[i  * nr + j],
+                    color='green')
+
+    # Radial lines
+
+    xlines_rad = [None] * np * nt
+    ylines_rad = [None] * np * nt
+    zlines_rad = [None] * np * nt
+    for i in range(0,nt):
+        for j in range(np):
+            xtemp = []
+            ytemp = []
+            ztemp = []
+            for k in range(nr):
+                xtemp.append(rad[i, j * nr + k] * math.cos(tor[i] * math.pi / 180))
+                ytemp.append(rad[i, j * nr + k] * math.sin(tor[i] * math.pi / 180))
+                ztemp.append(hor[i, j * nr + k])
+            xtemp.append(xtemp[0])
+            ytemp.append(ytemp[0])
+            ztemp.append(ztemp[0])
+            xlines_rad[i  * nr + j] = xtemp
+            ylines_rad[i  * nr + j] = ytemp
+            zlines_rad[i  * nr + j] = ztemp
+            ax.plot(xlines_rad[i * nr + j], ylines_rad[i * nr + j], zlines_rad[i * nr + j],
+                    color='green')
+
+    ax.set_xlabel('X [cm]')
+    ax.set_ylabel('Y [cm]')
+    ax.set_zlabel('Z [cm]')
+    ax.set_title('Mesh for Test Case EMC3')
+    if show:
+        plt.show()
+
+    return [xlines_pol,xlines_rad,xlines_tor,ylines_pol,ylines_rad,ylines_tor,zlines_pol,zlines_rad,zlines_tor]
+
+
+def plot2DContour(xnb,ynb,xmin,xmax,ymin,ymax,Rmajor, rminor, analy, xlabel,ylabel,title,show):
+    Bfield = [0] * xnb
+    xpoints = nump.linspace(xmin, xmax, xnb)
+    zpoints = nump.linspace(ymin, ymax, ynb)
+
+    for i in range(xnb):
+        Bfield[i] = [0] * ynb
+        for j in range(ynb):
+            Bfield[i][j] = analy(xpoints[j], zpoints[i], 0, Rmajor, rminor)
+    plt.contourf(xpoints, zpoints, Bfield, 1000, cmap='plasma')
+    plt.colorbar()
+    plt.xlabel(xlabel)
+    plt.ylabel(ylabel)
+    plt.title(title)
+    if show:
+        plt.show()
+    return [xpoints,zpoints,Bfield]
\ No newline at end of file

Scripts_Matthieu/PrepareGeometry.py

+# This file defines all the quantities and creates the files required to run a test case for the computation of
+# drift velocities in the EMC3 simulation code
+from MeshFunctions import *
+import os
+import shutil
+# The folder where the different cases should be stored is
+folder = "D:\Matthieu\Documents\\1.Master\Thesis\Data"
+# A specific folder for the case should be defined
+casefolder = "TestCase"
+# create the folder
+path = folder+"\\"+casefolder
+# Check if the folder should be removed => Control to not loose simulations
+if os.path.isdir(path):
+    Flag = int(input("Do you want to remove"+path+" Answer should be 1/0"))
+    if Flag:
+        shutil.rmtree(path)
+    else:
+        path = folder+"\\"+input("Give a new folder name")
+os.mkdir(path)
+path = path+"\\Inputs"
+os.mkdir(path)
+# Choose the parameters for the case to test
+# Geometry
+# Number of zones considered
+nz = 1
+# Number of toroidal, radial, poloidal surfaces
+nt = 20
+nr = 10*10
+np = 10*5
+# toroidal range (degrees)
+trange = 36
+# Major and minor radius (considering a torus) (in cm!)
+Rmajor = 500
+Rminor = 100
+# Create the file with the general info
+gi_path = path+"\\input.geo"
+CreateGeneralInfo(gi_path, nz,[nr],[np],[nt])
+
+
+# Create the GEOMETRY_3D_DATA file
+[tor,rad,hor] = create_toroidal_vertices(nt,nr,np,trange,Rmajor,Rminor)
+print(tor)
+writeGEOMETRY3DDATA(tor,rad,hor,nt,nr,np,path)
+
+# Create the magnetic field data
+# Create an analytical description
+def b_field_description(r,z,theta, Rmajor, Rminor):
+    # Radial B-field (linear and max at central surfaces of the torus
+    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
+    bval = 0.01/R*100  # R in cm
+    return bval
+bfieldpath = path+"\\BFIELD_STRENGTH"
+imposeAnalytic(tor,rad,hor,b_field_description,Rmajor, Rminor, bfieldpath)
+
+# Create the potential data
+# At the moment we use the same grid for magnetic and plasma cells
+def potential_description(r,z,theta,Rmajor,Rminor):
+    R = (r-Rmajor)
+    potval = math.exp(-R/100)*100
+    return potval
+potentialpath = path+"\\POTENTIAL"
+imposeAnalytic(tor,rad,hor,potential_description,Rmajor, Rminor,potentialpath)

Scripts_Matthieu/PrepareGeometryDescription1.py

+from MeshFunctions import *
+import os
+import shutil
+import AnalyticalPrescriptions
+# The folder where the different cases should be stored is
+basepath = "D:\Matthieu\Documents\\1.Master\Thesis\Data\\"
+folder = "D:\Matthieu\Documents\\1.Master\Thesis\Data\Description1"
+if os.path.isdir(folder):
+    Flag = int(input("Do you want to remove" + folder + " Answer should be 1/0"))
+    if Flag:
+        shutil.rmtree(folder)
+    else:
+        folder = basepath + "\\" + input("Give a new folder name")
+os.mkdir(folder)
+# Choose the parameters for the case to test
+# Geometry
+# Number of zones considered
+nz = 1
+# Loop over different cases to prepare
+nt = 10
+nrs = [10,50,100,200,400,1000]
+nps = nrs
+cases = [nrs,nps]
+caseNames = ['Case10', 'Case50', 'Case100','Case200', 'Case400', 'Case1000']
+# Number of toroidal, radial, poloidal surfaces
+
+for i in range(len(cases[0])):
+    # A specific folder for the case should be defined
+    path = folder + "\\" + caseNames[i]
+    if os.path.isdir(path):
+        Flag = int(input("Do you want to remove" + path + " Answer should be 1/0"))
+        if Flag:
+            shutil.rmtree(path)
+        else:
+            path = folder + "\\" + input("Give a new folder name")
+    os.mkdir(path)
+    path = path + "\\Inputs"
+    os.mkdir(path)
+    nr = cases[0][i]
+    np = cases[1][i]
+    potentialdescription = AnalyticalPrescriptions.potential_description1
+    bfielddescription = AnalyticalPrescriptions.bfieldstrength1
+    # toroidal range (degrees)
+    trange = 36
+    # Major and minor radius (considering a torus) (in cm!)
+    Rmajor = 500
+    Rminor = 100
+    # Create the file with the general info
+    gi_path = path+"\\input.geo"
+    CreateGeneralInfo(gi_path, nz,[nr],[np],[nt])
+
+
+    # Create the GEOMETRY_3D_DATA file
+    [tor,rad,hor] = create_toroidal_vertices(nt,nr,np,trange,Rmajor,Rminor)
+    print(tor)
+    writeGEOMETRY3DDATA(tor,rad,hor,nt,nr,np,path)
+
+    # Create the magnetic field data
+    # Create an analytical description
+    bfieldpath = path+"\\BFIELD_STRENGTH"
+    imposeAnalytic(tor,rad,hor,bfielddescription,Rmajor, Rminor, bfieldpath)
+
+    # Create the potential data
+    # At the moment we use the same grid for magnetic and plasma cells
+
+    potentialpath = path+"\\POTENTIAL"
+    imposeAnalytic(tor,rad,hor,potentialdescription,Rmajor, Rminor,potentialpath)

Scripts_Matthieu/ReadDriftResults.py

+# Read in the results computed in fortran
+# Read in the mesh points

Scripts_Matthieu/Testingfieldlines.py

+import numpy as np
+#import scipy as sp
+import matplotlib.pyplot as plt
+from MeshFunctions import *
+from TracingFieldLines import *
+import math
+
+
+# Test the functions
+# Create the 2D mesh to start with
+m2d = create_2D_mesh(10,10,'toroidal')
+# Create the 3D mesh using the field lines
+m3d = extend_mesh_3D(m2d,36,Bfieldvector,10,100000)
+
+# Plotting
+xdata = []; ydata = []; zdata = []
+rad = m3d[0]
+hor = m3d[1]
+tor = m3d[2]
+field_lines = m3d[3]
+# Plot grid points
+for i in range(len(tor)):
+    for j in range(len(rad[i])):
+            xdata.append(rad[i,j]*math.cos(tor[i]*math.pi/180))
+            ydata.append(rad[i,j]*math.sin(tor[i]*math.pi/180))
+            zdata.append(hor[i,j])
+
+
+ax = plt.axes(projection='3d')
+ax.scatter3D(xdata, ydata, zdata, c=zdata, cmap='Dark2')
+# Plot the field lines
+# Plot the field lines
+# Get the field lines in x,y,z
+lines = []
+for i in range(len(field_lines)):
+    field = []
+    x = []
+    y = []
+    z = []
+    for j in range(len(field_lines[i])):
+        x.append(field_lines[i][j][0] * math.cos(field_lines[i][j][2] * math.pi / 180))
+        y.append(field_lines[i][j][0] * math.sin(field_lines[i][j][2] * math.pi / 180))
+        z.append(field_lines[i][j][1])
+    field = [x,y,z]
+    lines.append(field)
+
+# Plot the field lines
+for i in range(len(field_lines)):
+    ax.plot(lines[i][0],lines[i][1],lines[i][2])
+plt.show()
+
+
+# plot only 2d original mesh
+xdata = []; ydata = []; zdata = []
+rad = m2d[0]
+hor = m2d[1]
+field_lines = m3d[3]
+# Plot grid points
+for j in range(len(rad)):
+        xdata.append(rad[j]*math.cos(0*math.pi/180))
+        ydata.append(rad[j]*math.sin(0*math.pi/180))
+        zdata.append(hor[j])
+
+ax = plt.axes(projection='3d')
+ax.scatter3D(xdata, ydata, zdata, c=zdata, cmap='Dark2')
+
+#plt.show()
\ No newline at end of file

Scripts_Matthieu/TracingFieldLines.py

+
+import numpy as np
+#import scipy as sp
+import matplotlib.pyplot as plt
+from MeshFunctions import *
+import math
+
+
+def create_2D_mesh(nb_rad, nb_pol,Rmajor,rminor,geometry):
+    if geometry == 'toroidal':
+        [toroidal_surfaces, radial_positions, vertical_positions] = create_toroidal_vertices(1,nb_rad, nb_pol, 0, Rmajor, rminor)
+
+        return [radial_positions,vertical_positions]
+    else:
+        print('not implemented yet')
+        return None
+
+
+def numerical_fieldline(B,res,nb_tor,Trange,r0,Rmajor,rminor):
+    points = []
+    points.append(r0)
+    r = []; z = []; theta = []
+    i = 0
+    while True:
+        r = points[i][0]
+        z = points[i] [1]
+        theta = points[i][2]
+        Bvector = B(r,theta,z,Rmajor,rminor)
+        Bnorm = np.linalg.norm(Bvector)
+        if Bnorm <0.001:
+            print('BNORM small')
+        n1 = points[i][0]+Bvector[0]*res/Bnorm
+        n2 = points[i][1] + Bvector[2] * res / Bnorm
+        n3 = points[i][2] + Bvector[1] * res / Bnorm
+        points.append([n1,n2,n3])
+        i = i+1
+        if theta>=Trange:
+            break
+    return points
+
+def extend_mesh_3D(mesh_2D,tor_range,B,nb_tor,precision, Rmajor, rminor):
+    field_lines = []
+    nb_2d = len(mesh_2D[0][0])
+    tor_res = tor_range/nb_tor
+    res = tor_res/precision
+    tor_surfaces = np.linspace(0, tor_range, nb_tor)  # Include first and last surface
+    # Initialize arrays for the radial and horizontal positions
+    # Radial position R = R0 (major radius) + a (minor radius)*cos(phi)
+    radial_positions = np.empty(shape=(nb_tor+1,nb_2d),dtype='float')
+    horizontal_positions = np.empty(shape=(nb_tor+1,nb_2d),dtype='float')
+    for i in range(nb_2d):
+        r0 = [mesh_2D[0][0][i],mesh_2D[1][0][i],0]
+        radial_positions[0,i] = r0[0]
+        horizontal_positions[0,i] = r0[1]
+        points = numerical_fieldline(B,res,nb_tor,tor_range,r0, Rmajor, rminor)
+        next_tor = tor_res
+        p = 1
+
+        for j in range(len(points)):
+            if points[j][2]>= next_tor:
+                radial_positions[p,i] = points[j][0]
+                horizontal_positions[p,i] = points[j][1]
+                p+=1
+                next_tor+=tor_res
+        field_lines.append(points)
+    return [radial_positions, horizontal_positions, tor_surfaces,field_lines]
\ No newline at end of file

Scripts_Matthieu/Utilities.py

+# General small functions
+import math
+def listdif(l1,l2):
+    l = []
+    for i in range(len(l1)):
+        l.append(l1[i]-l2[i])
+    return l
+def listnorm(l):
+    norm = math.sqrt(sum([x**2 for x in l]))
+    return norm
\ No newline at end of file