-
Notifications
You must be signed in to change notification settings - Fork 0
/
plotStages.py
157 lines (147 loc) · 5.25 KB
/
plotStages.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
import numpy
import tables
import scipy
import matplotlib
import matplotlib.pyplot as plt
import collections
import scipy.signal
import csv
import egan_vorpalUtil as egan
import os
from mpl_toolkits.axes_grid1 import host_subplot
import mpl_toolkits.axisartist as AA
import sys
sys.path.append('/scr_verus/wernerg/vrun/relRecon/relReconRepo')
#import calcStats
import vorpalUtil
widths=[15]
#width=1
for runNum in [2,3,4]:
print runNum
width=15
pathName="/scr_verus/wernerg/vrun/relRecon/plasmoids/m1-400-sig30-thd1p875-thbp01-res2-ppc3"+str(runNum)+"/"
resultsDir="./size_1600_"+str(runNum)+'/'
plotDir=resultsDir+"distplots/"
runName="relRecon2p_"
origfilename=pathName+runName+'yeeB_'+str(1)+".h5"
MULTPLOTS=False
HIST=False #2d histogram of wy vs psi
WIDTHCOMP=False #wy on y wx on x
AREA=False #wx*wy vs psi
NORMAL=True
DEBUG=False
STAGES=True
#Read in plasmoid info
if not os.path.exists(plotDir):
os.makedirs(plotDir)
execfile(pathName+"relReconVars.py")
print NX_TOTS[0]
filename=resultsDir+runName+"width_"+str(width)+"_extended.csv"
print filename
with open(filename,'r') as csvfile:
freader=csv.reader(csvfile,dialect='excel')
head1=next(freader)
#head2=next(freader)
columns=next(freader)
widths=[]
deltapsi=[]
dump=[]
xwidth=[]
time=[]
#print columns
#print columns[11]
for row in freader:
widths.append(float(row[11])/2+float(row[12])/2)
deltapsi.append(abs(float(row[4])-float(row[9])))
xwidth.append(abs(float(row[3])-float(row[10])))
dump.append(row[0])
time.append(row[13])
#print row
##NOW FOR TOP
filename=resultsDir+runName+"width_"+str(width)+"T_extended.csv"
print filename
with open(filename,'r') as csvfile:
freader=csv.reader(csvfile,dialect='excel')
head1=next(freader)
columns=next(freader)
#print columns
#print columns[11]
for row in freader:
#print row
widths.append(float(row[11])/2+float(row[12])/2)
deltapsi.append(abs(float(row[4])-float(row[9])))
xwidth.append(abs(float(row[3])-float(row[10])))
dump.append(row[0])
time.append(row[13])
#Get B_0 and system + cell size
execfile(pathName+"relReconVars.py")
print NX_TOTS[0]
##Dividing width in cell size by total number of cells should be enough
##wx/L=widths/400 (not 401 because the extra 1 does not count in system size)
widthNorm=[abs(float(i))/NX_TOTS[0] for i in widths]
fluxNorm=[abs(float(i))/(B_0*LX) for i in deltapsi]
xwidthNorm=[abs(float(i)/NX_TOTS[0]) for i in xwidth]
dumpint= numpy.array(dump,dtype=int)
time=numpy.array(time,dtype=float)
colors=["b","g","r","c","m","y"]
markerstyle=['.','o','v','^','<','>','8','s','p','*','+','x','d','D']
markerstyle=markerstyle*5
colors=colors*len(dump)
#######################
## SET UP PLOT #####
#######################
ax = host_subplot(111, axes_class=AA.Axes)
y2=ax.twinx() #For 2nd axes
ax.grid(true)
ax.set_xlim(10**-5,1)
ax.set_ylim(10**-5,1)
plotfile="plasmoidDist_"+str(width)
ylabel="Normalized y widths (wy/L)" #default val
xlabel="Normalized enclosed flux (deltaPsi/B_0*L)" #default val
ax.loglog([10**-5,1],[10**-5,1],linestyle="-",color='k',linewidth=1) #plot 1:1 line
ax.set_ylabel(ylabel)
ax.set_xlabel(xlabel)
x=numpy.arange(1,10**5,10**2)
x=x*(10**-5)
#SET RIGHT AXIS
if AREA==False:
y2.set_ylabel("Number of Cells")
y2.set_yscale("log")
y2.set_ylim(10**(-5)*NX_TOTS[0],1*NX_TOTS[0])
y2.axis["right"].label.set_color('k')
if NORMAL or HIST: #Plot lines
title="Distribution of Plasmoid fluxes wrt Y-width, PeakWidth="+str(width)
ax.loglog([.000005,1],[NX_TOTS[0]**-1,NX_TOTS[0]**-1],linestyle='--',color='k') ##Plot line of cell siz
ax.loglog([.000005,1],[DELTA/LX,DELTA/LX],linestyle='-',color='k') ##Plot line of reconnection layer width
ax.loglog(B_0*DELTA*numpy.log(numpy.cosh(x*NX_TOTS[0]/DELTA)),x,linestyle='-')
#####################
## Plot the values ##
#####################
######
# Get the stages ##
###################
simName=pathName+"relRecon2p"
fluxFn=calcStats.getUnreconnectedFluxVsTime(simName)
dnByName = "layerDnByLine"
byHistTimes = vorpalUtil.getHistoryTimes(simName, dnByName)[:,0]
returns=calcStats.fit2reconRate(byHistTimes, fluxFn(byHistTimes))
coords=returns[2]
stageTimes=1e6*coords[:,0]
##get times
numStages=len(stageTimes)
print stageTimes
print time
for i in range(0,len(stageTimes)-1):
if i==0:
stageWhere=numpy.where(time<stageTimes[i])
else:
stageWhere=numpy.where(numpy.logical_and(time<stageTimes[i],time>stageTimes[i-1]))
for j in stageWhere[0]:
ax.loglog(fluxNorm[j],widthNorm[j],color=colors[i],marker=markerstyle[i],linestyle="")
ax.plot(10,10, color=colors[i],marker=markerstyle[i],linestyle='',label="Stage "+str(i)) #plot pts for legend later
ax.set_title(title)
plt.legend(loc=4,numpoints=1,fontsize='small',title="Step #")
plt.draw()
plt.savefig(plotfile+"_bystages.png",dpi=500)
plt.show()
plt.close()