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

Filed line integration and other smaller python tools.

Browse Source
This commit is contained in:
Miikka Vaisala
2021-01-11 11:27:12 +08:00
parent 7878811820
commit 2dbf703c59
3 changed files with 473 additions and 26 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -21,6 +21,7 @@
import numpy as np
import os
import pandas as pd
#Optional YT interface
try:
@@ -302,9 +303,116 @@ class Mesh:
self.jj = curl_of_curl(self.aa, self.minfo)
if trim:
self.bb = ( self.bb[0][3:-3, 3:-3, 3:-3],self.bb[1][3:-3, 3:-3, 3:-3],self.bb[2][3:-3, 3:-3, 3:-3])
self.xx_trim = self.xx[3:-3]
self.yy_trim = self.yy[3:-3]
self.zz_trim = self.zz[3:-3]
if get_jj:
self.jj = (self.jj[0][3:-3, 3:-3, 3:-3],self.jj[1][3:-3, 3:-3, 3:-3],self.jj[2][3:-3, 3:-3, 3:-3])
def Bfieldlines(self, footloc = 'default', vartype = 'B', maxstep = 1000):
dx = self.minfo.contents['AC_dsx']
dy = self.minfo.contents['AC_dsy']
dz = self.minfo.contents['AC_dsz']
if vartype == 'U':
#Trim to match
self.uu = (self.uu[0][3:-3, 3:-3, 3:-3],self.uu[1][3:-3, 3:-3, 3:-3],self.uu[2][3:-3, 3:-3, 3:-3])
def field_line_step(self, coord, ds):
#TODO assume that grid is at a cell centre
ix = np.argmin(np.abs(self.xx_trim - coord[0]))
iy = np.argmin(np.abs(self.yy_trim - coord[1]))
iz = np.argmin(np.abs(self.zz_trim - coord[2]))
if vartype == 'U':
Bcell_vec = np.array([self.uu[0][ix, iy, iz],
self.uu[1][ix, iy, iz],
self.uu[2][ix, iy, iz]])
else:
Bcell_vec = np.array([self.bb[0][ix, iy, iz],
self.bb[1][ix, iy, iz],
self.bb[2][ix, iy, iz]])
Bcell_abs = np.sqrt(Bcell_vec[0]**2.0 + Bcell_vec[1]**2.0 + Bcell_vec[2]**2.0)
coord_new = coord + (Bcell_vec/Bcell_abs)*ds
return coord_new
self.df_lines = pd.DataFrame()
ds = np.amin([self.minfo.contents['AC_dsx'],
self.minfo.contents['AC_dsy'],
self.minfo.contents['AC_dsz']])
ii = 0
if footloc == 'middlez':
ixtot = 6
iytot = 6
iztot = 1
xfoots = np.linspace(self.xx_trim.min(), self.xx_trim.max(), num = ixtot)
yfoots = np.linspace(self.yy_trim.min(), self.yy_trim.max(), num = iytot)
zfoots = np.array([(self.zz_trim.max() - self.zz_trim.min())/2.0 + self.zz_trim.min()])
elif footloc == 'cube':
ixtot = 5
iytot = 5
iztot = 5
xfoots = np.linspace(self.xx_trim.min()+3.0*dx, self.xx_trim.max()-3.0*dx, num = ixtot)
yfoots = np.linspace(self.yy_trim.min()+3.0*dy, self.yy_trim.max()-3.0*dy, num = iytot)
zfoots = np.linspace(self.zz_trim.min()+3.0*dz, self.zz_trim.max()-3.0*dz, num = iztot)
else:
ixtot = 6
iytot = 6
iztot = 1
xfoots = np.linspace(self.xx_trim.min(), self.xx_trim.max(), num = ixtot)
yfoots = np.linspace(self.yy_trim.min(), self.yy_trim.max(), num = iytot)
zfoots = np.array([self.zz_trim.min()])
imax = ixtot * iytot * iztot
for zfoot in zfoots:
for yfoot in yfoots:
for xfoot in xfoots:
print(ii, "/", imax-1)
integrate = 1
counter = 0
dstot = 0.0
coord = np.array([xfoot, yfoot, zfoot])
self.df_lines = self.df_lines.append({"line_num":ii,
"dstot":dstot,
"coordx":coord[0],
"coordy":coord[1],
"coordz":coord[2]},
ignore_index=True)
while integrate:
coord = field_line_step(self, coord, ds)
dstot += ds
self.df_lines = self.df_lines.append({"line_num":ii,
"dstot":dstot,
"coordx":coord[0],
"coordy":coord[1],
"coordz":coord[2]},
ignore_index=True)
counter += 1
if counter >= maxstep:
integrate = 0
if ((coord[0] > self.xx_trim.max()) or
(coord[1] > self.yy_trim.max()) or
(coord[2] > self.zz_trim.max()) or
(coord[0] < self.xx_trim.min()) or
(coord[1] < self.yy_trim.min()) or
(coord[2] < self.zz_trim.min())):
#print("out of bounds")
integrate = 0
if (np.isnan(coord[0]) or
np.isnan(coord[1]) or
np.isnan(coord[2])):
integrate = 0
ii += 1
#print(self.df_lines)
def get_jj(self, trim=False):
self.jj = curl_of_curl(self.aa, minfo, trim=False)
if trim:

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

