From badc90d63f72e9a1ca2dec2149520006ff32acd9 Mon Sep 17 00:00:00 2001
From: Matthieu Jacobs <matthieu.jacobs@student.kuleuven.be>
Date: Wed, 2 Mar 2022 15:49:30 +0100
Subject: [PATCH] add gbdrifts
---
Scripts_Matthieu/AnalysisFunctions.py | 75 +++++++++++-----
Scripts_Matthieu/AnalyzeNumericalResults.py | 97 ++++++++++++++-------
2 files changed, 118 insertions(+), 54 deletions(-)
diff --git a/Scripts_Matthieu/AnalysisFunctions.py b/Scripts_Matthieu/AnalysisFunctions.py
index 748766a..b1a25f0 100644
--- a/Scripts_Matthieu/AnalysisFunctions.py
+++ b/Scripts_Matthieu/AnalysisFunctions.py
@@ -4,35 +4,51 @@ import Constants
import math
import Utilities
-def extractValues(positions, drifts, efield, edescr, bdescr, ExB):
+def extractValues(positions, drifts, efield, gradb, drgradb, edescr, bdescr, gradbdescr, ExB, gradBdrift):
pos = []
enum = []
- dnum = []
+ denum = []
+ gnum = []
+ dgnum = []
enumperp = []
- dnumperp = []
+ denumperp = []
+ gnumperp = []
+ dgnumperp = []
ea = []
ba = []
- da = []
+ dea = []
+ dga = []
+ ga = []
+
for i in range(len(drifts)):
# Convert to float values to allow computations and interpretation
- dn = [float(x) for x in drifts[i].split()]
+ den = [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()]
+ gbn = [float(x) for x in gradb[i].split()]
+ dbn = [float(x) for x in drgradb[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)
+ denum.append(dn)
pos.append(dpn)
enum.append(en)
+ gnum.append(gbn)
+ dgnum.append(dbn)
# 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)
+ gfa = gradbdescr(dpn[0], dpn[2], dpn[1], Constants.Rmajor, Constants.Rminor)
+ dga.append(gradBdrift(bga,gfa).tolist())
+ ga.append(gfa)
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]
+ ba.append(bfa)
+ dea.append(ExB(efa, bfa).tolist())
+ denumperp.append(transformCyltoRadPerp(denum[-1],pos[-1]))
+ enumperp.append(transformCyltoRadPerp(enum[-1],pos[-1]))
+ gnumperp.append(transformCyltoRadPerp(gnum[-1],pos[-1]))
+ dgnumperp.append(transformCyltoRadPerp(dgnum[-1],pos[-1]))
+ return [pos, enum, denum,gnum,dgnum, enumperp,denumperp, gnumperp ,dgnumperp,ea , ba, dea, ga,dga]
def extractVectorInfo(numeric,analytical):
res_mag = []
@@ -46,8 +62,8 @@ def extractVectorInfo(numeric,analytical):
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])
+ er_rel_vec.append([math.sqrt(x**2)/max(1,an_mag[i]) for x in er_vec[i]])
+ er_rel_mag.append(er_mag[i]/max(1,an_mag[i]))
return [res_mag,an_mag,er_mag,er_vec,er_rel_mag,er_rel_vec]
@@ -76,21 +92,26 @@ def getAllSlices(positions,allInfo,direction,coordinate):
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)
+def getAllInfo(positions,drifts,efield, gradb, drgradb,edescr,bdescr, gradbdescr, ExB, gradBdrift):
+ [pos, enum, denum,gnum,dgnum, enumperp,denumperp, gnumperp ,dgnumperp,ea , ba, dea, ga,dga] =
+ extractValues(positions, drifts, efield, gradb, drgradb, edescr, bdescr,gradbdescr, ExB, gradBdrift)
- [dres_mag, dan_mag, der_mag, der_vec, der_rel_mag, der_rel_vec] = extractVectorInfo(dnum,da)
+ [deres_mag, dean_mag, deer_mag, deer_vec, deer_rel_mag, deer_rel_vec] = extractVectorInfo(denum,dea)
[eres_mag, ean_mag, eer_mag, eer_vec, eer_rel_mag, eer_rel_vec] = extractVectorInfo(enum, ea)
+ [gres_mag, gan_mag, ger_mag, ger_vec, ger_rel_mag, ger_rel_vec] = extractVectorInfo(gnum, ga)
+ [dgres_mag, dgan_mag, dger_mag, dger_vec, dger_rel_mag, dger_rel_vec] = extractVectorInfo(dgnum, dga)
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']
+ keys = ['edrifts', 'efield','edrifts_perp','efield_perp' ,'ef_mag', 'ef_an_mag', 'edrift_mag', 'edrift_an_mag', 'ef_vec_er', 'ef_mag_er',
+ 'edrift_vec_er', 'edrift_er_mag','efield_er_rel_mag', 'efield_er_rel_vec', 'edrift_er_rel_mag', 'edrift_er_rel_vec',
+ 'gdrifts', 'gradb','gdrifts_perp','gradb_perp', 'gb_mag','gb_an_mag', 'gdrift_mag','gdrift_an_mag','gb_vec_er', 'gb_mag_er',
+ 'gdrift_vec_er','gdrift_er_mag','gb_er_rel_mag','gb_er_rel_vec','gdrift_er_rel_mag', 'gb_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]
+ values = [denum,enum,denumperp,enumperp,eres_mag,ean_mag,deres_mag,dean_mag,eer_vec,eer_mag,deer_vec,der_mag,eer_rel_mag,eer_rel_vec,deer_rel_mag,deer_rel_vec,
+ dgnum,enum,dgnumperp,enumperp,eres_mag,ean_mag,dgres_mag,dgan_mag,eer_vec,eer_mag,dger_vec,dgr_mag,eer_rel_mag,eer_rel_vec,dger_rel_mag,dger_rel_vec,pos]
for i in range(len(values)):
torSliceDict[keys[i]] = values[i]
@@ -106,10 +127,16 @@ def readText(filepath):
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)
+ result = []
+ phi = math.atan2(pos[2],pos[0]-500)
+ c = math.cos(phi)
+ s = math.sin(phi)
+ result.append(abs(v[0]*c+v[2]*s))
+ result.append(v[1])
+ result.append(abs(v[2]*c-v[0]*s))
+# print('phi',phi,s,c)
+ # print('v',v)
+ # print('res',result)
return result
diff --git a/Scripts_Matthieu/AnalyzeNumericalResults.py b/Scripts_Matthieu/AnalyzeNumericalResults.py
index ffc89d0..1742b0c 100644
--- a/Scripts_Matthieu/AnalyzeNumericalResults.py
+++ b/Scripts_Matthieu/AnalyzeNumericalResults.py
@@ -9,16 +9,24 @@ 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 = []
+path = '/u/mbj/Drifts/data-for-drift-computations/Description2/Case50/OutputGG/'
+edrifts_file = open(path+'ExB_RESULTS','r')
+driftpos_file = open(path+'DRIFT_POS','r')
+efield_file = open(path+'EFIELD','r')
+gradb_file = open(path+'GRADB', 'r')
+gdrifts_file = open(path+'GRADB_DRIFTS', 'r')
+edrifts = []
+edrift_pos = []
efield = []
+gradb = []
+gdrifts = []
# Read numeric drift values
-for line in drifts_file:
- drifts.append(line)
+for line in edrifts_file:
+ edrifts.append(line)
+
+for line in gdrifts_file:
+ gdrifts.append(line)
# Read drift positions
for line in driftpos_file:
@@ -26,32 +34,48 @@ for line in driftpos_file:
# Electric field
for line in efield_file:
- efield.append(line)
+ efield.append(line)
+
+# GradB
+for line in gradb_file:
+ gradb.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)
+[pos, enum, denum,gnum,dgnum, enumperp,denumperp, gnumperp ,dgnumperp,ea , ba, dea, ga,dga] = AnalysisFunctions.extractValues(drift_pos,drifts, efield,AnalyticalPrescriptions.e_descr3,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)
+[edrift_res_mag,edrift_an_mag,edrift_er_mag,edrift_er_vec,edrift_er_rel_mag,edrift_er_rel_vec] = AnalysisFunctions.extractVectorInfo(denum,dea)
+[efield_res_mag,efield_an_mag,efield_er_mag,efield_er_vec,efield_er_rel_mag,efield_er_rel_vec] = AnalysisFunctions.extractVectorInfo(enum,ea)
+[gdrift_res_mag,gdrift_an_mag,gdrift_er_mag,gdrift_er_vec,gdrift_er_rel_mag,gdrift_er_rel_vec] = AnalysisFunctions.extractVectorInfo(dgnum,dga)
+[gradb_res_mag,gradb_an_mag,gradb_er_mag,gradb_er_vec,gradb_er_rel_mag,gradb_er_rel_vec] = AnalysisFunctions.extractVectorInfo(gnum,ga)
+
# Take slices in different directions to allow for extra analysis
+#print(drift_pos_num)
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)
+keys = ['edrifts', 'efield','edrifts_perp','efield_perp' ,'ef_mag', 'ef_an_mag', 'edrift_mag', 'edrift_an_mag', 'ef_vec_er', 'ef_mag_er',
+ 'edrift_vec_er', 'edrift_er_mag','efield_er_rel_mag', 'efield_er_rel_vec', 'edrift_er_rel_mag', 'edrift_er_rel_vec',
+ 'gdrifts', 'gradb','gdrifts_perp','gradb_perp', 'gb_mag','gb_an_mag', 'gdrift_mag','gdrift_an_mag','gb_vec_er', 'gb_mag_er',
+ 'gdrift_vec_er','gdrift_er_mag','gb_er_rel_mag','gb_er_rel_vec','gdrift_er_rel_mag', 'gb_er_rel_vec', 'pos']
+allInfo = [enum, denum,gnum,dgnum, enumperp,denumperp, gnumperp ,dgnumperp,ea , ba, dea, ga,dga,
+ edrift_res_mag,edrift_an_mag,edrift_er_mag,edrift_er_vec,edrift_er_rel_mag,edrift_er_rel_vec,
+ efield_res_mag,efield_an_mag,efield_er_mag,efield_er_vec,efield_er_rel_mag,efield_er_rel_vec,
+ gdrift_res_mag,gdrift_an_mag,gdrift_er_mag,gdrift_er_vec,gdrift_er_rel_mag,gdrift_er_rel_vec,
+ gradb_res_mag,gradb_an_mag,gradb_er_mag,gradb_er_vec,gradb_er_rel_mag,gradb_er_rel_vec]
+
+torslices = AnalysisFunctions.getAllSlices(pos,allInfo,1,torsurf)
for i in range(len(torslices)):
torSliceDict[keys[i]] = torslices[i]
+print(torSliceDict.get('ef_mag'))
+print(torSliceDict.get('ef_an_mag'))
print(torSliceDict.get('efield_er_rel_mag'))
-print(sum(torSliceDict.get('efield_er_rel_mag')))
-print(len(torSliceDict.get('efield_er_rel_mag')))
+print('max',max(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')]
@@ -63,64 +87,77 @@ 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')
+title = 'Magnitude of Electic Field'
+plt.title(title)
sc1 = plt.scatter(x, y, s=20, c=torSliceDict.get('ef_mag'), cmap='plasma')
plt.colorbar(sc1)
-
+plt.savefig(path+title+'.png')
plot2 = plt.figure(2)
plt.xlabel('Radial Position [cm]')
plt.ylabel('Vertical Position [cm]')
-plt.title('Theta component of Electric Field')
+title = 'Theta component of Electric Field'
+plt.title(title)
sc2 = plt.scatter(x, y, s=20, c = [x[1] for x in torSliceDict.get('efield')], cmap='plasma')
plt.colorbar(sc2)
+plt.savefig(path+title+'.png')
plot3 = plt.figure(3)
plt.xlabel('Radial Position [cm]')
plt.ylabel('Vertical Position [cm]')
-plt.title('Relative magnitude error of Electric Field')
+title = 'Relative magnitude error of Electric Field'
+plt.title(title)
sc3 = plt.scatter(x, y, s=20, c=torSliceDict.get('efield_er_rel_mag'), cmap='plasma')
plt.colorbar(sc3)
+plt.savefig(path+title+'.png')
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)
+plt.savefig(path+title+'.png')
plot5 = plt.figure(5)
plt.xlabel('Radial Position [cm]')
plt.ylabel('Vertical Position [cm]')
-plt.title('Magnitude of ExB Drift (Numerical)')
+title = 'Magnitude of ExB Drift (Numerical)'
+plt.title(title)
sc5 = plt.scatter(x, y, s=20, c=torSliceDict.get('drift_mag'), cmap='plasma')
plt.colorbar(sc5)
+plt.savefig(path+title+'.png')
+
plot6 = plt.figure(6)
plt.xlabel('Radial Position [cm]')
plt.ylabel('Vertical Position [cm]')
-plt.title('Magnitude of ExB Drift (Analytical)')
+title = 'Magnitude of ExB Drift (Analytical)'
+plt.title(title)
sc6 = plt.scatter(x, y, s=20, c=torSliceDict.get('drift_an_mag'), cmap='plasma')
plt.colorbar(sc6)
+plt.savefig(path+title+'.png')
+
plot7 = plt.figure(7)
plt.xlabel('Radial Position [cm]')
plt.ylabel('Vertical Position [cm]')
-plt.title('Relative magnitude error of ExB Drift')
+title = 'Relative magnitude error of ExB Drift'
+plt.title(title)
sc7 = plt.scatter(x, y, s=20, c=torSliceDict.get('drift_er_rel_mag'), cmap='plasma')
plt.colorbar(sc7)
-
+plt.savefig(path+title+'.png')
plot8 = plt.figure(8)
plt.xlabel('Radial Position [cm]')
plt.ylabel('Vertical Position [cm]')
-plt.title('Magnitude of Efield (Rad)')
+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)
-
+plt.savefig(path+title+'.png')
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.savefig(path+title+'.png')
-
-plt.show()
+#plt.show()
--
GitLab