From 79436ac283a505c1f2d78d611746635530d28c59 Mon Sep 17 00:00:00 2001 From: Matthieu Jacobs <matthieu.jacobs@student.kuleuven.be> Date: Sun, 27 Feb 2022 12:00:15 +0100 Subject: [PATCH] python git --- Scripts_Matthieu/AnalysisFunctions.py | 118 +++++++++++ Scripts_Matthieu/AnalyticalPrescriptions.py | 51 +++++ Scripts_Matthieu/AnalyzeNumericalResults.py | 126 ++++++++++++ Scripts_Matthieu/AnalyzeResults.py | 102 +++++++++ Scripts_Matthieu/Constants.py | 4 + Scripts_Matthieu/CreateMesh.py | 10 + Scripts_Matthieu/ErrorConvergence.py | 68 ++++++ Scripts_Matthieu/MeshFunctions.py | 193 ++++++++++++++++++ .../NumericalAnalysisDescription1.py | 59 ++++++ Scripts_Matthieu/PlotCase1.py | 112 ++++++++++ Scripts_Matthieu/PlotCase3.py | 153 ++++++++++++++ Scripts_Matthieu/PlotMesh.py | 59 ++++++ .../PlotMeshFollowingFieldLines.py | 67 ++++++ Scripts_Matthieu/PlottingFunctions.py | 125 ++++++++++++ Scripts_Matthieu/PrepareGeometry.py | 62 ++++++ .../PrepareGeometryDescription1.py | 67 ++++++ Scripts_Matthieu/ReadDriftResults.py | 2 + Scripts_Matthieu/Testingfieldlines.py | 67 ++++++ Scripts_Matthieu/TracingFieldLines.py | 66 ++++++ Scripts_Matthieu/Utilities.py | 10 + Scripts_Matthieu/driftFunctions.py | 45 ++++ Scripts_Matthieu/plotSingleFluxTube.py | 96 +++++++++ Scripts_Matthieu/polarcontour.py | 28 +++ 23 files changed, 1690 insertions(+) create mode 100644 Scripts_Matthieu/AnalysisFunctions.py create mode 100644 Scripts_Matthieu/AnalyticalPrescriptions.py create mode 100644 Scripts_Matthieu/AnalyzeNumericalResults.py create mode 100644 Scripts_Matthieu/AnalyzeResults.py create mode 100644 Scripts_Matthieu/Constants.py create mode 100644 Scripts_Matthieu/CreateMesh.py create mode 100644 Scripts_Matthieu/ErrorConvergence.py create mode 100644 Scripts_Matthieu/MeshFunctions.py create mode 100644 Scripts_Matthieu/NumericalAnalysisDescription1.py create mode 100644 Scripts_Matthieu/PlotCase1.py create mode 100644 Scripts_Matthieu/PlotCase3.py create mode 100644 Scripts_Matthieu/PlotMesh.py create mode 100644 Scripts_Matthieu/PlotMeshFollowingFieldLines.py create mode 100644 Scripts_Matthieu/PlottingFunctions.py create mode 100644 Scripts_Matthieu/PrepareGeometry.py create mode 100644 Scripts_Matthieu/PrepareGeometryDescription1.py create mode 100644 Scripts_Matthieu/ReadDriftResults.py create mode 100644 Scripts_Matthieu/Testingfieldlines.py create mode 100644 Scripts_Matthieu/TracingFieldLines.py create mode 100644 Scripts_Matthieu/Utilities.py create mode 100644 Scripts_Matthieu/driftFunctions.py create mode 100644 Scripts_Matthieu/plotSingleFluxTube.py create mode 100644 Scripts_Matthieu/polarcontour.py diff --git a/Scripts_Matthieu/AnalysisFunctions.py b/Scripts_Matthieu/AnalysisFunctions.py new file mode 100644 index 0000000..7e7f002 --- /dev/null +++ b/Scripts_Matthieu/AnalysisFunctions.py @@ -0,0 +1,118 @@ +# 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 + + + + + diff --git a/Scripts_Matthieu/AnalyticalPrescriptions.py b/Scripts_Matthieu/AnalyticalPrescriptions.py new file mode 100644 index 0000000..d9d4fb7 --- /dev/null +++ b/Scripts_Matthieu/AnalyticalPrescriptions.py @@ -0,0 +1,51 @@ +# 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 diff --git a/Scripts_Matthieu/AnalyzeNumericalResults.py b/Scripts_Matthieu/AnalyzeNumericalResults.py new file mode 100644 index 0000000..ffc89d0 --- /dev/null +++ b/Scripts_Matthieu/AnalyzeNumericalResults.py @@ -0,0 +1,126 @@ +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() + diff --git a/Scripts_Matthieu/AnalyzeResults.py b/Scripts_Matthieu/AnalyzeResults.py new file mode 100644 index 0000000..f47c291 --- /dev/null +++ b/Scripts_Matthieu/AnalyzeResults.py @@ -0,0 +1,102 @@ +# 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 diff --git a/Scripts_Matthieu/Constants.py b/Scripts_Matthieu/Constants.py new file mode 100644 index 0000000..02dbe68 --- /dev/null +++ b/Scripts_Matthieu/Constants.py @@ -0,0 +1,4 @@ +# This file contains constant values, to be used in other files +Rmajor = 500 +Rminor = 100 +ns = 100 \ No newline at end of file diff --git a/Scripts_Matthieu/CreateMesh.py b/Scripts_Matthieu/CreateMesh.py new file mode 100644 index 0000000..07a86b6 --- /dev/null +++ b/Scripts_Matthieu/CreateMesh.py @@ -0,0 +1,10 @@ +# 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") + diff --git a/Scripts_Matthieu/ErrorConvergence.py b/Scripts_Matthieu/ErrorConvergence.py new file mode 100644 index 0000000..928e4f6 --- /dev/null +++ b/Scripts_Matthieu/ErrorConvergence.py @@ -0,0 +1,68 @@ +# 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() diff --git a/Scripts_Matthieu/MeshFunctions.py b/Scripts_Matthieu/MeshFunctions.py new file mode 100644 index 0000000..4f0ced4 --- /dev/null +++ b/Scripts_Matthieu/MeshFunctions.py @@ -0,0 +1,193 @@ +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 + + + + + + + diff --git a/Scripts_Matthieu/NumericalAnalysisDescription1.py b/Scripts_Matthieu/NumericalAnalysisDescription1.py new file mode 100644 index 0000000..0f54eaf --- /dev/null +++ b/Scripts_Matthieu/NumericalAnalysisDescription1.py @@ -0,0 +1,59 @@ +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() diff --git a/Scripts_Matthieu/PlotCase1.py b/Scripts_Matthieu/PlotCase1.py new file mode 100644 index 0000000..1535885 --- /dev/null +++ b/Scripts_Matthieu/PlotCase1.py @@ -0,0 +1,112 @@ +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 diff --git a/Scripts_Matthieu/PlotCase3.py b/Scripts_Matthieu/PlotCase3.py new file mode 100644 index 0000000..fc78245 --- /dev/null +++ b/Scripts_Matthieu/PlotCase3.py @@ -0,0 +1,153 @@ +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 diff --git a/Scripts_Matthieu/PlotMesh.py b/Scripts_Matthieu/PlotMesh.py new file mode 100644 index 0000000..8c68580 --- /dev/null +++ b/Scripts_Matthieu/PlotMesh.py @@ -0,0 +1,59 @@ +# 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 diff --git a/Scripts_Matthieu/PlotMeshFollowingFieldLines.py b/Scripts_Matthieu/PlotMeshFollowingFieldLines.py new file mode 100644 index 0000000..1a0008c --- /dev/null +++ b/Scripts_Matthieu/PlotMeshFollowingFieldLines.py @@ -0,0 +1,67 @@ +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) diff --git a/Scripts_Matthieu/PlottingFunctions.py b/Scripts_Matthieu/PlottingFunctions.py new file mode 100644 index 0000000..2c4d1bf --- /dev/null +++ b/Scripts_Matthieu/PlottingFunctions.py @@ -0,0 +1,125 @@ +# 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 diff --git a/Scripts_Matthieu/PrepareGeometry.py b/Scripts_Matthieu/PrepareGeometry.py new file mode 100644 index 0000000..392f68f --- /dev/null +++ b/Scripts_Matthieu/PrepareGeometry.py @@ -0,0 +1,62 @@ +# 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) diff --git a/Scripts_Matthieu/PrepareGeometryDescription1.py b/Scripts_Matthieu/PrepareGeometryDescription1.py new file mode 100644 index 0000000..6834292 --- /dev/null +++ b/Scripts_Matthieu/PrepareGeometryDescription1.py @@ -0,0 +1,67 @@ +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) diff --git a/Scripts_Matthieu/ReadDriftResults.py b/Scripts_Matthieu/ReadDriftResults.py new file mode 100644 index 0000000..94fee41 --- /dev/null +++ b/Scripts_Matthieu/ReadDriftResults.py @@ -0,0 +1,2 @@ +# Read in the results computed in fortran +# Read in the mesh points diff --git a/Scripts_Matthieu/Testingfieldlines.py b/Scripts_Matthieu/Testingfieldlines.py new file mode 100644 index 0000000..b98358f --- /dev/null +++ b/Scripts_Matthieu/Testingfieldlines.py @@ -0,0 +1,67 @@ +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 diff --git a/Scripts_Matthieu/TracingFieldLines.py b/Scripts_Matthieu/TracingFieldLines.py new file mode 100644 index 0000000..e8277fc --- /dev/null +++ b/Scripts_Matthieu/TracingFieldLines.py @@ -0,0 +1,66 @@ + +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 diff --git a/Scripts_Matthieu/Utilities.py b/Scripts_Matthieu/Utilities.py new file mode 100644 index 0000000..07f41f1 --- /dev/null +++ b/Scripts_Matthieu/Utilities.py @@ -0,0 +1,10 @@ +# 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 diff --git a/Scripts_Matthieu/driftFunctions.py b/Scripts_Matthieu/driftFunctions.py new file mode 100644 index 0000000..9391799 --- /dev/null +++ b/Scripts_Matthieu/driftFunctions.py @@ -0,0 +1,45 @@ +import numpy as np +import math +def ExBdrift(E,B): + E = np.array(E) + # print('E', E) + B = np.array(B) +# print('B',B) + normB = np.linalg.norm(B) + return np.cross(E,B)*1/normB**2 + +def gradBdrift(B,gradB): + Tq = 150 # keV => no e for electroncharge (should be 000) + return 2*Tq * np.cross(B,gradB)/np.linalg.norm(B)**3 + +def cyltocart(V,p): + alpha = p[1] + C = np.array([0.0,0.0,0.0]) + C[0] = -math.sin(alpha)*V[1] + math.cos(alpha)*V[0] + C[1] = math.cos(alpha)*V[1] + math.sin(alpha)*V[0] + C[2] = V[2] + return C + +def cylvec(V,p): + C = np.array([0.0,0.0,0.0]) + dif = np.array([0.0,0.0,0.0]) + pcart = np.array([0.0,0.0,0.0]) + V = np.array(V) + p = np.array(p)/100 + V[0] += p[0] + V[1] += p[1] + V[2] += p[2] + pcart[0] = p[0]*math.cos(p[1]) + pcart[1] = p[1]*math.sin(p[1]) + pcart[2] = p[2] + C[0]= V[0]*math.cos(V[1]) + C[1] = V[0]*math.sin(V[1]) + C[2] = V[2] + dif = C-pcart + return dif +def carttocyl(V): + C = np.array([0, 0, 0]) + C[0] = math.sqrt(V[0]**2+V[1]**2) + C[1] = math.atan2(V[1],V[0]) + C[2] = V[2] + return C \ No newline at end of file diff --git a/Scripts_Matthieu/plotSingleFluxTube.py b/Scripts_Matthieu/plotSingleFluxTube.py new file mode 100644 index 0000000..a185e88 --- /dev/null +++ b/Scripts_Matthieu/plotSingleFluxTube.py @@ -0,0 +1,96 @@ + +import math +import matplotlib +import matplotlib.pyplot as plt +import numpy as nump + +def plotfluxTube(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(0,2): + xtemp = [] + ytemp = [] + ztemp = [] + xtemp2 = [] + ytemp2= [] + ztemp2 = [] + for k in range(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]) + xtemp2.append(rad[k, i+nr] * math.cos(tor[k] * math.pi / 180)) + ytemp2.append(rad[k, i+nr] * math.sin(tor[k] * math.pi / 180)) + ztemp2.append(hor[k, i+nr]) + xlines_tor[i] = xtemp + ylines_tor[i] = ytemp + zlines_tor[i] = ztemp + xlines_tor[i+2] = xtemp2 + ylines_tor[i+2] = ytemp2 + zlines_tor[i+2] = ztemp2 + + if field_lines == None: + for i in range(4): + 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(nt): + for j in range(2): + xtemp = [] + ytemp = [] + ztemp = [] + for k in range(2): + 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(nt): + for j in range(2): + xtemp = [] + ytemp = [] + ztemp = [] + for k in range(2): + 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('Single flux tube for helical magnetic field') + if show: + plt.show() + + return [xlines_pol,xlines_rad,xlines_tor,ylines_pol,ylines_rad,ylines_tor,zlines_pol,zlines_rad,zlines_tor] diff --git a/Scripts_Matthieu/polarcontour.py b/Scripts_Matthieu/polarcontour.py new file mode 100644 index 0000000..3ca549f --- /dev/null +++ b/Scripts_Matthieu/polarcontour.py @@ -0,0 +1,28 @@ +import numpy as nump +import matplotlib.pyplot as plt +import math +azimuths = nump.radians(nump.linspace(0, 360, 20)) +zeniths = nump.arange(0, 100, 10) + + +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 + +r, theta = nump.meshgrid(zeniths, azimuths) +Bval = nump.empty((len(r),len(r[0]))) + +for i in range(len(r)): + for j in range(len(r[0])): + Bval[i][j] =b_field_description(r[i][j],0,theta[i][j],100,100) +values = nump.random.random((azimuths.size, zeniths.size)) +print(r) +print(theta) +#-- Plot... ------------------------------------------------ +fig, ax = plt.subplots(subplot_kw=dict(projection='polar')) +ax.contourf(theta, r, Bval) +ax.colorbar() +plt.show() -- GitLab