@@ -142,4 +142,103 @@ def plot_3(mesh, input_grid, title = '', fname = 'default', bitmap=False,
print('Saved %s_%s.png' % (fname, mesh.framenum))
plt.close(fig)
def volume_render(mesh, val1 = {"variable": None, "min": None, "max":None, "opacity":1.0}):
if val1["variable"] == "btot":
plt.figure()
bb_tot = np.sqrt(mesh.bb[0]**2.0 + mesh.bb[1]**2.0 + mesh.bb[2]**2.0)
array = bb_tot
varname = "btot"
meshxx = mesh.xx[3:-3]
meshyy = mesh.yy[3:-3]
meshzz = mesh.zz[3:-3]
if val1["variable"] == "utot":
plt.figure()
uu_tot = np.sqrt(mesh.uu[0]**2.0 + mesh.uu[1]**2.0 + mesh.uu[2]**2.0)
array = uu_tot
varname = "utot"
meshxx = mesh.xx
meshyy = mesh.yy
meshzz = mesh.zz
if val1["variable"] == "rho":
plt.figure()
array = np.exp(mesh.lnrho)
varname = "rho"
meshxx = mesh.xx
meshyy = mesh.yy
meshzz = mesh.zz
if val1["variable"] == "aa":
plt.figure()
aa_tot = np.sqrt(mesh.aa[0]**2.0 + mesh.aa[1]**2.0 + mesh.aa[2]**2.0)
array = aa_tot
varname = "aa"
meshxx = mesh.xx
meshyy = mesh.yy
meshzz = mesh.zz
#Histogram plot to find value ranges.
hist, bedges = np.histogram(array, bins=mesh.xx.size)
plt.plot(bedges[:-1], hist)
plt.yscale('log')
if val1["min"] != None or val1["max"] != None:
plt.plot([val1["min"],val1["min"]], [1,hist.max()], label=varname+" min")
plt.plot([val1["max"],val1["max"]], [1,hist.max()], label=varname+" max")
plt.legend()
plt.savefig('volrend_hist_%s_%s.png' % (varname, mesh.framenum))
plt.close()
if val1["min"] != None or val1["max"] != None:
#print(np.where(bb_tot < val1["min"]))
array[np.where(array < val1["min"])] = 0.0
array[np.where(array > val1["max"])] = 0.0
array[np.where(array > 0.0)] = val1["opacity"]
#plt.figure()
#plt.plot(bb_tot[:,64,64])
mapyz = array.sum(axis=0)
mapxz = array.sum(axis=1)
mapxy = array.sum(axis=2)
yy_yz, zz_yz = np.meshgrid(meshyy, meshzz, indexing='ij')
xx_xz, zz_xz = np.meshgrid(meshxx, meshzz, indexing='ij')
xx_xy, yy_xy = np.meshgrid(meshxx, meshyy, indexing='ij')
fig, ax = plt.subplots()
#plt.imshow(mapyz, vmin=0.0, vmax=1.0)
plt.pcolormesh(yy_yz, zz_yz, mapyz, vmin=0.0, vmax=1.0, shading='auto')
ax.set_aspect('equal')
ax.set_title(varname)
ax.set_xlabel('y')
ax.set_ylabel('z')
plt.savefig('volrend_%s_%s_%s.png' % (varname, "yz", mesh.framenum))
plt.close()
fig, ax = plt.subplots()
#plt.imshow(mapxz, vmin=0.0, vmax=1.0)
plt.pcolormesh(xx_xz, zz_xz, mapxz, vmin=0.0, vmax=1.0, shading='auto')
ax.set_aspect('equal')
ax.set_title(varname)
ax.set_xlabel('x')
ax.set_ylabel('z')
plt.savefig('volrend_%s_%s_%s.png' % (varname, "xz", mesh.framenum))
plt.close()
fig, ax = plt.subplots()
#plt.imshow(mapxy, vmin=0.0, vmax=1.0)
plt.pcolormesh(xx_xy, yy_xy, mapxy, vmin=0.0, vmax=1.0, shading='auto')
ax.set_aspect('equal')
ax.set_title(varname)
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.savefig('volrend_%s_%s_%s.png' % (varname, "xy", mesh.framenum))
plt.close()
#plt.show()

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

@@ -25,6 +25,16 @@ import sys
import os
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
#Optional YT interface
try:
import yt
yt_present = True
except ImportError:
yt_present = False
##mesh = ad.read.Mesh(500, fdir="/tiara/home/mvaisala/astaroth-code/astaroth_2.0/build/")
##
##print(np.shape(mesh.uu))
@@ -122,28 +132,110 @@ if "single" in sys.argv:
print( mesh.uu[2][100, 101, 00], "periodic")
if 'xline' in sys.argv:
mesh = ad.read.Mesh(0, fdir=meshdir)
plt.figure()
plt.plot(mesh.uu[0][100, 50, :] , label="z")
plt.plot(mesh.uu[0][100, :, 100], label="x")
plt.plot(mesh.uu[0][:, 50, 100] , label="y")
plt.legend()
mesh_file_numbers = ad.read.parse_directory(meshdir)
print(mesh_file_numbers)
maxfiles = np.amax(mesh_file_numbers)
plt.figure()
plt.plot(mesh.uu[0][197, 50, :] , label="z edge")
for i in mesh_file_numbers[-3:]:
mesh = ad.read.Mesh(i, fdir=meshdir)
mesh.Bfield(trim=True)
plt.figure()
plt.plot(mesh.uu[1][100, 50, :] , label="z")
plt.plot(mesh.uu[1][100, :, 100], label="x")
plt.plot(mesh.uu[1][:, 50, 100] , label="y")
plt.legend()
xhalf = int(mesh.uu[0].shape[0]/2.0)
yhalf = int(mesh.uu[0].shape[1]/2.0)
zhalf = int(mesh.uu[0].shape[2]/2.0)
print(xhalf, yhalf, zhalf)
plt.figure()
plt.plot(mesh.uu[0][xhalf, yhalf, :], label="z")
plt.plot(mesh.uu[0][xhalf, :, zhalf], label="y")
plt.plot(mesh.uu[0][ :, yhalf, zhalf], label="x")
plt.title("UUX")
plt.legend()
#plt.figure()
#plt.plot(mesh.uu[0][197, 50, :] , label="z edge")
plt.figure()
plt.plot(mesh.uu[1][xhalf, yhalf, :], label="z")
plt.plot(mesh.uu[1][xhalf, :, zhalf], label="y")
plt.plot(mesh.uu[1][ :, yhalf, zhalf], label="x")
plt.title("UUY")
plt.legend()
plt.figure()
plt.plot(mesh.uu[2][xhalf, yhalf, :], label="z")
plt.plot(mesh.uu[2][xhalf, :, zhalf], label="y")
plt.plot(mesh.uu[2][ :, yhalf, zhalf], label="x")
plt.legend()
plt.title("UUZ")
plt.figure()
plt.plot(mesh.bb[0][xhalf, yhalf, :], label="z")
plt.plot(mesh.bb[0][xhalf, :, zhalf], label="y")
plt.plot(mesh.bb[0][ :, yhalf, zhalf], label="x")
plt.title("BBX")
plt.legend()
plt.figure()
plt.plot(mesh.bb[1][xhalf, yhalf, :], label="z")
plt.plot(mesh.bb[1][xhalf, :, zhalf], label="y")
plt.plot(mesh.bb[1][ :, yhalf, zhalf], label="x")
plt.title("BBY")
plt.legend()
plt.figure()
plt.plot(mesh.bb[2][xhalf, yhalf, :], label="z")
plt.plot(mesh.bb[2][xhalf, :, zhalf], label="y")
plt.plot(mesh.bb[2][ :, yhalf, zhalf], label="x")
plt.legend()
plt.title("BBZ")
plt.figure()
plt.plot(mesh.aa[0][xhalf, yhalf, :], label="z")
plt.plot(mesh.aa[0][xhalf, :, zhalf], label="y")
plt.plot(mesh.aa[0][ :, yhalf, zhalf], label="x")
plt.title("AX")
plt.legend()
plt.figure()
plt.plot(mesh.aa[1][xhalf, yhalf, :], label="z")
plt.plot(mesh.aa[1][xhalf, :, zhalf], label="y")
plt.plot(mesh.aa[1][ :, yhalf, zhalf], label="x")
plt.title("AY")
plt.legend()
plt.figure()
plt.plot(mesh.aa[2][xhalf, yhalf, :], label="z")
plt.plot(mesh.aa[2][xhalf, :, zhalf], label="y")
plt.plot(mesh.aa[2][ :, yhalf, zhalf], label="x")
plt.legend()
plt.title("AZ")
uu_tot = np.sqrt(mesh.uu[0]**2.0 + mesh.uu[1]**2.0 + mesh.uu[2]**2.0)
bb_tot = np.sqrt(mesh.bb[0]**2.0 + mesh.bb[1]**2.0 + mesh.bb[2]**2.0)
plt.figure()
plt.plot(uu_tot[xhalf, yhalf, :], label="z")
plt.plot(uu_tot[xhalf, :, zhalf], label="y")
plt.plot(uu_tot[ :, yhalf, zhalf], label="x")
plt.legend()
plt.title("UTOT")
plt.figure()
plt.plot(bb_tot[xhalf, yhalf, :], label="z")
plt.plot(bb_tot[xhalf, :, zhalf], label="y")
plt.plot(bb_tot[ :, yhalf, zhalf], label="x")
plt.legend()
plt.title("BTOT")
plt.figure()
plt.plot(np.exp(mesh.lnrho[xhalf, yhalf, :]), label="z")
plt.plot(np.exp(mesh.lnrho[xhalf, :, zhalf]), label="y")
plt.plot(np.exp(mesh.lnrho[ :, yhalf, zhalf]), label="x")
plt.legend()
plt.title("RHO")
plt.figure()
plt.plot(mesh.uu[2][100, 50, :] , label="z")
plt.plot(mesh.uu[2][100, :, 100], label="x")
plt.plot(mesh.uu[2][:, 50, 100] , label="y")
plt.legend()
plt.show()
if 'check' in sys.argv:
mesh = ad.read.Mesh(0, fdir=meshdir)
@@ -180,22 +272,43 @@ if '1d' in sys.argv:
plt.legend()
plt.show()
if 'csv' in sys.argv:
filenum = sys.argv[1]
mesh = ad.read.Mesh(filenum, fdir=meshdir)
mesh.Bfield()
mesh.export_csv()
if 'raw' in sys.argv:
filenum = sys.argv[1]
mesh = ad.read.Mesh(filenum, fdir=meshdir)
mesh.Bfield()
mesh.export_raw()
if 'findnan' in sys.argv:
filenum = sys.argv[1]
mesh = ad.read.Mesh(filenum, fdir=meshdir)
print("nan uu", np.where(np.isnan(mesh.uu)))
print("nan aa", np.where(np.isnan(mesh.aa)))
print("nan lnrho", np.where(np.isnan(mesh.lnrho)))
print("inf uu", np.where(np.isinf(mesh.uu)))
print("inf aa", np.where(np.isinf(mesh.aa)))
print("inf lnrho", np.where(np.isinf(mesh.lnrho)))
if 'sl' in sys.argv:
#maxfiles = 200002
#stride = 10000
maxfiles = 500000001
stride = 1
for i in range(0, maxfiles, stride):
#mesh = ad.read.Mesh(i, fdir=meshdir)
mesh_file_numbers = ad.read.parse_directory(meshdir)
print(mesh_file_numbers)
maxfiles = np.amax(mesh_file_numbers)
for i in mesh_file_numbers:
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()
mesh.Bfield(trim=True)
bb_tot = np.sqrt(mesh.bb[0]**2.0 + mesh.bb[1]**2.0 + mesh.bb[2]**2.0)
if 'sym' in sys.argv:
@@ -248,6 +361,133 @@ if 'sl' in sys.argv:
vis.slices.plot_3(mesh, mesh.bb[1], title = r'$B_y$', bitmap = True, fname = 'bby', trimghost=3)#, bfieldlines=True)
vis.slices.plot_3(mesh, mesh.bb[2], title = r'$B_z$', bitmap = True, fname = 'bbz', trimghost=3)#, bfieldlines=True)
if 'yt' in sys.argv:
mesh.yt_conversion()
from mpl_toolkits.axes_grid1 import AxesGrid
coords = ['x', 'y','z']
for coord in coords:
fields = ['density', 'uux', 'uuy', 'uuz']
fig = plt.figure()
grid = AxesGrid(fig, (0.075,0.075,0.85,0.85),
nrows_ncols = (2, 2),
axes_pad = 1.0,
label_mode = "1",
share_all = True,
cbar_location="right",
cbar_mode="each",
cbar_size="3%",
cbar_pad="0%")
p = yt.SlicePlot(mesh.ytdata, coord, fields)
p.set_log('uux', False)
p.set_log('uuy', False)
p.set_log('uuz', False)
for i, field in enumerate(fields):
plot = p.plots[field]
plot.figure = fig
plot.axes = grid[i].axes
plot.cax = grid.cbar_axes[i]
p._setup_plots()
plt.savefig('yt_rho_uu_%s_%s.png' % (coord, mesh.framenum))
plt.close(fig=fig)
###
fields = ['density', 'bbx', 'bby', 'bbz']
fig = plt.figure()
grid = AxesGrid(fig, (0.075,0.075,0.85,0.85),
nrows_ncols = (2, 2),
axes_pad = 1.0,
label_mode = "1",
share_all = True,
cbar_location="right",
cbar_mode="each",
cbar_size="3%",
cbar_pad="0%")
p = yt.SlicePlot(mesh.ytdata, coord, fields)
p.set_log('bbx', False)
p.set_log('bby', False)
p.set_log('bbz', False)
for i, field in enumerate(fields):
plot = p.plots[field]
plot.figure = fig
plot.axes = grid[i].axes
plot.cax = grid.cbar_axes[i]
p._setup_plots()
plt.savefig('yt_rho_bb_%s_%s.png' % (coord, mesh.framenum))
plt.close(fig=fig)
elif 'csvall' in sys.argv:
mesh.export_csv()
elif 'rawall' in sys.argv:
mesh.export_raw()
if "vol" in sys.argv:
print("VOLUME RENDERING")
mesh_file_numbers = ad.read.parse_directory(meshdir)
print(mesh_file_numbers)
maxfiles = np.amax(mesh_file_numbers)
for i in mesh_file_numbers:
mesh = ad.read.Mesh(i, fdir=meshdir)
mesh.Bfield(trim=True)
print(" %i / %i" % (i, maxfiles))
if mesh.ok:
vis.slices.volume_render(mesh, val1 = {"variable": "btot", "min":0.5, "max":2.0, "opacity":0.05})
vis.slices.volume_render(mesh, val1 = {"variable": "utot", "min":0.5, "max":2.0, "opacity":0.05})
vis.slices.volume_render(mesh, val1 = {"variable": "rho", "min":10.0, "max":300.0, "opacity":0.05})
vis.slices.volume_render(mesh, val1 = {"variable": "aa", "min":0.1, "max":0.25, "opacity":0.05})
if (("bline" in sys.argv) or ("uline" in sys.argv)):
print("Field line computation")
mesh_file_numbers = ad.read.parse_directory(meshdir)
print(mesh_file_numbers)
maxfiles = np.amax(mesh_file_numbers)
for i in mesh_file_numbers:
mesh = ad.read.Mesh(i, fdir=meshdir)
mesh.Bfield(trim=True)
print(" %i / %i" % (i, maxfiles))
if mesh.ok:
if "uline" in sys.argv:
mesh.Bfieldlines(footloc = 'cube', vartype='U', maxstep = 200)
else:
mesh.Bfieldlines(footloc = 'default')
print(mesh.df_lines)
fig = plt.figure(figsize=(5.0,5.0))
ax = fig.gca(projection='3d')
for line_num in range(int(mesh.df_lines['line_num'].max()+1)):
df_myline = mesh.df_lines.loc[mesh.df_lines['line_num'] == line_num]
print(df_myline)
my_xscale = [mesh.xx_trim.min(), mesh.xx_trim.max()]
my_yscale = [mesh.yy_trim.min(), mesh.yy_trim.max()]
my_zscale = [mesh.zz_trim.min(), mesh.zz_trim.max()]
ax.plot(df_myline["coordx"], df_myline["coordy"], df_myline["coordz"], color="red")
ax.set_xlim3d(my_xscale)
ax.set_ylim3d(my_yscale)
ax.set_zlim3d(my_zscale)
if "uline" in sys.argv:
filename = 'Ugeometry_%s.png' % (mesh.framenum)
else:
filename = 'Bgeometry_%s.png' % (mesh.framenum)
plt.savefig(filename)
plt.close()
#plt.show()
if 'ts' in sys.argv:
ts = ad.read.TimeSeries(fdir=meshdir)