mirror of
https://github.com/JHUAPL/kaiju.git
synced 2026-01-09 15:17:56 -05:00
Adding new plotting tools to measure difference between two runs
This commit is contained in:
Jeffrey Garretson
@@ -23,11 +23,12 @@ eMax = 5.0 #Max current for contours
|
||||
|
||||
#Default pressure colorbar
|
||||
vP = kv.genNorm(vMin=1.0e-2,vMax=10.0,doLog=True)
|
||||
szStrs = ['small','std','big','dm']
|
||||
szStrs = ['small','std','big','bigger','fullD','dm']
|
||||
szBds = {}
|
||||
szBds["std"] = [-40.0 ,20.0,2.0]
|
||||
szBds["big"] = [-100.0,20.0,2.0]
|
||||
szBds["bigger"] = [-200.0,25.0,2.0]
|
||||
szBds["fullD"] = [-300.0,30.0,3.0] # full domain for double res
|
||||
szBds["small"] = [-10.0 , 5.0,2.0]
|
||||
szBds["dm"] = [-30.0 ,10.0,40.0/15.0]
|
||||
|
||||
@@ -48,6 +49,209 @@ def GetSizeBds(args):
|
||||
|
||||
return xyBds
|
||||
|
||||
#Plot absolute error in the requested, or given, equatorial field
|
||||
def PlotEqErrAbs(gsphP,gsphO,nStp,xyBds,Ax,fieldNames,AxCB=None,doClear=True,doDeco=True,vMin=1e-9,vMax=1e-4,doLog=True,doVerb=True):
|
||||
#specify two gsph objects as gsphPredicted (gsphP) and gsphObserved (gsphO)
|
||||
#
|
||||
normAbs = kv.genNorm(vMin,vMax,doLog=doLog)
|
||||
cmapAbs = "PRGn"
|
||||
if doLog:
|
||||
cmapAbs = "magma"
|
||||
|
||||
if (AxCB is not None):
|
||||
#Add the colorbar to AxCB
|
||||
AxCB.clear()
|
||||
|
||||
cbString = ""
|
||||
for fn in fieldNames:
|
||||
cbString = cbString + "'" + fn + "', "
|
||||
cbString = cbString + " Absolute Error"
|
||||
kv.genCB(AxCB,normAbs,cbString,cM=cmapAbs)
|
||||
#Now do main plotting
|
||||
if (doClear):
|
||||
Ax.clear()
|
||||
dataAbs = None
|
||||
for fn in fieldNames:
|
||||
dataP = gsphP.EggSlice(fn,nStp,doEq=True,doVerb=doVerb)
|
||||
dataO = gsphO.EggSlice(fn,nStp,doEq=True,doVerb=doVerb)
|
||||
if dataAbs is None:
|
||||
dataAbs = np.square(dataO - dataP)
|
||||
else:
|
||||
dataAbs = dataAbs + np.square(dataO - dataP)
|
||||
dataAbs = np.sqrt(dataAbs)
|
||||
Ax.pcolormesh(gsphP.xxi,gsphP.yyi,dataAbs,cmap=cmapAbs,norm=normAbs)
|
||||
|
||||
kv.SetAx(xyBds,Ax)
|
||||
|
||||
if (doDeco):
|
||||
kv.addEarth2D(ax=Ax)
|
||||
Ax.set_xlabel('SM-X [Re]')
|
||||
Ax.set_ylabel('SM-Y [Re]')
|
||||
return dataAbs
|
||||
|
||||
#Plot relative error in the requested, or given, equatorial field
|
||||
def PlotEqErrRel(gsphP,gsphO,nStp,xyBds,Ax,fieldNames,AxCB=None,doClear=True,doDeco=True,vMin=1e-16,vMax=1,doLog=True,doVerb=True):
|
||||
#specify two gsph objects as gsphPredicted (gsphP) and gsphObserved (gsphO)
|
||||
#
|
||||
normRel = kv.genNorm(vMin,vMax,doLog=doLog)
|
||||
cmapRel = "PRGn"
|
||||
if doLog:
|
||||
cmapRel = "magma"
|
||||
|
||||
if (AxCB is not None):
|
||||
#Add the colorbar to AxCB
|
||||
AxCB.clear()
|
||||
cbString = ""
|
||||
for fn in fieldNames:
|
||||
cbString = cbString + "'" + fn + "', "
|
||||
cbString = cbString + " Relative Error"
|
||||
kv.genCB(AxCB,normRel,cbString,cM=cmapRel)
|
||||
#Now do main plotting
|
||||
if (doClear):
|
||||
Ax.clear()
|
||||
dataAbs = None
|
||||
dataBase = None
|
||||
for fn in fieldNames:
|
||||
dataP = gsphP.EggSlice(fn,nStp,doEq=True,doVerb=doVerb)
|
||||
dataO = gsphO.EggSlice(fn,nStp,doEq=True,doVerb=doVerb)
|
||||
if dataAbs is None:
|
||||
dataBase = np.square(dataP)
|
||||
dataAbs = np.square(dataO - dataP)
|
||||
else:
|
||||
dataBase = dataBase + np.square(dataP)
|
||||
dataAbs = dataAbs + np.square(dataO - dataP)
|
||||
dataBase = np.sqrt(dataBase)
|
||||
dataAbs = np.sqrt(dataAbs)
|
||||
# math breaks when Base field is exactly 0
|
||||
dataBase[dataBase == 0] = np.finfo(np.float32).tiny
|
||||
dataRel = np.absolute(dataAbs/dataBase)
|
||||
Ax.pcolormesh(gsphP.xxi,gsphP.yyi,dataRel,cmap=cmapRel,norm=normRel)
|
||||
|
||||
kv.SetAx(xyBds,Ax)
|
||||
|
||||
if (doDeco):
|
||||
kv.addEarth2D(ax=Ax)
|
||||
Ax.set_xlabel('SM-X [Re]')
|
||||
Ax.set_ylabel('SM-Y [Re]')
|
||||
return dataRel
|
||||
|
||||
#Plot absolute error along the requested logical axis
|
||||
def PlotLogicalErrAbs(gsphP,gsphO,nStp,Ax,fieldNames,meanAxis,AxCB=None,doClear=True,doDeco=True,vMin=1e-16,vMax=1,doLog=True,doVerb=True):
|
||||
#specify two gsph objects as gsphPredicted (gsphP) and gsphObserved (gsphO)
|
||||
#
|
||||
normAbs = kv.genNorm(vMin,vMax,doLog=doLog)
|
||||
cmapAbs = "PRGn"
|
||||
if doLog:
|
||||
cmapAbs = "magma"
|
||||
|
||||
if (AxCB is not None):
|
||||
#Add the colorbar to AxCB
|
||||
AxCB.clear()
|
||||
cbString = ""
|
||||
for fn in fieldNames:
|
||||
cbString = cbString + "'" + fn + "', "
|
||||
cbString = cbString + " Absolute Error along "
|
||||
if meanAxis == 0:
|
||||
cbString = cbString + "I axis"
|
||||
elif meanAxis == 1:
|
||||
cbString = cbString + "J axis"
|
||||
elif meanAxis == 2:
|
||||
cbString = cbString + "K axis"
|
||||
kv.genCB(AxCB,normAbs,cbString,cM=cmapAbs)
|
||||
#Now do main plotting
|
||||
if (doClear):
|
||||
Ax.clear()
|
||||
dataAbs = CalcTotalErrAbs(gsphP,gsphO,nStp,fieldNames,doVerb=doVerb,meanAxis=meanAxis)
|
||||
dataAbs = np.transpose(dataAbs) # transpose to put I/J on the horizontal axis
|
||||
Ax.pcolormesh(dataAbs,cmap=cmapAbs,norm=normAbs)
|
||||
|
||||
if (doDeco):
|
||||
if meanAxis == 0:
|
||||
Ax.set_xlabel('J Indices')
|
||||
Ax.set_ylabel('K Indices')
|
||||
elif meanAxis == 1:
|
||||
Ax.set_xlabel('I Indices')
|
||||
Ax.set_ylabel('K Indices')
|
||||
elif meanAxis == 2:
|
||||
Ax.set_xlabel('I Indices')
|
||||
Ax.set_ylabel('J Indices')
|
||||
return dataAbs
|
||||
|
||||
#Plot relative error along the requested logical axis
|
||||
def PlotLogicalErrRel(gsphP,gsphO,nStp,Ax,fieldNames,meanAxis,AxCB=None,doClear=True,doDeco=True,vMin=1e-16,vMax=1,doLog=True,doVerb=True):
|
||||
#specify two gsph objects as gsphPredicted (gsphP) and gsphObserved (gsphO)
|
||||
#
|
||||
normRel = kv.genNorm(vMin,vMax,doLog=doLog)
|
||||
cmapRel = "PRGn"
|
||||
if doLog:
|
||||
cmapRel = "magma"
|
||||
|
||||
if (AxCB is not None):
|
||||
#Add the colorbar to AxCB
|
||||
AxCB.clear()
|
||||
cbString = ""
|
||||
for fn in fieldNames:
|
||||
cbString = cbString + "'" + fn + "', "
|
||||
cbString = cbString + " Relative Error along "
|
||||
if meanAxis == 0:
|
||||
cbString = cbString + "I axis"
|
||||
elif meanAxis == 1:
|
||||
cbString = cbString + "J axis"
|
||||
elif meanAxis == 2:
|
||||
cbString = cbString + "K axis"
|
||||
kv.genCB(AxCB,normAbs,cbString,cM=cmapAbs)
|
||||
#Now do main plotting
|
||||
if (doClear):
|
||||
Ax.clear()
|
||||
dataRel = CalcTotalErrRel(gsphP,gsphO,nStp,fieldNames,doVerb=doVerb,meanAxis=meanAxis)
|
||||
dataRel = np.transpose(dataRel) # transpose to put I/J on the horizontal axis
|
||||
Ax.pcolormesh(dataRel,cmap=cmapRel,norm=normRel)
|
||||
|
||||
if (doDeco):
|
||||
if meanAxis == 0:
|
||||
Ax.set_xlabel('J Indices')
|
||||
Ax.set_ylabel('K Indices')
|
||||
elif meanAxis == 1:
|
||||
Ax.set_xlabel('I Indices')
|
||||
Ax.set_ylabel('K Indices')
|
||||
elif meanAxis == 2:
|
||||
Ax.set_xlabel('I Indices')
|
||||
Ax.set_ylabel('J Indices')
|
||||
return dataRel
|
||||
|
||||
#Calculate total cumulative absolute error between two cases
|
||||
def CalcTotalErrAbs(gsphP,gsphO,nStp,fieldNames,doVerb=True,meanAxis=None):
|
||||
dataAbs = None
|
||||
for fn in fieldNames:
|
||||
dataP = gsphP.GetVar(fn,nStp,doVerb=doVerb)
|
||||
dataO = gsphO.GetVar(fn,nStp,doVerb=doVerb)
|
||||
if dataAbs is None:
|
||||
dataAbs = np.square(dataO - dataP)
|
||||
else:
|
||||
dataAbs = dataAbs + np.square(dataO - dataP)
|
||||
dataAbs = np.sqrt(dataAbs)
|
||||
return np.mean(dataAbs,axis=meanAxis)
|
||||
|
||||
#Calculate total cumulative relative error between two cases
|
||||
def CalcTotalErrRel(gsphP,gsphO,nStp,fieldNames,doVerb=True,meanAxis=None):
|
||||
dataAbs = None
|
||||
dataBase = None
|
||||
for fn in fieldNames:
|
||||
dataP = gsphP.GetVar(fn,nStp,doVerb=doVerb)
|
||||
dataO = gsphO.GetVar(fn,nStp,doVerb=doVerb)
|
||||
if dataAbs is None:
|
||||
dataBase = np.square(dataP)
|
||||
dataAbs = np.square(dataO - dataP)
|
||||
else:
|
||||
dataBase = dataBase + np.square(dataP)
|
||||
dataAbs = dataAbs + np.square(dataO - dataP)
|
||||
dataBase = np.sqrt(dataBase)
|
||||
dataAbs = np.sqrt(dataAbs)
|
||||
# math breaks when Base field is exactly 0
|
||||
dataBase[dataBase == 0] = np.finfo(np.float32).tiny
|
||||
dataRel = np.absolute(dataAbs/dataBase)
|
||||
return np.mean(dataRel,axis=meanAxis)
|
||||
|
||||
#Plot equatorial field
|
||||
def PlotEqB(gsph,nStp,xyBds,Ax,AxCB=None,doClear=True,doDeco=True,doBz=False):
|
||||
vBZ = kv.genNorm(dbMax)
|
||||
|
||||
211
scripts/quicklook/gamerrVid.py
Executable file
211
scripts/quicklook/gamerrVid.py
Executable file
@@ -0,0 +1,211 @@
|
||||
#!/usr/bin/env python
|
||||
#Make video of error between two Gamera cases
|
||||
import argparse
|
||||
from argparse import RawTextHelpFormatter
|
||||
import matplotlib as mpl
|
||||
mpl.use('Agg')
|
||||
import matplotlib.pyplot as plt
|
||||
import kaipy.kaiViz as kv
|
||||
import matplotlib.gridspec as gridspec
|
||||
import numpy as np
|
||||
import numpy as np
|
||||
import kaipy.gamera.msphViz as mviz
|
||||
import kaipy.gamera.magsphere as msph
|
||||
import kaipy.gamera.rcmpp as rcmpp
|
||||
from alive_progress import alive_bar
|
||||
import kaipy.kdefs as kdefs
|
||||
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
|
||||
import os
|
||||
import errno
|
||||
|
||||
cLW = 0.25
|
||||
|
||||
if __name__ == "__main__":
|
||||
#Defaults
|
||||
fdir1 = os.getcwd()
|
||||
ftag1 = "msphere"
|
||||
fdir2 = os.getcwd()
|
||||
ftag2 = "msphere"
|
||||
oDir = "vid2D"
|
||||
ts = 0 #[min]
|
||||
te = 200 #[min]
|
||||
dt = 60.0 #[sec]
|
||||
Nblk = 1 #Number of blocks
|
||||
nID = 1 #Block ID of this job
|
||||
noMPI = False # Don't add MPI tiling
|
||||
noLog = False
|
||||
fieldNames = "Bx, By, Bz"
|
||||
doVerb = False
|
||||
|
||||
MainS = """Creates simple multi-panel figure for Gamera magnetosphere run
|
||||
Left Panel - Residual vertical magnetic field
|
||||
Right Panel - Pressure (or density) and hemispherical insets
|
||||
"""
|
||||
|
||||
parser = argparse.ArgumentParser(description=MainS, formatter_class=RawTextHelpFormatter)
|
||||
parser.add_argument('-d1',type=str,metavar="directory",default=fdir1,help="Directory to read first dataset from (default: %(default)s)")
|
||||
parser.add_argument('-id1',type=str,metavar="runid",default=ftag1,help="RunID of first dataset (default: %(default)s)")
|
||||
parser.add_argument('-d2',type=str,metavar="directory",default=fdir2,help="Directory to read second dataset from (default: %(default)s)")
|
||||
parser.add_argument('-id2',type=str,metavar="runid",default=ftag2,help="RunID of second dataset (default: %(default)s)")
|
||||
parser.add_argument('-o',type=str,metavar="directory",default=oDir,help="Subdirectory to write to (default: %(default)s)")
|
||||
parser.add_argument('-ts' ,type=int,metavar="tStart",default=ts,help="Starting time [min] (default: %(default)s)")
|
||||
parser.add_argument('-te' ,type=int,metavar="tEnd" ,default=te,help="Ending time [min] (default: %(default)s)")
|
||||
parser.add_argument('-dt' ,type=int,metavar="dt" ,default=dt,help="Cadence [sec] (default: %(default)s)")
|
||||
parser.add_argument('-Nblk' ,type=int,metavar="Nblk",default=Nblk,help="Number of job blocks (default: %(default)s)")
|
||||
parser.add_argument('-nID' ,type=int,metavar="nID" ,default=nID,help="Block ID of this job [1-Nblk] (default: %(default)s)")
|
||||
parser.add_argument('-f',type=str,metavar="fieldnames",default=fieldNames,help="Comma-separated fields to plot (default: %(default)s)")
|
||||
parser.add_argument('-linear',action='store_true', default=noLog,help="Plot linear line plot instead of logarithmic (default: %(default)s)")
|
||||
parser.add_argument('-v',action='store_true', default=doVerb,help="Do verbose output (default: %(default)s)")
|
||||
#parser.add_argument('-nompi', action='store_true', default=noMPI,help="Don't show MPI boundaries (default: %(default)s)")
|
||||
|
||||
mviz.AddSizeArgs(parser)
|
||||
|
||||
#Finalize parsing
|
||||
args = parser.parse_args()
|
||||
fdir1 = args.d1
|
||||
ftag1 = args.id1
|
||||
fdir2 = args.d2
|
||||
ftag2 = args.id2
|
||||
ts = args.ts
|
||||
te = args.te
|
||||
dt = args.dt
|
||||
oSub = args.o
|
||||
Nblk = args.Nblk
|
||||
nID = args.nID
|
||||
fieldNames = args.f
|
||||
noLog = args.linear
|
||||
doVerb = args.v
|
||||
#noMPI = args.noMPI
|
||||
|
||||
fnList = [item.strip() for item in fieldNames.split(',')]
|
||||
|
||||
#Setup timing info
|
||||
tOut = np.arange(ts*60.0,te*60.0,dt)
|
||||
Nt = len(tOut)
|
||||
vO = np.arange(0,Nt)
|
||||
|
||||
print("Writing %d outputs between minutes %d and %d"%(Nt,ts,te))
|
||||
if (Nblk>1):
|
||||
#Figure out work bounds
|
||||
dI = (Nt//Nblk)
|
||||
i0 = (nID-1)*dI
|
||||
i1 = i0+dI
|
||||
if (nID == Nblk):
|
||||
i1 = Nt #Make sure we get last bit
|
||||
print("\tBlock #%d: %d to %d"%(nID,i0,i1))
|
||||
else:
|
||||
i0 = 0
|
||||
i1 = Nt
|
||||
|
||||
#Setup output directory
|
||||
oDir = os.getcwd() + "/" + oSub
|
||||
print("Writing output to %s"%(oDir))
|
||||
|
||||
#Check/create directory if necessary
|
||||
if (not os.path.exists(oDir)):
|
||||
try:
|
||||
print("Creating directory %s"%(oDir))
|
||||
os.makedirs(oDir)
|
||||
except OSError as exc:
|
||||
if exc.errno == errno.EEXIST and os.path.isdir(oDir):
|
||||
pass
|
||||
else:
|
||||
raise
|
||||
|
||||
#Get domain size
|
||||
xyBds = mviz.GetSizeBds(args)
|
||||
|
||||
#---------
|
||||
#Figure parameters
|
||||
figSz = (12,7.5)
|
||||
|
||||
|
||||
#======
|
||||
#Init data
|
||||
gsph1 = msph.GamsphPipe(fdir1,ftag1)
|
||||
gsph2 = msph.GamsphPipe(fdir2,ftag2)
|
||||
|
||||
#======
|
||||
#Setup figure
|
||||
fig = plt.figure(figsize=figSz)
|
||||
gs = gridspec.GridSpec(5,2,height_ratios=[20,5,1,5,9],hspace=0.025)
|
||||
|
||||
|
||||
AxTL = fig.add_subplot(gs[0,0])
|
||||
AxTR = fig.add_subplot(gs[0,1])
|
||||
AxB = fig.add_subplot(gs[-1,0:2])
|
||||
AxB2 = AxB.twinx() # second plot on bottom axis
|
||||
|
||||
AxCT = fig.add_subplot(gs[2,0:2])
|
||||
|
||||
print(fig.axes)
|
||||
|
||||
errTimes = []
|
||||
errListRel = []
|
||||
errListAbs = []
|
||||
relColor = "tab:blue"
|
||||
absColor = "tab:orange"
|
||||
|
||||
#Loop over sub-range
|
||||
titstr = "Comparing '%s' to '%s'"%(fdir1,fdir2)
|
||||
with alive_bar(i1-i0,title=titstr.ljust(kdefs.barLab),length=kdefs.barLen,disable=doVerb) as bar:
|
||||
for i in range(i0,i1):
|
||||
#Convert time (in seconds) to Step #
|
||||
nStp = np.abs(gsph1.T-tOut[i]).argmin()+gsph1.s0
|
||||
if doVerb:
|
||||
print("Minute = %5.2f / Step = %d"%(tOut[i]/60.0,nStp))
|
||||
npl = vO[i]
|
||||
|
||||
AxTL.clear()
|
||||
AxTR.clear()
|
||||
AxB.clear()
|
||||
AxB2.clear()
|
||||
|
||||
#plot upper left msph error
|
||||
mviz.PlotEqErrRel(gsph1,gsph2,nStp,xyBds,AxTL,fnList,AxCB=AxCT,doVerb=doVerb)
|
||||
AxTL.set_title("Equatorial Slice of Relative Error")
|
||||
|
||||
#plot upper right k-axis error
|
||||
mviz.PlotLogicalErrRel(gsph1,gsph2,nStp,AxTR,fnList,2,doVerb=doVerb)
|
||||
AxTR.set_title("Per-Cell Relative Error along K-Axis")
|
||||
if (not noMPI):
|
||||
#plot I-MPI decomp on logical plot
|
||||
if(gsph2.Ri > 1):
|
||||
for im in range(gsph2.Ri):
|
||||
i0 = im*gsph2.dNi
|
||||
AxTR.plot([i0, i0],[0, gsph2.Nj],"deepskyblue",linewidth=0.25,alpha=0.5)
|
||||
#plot J-MPI decomp on logical plot
|
||||
if (gsph2.Rj>1):
|
||||
for jm in range(1,gsph2.Rj):
|
||||
j0 = jm*gsph2.dNj
|
||||
AxTR.plot([0, gsph2.Ni],[j0, j0],"deepskyblue",linewidth=0.25,alpha=0.5)
|
||||
|
||||
#plot bottom line plot
|
||||
errTimes.append(tOut[i]/60.0)
|
||||
errListRel.append(mviz.CalcTotalErrRel(gsph1,gsph2,nStp,fnList,doVerb=doVerb))
|
||||
errListAbs.append(mviz.CalcTotalErrAbs(gsph1,gsph2,nStp,fnList,doVerb=doVerb))
|
||||
if noLog:
|
||||
AxB.plot(errTimes, errListRel,color=relColor)
|
||||
AxB2.plot(errTimes, errListAbs,color=absColor)
|
||||
else:
|
||||
AxB.semilogy(errTimes, errListRel,color=relColor)
|
||||
AxB2.semilogy(errTimes, errListAbs,color=absColor)
|
||||
AxB.set_xlabel('Time (min)')
|
||||
AxB.set_ylabel('Per-Cell Mean Relative Error',color=relColor)
|
||||
AxB.tick_params(axis='y',which='both',colors=relColor,left=True,right=True,labelleft=True,labelright=False)
|
||||
AxB2.set_ylabel('Per-Cell Mean Absolute Error',color=absColor)
|
||||
#AxB2.yaxis.tick_right()
|
||||
AxB2.tick_params(axis='y',which='both',colors=absColor,left=True,right=True,labelleft=False,labelright=True)
|
||||
AxB.set_title("'" + fieldNames + "' Per-Cell Error Over Time")
|
||||
|
||||
gsph1.AddTime(nStp,AxTL,xy=[0.025,0.84],fs="x-large")
|
||||
|
||||
#Add MPI decomp
|
||||
if (not noMPI):
|
||||
mviz.PlotMPI(gsph2,AxTL)
|
||||
|
||||
fOut = oDir+"/vid.%04d.png"%(npl)
|
||||
kv.savePic(fOut,bLenX=45)
|
||||
|
||||
bar()
|
||||
|
||||
99
scripts/quicklook/gamerrpic.py
Executable file
99
scripts/quicklook/gamerrpic.py
Executable file
@@ -0,0 +1,99 @@
|
||||
#!/usr/bin/env python
|
||||
#Make video of error between two Gamera cases
|
||||
import argparse
|
||||
from argparse import RawTextHelpFormatter
|
||||
import matplotlib as mpl
|
||||
mpl.use('Agg')
|
||||
import matplotlib.pyplot as plt
|
||||
import kaipy.kaiViz as kv
|
||||
import matplotlib.gridspec as gridspec
|
||||
import numpy as np
|
||||
import numpy as np
|
||||
import kaipy.gamera.msphViz as mviz
|
||||
import kaipy.gamera.magsphere as msph
|
||||
import kaipy.gamera.rcmpp as rcmpp
|
||||
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
|
||||
import os
|
||||
import errno
|
||||
|
||||
cLW = 0.25
|
||||
|
||||
if __name__ == "__main__":
|
||||
#Defaults
|
||||
fdir1 = os.getcwd()
|
||||
ftag1 = "msphere"
|
||||
fdir2 = os.getcwd()
|
||||
ftag2 = "msphere"
|
||||
nStp=1
|
||||
fieldNames = "Bx, By, Bz"
|
||||
doMPI = False #[Add MPI tiling]
|
||||
noMPI = False
|
||||
|
||||
MainS = """Creates simple multi-panel figure for Gamera magnetosphere run
|
||||
Left Panel - Residual vertical magnetic field
|
||||
Right Panel - Pressure (or density) and hemispherical insets
|
||||
"""
|
||||
|
||||
parser = argparse.ArgumentParser(description=MainS, formatter_class=RawTextHelpFormatter)
|
||||
parser.add_argument('-d1',type=str,metavar="directory1",default=fdir1,help="Directory to read first dataset from (default: %(default)s)")
|
||||
parser.add_argument('-id1',type=str,metavar="runid1",default=ftag1,help="RunID of first dataset (default: %(default)s)")
|
||||
parser.add_argument('-d2',type=str,metavar="directory2",default=fdir2,help="Directory to read second dataset from (default: %(default)s)")
|
||||
parser.add_argument('-id2',type=str,metavar="runid2",default=ftag2,help="RunID of second dataset (default: %(default)s)")
|
||||
parser.add_argument('-n',type=int,metavar="nStp",default=nStp,help="Step number to plot (default: %(default)s)")
|
||||
parser.add_argument('-f',type=str,metavar="fieldnames",default=fieldNames,help="Comma-separated fields to plot (default: %(default)s)")
|
||||
#parser.add_argument('-nompi', action='store_true', default=noMPI,help="Don't show MPI boundaries (default: %(default)s)")
|
||||
|
||||
mviz.AddSizeArgs(parser)
|
||||
|
||||
#Finalize parsing
|
||||
args = parser.parse_args()
|
||||
fdir1 = args.d1
|
||||
ftag1 = args.id1
|
||||
fdir2 = args.d2
|
||||
ftag2 = args.id2
|
||||
nStp = args.n
|
||||
fieldNames = args.f
|
||||
oName = "gamErrPic.png"
|
||||
|
||||
#Get domain size
|
||||
xyBds = mviz.GetSizeBds(args)
|
||||
|
||||
#---------
|
||||
#Figure parameters
|
||||
figSz = (12,7.5)
|
||||
|
||||
|
||||
#======
|
||||
#Init data
|
||||
gsph1 = msph.GamsphPipe(fdir1,ftag1)
|
||||
gsph2 = msph.GamsphPipe(fdir2,ftag2)
|
||||
|
||||
#======
|
||||
#Setup figure
|
||||
fig = plt.figure(figsize=figSz)
|
||||
gs = gridspec.GridSpec(2,2,height_ratios=[20,1],hspace=0.025)
|
||||
|
||||
|
||||
AxL = fig.add_subplot(gs[0,0])
|
||||
AxR = fig.add_subplot(gs[0,1])
|
||||
|
||||
AxCL = fig.add_subplot(gs[-1,0])
|
||||
AxCR = fig.add_subplot(gs[-1,1])
|
||||
|
||||
fnList = [item.strip() for item in fieldNames.split(',')]
|
||||
|
||||
AxL.clear()
|
||||
AxR.clear()
|
||||
|
||||
mviz.PlotEqErrRel(gsph1,gsph2,nStp,xyBds,AxL,fnList,AxCB=AxCL)
|
||||
mviz.PlotEqErrAbs(gsph1,gsph2,nStp,xyBds,AxR,fnList,AxCB=AxCR)
|
||||
|
||||
gsph1.AddTime(nStp,AxL,xy=[0.025,0.89],fs="x-large")
|
||||
|
||||
#Add MPI decomp
|
||||
if (doMPI):
|
||||
mviz.PlotMPI(gsph2,AxL)
|
||||
mviz.PlotMPI(gsph2,AxR)
|
||||
|
||||
kv.savePic(oName,bLenX=45)
|
||||
|
||||
@@ -48,15 +48,28 @@ if __name__ == "__main__":
|
||||
fOut = swtag.split('.')[0]+'.'+imgtype
|
||||
|
||||
# pulling UT variable for plotting
|
||||
t0Fmt = "%Y-%m-%d %H:%M:%S"
|
||||
t0Fmts = ["%Y-%m-%d %H:%M:%S","%Y-%m-%dT%H:%M:%S.%f"]
|
||||
utfmt='%H:%M \n%Y-%m-%d'
|
||||
|
||||
UTall = kh5.PullVar(swIn,"UT")
|
||||
|
||||
#Identify the correct time format
|
||||
t0Fmt = None
|
||||
for tfmt in t0Fmts:
|
||||
try:
|
||||
datetime.datetime.strptime(UTall[1].decode('utf-8'),tfmt)
|
||||
t0Fmt = tfmt
|
||||
break # datetime parse succeeded
|
||||
except ValueError:
|
||||
pass # datetime parse failed
|
||||
if t0Fmt is None:
|
||||
print("Time format in bcwind.h5 did not match any expected format.")
|
||||
sys.exit()
|
||||
|
||||
utall = []
|
||||
for n in range(len(UTall)):
|
||||
utall.append(datetime.datetime.strptime(UTall[n].decode('utf-8'),t0Fmt))
|
||||
|
||||
|
||||
# pulling the solar wind values from the table
|
||||
varlist = kh5.getRootVars(swIn)
|
||||
D = kh5.PullVar(swIn,"D")
|
||||
@@ -74,6 +87,6 @@ if __name__ == "__main__":
|
||||
else:
|
||||
pltInterp = 0*D
|
||||
doEps = False
|
||||
swBCplots.swQuickPlot(UTall,D,Temp,Vx,Vy,Vz,Bx,By,Bz,SYMH,pltInterp,fOut,doEps=doEps,doTrim=doTrim)
|
||||
swBCplots.swQuickPlot(UTall,D,Temp,Vx,Vy,Vz,Bx,By,Bz,SYMH,pltInterp,fOut,doEps=doEps,doTrim=doTrim,t0fmt=t0Fmt)
|
||||
|
||||
|
||||
|
||||
@@ -706,7 +706,9 @@ module streamline
|
||||
if (inDom) Np = Np + 1
|
||||
enddo
|
||||
|
||||
if (MaxFL == Np) then
|
||||
! check if exceeded tube bounds
|
||||
if(Np > MaxFL) then
|
||||
Np = MaxFL
|
||||
!$OMP CRITICAL
|
||||
write(*,*) ANSIRED
|
||||
write(*,*) "<WARNING! genTrace hit max tube size!>"
|
||||
|
||||
Reference in New Issue
Block a user