cwpearson/astaroth
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Can now picture the magnetic field and streamlines. And some other minor improvements.

Browse Source
This commit is contained in:
Miikka Vaisala
2019-11-04 11:27:53 +08:00
parent 8894b7c7d6
commit d3260edd2a
8 changed files with 157 additions and 1148 deletions
Download Patch File Download Diff File Expand all files Collapse all files

49
analysis/python/astar/data/read.py
View File

@@ -21,6 +21,7 @@
import numpy as np
def set_dtype(endian, AcRealSize):
if endian == 0:
en = '>'
@@ -83,6 +84,43 @@ def read_meshtxt(fdir, fname):
return contents
#Now just 2nd order
def DERX(array, dx):
output = np.zeros_like(array)
for i in range(3, array.shape[0]-3): #Keep boundary poits as 0
#output[i,:,:] = (- array[i-1,:,:] + array[i+1,:,:])/(2.0*dx)
output[i,:,:] =( -45.0*array[i-1,:,:] + 45.0*array[i+1,:,:]
+ 9.0*array[i-2,:,:] - 9.0*array[i+2,:,:]
- array[i-3,:,:] + array[i+3,:,:] )/(60.0*dx)
return output
def DERY(array, dy):
output = np.zeros_like(array)
for i in range(3,array.shape[1]-3):
#output[:,i,:] = (- array[:,i-1,:] + array[:,i+1,:])/(2.0*dy)
output[:,i,:] =( -45.0*array[:,i-1,:] + 45.0*array[:,i+1,:]
+ 9.0*array[:,i-2,:] - 9.0*array[:,i+2,:]
- array[:,i-3,:] + array[:,i+3,:] )/(60.0*dy)
return output
def DERZ(array, dz):
output = np.zeros_like(array)
for i in range(3, array.shape[2]-3):
#output[:,:,i] = (- array[:,:,i-1] + array[:,:,i+1])/(2.0*dz)
output[:,:,i] =( -45.0*array[:,:,i-1] + 45.0*array[:,:,i+1]
+ 9.0*array[:,:,i-2] - 9.0*array[:,:,i+2]
- array[:,:,i-3] + array[:,:,i+3] )/(60.0*dz)
return output
def curl(aa, minfo):
dx = minfo.contents['AC_dsx']
dy = minfo.contents['AC_dsy']
dz = minfo.contents['AC_dsz']
return (DERY(aa[2], dy)-DERZ(aa[1], dz),
DERZ(aa[0], dz)-DERX(aa[2], dx),
DERX(aa[1], dx)-DERY(aa[0], dy))
class MeshInfo():
'''Object that contains all mesh info'''
@@ -123,6 +161,14 @@ class Mesh:
aay = []
aaz = []
#self.aa[0][:,:,:] = 0.0
#self.aa[1][:,:,:] = 0.0
#self.aa[2][:,:,:] = 0.0
#for i in range(0, self.aa[0].shape[0]):
# self.aa[0][:,i,:] = float(i)
self.xx = np.arange(self.minfo.contents['AC_mx']) * self.minfo.contents['AC_dsx']
self.yy = np.arange(self.minfo.contents['AC_my']) * self.minfo.contents['AC_dsy']
self.zz = np.arange(self.minfo.contents['AC_mz']) * self.minfo.contents['AC_dsz']
@@ -130,6 +176,9 @@ class Mesh:
self.xmid = int(self.minfo.contents['AC_mx']/2)
self.ymid = int(self.minfo.contents['AC_my']/2)
self.zmid = int(self.minfo.contents['AC_mz']/2)
def Bfield(self):
self.bb = curl(self.aa, self.minfo)
print(self.bb[2])
def parse_ts(fdir, fname):

37
analysis/python/astar/visual/lineplot.py
View File

