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

Merge branch 'master' of https://bitbucket.org/jpekkila/astaroth

Browse Source
This commit is contained in:
Johannes Pekkila
2019-11-27 08:55:23 +01:00
parent 981331e7d7 5e3caf086e
commit 3eabf94f92
10 changed files with 186 additions and 1153 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):

27
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):
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'
@@ -52,6 +69,14 @@ def plot_ts(ts, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False,
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()
xaxis = 't_step'

42
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] )
@@ -53,6 +55,16 @@ def plot_3(mesh, input_grid, title = '', fname = 'default', bitmap=False, slicet
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:
@@ -64,6 +76,8 @@ def plot_3(mesh, input_grid, title = '', fname = 'default', bitmap=False, slicet
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:
ax10.contourf(xx, zz, xz_slice, norm=plotnorm, cmap=colormap, nlev=10)
@@ -73,6 +87,8 @@ def plot_3(mesh, input_grid, title = '', fname = 'default', bitmap=False, slicet
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:
ax11.contourf(xx, yy, xy_slice, norm=plotnorm, cmap=colormap, nlev=10)
@@ -82,6 +98,30 @@ def plot_3(mesh, input_grid, title = '', fname = 'default', bitmap=False, slicet
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)
if bitmap:

9
analysis/python/calc/convert.sh
View File

@@ -1,9 +0,0 @@
#!/bin/bash
#gm convert -delay 40 colden_*.png colden.gif
DATE=`date '+%Y_%m_%d_%H_%M'`
echo $DATE
gm convert -delay 15 $1_*.png $1_$DATE.gif

835
analysis/python/calc/galli_shu_plotter.py
View File

@@ -1,835 +0,0 @@
'''
Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
This file is part of Astaroth.
Astaroth is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
Astaroth is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with Astaroth. If not, see <http://www.gnu.org/licenses/>.
'''
import numpy as np
import pylab as plt
import scipy as scp
import matplotlib.colors as colors
G_newton = 6.674e-8 #cm**3 g**-1 s**-2
# Time to convert to physical quantities
yr = 3.154e+7 #s
kyr = 1000.0*yr
km = 1e5 #cm
AU = 1.496e+13 #cm
Msun = 1.98847e33 #g
#cs0 = 20000.0 #cs cm/s "a" in Shu notation
cs0 = 35000.0 #cs cm/s "a" in Shu notation
B0 = 30e-6 #G
ksii = 11.3 #
#GS Eq. 10
ttm = 9.03e12*(cs0/35000.0)/(B0/30e-6)
CM_INFERNO = plt.cm.get_cmap('inferno')
def P_harmonics(theta, J=666):
#Vector spherical harmonics in e_r direction
if J == 0:
P = np.ones_like(theta) # 1.0
elif J == 2:
cos_theta = np.cos(theta)
P = (1.0/2.0)*(3.0*(cos_theta**2.0) - 1.0)
else:
P = 0.0
#print("P_2", P)
return P
def B_harmonics(theta, J=666):
#Vector spherical harmonics in e_theta direction
#print("B_harmonics theta", theta)
if J == 2:
sin_theta = np.abs(np.sin(theta))
cos_theta = np.cos(theta)
#B = -(3.0/np.sqrt(6.0))*cos_theta*sin_theta #Morse & Feshbach 1953 book
B = -3.0*cos_theta*sin_theta #GS93 Appendix B
else:
B = 0.0*theta
#print("B_harmonics", B)
return B
def get_tau(tt):
return tt/ttm
def get_SHU77_potential(xx_point):
#Copied here again for convenience
m0 = 0.975 #Shu 77 core reduced mass
xx_SHU_table = np.array([ 0.05, 0.10, 0.15, 0.20, 0.25,
0.30, 0.35, 0.40, 0.45, 0.50,
0.55, 0.60, 0.65, 0.70, 0.75,
0.80, 0.85, 0.90, 0.95, 1.00])
mm_SHU77_table = np.array([0.981, 0.993, 1.01, 1.03, 1.05,
1.08, 1.12, 1.16, 1.20, 1.25,
1.30, 1.36, 1.42, 1.49, 1.56,
1.64, 1.72, 1.81, 1.90, 2.00])
xx = xx_SHU_table[ np.where(xx_SHU <= xx_point)]
mm = mm_SHU77_table[np.where(xx_SHU <= xx_point)]
psi = - m0/xx_point + np.trapz(mm/(xx**2.0), xx)
return psi
def psi2(xx_SHU, mm_term, pp_term, J=666):
#GS93 Eq. 113
if J == 0:
psi2 = - mm_term/xx_SHU + pp_term
elif J == 2:
psi2 = - mm_term/(xx_SHU**3.0) + (xx_SHU**2.0)*pp_term
else:
psi2 = 0.0
#print('psi2', psi2, 'J', J, 'mm_term', mm_term, 'xx_SHU', xx_SHU, 'pp_term', pp_term)
return psi2
# Calculate the directional parameter
def dv_dx(xx,vv, alpha):
EE = alpha*(xx-vv) - 2.0/xx
HH = (xx-vv)**2.0 - 1.0
return (EE/HH)*(xx-vv)
def dalpha_dx(xx,vv, alpha):
EE = alpha*(alpha - (2.0/xx)*(xx-vv))
HH = (xx-vv)**2.0 - 1.0
return (EE/HH)*(xx-vv)
def dpsi_dx(xx, mm):
return mm/(xx**2.0)
def dmm_dx(xx, alpha):
return (xx**2.0)*alpha
def dphi_dx(xx, alpha, mm, theta):
ff_zero_der = 0.5*mm*dmm_dx(xx, alpha)
sin_theta = np.sin(theta)
return ff_zero_der*(sin_theta*2.0)
def deltaspace(theta, tau):
#Assuming J= 0, 2 only
v0 = -2.222e-1
v2 = 2.177e-1
deltaJ2 = -(1.0/3.0)*((v0+2.0/3.0)*P_harmonics(theta, J=0) + (v2 - 2.0/3.0)*P_harmonics(theta, J=2))
delta = 1 + (tau**2.0)*deltaJ2
return delta
def delta2(theta, tau):
#Assuming J= 0, 2 only
return deltaspace(theta, tau)**2.0
def yy_transform(xx_SHU, alpha_SHU77, alpha_mono_GS93, alpha_quad_GS93):
return alpha_mono_GS93, alpha_quad_GS93
# Calculating the perturbation stage
def alpha_perturb(tau, xx_SHU, vv_SHU77, alpha_SHU77, alpha_mono_GS93, alpha_quad_GS93, theta):
#Assuming J= 0, 2 only
directional = xx_SHU*dalpha_dx(xx_SHU, vv_SHU77, alpha_SHU77)*delta2(theta, tau)
directional = 0.0 #
alpha = alpha_mono_GS93*P_harmonics(theta, J=0) + alpha_quad_GS93*P_harmonics(theta, J=2) + directional
return alpha
def vv_perturb(tau, xx_SHU, vv_SHU77, alpha_SHU77, vv_ww_mono_GS93, vv_ww_quad_GS93, theta):
#Assuming J= 0, 2 only
directional = xx_SHU*dv_dx(xx_SHU, vv_SHU77, alpha_SHU77)*delta2(theta, tau)
directional = 0.0 #
vv_mono = vv_ww_mono_GS93[0]
vv_quad = vv_ww_quad_GS93[0]
ww_mono = vv_ww_mono_GS93[1]
ww_quad = vv_ww_quad_GS93[1]
#print('vv_mono, vv_quad, ww_mono, ww_quad', vv_mono, vv_quad, ww_mono, ww_quad)
vv_r = vv_mono*P_harmonics(theta, J=0) + vv_quad*P_harmonics(theta, J=2) + directional ## vv
vv_theta = ww_mono*B_harmonics(theta, J=0) + ww_quad*B_harmonics(theta, J=2) + directional ## ww
#print("vv_r, vv_theta", vv_r, vv_theta)
vv = np.array([vv_r, vv_theta])
return vv
def psi_perturb(tau, xx_SHU, mm_SHU77, mm_pp_mono_GS93, mm_pp_quad_GS93, theta):
#Assuming J= 0, 2 only
directional = xx_SHU*dpsi_dx(xx_SHU, mm_SHU77)*delta2(theta, tau)
directional = 0.0 #
mm_mono = mm_pp_mono_GS93[0]
mm_quad = mm_pp_quad_GS93[0]
pp_mono = mm_pp_mono_GS93[1]
pp_quad = mm_pp_quad_GS93[1]
#print('mm_pp_mono_GS93', mm_pp_mono_GS93)
#print('mm_mono', mm_mono)
psi = psi2(xx_SHU, mm_mono, pp_mono, J=0)*P_harmonics(theta, J=0) \
+ psi2(xx_SHU, mm_quad, pp_quad, J=0)*P_harmonics(theta, J=2) \
+ directional
#print('psi_perturb', psi)
return psi
def phi_vecpot_second_order(tau, xx_SHU, mm_SHU77, alpha_SHU77, FF_DD_mono_GS93, FF_DD_quad_GS93, theta):
directional = xx_SHU*dphi_dx(xx_SHU, alpha_SHU77, mm_SHU77, theta)*delta2(theta, tau)
directional = 0.0 #
sin_theta = np.sin(theta)
#print(FF_DD_mono_GS93)
#print(FF_DD_quad_GS93)
#print(ksii, P_harmonics(theta, J=0), P_harmonics(theta, J=2))
mono_term = (FF_DD_mono_GS93[0] + (1.0/ksii)*FF_DD_mono_GS93[1])
quad_term = (FF_DD_quad_GS93[0] + (1.0/ksii)*FF_DD_quad_GS93[1])
phi_vecpot_second = (sin_theta**2.0)*( mono_term*P_harmonics(theta, J=0) \
+ quad_term*P_harmonics(theta, J=2) ) \
+ directional
return phi_vecpot_second
def phi_vecpot_zero_order(xx_SHU, mm_SHU77, theta):
ff_zero = 0.25*(mm_SHU77**2.0)
sin_theta = np.sin(theta)
phi_vecpot_zero = ff_zero*(sin_theta*2.0)
return phi_vecpot_zero
# Combining the perturbation stage.
def alpha_xvec_tau(tau, xx_SHU, vv_SHU77, alpha_SHU77, alpha_mono_GS93, alpha_quad_GS93, theta):
alpha = alpha_SHU77 + (tau**2.0)*alpha_perturb(tau, xx_SHU, vv_SHU77, alpha_SHU77, alpha_mono_GS93, alpha_quad_GS93, theta)
return alpha
def vv_xvec_tau(tau, xx_SHU, vv_SHU77, alpha_SHU77, vv_ww_mono_GS93, vv_ww_quad_GS93, theta):
vv = (tau**2.0)*vv_perturb(tau, xx_SHU, vv_SHU77, alpha_SHU77, vv_ww_mono_GS93, vv_ww_quad_GS93, theta)
#print("BF",vv, vv_ww_mono_GS93, vv_ww_quad_GS93)
vv[0] = vv_SHU77 + vv[0]
vv[1] = 0.0 + vv[1] #No poloidal velocity in Shu77
#print("AF",vv)
return vv
def psi_xvec_tau(tau, xx_SHU, mm_SHU77, mm_pp_mono, mm_pp_quad, theta):
#print("psi_xvec_tau --- tau, xx_SHU, mm_SHU7, mm_pp_mono, mm_pp_quad, theta", tau, xx_SHU, mm_SHU77, mm_pp_mono, mm_pp_quad, theta)
psi = (tau**2.0)*psi_perturb(tau, xx_SHU, mm_SHU77, mm_pp_mono, mm_pp_quad, theta)
psi77 = get_SHU77_potential(xx_SHU)
#print('psi77', psi77)
psi = psi77 + psi
#print('psi_xvec_tau', psi)
return psi
def phi_vecpot_xvec_tau(tau, xx_SHU, mm_SHU77, alpha_SHU77, FF_DD_mono_GS93, FF_DD_quad_GS93, theta):
phi_vecpot_second = (tau**2.0)*phi_vecpot_second_order(tau, xx_SHU, mm_SHU77, alpha_SHU77, FF_DD_mono_GS93, FF_DD_quad_GS93, theta)
phi_vecpot_zero = phi_vecpot_zero_order(xx_SHU, mm_SHU77, theta)
phi_vecpot = phi_vecpot_zero + phi_vecpot_second
return phi_vecpot
#Physical unit converion stage
def rho_rt(tt, xx_SHU, vv_SHU77, alpha_SHU77, alpha_mono_GS93, alpha_quad_GS93, theta):
tau = get_tau(tt)
alpha_xvec = alpha_xvec_tau(tau, xx_SHU, vv_SHU77, alpha_SHU77, alpha_mono_GS93, alpha_quad_GS93, theta)
rho = (1.0/(4.0*np.pi*G_newton*(tt**2.0))) * alpha_xvec
return rho, alpha_xvec
def uu_rt(tt, xx_SHU, vv_SHU77, alpha_SHU77, vv_ww_mono_GS93, vv_ww_quad_GS93, theta):
tau = get_tau(tt)
vv_xvec = vv_xvec_tau(tau, xx_SHU, vv_SHU77, alpha_SHU77, vv_ww_mono_GS93, vv_ww_quad_GS93, theta)
uu = cs0*vv_xvec
return uu, vv_xvec
def grav_psi_rt(tt, xx_SHU, mm_SHU77, mm_pp_mono, mm_pp_quad, theta):
tau = get_tau(tt)
#print("tt , xx_SHU, mm_SHU77, mm_pp_mono, mm_pp_quad, theta", tt, xx_SHU, mm_SHU77, mm_pp_mono, mm_pp_quad, theta)
psi_xvec = psi_xvec_tau(tau, xx_SHU, mm_SHU77, mm_pp_mono, mm_pp_quad, theta)
Vpot = (cs0**2.0)*psi_xvec
return Vpot, psi_xvec
def vectorpot_rt(tt, xx_SHU, mm_SHU77, alpha_SHU77, FF_DD_mono_GS93, FF_DD_quad_GS93, theta):
tau = get_tau(tt)
phi_vecpot_xvec = phi_vecpot_xvec_tau(tau, xx_SHU, mm_SHU77, alpha_SHU77, FF_DD_mono_GS93, FF_DD_quad_GS93, theta)
Phi_flux = np.pi*B0*((cs0*tt)**2.0)*phi_vecpot_xvec
return Phi_flux, phi_vecpot_xvec
###def match_xx(xx_rad, xx_SHU):
### xx_buffer = np.empty_like(xx_rad)
### stride = np.abs(xx_SHU[1] - xx_SHU[0])
### for xx in xx_SHU:
### #where xx - stride < xx_rad < xx + stride -> xx_rad[i] = xx
### #loc = np.where((xx_rad <= (xx + stride) and xx_rad > (xx - stride) ))
### loc = np.where(xx_rad <= (xx + stride) )
### print(loc)
def get_shu_index(xx, xx_SHU):
stride = np.abs(xx_SHU[1] - xx_SHU[0])/2.0
#ishu = np.where((xx_SHU <= (xx + stride)) & (xx_SHU > (xx - stride)))[0]
#TODO Now a purkka version. Do better.
# Can be improve by taking the treatment of the actual low and high x cases.
if (xx > xx_SHU[xx_SHU.size-1]):
ishu = xx_SHU.size-1
elif (xx < xx_SHU[0]):
ishu = 0
else:
ishu = np.where((xx_SHU <= (xx + stride)) & (xx_SHU > (xx - stride)))[0]
#print("get_shu_index", ishu, ishu.size)
ishu = ishu[0]
#print("get_shu_index", ishu, ishu.size)
#print(ishu, xx_SHU[ishu], xx)
return ishu
def plot_figure(tt, xx_horizontal_corners, xx_vertical_corners, xx_horizontal, xx_vertical, xxvar, physvar,
vv_hor=np.array(None), vv_ver=np.array(None), uu_hor=np.array(None), uu_ver=np.array(None),
title1=r"\alpha", title2=r"\rho", filetitle='density',
var_min=[None, None], var_max=[None, None], colmap=CM_INFERNO, normtype='log',
streamlines = 0, contourplot = 0):
if var_min[0] != None:
if normtype == 'log':
mynorm1 = colors.LogNorm( vmin=var_min[0], vmax=var_max[0] )
mynorm2 = colors.LogNorm( vmin=var_min[1], vmax=var_max[1] )
else:
mynorm1 = colors.Normalize( vmin=var_min[0], vmax=var_max[0] )
mynorm2 = colors.Normalize( vmin=var_min[1], vmax=var_max[1] )
else:
mynorm1 = colors.Normalize( )
mynorm2 = colors.Normalize( )
if contourplot:
if normtype =='cdensity':
numbers = np.arange(0, 20, dtype=np.float64)
contourlevs = 1e-20*(np.sqrt(2.0)**numbers)
contournorm = colors.LogNorm( vmin=contourlevs.min(), vmax=contourlevs.max() )
elif normtype =='cflux':
contourlevs = np.linspace(1.0, 1e31, num=20)
contournorm = colors.Normalize( vmin=contourlevs.min(), vmax=contourlevs.max() )
else:
contourlevs = np.linspace(physvar.min(), physvar.max(), num=10)
contournorm = colors.Normalize( vmin=contourlevs.min(), vmax=contourlevs.max() )
##rr_horizontal_corners = xx_horizontal_corners*(cs0*tt)/AU
##rr_vertical_corners = xx_vertical_corners* (cs0*tt)/AU
##rr_horizontal = xx_horizontal*(cs0*tt)/AU
##rr_vertical = xx_vertical* (cs0*tt)/AU
rr_horizontal_corners = xx_horizontal_corners*(cs0*tt)/1e17
rr_vertical_corners = xx_vertical_corners* (cs0*tt)/1e17
rr_horizontal = xx_horizontal*(cs0*tt)/1e17
rr_vertical = xx_vertical* (cs0*tt)/1e17
figa, axa = plt.subplots(nrows=1, ncols=2, figsize=(16,6))
if contourplot:
mapa = axa[0].contourf(xx_horizontal, xx_vertical, xxvar, cmap=colmap, norm=mynorm1)
maprho = axa[1].contourf(rr_horizontal, rr_vertical, physvar, contourlevs, cmap=colmap, norm=contournorm)
else:
mapa = axa[0].pcolormesh(xx_horizontal_corners, xx_vertical_corners, xxvar, cmap=colmap, norm=mynorm1 )
maprho = axa[1].pcolormesh(rr_horizontal_corners, rr_vertical_corners, physvar, cmap=colmap, norm=mynorm2)
#mapa = axa[0].contourf(xx_horizontal, xx_vertical, alpha, cmap=CM_INFERNO, norm=colors.LogNorm(vmin=0.1, vmax=50.0))
#maprho = axa[1].contourf(xx_horizontal*(cs0*tt)/AU, xx_vertical*(cs0*tt)/AU, rho, cmap=CM_INFERNO, norm=colors.LogNorm(vmin=1e15, vmax=1e20))
if vv_hor.any() != None:
if streamlines:
#vv_tot = np.sqrt(vv_hor**2.0 + vv_ver**2.0)
#vv_tot = np.log(vv_tot/vv_tot.max())
axa[0].streamplot(xx_horizontal, xx_vertical, vv_hor, vv_ver, color = 'k')
axa[1].streamplot(rr_horizontal, rr_vertical, uu_hor, uu_ver, color = 'k' )
else:
axa[0].quiver(xx_horizontal, xx_vertical, vv_hor, vv_ver, pivot = 'middle')
axa[1].quiver(rr_horizontal, rr_vertical, uu_hor, uu_ver, pivot = 'middle')
fig.colorbar(mapa, ax=axa[0])
fig.colorbar(maprho, ax=axa[1])
tau = get_tau(tt)
tt_kyr = tt/kyr
axa[0].set_title(r'$%s(x, \tau = %.3f)$ ' % (title1, tau))
axa[1].set_title(r'$%s(r, t = %.3f \mathrm{kyr})$ ' % (title2, tt_kyr))
axa[0].set_xlabel('x')
axa[0].set_ylabel('x')
#axa[1].set_xlabel('r (AU)')
#axa[1].set_ylabel('r (AU)')
axa[1].set_xlabel(r'r ($10^{17}$ cm)')
axa[1].set_ylabel(r'r ($10^{17}$ cm)' )
##axa[1].set_xlim(0.0, 3e17/AU)
##axa[1].set_ylim(0.0, 3e17/AU)
axa[1].set_xlim(0.0, 3.0)
axa[1].set_ylim(0.0, 3.0)
axa[0].set_aspect('equal', 'datalim')
#axa[1].set_aspect('equal', 'datalim')
figfile = '%s_%s.png' % (filetitle, str(numslice).zfill(6))
print(figfile)
figa.savefig(figfile)
plt.close(figa)
xx_SHU = np.array([ 0.05, 0.10, 0.15, 0.20, 0.25,
0.30, 0.35, 0.40, 0.45, 0.50,
0.55, 0.60, 0.65, 0.70, 0.75,
0.80, 0.85, 0.90, 0.95, 1.00])
alpha_SHU77 = np.array([ 71.5, 27.8, 16.4, 11.5, 8.76,
7.09, 5.95, 5.14, 4.52, 4.04,
3.66, 3.35, 3.08, 2.86, 2.67,
2.50, 2.35, 2.22, 2.10, 2.00])
vv_SHU77 = -np.array([ 5.44, 3.47, 2.58, 2.05, 1.68,
1.40, 1.18, 1.01, 0.861, 0.735,
0.625, 0.528, 0.442, 0.363, 0.291,
0.225, 0.163, 0.106, 0.051, 0.00])
mm_SHU77 = np.array([0.981, 0.993, 1.01, 1.03, 1.05,
1.08, 1.12, 1.16, 1.20, 1.25,
1.30, 1.36, 1.42, 1.49, 1.56,
1.64, 1.72, 1.81, 1.90, 2.00])
#GS Table 1
alpha_mono_GS93 = np.array([ 6.304, 2.600, 1.652, 1.156, 9.005e-1,
7.314e-1, 6.084e-1, 5.084e-1, 4.256e-1, 3.517e-1,
2.829e-1, 2.172e-1, 1.488e-1, 8.091e-2, 8.360e-3,
-6.826e-2, -1.512e-1, -2.406e-1, -3.382e-1, -4.444e-1])
vv_ww_mono_GS93 = np.array([[4.372e-1, 3.335e-1, 2.390e-1, 1.918e-1, 1.522e-1,
1.226e-1, 9.579e-2, 7.103e-2, 4.828e-2, 2.640e-2,
5.058e-3, -1.588e-2, -3.791e-2, -5.975e-2, -8.293e-2,
-1.071e-1, -1.330e-1, -1.605e-1, -1.902e-1, -2.222e-1],
[ 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0]])
mm_pp_mono_GS93 = np.array([[8.634e-4, 1.959e-3, 3.560e-3, 5.661e-3, 8.235e-3,
1.130e-2, 1.482e-2, 1.873e-2, 2.293e-2, 2.730e-2,
3.166e-2, 3.579e-2, 3.935e-2, 4.196e-2, 4.312e-2,
4.221e-2, 3.847e-2, 3.097e-2, 1.859e-2, 0.0],
[ 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0]])
FF_DD_mono_GS93 = np.array([[ -1.130, -3.275e-1, -1.355e-1, -6.415e-2, -2.889e-2, #F
-8.387e-3, 5.358e-3, 1.534e-2, 2.303e-2, 2.931e-2,
3.454e-2, 3.888e-2, 4.225e-2, 4.442e-2, 4.504e-2,
4.358e-2, 3.935e-2, 3.146e-2, 1.881e-2, 0.0],
[ -1.246e1, -3.168, -1.141, -5.740e-1, -3.178e-1, #D
-1.878e-1, -1.049e-1, -4.547e-2, 3.393e-4, 3.924e-2,
7.431e-2, 1.070e-1, 1.376e-1, 1.650e-1, 1.867e-1,
1.992e-1, 1.966e-1, 1.708e-1, 1.103e-1, 0.0]])
#GS Table 2
alpha_quad_GS93 = np.array([ -1.096e3, -1.191e2, -3.148e1, -1.158e1, -5.105,
-2.456, -1.217, -5.889e-1, -2.569e-1, -7.024e-2,
3.790e-2, 1.042e-1, 1.505e-1, 1.845e-1, 2.163e-1,
2.492e-1, 2.865e-1, 3.302e-1, 3.823e-1, 4.437e-1])
vv_ww_quad_GS93 = np.array([[ -2.581, -1.533, -8.072e-1, -5.666e-1, -3.905e-1, #v
-2.790e-1, -1.928e-1, -1.254e-1, -7.156e-2, -2.614e-2,
1.267e-2, 4.650e-2, 7.724e-2, 1.042e-1, 1.288e-1,
1.510e-1, 1.711e-1, 1.889e-1, 2.045e-1, 2.177e-1],
[ -2.085, -4.890, -1.811, -8.842e-1, -4.816e-1, #w
-2.807e-1, -1.628e-1, -8.779e-2, -3.852e-2, -4.481e-3,
1.928e-2, 3.578e-2, 4.683e-2, 5.306e-2, 5.512e-2,
5.312e-2, 4.704e-2, 3.670e-2, 2.179e-2, 1.898e-3]])
mm_pp_quad_GS93 = np.array([[-3.860e-5, -1.541e-4, -3.044e-4, -4.847e-4, -6.831e-4, #m
-8.874e-4, -1.083e-3, -1.253e-3, -1.385e-3, -1.462e-3,
-1.470e-3, -1.389e-3, -1.191e-3, -8.405e-4, -2.841e-4,
5.579e-4, 1.800e-3, 3.609e-3, 6.218e-3, 9.951e-3],
[ -7.539e1, -7.275, -1.730, -5.586e-1, -1.999e-1, #p
-6.591e-1, -1.062e-2, 1.294e-2, 2.267e-2, 2.600e-2,
2.625e-2, 2.500e-2, 2.294e-2, 2.046e-2, 1.769e-2,
1.469e-2, 1.146e-2, 7.941e-3, 4.102e-3, -1.214e-4]])
FF_DD_quad_GS93 = np.array([[ -2.253, -6.517e-1, -2.722e-1, -1.345e-1, -6.993e-2, #F
-3.593e-2, -1.660e-2, -5.864e-3, -6.809e-4, 8.213e-4,
-3.086e-4, -3.338e-3, -7.681e-3, -1.272e-2, -1.778e-2,
-2.191e-2, -2.392e-2, -2.219e-2, -1.457e-2, 1.729e-3],
[ -2.484e1, -6.258, -2.221, -1.102, -6.127e-1, #D
-3.645e-1, -2.213e-1, -1.297e-1, -7.020e-2, -1.112e-2,
-2.139e-3, -1.615e-2, 2.744e-2, 3.252e-2, 3.269e-2,
2.839e-2, 2.104e-2, 1.199e-2, 3.732e-3, 0.0]])
tt = 0.3*ttm
theta = 0.5*np.pi
xx_SHU = xx_SHU[:-1]
vv_SHU77 = vv_SHU77[:-1]
alpha_SHU77 = alpha_SHU77[:-1]
alpha_mono_GS93 = alpha_mono_GS93[:-1]
alpha_quad_GS93 = alpha_quad_GS93[:-1]
vv_ww_mono_GS93 = np.array([vv_ww_mono_GS93[0][:-1], vv_ww_mono_GS93[1][:-1]])
vv_ww_quad_GS93 = np.array([vv_ww_quad_GS93[0][:-1], vv_ww_quad_GS93[1][:-1]])
rho, alpha_xvec = rho_rt(tt, xx_SHU, vv_SHU77, alpha_SHU77, alpha_mono_GS93, alpha_quad_GS93, theta)
rr = xx_SHU*cs0*tt
np.set_printoptions(linewidth=200)
print(rho.shape)
print(xx_SHU.shape)
print(rho)
print(xx_SHU)
print(vv_ww_mono_GS93)
print(vv_ww_quad_GS93)
print(vv_ww_quad_GS93[0])
print(vv_ww_quad_GS93[1])
#plt.figure()
#plt.plot(rr, rho)
#
#plt.figure()
#plt.plot(xx_SHU, alpha_xvec, label = "GS93")
#plt.plot(xx_SHU, alpha_SHU77, label = "Shu77")
#plt.legend()
#alpha_mono_yy, alpha_quad_yy, alpha_mono_yy = yy_transform(xx_SHU, alpha_SHU77, alpha_mono_GS93, alpha_quad_GS93)
plt.figure()
plt.plot(xx_SHU, alpha_SHU77, label=r"$\alpha^{(0)}$")
plt.plot(xx_SHU, alpha_mono_GS93, label=r"$\alpha^{(2)}_0$")
plt.plot(xx_SHU, alpha_quad_GS93, label=r"$\alpha^{(2)}_2$")
plt.ylim([-5.0,5.0])
plt.legend()
plt.show()
'''
ii = 0
theta_axis = np.linspace(0.0, np.pi)
xx_theta = np.array([])
print("PIIP")
plt.figure()
for ii in range(0,xx_SHU.size):
alpha_theta = np.array([])
alpha_shuref = np.array([])
for theta in theta_axis:
rho, alpha_xvec = rho_rt(tt, xx_SHU[ii], vv_SHU77[ii], alpha_SHU77[ii], alpha_mono_GS93[ii], alpha_quad_GS93[ii])
alpha_theta = np.append(alpha_theta, alpha_xvec)
alpha_shuref = np.append(alpha_shuref, alpha_SHU77[ii])
plt.plot(alpha_theta, theta_axis, label = "GS93")
#plt.plot(alpha_shuref, theta_axis, label = "GS93")
'''
#Interpolate a mesh.
xx_SHU_GRID = np.insert(xx_SHU, 0, 0.0)
print(xx_SHU_GRID)
xx_horizontal, xx_vertical = np.meshgrid(xx_SHU_GRID, xx_SHU_GRID, indexing='xy')
theta = np.arctan2(xx_horizontal, xx_vertical)
#Take pcolormesh coordinate system into account, which marks corners instead of centre points.
dxx = np.abs(xx_horizontal[0,1] - xx_horizontal[0,0])
print(dxx)
xx_horizontal_corners = xx_horizontal - dxx/2.0
xx_vertical_corners = xx_vertical - dxx/2.0
xx_rad = np.sqrt(xx_horizontal**2.0 + xx_vertical**2.0)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(16,4))
map1 = ax[0].pcolormesh(xx_horizontal_corners, xx_vertical_corners, theta)
map2 = ax[1].pcolormesh(xx_horizontal_corners, xx_vertical_corners, xx_rad)
ax[0].set_title(r"$\theta$")
ax[1].set_title(r"$x_\mathrm{rad}$")
fig.colorbar(map1, ax=ax[0])
fig.colorbar(map2, ax=ax[1])
ax[0].set_aspect('equal', 'datalim')
ax[1].set_aspect('equal', 'datalim')
Pfig, Pax = plt.subplots(nrows=1, ncols=3, figsize=(16,4))
print("P_harmonics(theta, J=0)", P_harmonics(theta, J=0))
Pmap1 = Pax[0].pcolormesh(xx_horizontal_corners, xx_vertical_corners, P_harmonics(theta, J=0))
Pmap2 = Pax[1].pcolormesh(xx_horizontal_corners, xx_vertical_corners, P_harmonics(theta, J=2))
Pmap3 = Pax[2].pcolormesh(xx_horizontal_corners, xx_vertical_corners, deltaspace(theta, 0.5))
Pax[0].set_title(r"$P_0(\theta)$")
Pax[1].set_title(r"$P_2(\theta)$")
Pax[2].set_title(r"$\Delta(\theta, \tau = 0.5)$")
Pfig.colorbar(Pmap1, ax=Pax[0])
Pfig.colorbar(Pmap2, ax=Pax[1])
Pfig.colorbar(Pmap3, ax=Pax[2])
Pax[0].set_aspect('equal', 'datalim')
Pax[1].set_aspect('equal', 'datalim')
Pax[2].set_aspect('equal', 'datalim')
Bfig, Bax = plt.subplots(nrows=1, ncols=2, figsize=(16,4))
print("B_harmonics(theta, J=0)", B_harmonics(theta, J=0))
Bmap1 = Bax[0].pcolormesh(xx_horizontal_corners, xx_vertical_corners, B_harmonics(theta, J=0))
Bmap2 = Bax[1].pcolormesh(xx_horizontal_corners, xx_vertical_corners, B_harmonics(theta, J=2))
Bax[0].set_title(r"$B_0(\theta)$")
Bax[1].set_title(r"$B_2(\theta)$")
Bfig.colorbar(Bmap1, ax=Bax[0])
Bfig.colorbar(Bmap2, ax=Bax[1])
Bax[0].set_aspect('equal', 'datalim')
Bax[1].set_aspect('equal', 'datalim')
plt.show()
##xx_horizontal_corners = np.append(xx_horizontal_corners, (np.amax(xx_horizontal_corners)+dxx)*np.ones((xx_horizontal_corners.shape[1],1)), axis=1)
print(xx_horizontal_corners[-1,:])
print(xx_horizontal_corners)
##xx_horizontal_corners = np.vstack((xx_horizontal_corners, xx_horizontal_corners[-1,:]))
##print(xx_horizontal_corners)
##xx_vertical_corners = np.append(xx_vertical_corners, (np.amax(xx_vertical_corners)+dxx)*np.ones((1,xx_vertical_corners.shape[0])), axis=0)
print(xx_vertical_corners[:, -1])
print(xx_vertical_corners)
##xx_vertical_corners = np.hstack((xx_vertical_corners, xx_vertical_corners[:,-1]))
print(xx_vertical_corners)
numslice = 0
frametot = 201
#frametot = 101
#frametot = 11
for tt in np.linspace(0.1, ttm, num=frametot):
alpha = np.empty_like(xx_rad)
alpha77 = np.empty_like(xx_rad)
rho = np.empty_like(xx_rad)
vv_rad = np.empty_like(xx_rad)
vv_pol = np.empty_like(xx_rad)
uu_rad = np.empty_like(xx_rad)
uu_pol = np.empty_like(xx_rad)
psi = np.empty_like(xx_rad)
Vpot = np.empty_like(xx_rad)
Delta = np.empty_like(xx_rad)
Phi_flux = np.empty_like(xx_rad)
phi_vecpot = np.empty_like(xx_rad)
alpha_2_J = np.empty_like(xx_rad)
for ii in range(xx_SHU_GRID.size):
for kk in range(xx_SHU_GRID.size):
xx = xx_rad[ii,kk]
th = theta[ii,kk]
ishu = get_shu_index(xx, xx_SHU)
rho[ii, kk], alpha[ii, kk] = rho_rt(tt, xx_SHU[ishu],
vv_SHU77[ishu],
alpha_SHU77[ishu],
alpha_mono_GS93[ishu],
alpha_quad_GS93[ishu], th)
alpha77[ii, kk] = alpha_SHU77[ishu]
vv_ww_mono_point = vv_ww_mono_GS93[:, ishu]
vv_ww_quad_point = vv_ww_quad_GS93[:, ishu]
uu_dump, vv_dump = uu_rt(tt, xx_SHU[ishu], vv_SHU77[ishu], alpha_SHU77[ishu], vv_ww_mono_point, vv_ww_quad_point, th)
vv_rad[ii, kk] = vv_dump[0]
vv_pol[ii, kk] = vv_dump[1]
uu_rad[ii, kk] = uu_dump[0]
uu_pol[ii, kk] = uu_dump[1]
mm_pp_mono_point = mm_pp_mono_GS93[:, ishu]
mm_pp_quad_point = mm_pp_quad_GS93[:, ishu]
Vpot[ii, kk], psi[ii, kk] = grav_psi_rt(tt, xx_SHU[ishu], mm_SHU77[ishu], mm_pp_mono_point, mm_pp_quad_point, th)
Phi_flux[ii, kk], phi_vecpot[ii, kk] = vectorpot_rt(tt, xx_SHU[ishu], mm_SHU77[ishu], alpha_SHU77[ishu],
FF_DD_mono_GS93[:, ishu],
FF_DD_quad_GS93[:, ishu], th)
Delta[ii, kk] = deltaspace(th, get_tau(tt))
alpha_2_J[ii, kk] = alpha_mono_GS93[ishu]*P_harmonics(th, J=0) + alpha_quad_GS93[ishu]*P_harmonics(th, J=2)
vv_hor = vv_pol*np.cos(theta) + vv_rad*np.sin(theta)
vv_ver = - vv_pol*np.sin(theta) + vv_rad*np.cos(theta)
uu_hor = uu_pol*np.cos(theta) + uu_rad*np.sin(theta)
uu_ver = - uu_pol*np.sin(theta) + uu_rad*np.cos(theta)
rho77 = alpha77 * (1.0/(4.0*np.pi*G_newton)*tt) #TODO WRONG COEFFS!!!
#Apply mask
rad_mask = 0.2
alpha = np.ma.masked_where(xx_rad < rad_mask, alpha)
rho = np.ma.masked_where(xx_rad < rad_mask, rho)
vv_rad = np.ma.masked_where(xx_rad < rad_mask, vv_rad)
uu_rad = np.ma.masked_where(xx_rad < rad_mask, uu_rad)
vv_pol = np.ma.masked_where(xx_rad < rad_mask, vv_pol)
uu_pol = np.ma.masked_where(xx_rad < rad_mask, uu_pol)
vv_hor = np.ma.masked_where(xx_rad < rad_mask, vv_hor)
vv_ver = np.ma.masked_where(xx_rad < rad_mask, vv_ver)
uu_hor = np.ma.masked_where(xx_rad < rad_mask, uu_hor)
uu_ver = np.ma.masked_where(xx_rad < rad_mask, uu_ver)
psi = np.ma.masked_where(xx_rad < rad_mask, psi )
Vpot = np.ma.masked_where(xx_rad < rad_mask, Vpot)
phi_vecpot = np.ma.masked_where(xx_rad < rad_mask, phi_vecpot)
Phi_flux = np.ma.masked_where(xx_rad < rad_mask, Phi_flux )
alpha_2_J = np.ma.masked_where(xx_rad < rad_mask, alpha_2_J)
Delta = np.ma.masked_where(xx_rad < rad_mask, Delta )
plot_figure(tt, xx_horizontal_corners, xx_vertical_corners, xx_horizontal, xx_vertical, alpha, rho,
vv_hor=vv_hor, vv_ver=vv_ver, uu_hor=uu_hor, uu_ver=uu_ver,
title1=r"\alpha", title2=r"\rho", filetitle='GS93density',
streamlines = 1, contourplot=1,
var_min=[0.00, 1e15], var_max=[16, 1e21],
normtype = 'cdensity')
plot_figure(tt, xx_horizontal_corners, xx_vertical_corners, xx_horizontal, xx_vertical, alpha77, rho77,
#var_min=[0.00, 0], var_max=[16, 1e20],
title1=r"\alpha", title2=r"\rho", filetitle='S77density')
plot_figure(tt, xx_horizontal_corners, xx_vertical_corners, xx_horizontal, xx_vertical, vv_rad, uu_rad,
vv_hor=vv_hor, vv_ver=vv_ver, uu_hor=uu_hor, uu_ver=uu_ver,
title1=r"v_r", title2=r"u_r", filetitle='GS93velocity_rad',
var_min=[-2.5, -2.5*cs0], var_max=[0.0, 0.0*cs0],
normtype = 'lin')
plot_figure(tt, xx_horizontal_corners, xx_vertical_corners, xx_horizontal, xx_vertical, vv_pol, uu_pol,
vv_hor=vv_hor, vv_ver=vv_ver, uu_hor=uu_hor, uu_ver=uu_ver,
title1=r"v_\theta", title2=r"u_\theta", filetitle='GS93velocity_pol',
var_min=[0.0, 0.0*cs0], var_max=[0.5, 0.5*cs0],
normtype = 'lin')
plot_figure(tt, xx_horizontal_corners, xx_vertical_corners, xx_horizontal, xx_vertical, psi, Vpot,
vv_hor=vv_hor, vv_ver=vv_ver, uu_hor=uu_hor, uu_ver=uu_ver,
title1=r"\psi", title2=r"V_\mathrm{pot}", filetitle='GS93gravpot',
var_min=[12.0, 12.0*(cs0**2.0)], var_max=[21.0, 21.0*(cs0**2.0)],
normtype = 'lin')
plot_figure(tt, xx_horizontal_corners, xx_vertical_corners, xx_horizontal, xx_vertical, phi_vecpot, Phi_flux,
title1=r"\phi", title2=r"\Phi_\mathrm{flux}", filetitle='GS93vecpot',
vv_hor=vv_hor, vv_ver=vv_ver, uu_hor=uu_hor, uu_ver=uu_ver,
streamlines = 1, contourplot=1,
normtype = 'cflux')
plot_figure(tt, xx_horizontal_corners, xx_vertical_corners, xx_horizontal, xx_vertical, np.sqrt(vv_hor**2.0 + vv_ver**2.0), np.sqrt(uu_hor**2.0 + uu_ver**2.0),
title1=r"|v|", title2=r"|u| (cm/s)", filetitle='GS93vel2',
var_min=[0.0, 0.0*cs0], var_max=[2.5, 2.5*cs0],
vv_hor=vv_hor, vv_ver=vv_ver, uu_hor=uu_hor, uu_ver=uu_ver,
streamlines = 1,
normtype = 'lin')
##plot_figure(tt, xx_horizontal_corners, xx_vertical_corners, xx_horizontal, xx_vertical, Delta, Delta,
## title1=r"\Delta", title2=r"\Delta", filetitle='Delta',
## normtype = 'lin')
##plot_figure(tt, xx_horizontal_corners, xx_vertical_corners, xx_horizontal, xx_vertical, alpha_2_J, alpha_2_J,
## title1=r"\sum \alpha^{(2)}_J", title2=r"\sum \alpha^{(2)}_J", filetitle='alpha_2_J',
## normtype = 'lin')
numslice += 1

1
analysis/python/calc/purge.sh
View File

@@ -1 +0,0 @@
rm *.png

279
analysis/python/calc/shu_selfsim.py
View File

@@ -1,279 +0,0 @@
'''
Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
This file is part of Astaroth.
Astaroth is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
Astaroth is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with Astaroth. If not, see <http://www.gnu.org/licenses/>.
'''
import numpy as np
import pylab as plt
G_newton = 6.674e-8 #cm**3 g**-1 s**-2
def dv_dx(xx,vv, alpha):
EE = alpha*(xx-vv) - 2.0/xx
HH = (xx-vv)**2.0 - 1.0
return (EE/HH)*(xx-vv)
def dalpha_dx(xx,vv, alpha):
EE = alpha*(alpha - (2.0/xx)*(xx-vv))
HH = (xx-vv)**2.0 - 1.0
return (EE/HH)*(xx-vv)
###def dv_dx(xx,vv, alpha):
### return 2.0*(xx-vv)
###
###def dalpha_dx(xx,vv, alpha):
### return -1.0*(xx-vv)
def get_m(xx, vv, alpha):
mm = xx**2.0 * alpha * (xx-vv)
return mm
def alpha_to_rho(alpha, tt):
rho = alpha/(4.0*np.pi*G_newton*(tt**2.0))
return rho
def vv_to_uu(vv, cs0):
uu = cs0*vv
return uu
def mm_to_MM(mm, tt, cs0):
MM = (((cs0**3.0)*tt)/G_newton)*mm
return MM
def euler(xx_step, xx, vv, alpha, mm, target):
diff = target - xx[-1]
if diff >= 0:
while xx[-1] <= target:
vv_step = vv[-1] + xx_step*dv_dx(xx[-1], vv[-1], alpha[-1])
alpha_step = alpha[-1] + xx_step*dalpha_dx(xx[-1], vv[-1], alpha[-1])
xx = np.append(xx, xx[-1]+xx_step)
alpha = np.append(alpha, alpha_step)
vv = np.append(vv, vv_step)
mm_step = get_m(xx[-1], vv[-1], alpha[-1])
mm = np.append(mm, mm_step)
else:
while xx[-1] <= target:
vv_step = vv[-1] + xx_step*dv_dx(xx[-1], vv[-1], alpha[-1])
alpha_step = alpha[-1] + xx_step*dalpha_dx(xx[-1], vv[-1], alpha[-1])
xx = np.append(xx, xx[-1]+xx_step)
alpha = np.append(alpha, alpha_step)
vv = np.append(vv, vv_step)
mm_step = get_m(xx[-1], vv[-1], alpha[-1])
mm = np.append(mm, mm_step)
return xx, vv, alpha, mm
def RK4_step(vv, xx, alpha, xx_step):
vv1 = xx_step*dv_dx(xx[-1], vv[-1], alpha[-1])
alpha1 = xx_step*dalpha_dx(xx[-1], vv[-1], alpha[-1])
vv2 = xx_step*dv_dx(xx[-1]+xx_step/2.0, vv[-1]+vv1/2.0, alpha[-1]+alpha1/2.0)
alpha2 = xx_step*dalpha_dx(xx[-1]+xx_step/2.0, vv[-1]+vv1/2.0, alpha[-1]+alpha1/2.0)
vv3 = xx_step*dv_dx(xx[-1]+xx_step/2.0, vv[-1]+vv2/2.0, alpha[-1]+alpha2/2.0)
alpha3 = xx_step*dalpha_dx(xx[-1]+xx_step/2.0, vv[-1]+vv2/2.0, alpha[-1]+alpha2/2.0)
vv4 = xx_step*dv_dx(xx[-1]+xx_step, vv[-1]+vv3, alpha[-1]+alpha3)
alpha4 = xx_step*dalpha_dx(xx[-1]+xx_step, vv[-1]+vv3, alpha[-1]+alpha3)
vv_step = vv[-1] + (1.0/6.0)*(vv1 + 2.0*vv2 + 2.0*vv3 + vv4)
alpha_step = alpha[-1] + (1.0/6.0)*(alpha1 + 2.0*alpha2 + 2.0*alpha3 + alpha4)
return vv_step, alpha_step
def RK4(xx_step, xx, vv, alpha, mm, target, epsilon):
#Runge-Kutta RK4
diff = target - xx[-1]
#if diff < 0:
if diff >= 0:
while xx[-1] <= target:
if (np.abs(xx[-1] - vv[-1] - 1.0) > epsilon):
vv_step, alpha_step = RK4_step(vv, xx, alpha, xx_step)
print( vv_step, alpha_step)
else:
vv_step = vv[-1]
alpha_step = alpha[-1]
print("PIIP")
#print(np.abs(xx[-1] - vv[-1]), epsilon)
xx = np.append(xx, xx[-1]+xx_step)
alpha = np.append(alpha, alpha_step)
vv = np.append(vv, vv_step)
mm_step = get_m(xx[-1], vv[-1], alpha[-1])
mm = np.append(mm, mm_step)
else:
while xx[-1] >= target:
if (np.abs(xx[-1] - vv[-1] - 1.0) > epsilon):
vv_step, alpha_step = RK4_step(vv, xx, alpha, xx_step)
print( vv_step, alpha_step)
else:
vv_step = vv[-1]
alpha_step = alpha[-1]
print("PIIP")
#print(np.abs(xx[-1] - vv[-1]), epsilon)
xx = np.append(xx, xx[-1]+xx_step)
alpha = np.append(alpha, alpha_step)
vv = np.append(vv, vv_step)
mm_step = get_m(xx[-1], vv[-1], alpha[-1])
mm = np.append(mm, mm_step)
return xx, vv, alpha, mm
# From Shu 1977 TABLE II
xx_SHU = np.array([0.05 , 0.10 , 0.15 , 0.20 , 0.25 , 0.30 , 0.35 , 0.40 , 0.45 ,
0.50 , 0.55 , 0.60 , 0.65 , 0.70 , 0.75 , 0.80 , 0.85 ,
0.90 , 0.95 , 1.00])
alpha_SHU = np.array([71.5 , 27.8 , 16.4 , 11.5 , 8.76 , 7.09 , 5.95 , 5.14 , 4.52 ,
4.04 , 3.66 , 3.35 , 3.08 , 2.86 , 2.67 , 2.50 , 2.35 ,
2.22 , 2.10 , 2.00])
vv_SHU = -np.array([5.44 , 3.47 , 2.58 , 2.05 , 1.68 , 1.40 , 1.18 , 1.01 , 0.861,
0.735, 0.625, 0.528, 0.442, 0.363, 0.291, 0.225, 0.163,
0.106, 0.051, 0.00])
mm_SHU = np.array([0.981, 0.993, 1.01 , 1.03 , 1.05 , 1.08 , 1.12 , 1.16 , 1.20 ,
1.25 , 1.30 , 1.36 , 1.42 , 1.49 , 1.56 , 1.64 , 1.72 ,
1.81 , 1.90 , 2.00])
##From Shu (1977)
#AA = [ 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0]
#m0 = [0.975, 1.45, 1.88, 2.31, 2.74, 3.18, 3.63, 4.10, 4.58, 5.08, 5.58]
#AA = np.array(AA)
#m0 = np.array(m0)
#xx0 = xx_SHU[1]
#alpha0 = alpha_SHU[1]
#vv0 = vv_SHU[1]
#xx_step = 0.005
#target = 1.0
xx0 = xx_SHU[-3]
alpha0 = alpha_SHU[-3]
vv0 = vv_SHU[-3]
target = 0.05
xx_step = -0.005
xx_step = -0.001
print(get_m(xx0, alpha0, vv0))
xx = np.array([])
alpha = np.array([])
vv = np.array([])
mm = np.array([])
xx = np.append(xx, xx0)
alpha = np.append(alpha, alpha0)
vv = np.append(vv, vv0)
mm = np.append(mm, get_m(xx0, alpha0, vv0))
print(xx, alpha, vv, mm)
xx_EUL, vv_EUL, alpha_EUL, mm_EUL = euler(xx_step, xx, vv, alpha, mm, target)
xx_RK , vv_RK , alpha_RK , mm_RK = RK4(xx_step, xx, vv, alpha, mm, target, epsilon = 0.000001)
mm_EUL = get_m(xx_EUL, alpha_EUL, vv_EUL)
mm_RK = get_m(xx_RK , alpha_RK , vv_RK )
mm_SHU = get_m(xx_SHU, alpha_SHU, vv_SHU)
# Plotting time
figQ, axQ = plt.subplots(nrows=2, ncols=2, sharex=True)
axQ[0,0].plot(xx_EUL, alpha_EUL, label=r'$\alpha$ (Euler)', linewidth = 3.0)
axQ[0,0].plot(xx_RK , alpha_RK , label=r'$\alpha$ (RK4)', linewidth = 3.0)
axQ[0,0].plot(xx_SHU, alpha_SHU, 'd', label=r'$\alpha$ (Shu)', linewidth = 3.0)
axQ[0,0].set_xlabel(r'x')
axQ[0,0].set_ylabel(r'$\alpha$')
axQ[0,0].legend()
axQ[0,1].plot(xx_EUL, np.abs(vv_EUL), label='v (Euler)', linewidth = 3.0)
axQ[0,1].plot(xx_RK , np.abs(vv_RK ), label='v (RK4)', linewidth = 3.0)
axQ[0,1].plot(xx_SHU, np.abs(vv_SHU),'d', label='v (Shu)', linewidth = 3.0)
axQ[0,1].set_xlabel(r'x')
axQ[0,1].set_ylabel(r'-v')
axQ[0,1].legend()
axQ[1,0].plot(xx_EUL, mm_EUL, label='m (Euler)', linewidth = 3.0)
axQ[1,0].plot(xx_RK , mm_RK , label='m (RK4)', linewidth = 3.0)
axQ[1,0].plot(xx_SHU , mm_SHU , 'd', label='m (Shu)', linewidth = 3.0)
axQ[1,0].set_xlabel(r'x')
axQ[1,0].set_ylabel(r'm')
axQ[1,0].legend()
axQ[1,1].plot(xx_EUL, xx_EUL-vv_EUL, label='x-v (Euler)', linewidth = 3.0)
axQ[1,1].plot(xx_RK , xx_RK -vv_RK , label='x-v (RK4)', linewidth = 3.0)
axQ[1,1].plot(xx_SHU, xx_SHU-vv_SHU, 'd', label='x-v (Shu)', linewidth = 3.0)
axQ[1,1].set_xlabel(r'x')
axQ[1,1].set_ylabel(r'x-v')
axQ[1,1].legend()
# Time to convert to physical quantities
yr = 3.154e+7 #s
kyr = 1000.0*yr
km = 1e5 #cm
AU = 1.496e+13 #cm
Msun = 1.98847e33 #g
cs0 = 20000 #cs cm/s "a" in Shu notation
tt_list = np.linspace(10*kyr, 20.0*kyr, num=4)
mm = get_m(xx_RK, vv_RK, alpha_RK)
fig, ax = plt.subplots(nrows=1, ncols=3, sharex=True)
for tt in tt_list:
rho = alpha_to_rho(alpha_RK, tt)
RR = xx_RK*(cs0*tt)
time = r'%.2f $\mathrm{kyr}$' % (tt/kyr)
ax[0].plot(RR/AU, rho, label= r'$\rho$, t = ' + time, linewidth = 3.0)
ax[0].set_xlabel(r'R (AU)')
ax[0].set_ylabel(r'$\rho$ (g/cm$^3$)')
ax[0].set_xscale('log')
ax[0].set_yscale('log')
ax[0].legend()
uu = vv_to_uu(vv_RK, cs0)
ax[1].plot(RR/AU, -uu/km, label= r'$u$, t = ' + time, linewidth = 3.0)
ax[1].set_xlabel(r'R (AU)')
ax[1].set_ylabel(r'-$u$ (km/s)')
ax[1].set_yscale('log')
ax[1].legend()
MM = mm_to_MM(mm, tt, cs0)
ax[2].plot(RR/AU, MM/Msun, label= r'$M$, t = ' + time, linewidth = 3.0)
ax[2].set_xlabel(r'R (AU)')
ax[2].set_ylabel(r'$M$ ($M_\odot}$)')
ax[2].legend()
plt.show()

49
analysis/python/samples/readtest.py
View File

@@ -42,10 +42,7 @@ AC_unit_length = 1.496e+13
print("sys.argv", sys.argv)
#meshdir = "/tiara/home/mvaisala/astaroth-code/astaroth_2.0/build/"
#meshdir = "/tiara/ara/data/mvaisala/tmp/astaroth-code/astaroth_2.0/build/"
#meshdir = "/tiara/ara/data/mvaisala/asth_testbed_double/"
meshdir = "/home/mvaisala/astaroth/build/"
meshdir = "$HOME/astaroth/build/"
if "xtopbound" in sys.argv:
for i in range(0, 171):
@@ -169,36 +166,58 @@ if '1d' in sys.argv:
if 'sl' in sys.argv:
#maxfiles = 200002
#stride = 10000
maxfiles = 50001
stride = 500
maxfiles = 500000001
stride = 1
for i in range(0, maxfiles, stride):
#mesh = ad.read.Mesh(i, fdir=meshdir)
mesh = ad.read.Mesh(i, fdir=meshdir)
print(" %i / %i" % (i, maxfiles))
if mesh.ok:
uu_tot = np.sqrt(mesh.uu[0]**2.0 + mesh.uu[1]**2.0 + mesh.uu[2]**2.0)
aa_tot = np.sqrt(mesh.aa[0]**2.0 + mesh.aa[1]**2.0 + mesh.aa[2]**2.0)
mesh.Bfield()
bb_tot = np.sqrt(mesh.bb[0]**2.0 + mesh.bb[1]**2.0 + mesh.bb[2]**2.0)
if 'lim' in sys.argv:
vis.slices.plot_3(mesh, mesh.lnrho, title = r'$\ln \rho$', bitmap = True, fname = 'lnrho', colrange=[0.7, 1.15])
vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$\rho$', bitmap = True, fname = 'rho', colrange=[2.0, 4.0])
vis.slices.plot_3(mesh, mesh.uu[0], title = r'$u_x$', bitmap = True, fname = 'uux', colrange=[-1.00, 1.00])
vis.slices.plot_3(mesh, mesh.uu[1], title = r'$u_y$', bitmap = True, fname = 'uuy', colrange=[-1.00, 1.00])
vis.slices.plot_3(mesh, mesh.uu[2], title = r'$u_z$', bitmap = True, fname = 'uuz', colrange=[-1.00, 1.00])
vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$N_\mathrm{col}$', bitmap = True, fname = 'colden', slicetype = 'sum', colrange=[330.0, 395.0])
vis.slices.plot_3(mesh, uu_tot, title = r'$|u|$', bitmap = True, fname = 'uutot', colrange=[0.00, 1.2])
vis.slices.plot_3(mesh, mesh.lnrho, title = r'$\ln \rho$', bitmap = True, fname = 'lnrho', colrange=[-0.0, 3.5])
vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$\rho$', bitmap = True, fname = 'rho', colrange=[1.0, 5.0 ])
vis.slices.plot_3(mesh, mesh.uu[0], title = r'$u_x$', bitmap = True, fname = 'uux', colrange=[-5.0, 5.0])
vis.slices.plot_3(mesh, mesh.uu[1], title = r'$u_y$', bitmap = True, fname = 'uuy', colrange=[-5.0, 5.0])
vis.slices.plot_3(mesh, mesh.uu[2], title = r'$u_z$', bitmap = True, fname = 'uuz', colrange=[-5.0, 5.0])
vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$N_\mathrm{col}$', bitmap = True, fname = 'colden', slicetype = 'sum', colrange=[0.0, 200.0])
vis.slices.plot_3(mesh, uu_tot, title = r'$|u|$', bitmap = True, fname = 'uutot', colrange=[0.00, 5.0])
vis.slices.plot_3(mesh, aa_tot, title = r'$\|A\|$', bitmap = True, fname = 'aatot', colrange=[0.0,0.01])
vis.slices.plot_3(mesh, mesh.aa[0], title = r'$A_x$', bitmap = True, fname = 'aax', colrange=[-0.01,0.01])
vis.slices.plot_3(mesh, mesh.aa[1], title = r'$A_y$', bitmap = True, fname = 'aay', colrange=[-0.01,0.01])
vis.slices.plot_3(mesh, mesh.aa[2], title = r'$A_z$', bitmap = True, fname = 'aaz', colrange=[-0.01,0.01])
vis.slices.plot_3(mesh, mesh.accretion, title = r'$Accretion$', bitmap = True, fname = 'accretion', colrange=[0.0,0.000001])
vis.slices.plot_3(mesh, bb_tot, title = r'$\|B\|$', bitmap = True, fname = 'bbtot', colrange=[0.0,1.0e-4])
vis.slices.plot_3(mesh, mesh.bb[0], title = r'$B_x$', bitmap = True, fname = 'bbx', colrange=[-1.0e-4, 1.0e-4])#, bfieldlines=True)
vis.slices.plot_3(mesh, mesh.bb[1], title = r'$B_y$', bitmap = True, fname = 'bby', colrange=[-1.0e-4, 1.0e-4])#, bfieldlines=True)
vis.slices.plot_3(mesh, mesh.bb[2], title = r'$B_z$', bitmap = True, fname = 'bbz', colrange=[-1.0e-4, 1.0e-4])#, bfieldlines=True)
else:
vis.slices.plot_3(mesh, mesh.lnrho, title = r'$\ln \rho$', bitmap = True, fname = 'lnrho')
vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$\rho$', bitmap = True, fname = 'rho')
#vis.slices.plot_3(mesh, mesh.ss, title = r'$s$', bitmap = True, fname = 'ss')
vis.slices.plot_3(mesh, mesh.uu[0], title = r'$u_x$', bitmap = True, fname = 'uux')
vis.slices.plot_3(mesh, mesh.uu[0], title = r'$u_x$', bitmap = True, fname = 'uux')#, velfieldlines=True)
vis.slices.plot_3(mesh, mesh.uu[1], title = r'$u_y$', bitmap = True, fname = 'uuy')
vis.slices.plot_3(mesh, mesh.uu[2], title = r'$u_z$', bitmap = True, fname = 'uuz')
vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$N_\mathrm{col}$', bitmap = True, fname = 'colden', slicetype = 'sum')
vis.slices.plot_3(mesh, uu_tot, title = r'$|u|$', bitmap = True, fname = 'uutot')
vis.slices.plot_3(mesh, mesh.accretion, title = r'$Accretion$', bitmap = True, fname = 'accretion')
vis.slices.plot_3(mesh, aa_tot, title = r'$\|A\|$', bitmap = True, fname = 'aatot')
vis.slices.plot_3(mesh, mesh.aa[0], title = r'$A_x$', bitmap = True, fname = 'aax')
vis.slices.plot_3(mesh, mesh.aa[1], title = r'$A_y$', bitmap = True, fname = 'aay')
vis.slices.plot_3(mesh, mesh.aa[2], title = r'$A_z$', bitmap = True, fname = 'aaz')
vis.slices.plot_3(mesh, bb_tot, title = r'$\|B\|$', bitmap = True, fname = 'bbtot')
vis.slices.plot_3(mesh, mesh.bb[0], title = r'$B_x$', bitmap = True, fname = 'bbx')#, bfieldlines=True)
vis.slices.plot_3(mesh, mesh.bb[1], title = r'$B_y$', bitmap = True, fname = 'bby')#, bfieldlines=True)
vis.slices.plot_3(mesh, mesh.bb[2], title = r'$B_z$', bitmap = True, fname = 'bbz')#, bfieldlines=True)
if 'ts' in sys.argv:
ts = ad.read.TimeSeries(fdir=meshdir)
vis.lineplot.plot_ts(ts, show_all=True)
#vis.lineplot.plot_ts(ts, isotherm=True)

6
src/core/device.cu
View File

@@ -799,7 +799,7 @@ mod(const int a, const int b)
return r < 0 ? r + b : r;
}
static int
static inline int
get_neighbor(const int3 offset)
{
// The number of nodes is n^3 = m = num_processes
@@ -1246,7 +1246,9 @@ acDeviceRunMPITest(void)
////////////////////////////////////////////////////////////////////////////////////////////////
Device device;
acDeviceCreate(0, submesh_info, &device);
int devices_per_node = -1;
cudaGetDeviceCount(&devices_per_node);
acDeviceCreate(pid % devices_per_node, submesh_info, &device);
acDeviceLoadMesh(device, STREAM_DEFAULT, submesh);
// Warmup

26
src/standalone/benchmark.cc
View File

@@ -40,6 +40,12 @@
#include <math.h>
#include <vector>
static bool
smaller_than(const double& a, const double& b)
{
return a < b;
}
int
run_benchmark(const char* config_path)
{
@@ -59,21 +65,37 @@ run_benchmark(const char* config_path)
acInit(mesh_info);
acLoad(*mesh);
std::vector<double> results;
results.reserve(num_iters);
// Warmup
for (int i = 0; i < 10; ++i) {
acIntegrate(0);
}
acSynchronize();
const AcReal dt = FLT_EPSILON;
Timer total_time;
timer_reset(&total_time);
Timer step_time;
for (int i = 0; i < num_iters; ++i) {
const AcReal dt = FLT_EPSILON;
timer_reset(&step_time);
acIntegrate(dt);
acSynchronize();
results.push_back(timer_diff_nsec(step_time) / 1e6);
}
acSynchronizeStream(STREAM_ALL);
acSynchronize();
const double ms_elapsed = timer_diff_nsec(total_time) / 1e6;
const double nth_percentile = 0.95;
std::sort(results.begin(), results.end(), smaller_than);
printf("vertices: %d^3, iterations: %d\n", nn, num_iters);
printf("Total time: %f ms\n", ms_elapsed);
printf("%dth percentile per step: %f ms\n", int(100 * nth_percentile),
results[int(nth_percentile * num_iters)]);
acQuit();
acmesh_destroy(mesh);