Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add support for plotting isocontours of vector potential in snap.py #7

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Add support for plotting isocontours of vector potential in snap.py #7

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f2f57d8
Add support for plotting isocontours of vector potential in snap.py
soumide1102 Oct 20, 2020
70a66c0
Adjust zorder
soumide1102 Oct 20, 2020
fecb200
Make zorder an input param in snap.py and pass it to overlay_field
soumide1102 Oct 20, 2020
6e395d6
correct type and default for zorder
soumide1102 Oct 20, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 22 additions & 12 deletions script/analysis/plot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -123,7 +123,8 @@ def plot_X1X3(ax, geom, var, dump, cmap='jet', vmin=None, vmax=None, cbar=True,
if not show_axes:
ax.set_xticklabels([]); ax.set_yticklabels([])

def plot_xz(ax, geom, var, dump, cmap='jet', vmin=None, vmax=None, cbar=True,

def plot_xz(ax, geom, var, dump, cmap='jet', vmin=None, vmax=None, cbar=True,
label=None, ticks=None, shading='gouraud',
l_scale = None,
reverse_x = False,
Expand Down Expand Up @@ -172,7 +173,6 @@ def plot_xz(ax, geom, var, dump, cmap='jet', vmin=None, vmax=None, cbar=True,
ax.set_aspect('equal')
ax.set_xlabel('x/M'); ax.set_ylabel('z/M')
return mesh
#ax.grid(True, linestyle=':', color='k', alpha=0.5, linewidth=0.5)

def contour_xz(ax, geom, var, dump,
levels = None,
Expand Down Expand Up @@ -329,29 +329,39 @@ def overlay_field(ax, geom, dump, NLEV=20, linestyle='-', linewidth=1,
linecolor='k'):
from scipy.integrate import trapz
hdr = dump['hdr']
N1 = hdr['N1']; N2 = hdr['N2']
x = flatten_xz(geom['x'], hdr).transpose()
z = flatten_xz(geom['z'], hdr).transpose()
N1 = hdr['N1']//2; N2 = hdr['N2']
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are you dividing hdr['N1'] by 2?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, this is the line that I wanted to check with you.

Is A_phi[j,N1-1-i] calculating A in the -ve direction and A_phi[j,i+N1] in the +ve direction? So if I don't divide by 2, ultimately, in line 363, rcyl and z ends up having dimensions 256 X 256, and A_phi is 256 X 512. So I divided by 2 to make N1 represent cells in either the +ve/-ve direction only.

x = geom['x']
y = geom['y']
z = geom['z']
if dump['hdr']['N3'] > 1. and dump['hdr']['stopx'][3] >= np.pi:
x = flatten_xz(x, dump['hdr'], flip=True)
y = flatten_xz(y, dump['hdr'], flip=True)
z = flatten_xz(z, dump['hdr'])
rcyl = np.sqrt(x**2 + y**2)
rcyl[np.where(x<0)] *= -1
else:
x = x[:,:,0]
y = y[:,:,0]
z = z[:,:,0]
rcyl = np.sqrt(x**2 + y**2)
A_phi = np.zeros([N2, 2*N1])
gdet = geom['gdet'].transpose()
B1 = dump['B1'].mean(axis=-1).transpose()
B2 = dump['B2'].mean(axis=-1).transpose()
print(gdet.shape)
for j in range(N2):
for i in range(N1):
A_phi[j,N1-1-i] = (trapz(gdet[j,:i]*B2[j,:i], dx=hdr['dx'][1]) -
A_phi[j,N1-1-i] = (trapz(gdet[j,:i]*B2[j,:i], dx=hdr['dx'][1]) -
trapz(gdet[:j, i]*B1[:j, i], dx=hdr['dx'][2]))
A_phi[j,i+N1] = (trapz(gdet[j,:i]*B2[j,:i], dx=hdr['dx'][1]) -
A_phi[j,i+N1] = (trapz(gdet[j,:i]*B2[j,:i], dx=hdr['dx'][1]) -
trapz(gdet[:j, i]*B1[:j, i], dx=hdr['dx'][2]))
A_phi -= (A_phi[N2//2-1,-1] + A_phi[N2//2,-1])/2.
Apm = np.fabs(A_phi).max()
if np.fabs(A_phi.min()) > A_phi.max():
A_phi *= -1.
#NLEV = 20
levels = np.concatenate((np.linspace(-Apm,0,NLEV)[:-1],
levels = np.concatenate((np.linspace(-Apm,0,NLEV)[:-1],
np.linspace(0,Apm,NLEV)[1:]))
ax.contour(x, z, A_phi, levels=levels, colors=linecolor, linestyles=linestyle,
linewidths=linewidth)
ax.contour(rcyl, z, A_phi, levels=levels, colors=linecolor, linestyles=linestyle,
linewidths=linewidth, zorder=2)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you may want to make zorder an optional parameter that can be passed in to the function, in case the user wants to layer them some different way.


def plot_xy(ax, geom, var, dump, cmap='jet', vmin=None, vmax=None, cbar=True,
label=None, ticks=None, shading='gouraud', fix_bounds=True):
Expand Down
23 changes: 23 additions & 0 deletions script/analysis/snap.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,21 @@
parser.add_argument('-c','--cmap',
type=str,default='jet',
help='Colormap used')
parser.add_argument('--A-contours', action="store_true",
default=False,
help='Flag to overlay contours of the vector potential')
parser.add_argument('--nlev',
type=int, default=20,
help='Number of positive/negative A-contour levels. Needed only if using --A-contours')
parser.add_argument('--linestyle',
type=str, default='-',
help='Linestyle of A-contours. Needed only if using --A-contours')
parser.add_argument('--linewidth',
type=int, default=1,
help='Linewidth of A-contours. Needed only if using --A-contours')
parser.add_argument('--linecolor',
type=str, default='k',
help='Linecolor of A-contours. Needed only if using --A-contours')
parser.add_argument('--save',
type=str,default=None,
help='Figure filename if you want to save the figure')
Expand Down Expand Up @@ -101,6 +116,10 @@ def make_snap(dfnam,vnam,coords,size,cmap,logplot,
bplt.plot_xz(ax, geom, var, dump, cmap=cmap, vmin=vmin, vmax=vmax,
cbar=False, label=label, ticks=None, shading='gouraud')
ax.set_xlim([-size,size]); ax.set_ylim([-size,size])
if args.A_contours:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I just noticed that args is referenced in make_snap now. This is going to break movie.py which calls this method. Can you pass these through as optional arguments to make_snap so that the args object isn't used in scope it might not be available?

bplt.overlay_field(ax, geom, dump, NLEV=args.nlev,
linestyle=args.linestyle, linewidth=args.linewidth,
linecolor=args.linecolor)
ax = a1
if coords == 'mks':
bplt.plot_X1X3(ax, geom, var, dump, cmap=cmap, vmin=vmin, vmax=vmax,
Expand All @@ -120,6 +139,10 @@ def make_snap(dfnam,vnam,coords,size,cmap,logplot,
bplt.plot_xz(ax, geom, var, dump, cmap=cmap, vmin=vmin, vmax=vmax,
cbar=True, label=label, ticks=None, shading='gouraud')
ax.set_xlim([0,size]); ax.set_ylim([-size,size])
if args.A_contours:
bplt.overlay_field(ax, geom, dump, NLEV=args.nlev,
linestyle=args.linestyle, linewidth=args.linewidth,
linecolor=args.linecolor)

if savefig == False:
plt.show()
Expand Down