@@ -26,16 +26,22 @@ CM_INFERNO = plt.get_cmap('inferno')
end_rm = -1 #-35#-40
def plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3):
plt.plot(ts.var[xaxis][:end_rm], ts.var[yaxis1][:end_rm], label=yaxis1)
plt.plot(ts.var[xaxis][:end_rm], ts.var[yaxis2][:end_rm], label=yaxis2)
plt.plot(ts.var[xaxis][:end_rm], ts.var[yaxis3][:end_rm], label=yaxis3)
plt.xlabel(xaxis)
plt.legend()
if yaxis1 in ts.var:
plt.plot(ts.var[xaxis][:end_rm], ts.var[yaxis1][:end_rm], label=yaxis1)
plt.plot(ts.var[xaxis][:end_rm], ts.var[yaxis2][:end_rm], label=yaxis2)
plt.plot(ts.var[xaxis][:end_rm], ts.var[yaxis3][:end_rm], label=yaxis3)
plt.xlabel(xaxis)
plt.legend()
else:
print("%s %s and %s not found! Skipping...", yaxis1, yaxis2, yaxis3)
plt.close()
def plot_ts(ts, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False, uuz=False, ss=False, acc=False, sink=False):
def plot_ts(ts, isotherm=False, show_all=False, lnrho=False, uutot=False,
uux=False, uuy=False, uuz=False, ss=False, acc=False, sink=False):
if show_all:
lnrho=True
rho=True
uutot=True
uux=True
uuy=True
@@ -44,6 +50,17 @@ def plot_ts(ts, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False,
acc=True
sink=True
if isotherm:
lnrho=True
rho=True
uutot=True
uux=True
uuy=True
uuz=True
ss=False
acc=True
sink=True
if lnrho:
plt.figure()
xaxis = 't_step'
@@ -51,6 +68,14 @@ def plot_ts(ts, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False,
yaxis2 = 'lnrho_min'
yaxis3 = 'lnrho_max'
plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3)
if rho:
plt.figure()
xaxis = 't_step'
yaxis1 = 'rho_rms'
yaxis2 = 'rho_min'
yaxis3 = 'rho_max'
plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3)
if uutot:
plt.figure()

44
analysis/python/astar/visual/slices.py
View File

