Skip to content

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

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

Open
soumide1102 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 all 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
36 changes: 23 additions & 13 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 @@ -326,32 +326,42 @@ def quiver_xy(ax, geom, dump, varx, vary, C=None,
coordinates='figure')

def overlay_field(ax, geom, dump, NLEV=20, linestyle='-', linewidth=1,
linecolor='k'):
linecolor='k', zorder=2):
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

@Yurlungur Yurlungur Oct 20, 2020

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

@soumide1102 soumide1102 Oct 20, 2020

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=zorder)

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
26 changes: 26 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,24 @@
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('--zorder',
type=int, default=2,
help='zorder for 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 +119,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

@Yurlungur Yurlungur Oct 20, 2020

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, zorder=args.zorder)
ax = a1
if coords == 'mks':
bplt.plot_X1X3(ax, geom, var, dump, cmap=cmap, vmin=vmin, vmax=vmax,
Expand All @@ -120,6 +142,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, zorder=args.zorder)

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