From 2dbf703c59ec6a3167b73242a943e511082160d9 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Mon, 11 Jan 2021 11:27:12 +0800 Subject: [PATCH] Filed line integration and other smaller python tools. --- analysis/python/astar/data/read.py | 108 +++++++++ analysis/python/astar/visual/slices.py | 99 +++++++++ analysis/python/samples/readtest.py | 292 ++++++++++++++++++++++--- 3 files changed, 473 insertions(+), 26 deletions(-) diff --git a/analysis/python/astar/data/read.py b/analysis/python/astar/data/read.py index 3ae6ff8..4fa8552 100644 --- a/analysis/python/astar/data/read.py +++ b/analysis/python/astar/data/read.py @@ -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: diff --git a/analysis/python/astar/visual/slices.py b/analysis/python/astar/visual/slices.py index 4d95ce8..5c151e5 100644 --- a/analysis/python/astar/visual/slices.py +++ b/analysis/python/astar/visual/slices.py @@ -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() diff --git a/analysis/python/samples/readtest.py b/analysis/python/samples/readtest.py index ddede1a..6392f07 100644 --- a/analysis/python/samples/readtest.py +++ b/analysis/python/samples/readtest.py @@ -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)