@@ -24,7 +24,9 @@ import matplotlib.colors as colors
CM_INFERNO = plt.get_cmap('inferno')
def plot_3(mesh, input_grid, title = '', fname = 'default', bitmap=False, slicetype = 'middle', colrange=None, colormap=CM_INFERNO , contourplot=False):
def plot_3(mesh, input_grid, title = '', fname = 'default', bitmap=False,
slicetype = 'middle', colrange=None, colormap=CM_INFERNO ,
contourplot=False, points_from_centre = -1, bfieldlines=False, velfieldlines=False):
fig = plt.figure(figsize=(8, 8))
grid = gridspec.GridSpec(2, 3, wspace=0.4, hspace=0.4, width_ratios=[1,1, 0.15])
ax00 = fig.add_subplot( grid[0,0] )
@@ -52,7 +54,17 @@ def plot_3(mesh, input_grid, title = '', fname = 'default', bitmap=False, slicet
plotnorm = colors.Normalize(vmin=cmin,vmax=cmax)
else:
plotnorm = colors.Normalize(vmin=colrange[0],vmax=colrange[1])
if points_from_centre > 0:
yz_slice = yz_slice[int(yz_slice.shape[0]/2)-points_from_centre : int(yz_slice.shape[0]/2)+points_from_centre,
int(yz_slice.shape[1]/2)-points_from_centre : int(yz_slice.shape[1]/2)+points_from_centre]
xz_slice = xz_slice[int(xz_slice.shape[0]/2)-points_from_centre : int(xz_slice.shape[0]/2)+points_from_centre,
int(xz_slice.shape[1]/2)-points_from_centre : int(xz_slice.shape[1]/2)+points_from_centre]
xy_slice = xy_slice[int(xy_slice.shape[0]/2)-points_from_centre : int(xy_slice.shape[0]/2)+points_from_centre,
int(xy_slice.shape[1]/2)-points_from_centre : int(xy_slice.shape[1]/2)+points_from_centre]
mesh.xx = mesh.xx[int(mesh.xx.shape[0]/2)-points_from_centre : int(mesh.xx.shape[0]/2)+points_from_centre]
mesh.yy = mesh.yy[int(mesh.yy.shape[0]/2)-points_from_centre : int(mesh.yy.shape[0]/2)+points_from_centre]
mesh.zz = mesh.zz[int(mesh.zz.shape[0]/2)-points_from_centre : int(mesh.zz.shape[0]/2)+points_from_centre]
yy, zz = np.meshgrid(mesh.yy, mesh.zz, indexing='ij')
if contourplot:
@@ -63,6 +75,8 @@ def plot_3(mesh, input_grid, title = '', fname = 'default', bitmap=False, slicet
ax00.set_ylabel('z')
ax00.set_title('%s t = %.4e' % (title, mesh.timestamp) )
ax00.set_aspect('equal')
ax00.contour(yy, zz, np.sqrt((yy-yy.max()/2.0)**2.0 + (zz-zz.max()/2.0)**2.0), [mesh.minfo.contents["AC_accretion_range"]])
xx, zz = np.meshgrid(mesh.xx, mesh.zz, indexing='ij')
if contourplot:
@@ -72,6 +86,8 @@ def plot_3(mesh, input_grid, title = '', fname = 'default', bitmap=False, slicet
ax10.set_xlabel('x')
ax10.set_ylabel('z')
ax10.set_aspect('equal')
ax10.contour(xx, zz, np.sqrt((xx-xx.max()/2.0)**2.0 + (zz-zz.max()/2.0)**2.0), [mesh.minfo.contents["AC_accretion_range"]])
xx, yy = np.meshgrid(mesh.xx, mesh.yy, indexing='ij')
if contourplot:
@@ -81,6 +97,30 @@ def plot_3(mesh, input_grid, title = '', fname = 'default', bitmap=False, slicet
ax11.set_xlabel('x')
ax11.set_ylabel('y')
ax11.set_aspect('equal')
ax11.contour(xx, yy, np.sqrt((xx-xx.max()/2.0)**2.0 + (yy-yy.max()/2.0)**2.0), [mesh.minfo.contents["AC_accretion_range"]])
if bfieldlines:
ax00.streamplot(mesh.yy, mesh.zz, np.mean(mesh.bb[1], axis=0), np.mean(mesh.bb[2], axis=0))
ax10.streamplot(mesh.xx, mesh.zz, np.mean(mesh.bb[0], axis=1), np.mean(mesh.bb[2], axis=1))
ax11.streamplot(mesh.xx, mesh.yy, np.mean(mesh.bb[0], axis=2), np.mean(mesh.bb[1], axis=2))
#ax00.streamplot(mesh.yy, mesh.zz, mesh.bb[1][mesh.xmid, :, :], mesh.bb[2][mesh.xmid, :, :])
#ax10.streamplot(mesh.xx, mesh.zz, mesh.bb[0][:, mesh.ymid, :], mesh.bb[2][:, mesh.ymid, :])
#ax11.streamplot(mesh.xx, mesh.yy, mesh.bb[0][:, : ,mesh.zmid], mesh.bb[1][:, :, mesh.zmid])
#ax00.quiver(mesh.bb[2][mesh.xmid, ::10, ::10], mesh.bb[1][mesh.xmid, ::10, ::10], pivot='middle')
#ax10.quiver(mesh.bb[2][::10, mesh.ymid, ::10], mesh.bb[0][::10, mesh.ymid, ::10], pivot='middle')
#ax11.quiver(mesh.bb[1][::10, ::10, mesh.zmid], mesh.bb[0][::10, ::10, mesh.zmid], pivot='middle')
#ax00.quiver(mesh.yy, mesh.zz, mesh.bb[2][mesh.xmid, :, :], mesh.bb[1][mesh.xmid, :, :], pivot='middle')
#ax10.quiver(mesh.xx, mesh.zz, mesh.bb[2][:, mesh.ymid, :], mesh.bb[0][:, mesh.ymid, :], pivot='middle')
#ax11.quiver(mesh.xx, mesh.yy, mesh.bb[1][:, :, mesh.zmid], mesh.bb[0][:, :, mesh.zmid], pivot='middle')
if velfieldlines:
ax00.streamplot(mesh.yy, mesh.zz, mesh.uu[2][mesh.xmid, :, :], mesh.uu[1][mesh.xmid, :, :])
ax10.streamplot(mesh.xx, mesh.zz, mesh.uu[2][:, mesh.ymid, :], mesh.uu[0][:, mesh.ymid, :])
ax11.streamplot(mesh.xx, mesh.yy, mesh.uu[1][:, :, mesh.zmid], mesh.uu[0][:, : ,mesh.zmid])
cbar = plt.colorbar(map1, cax=axcbar)