From b2497441292a719e6cb38d78e325bfb54ea0d9a1 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Mon, 11 Dec 2023 10:55:20 -0700 Subject: [PATCH 01/11] Starting in on this --- scripts/python/get_torus_fluxes.py | 84 ++++++++++++------------- scripts/python/grmhd_radial_profiles.py | 57 +++++++++++++++++ scripts/python/phoedf.py | 33 +++++++++- scripts/python/plot_torus.py | 2 +- scripts/python/plot_torus_3d.py | 8 +-- scripts/python/torus_analyzer.py | 15 ++++- 6 files changed, 149 insertions(+), 50 deletions(-) create mode 100644 scripts/python/grmhd_radial_profiles.py diff --git a/scripts/python/get_torus_fluxes.py b/scripts/python/get_torus_fluxes.py index ffd27db30..cd4d3486d 100644 --- a/scripts/python/get_torus_fluxes.py +++ b/scripts/python/get_torus_fluxes.py @@ -1,6 +1,6 @@ #!/usr/bin/env python -# © 2022. Triad National Security, LLC. All rights reserved. This +# © 2022-2023. Triad National Security, LLC. All rights reserved. This # program was produced under U.S. Government contract # 89233218CNA000001 for Los Alamos National Laboratory (LANL), which # is operated by Triad National Security, LLC for the U.S. Department @@ -30,7 +30,6 @@ import phoebus_utils from phoedf import phoedf - def get_torus_fluxes(dfile): fluxes = {} @@ -72,10 +71,8 @@ def get_torus_fluxes(dfile): Pjm = 0.0 Pjr = 0.0 for b in range(dfile.NumBlocks): - blockBounds = dfile.BlockBounds[b] - - dx2 = (blockBounds[3] - blockBounds[2]) / dfile.MeshBlockSize[1] - dx3 = (blockBounds[5] - blockBounds[4]) / dfile.MeshBlockSize[2] + dx2 = dfile.Dx2[b] + dx3 = dfile.Dx3[b] # Block contains event horizon if blockBounds[0] < xh and blockBounds[1] > xh: @@ -94,25 +91,26 @@ def get_torus_fluxes(dfile): Eg_in += -dx2 * dx3 * gdet[b, k, j, i_eh] * Tmunu_concov[1, 0] Lg_in += dx2 * dx3 * gdet[b, k, j, i_eh] * Tmunu_concov[1, 3] - for j in range(Nx2): - for k in range(Nx3): - Rmunu_concov = dfile.GetRmunu_concov(b, k, j, i_eh) - for ispec in range(dfile.NumSpecies): - Er_in += ( - -dx2 * dx3 * gdet[b, k, j, i_eh] * Rmunu_concov[1, 0, ispec] - ) - Lr_in += ( - dx2 * dx3 * gdet[b, k, j, i_eh] * Rmunu_concov[1, 3, ispec] - ) + if dfile.RadiationActive: + for j in range(Nx2): + for k in range(Nx3): + Rmunu_concov = dfile.GetRmunu_concov(b, k, j, i_eh) + for ispec in range(dfile.NumSpecies): + Er_in += ( + -dx2 * dx3 * gdet[b, k, j, i_eh] * Rmunu_concov[1, 0, ispec] + ) + Lr_in += ( + dx2 * dx3 * gdet[b, k, j, i_eh] * Rmunu_concov[1, 3, ispec] + ) - for j in range(Nx2): - for k in range(Nx3): - Phi += ( - dx2 - * dx3 - * gdet[b, k, j, i_eh] - * np.fabs(Bcon[b, 0, k, j, i_eh] / alpha[b, k, j, i_eh]) - ) + for j in range(Nx2): + for k in range(Nx3): + Phi += ( + dx2 + * dx3 + * gdet[b, k, j, i_eh] + * np.fabs(Bcon[b, 0, k, j, i_eh] / alpha[b, k, j, i_eh]) + ) if blockBounds[0] < 1.6 and blockBounds[1] > 1.6: # Nx1 // 2 is kind of cheating... @@ -121,7 +119,6 @@ def get_torus_fluxes(dfile): for k in range(Nx3): if sigma[b, k, j, i] > 1.0: Tmunu_concov = dfile.GetTmunu_concov(b, k, j, i) - Rmunu_concov = dfile.GetRmunu_concov(b, k, j, i) Pj -= dx2 * dx3 * gdet[b, k, j, i] * Tmunu_concov[1, 0] Pjm -= ( dx2 @@ -130,10 +127,12 @@ def get_torus_fluxes(dfile): * rho[b, k, j, i] * ucon[b, 1, k, j, i] ) - for ispec in range(dfile.NumSpecies): - Pjr -= ( - dx2 * dx3 * gdet[b, k, j, i] * Rmunu_concov[1, 0, ispec] - ) + if dfile.RadiationActive: + Rmunu_concov = dfile.GetRmunu_concov(b, k, j, i) + for ispec in range(dfile.NumSpecies): + Pjr -= ( + dx2 * dx3 * gdet[b, k, j, i] * Rmunu_concov[1, 0, ispec] + ) # Block contains outer boundary if np.fabs(blockBounds[1] - x1max) / x1max < 1.0e-10: @@ -147,16 +146,17 @@ def get_torus_fluxes(dfile): Eg_out += -dx2 * dx3 * gdet[b, k, j, -1] * Tmunu_concov[1, 0] Lg_out += dx2 * dx3 * gdet[b, k, j, -1] * Tmunu_concov[1, 3] - for j in range(Nx2): - for k in range(Nx3): - Rmunu_concov = dfile.GetRmunu_concov(b, k, j, -1) - for ispec in range(dfile.NumSpecies): - Er_out += ( - -dx2 * dx3 * gdet[b, k, j, -1] * Rmunu_concov[1, 0, ispec] - ) - Lr_out += ( - dx2 * dx3 * gdet[b, k, j, -1] * Rmunu_concov[1, 3, ispec] - ) + if dfile.RadiationActive: + for j in range(Nx2): + for k in range(Nx3): + Rmunu_concov = dfile.GetRmunu_concov(b, k, j, -1) + for ispec in range(dfile.NumSpecies): + Er_out += ( + -dx2 * dx3 * gdet[b, k, j, -1] * Rmunu_concov[1, 0, ispec] + ) + Lr_out += ( + dx2 * dx3 * gdet[b, k, j, -1] * Rmunu_concov[1, 3, ispec] + ) fluxes["mdot_in"] = mdot_in fluxes["mdot_edd"] = ( @@ -191,14 +191,14 @@ def get_torus_fluxes(dfile): parser = argparse.ArgumentParser(description="Get fluxes from torus dump") parser.add_argument( - "files", type=str, nargs="+", help="Files to take a snapshot of" + "filenames", type=str, nargs="+", help="Files to derive fluxes for" ) args = parser.parse_args() logfile = open("fluxes_logfile.txt", "a") - for n, fname in enumerate(args.files): - print(f"Opening file {fname}... ", end="") + for n, fname in enumerate(args.filenames): + print(f"Opening file {fname}... ", end="", flush=True) dfile = phoedf(fname) print("done") diff --git a/scripts/python/grmhd_radial_profiles.py b/scripts/python/grmhd_radial_profiles.py new file mode 100644 index 000000000..de2670722 --- /dev/null +++ b/scripts/python/grmhd_radial_profiles.py @@ -0,0 +1,57 @@ +#!/usr/bin/env python + +# © 2023. Triad National Security, LLC. All rights reserved. This +# program was produced under U.S. Government contract +# 89233218CNA000001 for Los Alamos National Laboratory (LANL), which +# is operated by Triad National Security, LLC for the U.S. Department +# of Energy/National Nuclear Security Administration. All rights in +# the program are reserved by Triad National Security, LLC, and the +# U.S. Department of Energy/National Nuclear Security +# Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +# license in this material to reproduce, prepare derivative works, +# distribute copies to the public, perform publicly and display +# publicly, and to permit others to do so. + +import argparse +import numpy as np +import sys +import os +import phoebus_utils +from phoedf import phoedf + +def get_torus_radial_profiles(dfile): + + profiles = {} + + gdet = dfile.gdet + dx1 = dfile.dx1 + + print(dx1) + + for b in range(dfile.NumBlocks): + dx2 = dfile.Dx2[b] + dx3 = dfile.Dx3[b] + print(f"b: {b} dx2: {dx2} dx3: {dx3})") + + + return profiles + +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description="Get radial profiles from torus dump") + parser.add_argument( + "filenames", type=str, nargs="+", help="Files to derive radial profiles for" + ) + args = parser.parse_args() + + print("Phoebus GRMHD analysis script for generating radial profiles from dumps") + print(f" Number of files: {len(args.filenames)}") + + for filename in args.filenames: + print(f"Opening file {os.path.basename(filename)}... ", end="", flush=True) + dfile = phoedf(filename) + print("done") + + profiles = get_torus_radial_profiles(dfile) + diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py index 2d627f5cd..6cbb59d32 100644 --- a/scripts/python/phoedf.py +++ b/scripts/python/phoedf.py @@ -53,6 +53,17 @@ def __init__(self, filename, geomfile=None): self.Nx2 = self.MeshBlockSize[1] self.Nx3 = self.MeshBlockSize[2] + self.Dx1 = np.zeros([self.NumBlocks]) + self.Dx2 = np.zeros([self.NumBlocks]) + self.Dx3 = np.zeros([self.NumBlocks]) + for b in range(self.NumBlocks): + bounds = self.BlockBounds[b] + self.Dx1[b] = (bounds[1] - bounds[0])/self.MeshBlockSize[0] + self.Dx2[b] = (bounds[3] - bounds[2])/self.MeshBlockSize[1] + self.Dx3[b] = (bounds[5] - bounds[4])/self.MeshBlockSize[2] + + assert self.MaxLevel == 0, "phoedf does not currently support mesh refinement!" + self.ScalarField = [self.NumBlocks, self.Nx3, self.Nx2, self.Nx1] self.ThreeVectorField = [self.NumBlocks, 3, self.Nx3, self.Nx2, self.Nx1] self.FourVectorField = [self.NumBlocks, 4, self.Nx3, self.Nx2, self.Nx1] @@ -106,6 +117,8 @@ def flatten_indices(mu, nu): else: kstart = None kend = None + kstart = None + kend = None # Unroll gcov for convenience for mu in range(4): @@ -399,6 +412,9 @@ def Getucon(self): return self.ucon def GetXi(self): + if not self.RadiationActive: + return None + if self.xi is None: self.xi = np.zeros( [self.NumBlocks, self.NumSpecies, self.Nx3, self.Nx2, self.Nx1] @@ -425,6 +441,9 @@ def GetXi(self): return self.xi def GetE(self): + if not self.RadiationActive: + return None + if self.E is None: self.E = np.zeros( [self.NumBlocks, self.NumSpecies, self.Nx3, self.Nx2, self.Nx1] @@ -452,6 +471,9 @@ def GetE(self): return self.E def GetF(self): + if not self.RadiationActive: + return None + if self.F is None: self.F = np.zeros( [self.NumBlocks, 3, self.NumSpecies, self.Nx3, self.Nx2, self.Nx1] @@ -492,6 +514,9 @@ def GetF(self): return self.F def GetP(self): + if not self.RadiationActive: + return None + if self.P is None: self.P = np.zeros( [self.NumBlocks, 3, 3, self.NumSpecies, self.Nx3, self.Nx2, self.Nx1] @@ -574,6 +599,9 @@ def GetTmunu_concov(self, b, k, j, i): return Tmunu_concov def GetRmunu_concon(self, b, k, j, i): + if not self.RadiationActive: + return None + Rmunu = np.zeros([4, 4, self.NumSpecies]) ncon = np.zeros(4) @@ -603,6 +631,9 @@ def GetRmunu_concon(self, b, k, j, i): return Rmunu def GetRmunu_concov(self, b, k, j, i): + if not self.RadiationActive: + return None + Rmunu_concon = self.GetRmunu_concon(b, k, j, i) gcov = self.gcov[b, :, :, k, j, i] @@ -619,4 +650,4 @@ def GetRmunu_concov(self, b, k, j, i): def GetMdotEddington(self, eff=0.1): Mbh = self.LengthCodeToCGS * cgs["CL"] ** 2 / cgs["GNEWT"] - return 1.4e18 * Mbh / cgs["MSOLAR"] # Nominal eff = 0.1 + return 1.4e19 * eff * Mbh / cgs["MSOLAR"] diff --git a/scripts/python/plot_torus.py b/scripts/python/plot_torus.py index efc596ceb..e08d50847 100644 --- a/scripts/python/plot_torus.py +++ b/scripts/python/plot_torus.py @@ -24,7 +24,7 @@ import os from subprocess import call, DEVNULL import glob -from parthenon_tools import phdf +#from parthenon_tools import phdf import time from enum import Enum from phoebus_constants import cgs, scalefree diff --git a/scripts/python/plot_torus_3d.py b/scripts/python/plot_torus_3d.py index 0e6f680dc..1e7a2e1f2 100644 --- a/scripts/python/plot_torus_3d.py +++ b/scripts/python/plot_torus_3d.py @@ -120,10 +120,10 @@ def plot_frame_from_phoedf( ldensity = np.log10(dfile.GetRho()) density = dfile.Get("p.density", flatten=False) - J = dfile.Get("r.p.J", flatten=False) - # print(f'density min: {density.min()} max: {density.max()}') - # if J is not None: - # print(f'J min: {J.min()} max: {J.max()}') + if dfile.RadiationActive: + J = dfile.Get("r.p.J", flatten=False) + else: + J = None if ndim == 2: k = 0 diff --git a/scripts/python/torus_analyzer.py b/scripts/python/torus_analyzer.py index fb6243d92..c7ed4d61e 100644 --- a/scripts/python/torus_analyzer.py +++ b/scripts/python/torus_analyzer.py @@ -1,6 +1,6 @@ #!/usr/bin/env python -# © 2022. Triad National Security, LLC. All rights reserved. This +# © 2022-2023. Triad National Security, LLC. All rights reserved. This # program was produced under U.S. Government contract # 89233218CNA000001 for Los Alamos National Laboratory (LANL), which # is operated by Triad National Security, LLC for the U.S. Department @@ -38,7 +38,7 @@ def process_file(in_filename, fig_folder, log_folder, avg_folder): - print(f"Processing {in_filename}") + print(f"Processing {os.path.basename(in_filename)}") log_filename = ( os.path.join(log_folder, in_filename.split(os.sep)[-1]).rstrip(".phdf") + ".log" @@ -53,6 +53,7 @@ def process_file(in_filename, fig_folder, log_folder, avg_folder): and os.path.exists(log_filename) and os.path.exists(fig_filename) ): + print(" Output already exists! Skipping...") return # Load data @@ -142,6 +143,16 @@ def process_file(in_filename, fig_folder, log_folder, avg_folder): base_filename = args.problem_name + "*.phdf" filenames = np.sort(glob.glob(os.path.join(args.directory, base_filename))) + print("Phoebus GRMHD analysis script for generating derived quantities from dumps") + print(f" Directory: {args.directory}") + print(f" Number of files: {len(filenames)}") + print(f" Generating figures: {not args.no_figs}") + print(f" Generating logs: {not args.no_logs}") + print(f" Generating averages: {not args.no_avgs}") + print(f" Serial? {args.serial}") + print(f" Overwrite? {args.overwrite}") + print("") + fig_folder = os.path.join(args.directory, "figures") log_folder = os.path.join(args.directory, "logs") avg_folder = os.path.join(args.directory, "averages") From 652d6fb05e9bfd27da700e6433e5665bb09e4c14 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Mon, 11 Dec 2023 14:03:45 -0700 Subject: [PATCH 02/11] Still working --- scripts/python/get_torus_fluxes.py | 7 +- scripts/python/get_torus_radial_profiles.py | 145 ++++++++++++++++++++ scripts/python/grmhd_radial_profiles.py | 57 -------- scripts/python/phoedf.py | 27 +++- 4 files changed, 175 insertions(+), 61 deletions(-) create mode 100644 scripts/python/get_torus_radial_profiles.py delete mode 100644 scripts/python/grmhd_radial_profiles.py diff --git a/scripts/python/get_torus_fluxes.py b/scripts/python/get_torus_fluxes.py index cd4d3486d..0aaa5937a 100644 --- a/scripts/python/get_torus_fluxes.py +++ b/scripts/python/get_torus_fluxes.py @@ -43,9 +43,10 @@ def get_torus_fluxes(dfile): bsq = dfile.GetPm() * 2 sigma = bsq / rho - E = dfile.GetE() - F = dfile.GetF() - P = dfile.GetP() + if dfile.RadiationActive: + E = dfile.GetE() + F = dfile.GetF() + P = dfile.GetP() Nx1 = dfile.MeshBlockSize[0] Nx2 = dfile.MeshBlockSize[1] diff --git a/scripts/python/get_torus_radial_profiles.py b/scripts/python/get_torus_radial_profiles.py new file mode 100644 index 000000000..057378d2c --- /dev/null +++ b/scripts/python/get_torus_radial_profiles.py @@ -0,0 +1,145 @@ +#!/usr/bin/env python + +# © 2023. Triad National Security, LLC. All rights reserved. This +# program was produced under U.S. Government contract +# 89233218CNA000001 for Los Alamos National Laboratory (LANL), which +# is operated by Triad National Security, LLC for the U.S. Department +# of Energy/National Nuclear Security Administration. All rights in +# the program are reserved by Triad National Security, LLC, and the +# U.S. Department of Energy/National Nuclear Security +# Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +# license in this material to reproduce, prepare derivative works, +# distribute copies to the public, perform publicly and display +# publicly, and to permit others to do so. + +import argparse +import numpy as np +import sys +import os +import phoebus_utils +from phoedf import phoedf + +def get_torus_radial_profiles(dfile): + + profiles = {} + + # Geometry + gdet = dfile.gdet + gcov = dfile.gcov + + # Set up size of 1D profiles in memory. Use max and min of global domain and minimum + # dx1 across all blocks. + x1Max = max(np.array(dfile.BlockBounds)[:,1]) + x1Min = min(np.array(dfile.BlockBounds)[:,0]) + dx1_profile = min(dfile.Dx1) + Nx1_profile = round((x1Max - x1Min) / dx1_profile) + + # Simulation state + rho = dfile.GetRho() + ug = dfile.GetUg() + Pg = dfile.GetPg() + ucon = dfile.Getucon() + ucov = dfile.Getucov() + vpcon = dfile.GetVpcon() + bsq = dfile.GetPm() * 2 + sigma = bsq / rho + alpha = dfile.alpha + Gamma = dfile.GetGamma() + + # Initialize profiles + Volume = np.zeros(Nx1_profile) # For volume averages + Mass = np.zeros(Nx1_profile) # For density-weighted volume averages + #Mdot = np.zeros(Nx1_profile) # Total accretion rate + F_M_in = np.zeros(Nx1_profile) # Inflow mass flux + F_M_out = np.zeros(Nx1_profile) # Outflow mass flux + F_Eg = np.zeros(Nx1_profile) # Gas (MHD) energy flux + F_Pg = np.zeros(Nx1_profile) # Gas radial momentum flux + F_Lg = np.zeros(Nx1_profile) # Gas angular momentum flux + Pg_sadw = np.zeros(Nx1_profile) # SADW gas pressure + Pm_sadw = np.zeros(Nx1_profile) # SADW magnetic pressure + JetP_m = np.zeros(Nx1_profile) # MHD jet power (sigma > 1) + JetP_g = np.zeros(Nx1_profile) # hydro jet power (sigma > 1) + beta = np.zeros(Nx1_profile) # Plasma beta + vr = np.zeros(Nx1_profile) # SADW radial velocity + # TODO(BRR) Precalculate some azimuthal averages like for Bernoulli. Store + # these 2D arrays as well? + + for b in range(dfile.NumBlocks): + dx1 = dfile.Dx1[b] + dx2 = dfile.Dx2[b] + dx3 = dfile.Dx3[b] + dA = dx2*dx3*gdet + dV = dx1*dA + + for i in range(dfile.Nx1): + # Get min and max indices for which this value fits + assert dfile.MaxLevel == 0, "This does not support mesh refinement!" + fine_i = 1 # TODO(BRR) placeholder + x1Min_block = dfile.BlockBounds[b][0] + ix1_min = round(x1Min_block / dx1_profile) + + # Bernoulli number to distinguish inflows and outflows + Be = - ((rho[b,:,:,i] - Gamma[b,:,:,i]*ug[b,:,:,i])* ucov[b,0,:,:,i]) / rho[b,:,:,i] - 1. + + for ii in range(fine_i): + ip = i + ix1_min + Volume[ip] += np.sum(dV[b,:,:,i]) + Mass[ip] += np.sum(dV[b,:,:,i]*rho[b,:,:,i]) + + # SADW, rather than / + # TODO(BRR) provide / + Pg_sadw[ip] += np.sum(dV[b,:,:,i]*rho[b,:,:,i] * Pg[b,:,:,i]) + Pm_sadw[ip] += np.sum(dV[b,:,:,i]*rho[b,:,:,i] * bsq[b,:,:,i] / 2.) + + for j in range(dfile.Nx2): + for k in range(dfile.Nx3): + Tmunu_concov = dfile.GetTmunu_concov(b, k, j, i) + F_Eg[ip] += dA[b,k,j,i] * Tmunu_concov[1, 0] + F_Pg[ip] += dA[b,k,j,i] * Tmunu_concov[1, 1] + F_Lg[ip] += dA[b,k,j,i] * Tmunu_concov[1, 3] + if Be[k,j] > 0: + F_M_out[ip] += dA[b,k,j,i]*rho[b,k,j,i]*ucon[b,1,k,j,i] + else: + F_M_in[ip] += dA[b,k,j,i]*rho[b,k,j,i]*ucon[b,1,k,j,i] + if sigma[b,k,j,i] > 1.: + JetP_m[ip] += dA[b,k,j,i]*Tmunu_concov[1,0] + JetP_g[ip] += dA[b,k,j,i]*rho[b,k,j,i]*ucon[b,1,k,j,i] + + vr[ip] += np.sum(dV[b,:,:,i] * vpcon[b,0,:,:,i] / Gamma[b,:,:,i]) + + beta = Pg_sadw / Pm_sadw + + profiles['Volume'] = Volume + profiles['Mass'] = Mass + profiles['F_M'] = F_M_in + F_M_out + profiles['F_M_in'] = F_M_in + profiles['F_M_out'] = F_M_out + profiles['F_Eg'] = F_Eg + profiles['F_Pg'] = F_Pg + profiles['F_Lg'] = F_Lg + profiles['beta'] = beta + profiles['JetP_m'] = JetP_m + profiles['JetP_g'] = JetP_g + profiles['vr'] = vr + + return profiles + +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description="Get radial profiles from torus dump") + parser.add_argument( + "filenames", type=str, nargs="+", help="Files to derive radial profiles for" + ) + args = parser.parse_args() + + print("Phoebus GRMHD analysis script for generating radial profiles from dumps") + print(f" Number of files: {len(args.filenames)}") + + for filename in args.filenames: + print(f"Opening file {os.path.basename(filename)}... ") + dfile = phoedf(filename) + print(f"File {os.path.basename(filename)} opened.") + + profiles = get_torus_radial_profiles(dfile) + diff --git a/scripts/python/grmhd_radial_profiles.py b/scripts/python/grmhd_radial_profiles.py deleted file mode 100644 index de2670722..000000000 --- a/scripts/python/grmhd_radial_profiles.py +++ /dev/null @@ -1,57 +0,0 @@ -#!/usr/bin/env python - -# © 2023. Triad National Security, LLC. All rights reserved. This -# program was produced under U.S. Government contract -# 89233218CNA000001 for Los Alamos National Laboratory (LANL), which -# is operated by Triad National Security, LLC for the U.S. Department -# of Energy/National Nuclear Security Administration. All rights in -# the program are reserved by Triad National Security, LLC, and the -# U.S. Department of Energy/National Nuclear Security -# Administration. The Government is granted for itself and others -# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide -# license in this material to reproduce, prepare derivative works, -# distribute copies to the public, perform publicly and display -# publicly, and to permit others to do so. - -import argparse -import numpy as np -import sys -import os -import phoebus_utils -from phoedf import phoedf - -def get_torus_radial_profiles(dfile): - - profiles = {} - - gdet = dfile.gdet - dx1 = dfile.dx1 - - print(dx1) - - for b in range(dfile.NumBlocks): - dx2 = dfile.Dx2[b] - dx3 = dfile.Dx3[b] - print(f"b: {b} dx2: {dx2} dx3: {dx3})") - - - return profiles - -if __name__ == "__main__": - - parser = argparse.ArgumentParser(description="Get radial profiles from torus dump") - parser.add_argument( - "filenames", type=str, nargs="+", help="Files to derive radial profiles for" - ) - args = parser.parse_args() - - print("Phoebus GRMHD analysis script for generating radial profiles from dumps") - print(f" Number of files: {len(args.filenames)}") - - for filename in args.filenames: - print(f"Opening file {os.path.basename(filename)}... ", end="", flush=True) - dfile = phoedf(filename) - print("done") - - profiles = get_torus_radial_profiles(dfile) - diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py index 6cbb59d32..2623eb857 100644 --- a/scripts/python/phoedf.py +++ b/scripts/python/phoedf.py @@ -62,8 +62,9 @@ def __init__(self, filename, geomfile=None): self.Dx2[b] = (bounds[3] - bounds[2])/self.MeshBlockSize[1] self.Dx3[b] = (bounds[5] - bounds[4])/self.MeshBlockSize[2] - assert self.MaxLevel == 0, "phoedf does not currently support mesh refinement!" + assert self.MaxLevel == 0, "phoedf and its derivatives do not currently support mesh refinement!" + # Common array dimensions self.ScalarField = [self.NumBlocks, self.Nx3, self.Nx2, self.Nx1] self.ThreeVectorField = [self.NumBlocks, 3, self.Nx3, self.Nx2, self.Nx1] self.FourVectorField = [self.NumBlocks, 4, self.Nx3, self.Nx2, self.Nx1] @@ -173,7 +174,9 @@ def flatten_indices(mu, nu): self.tau = None self.Tg = None self.ucon = None + self.ucov = None self.bcon = None + self.bcov = None self.J = None self.Hcov = None self.E = None @@ -253,6 +256,17 @@ def Getbcon(self): return self.bcon + def Getbcov(self): + if self.bcov is None: + self.bcov = np.zeros(self.FourVectorField) + + bcon = self.Getbcon() + + for mu in range(4): + self.bcov[:,:,:,:,:] = self.gcov[:,0,:,:,:,:] * bcon[:,:,:,:,:] + + return self.bcov + def GetEOS(self): return eos_type_dict[self.eos_type](self.Params) @@ -411,6 +425,17 @@ def Getucon(self): return self.ucon + def Getucov(self): + if self.ucov is None: + self.ucov = np.zeros(self.FourVectorField) + + ucon = self.Getucon() + + for mu in range(4): + self.ucov[:,:,:,:,:] = self.gcov[:,0,:,:,:,:] * ucon[:,:,:,:,:] + + return self.ucov + def GetXi(self): if not self.RadiationActive: return None From ccd84a5c2e741b282f0f8e31f9140ceefa913315 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Mon, 11 Dec 2023 14:32:05 -0700 Subject: [PATCH 03/11] More progress --- scripts/python/get_torus_radial_profiles.py | 66 +++++++++++++++++++-- scripts/python/phoedf.py | 6 +- scripts/python/torus_analyzer.py | 2 +- 3 files changed, 66 insertions(+), 8 deletions(-) diff --git a/scripts/python/get_torus_radial_profiles.py b/scripts/python/get_torus_radial_profiles.py index 057378d2c..5d0fa990c 100644 --- a/scripts/python/get_torus_radial_profiles.py +++ b/scripts/python/get_torus_radial_profiles.py @@ -17,8 +17,8 @@ import numpy as np import sys import os -import phoebus_utils from phoedf import phoedf +from multiprocessing import Pool def get_torus_radial_profiles(dfile): @@ -34,6 +34,9 @@ def get_torus_radial_profiles(dfile): x1Min = min(np.array(dfile.BlockBounds)[:,0]) dx1_profile = min(dfile.Dx1) Nx1_profile = round((x1Max - x1Min) / dx1_profile) + r = np.zeros(Nx1_profile) + for i in range(Nx1_profile): + r[i] = x1Min + (0.5 + i) * dx1_profile # Simulation state rho = dfile.GetRho() @@ -110,6 +113,7 @@ def get_torus_radial_profiles(dfile): beta = Pg_sadw / Pm_sadw + profiles['r'] = r profiles['Volume'] = Volume profiles['Mass'] = Mass profiles['F_M'] = F_M_in + F_M_out @@ -125,21 +129,71 @@ def get_torus_radial_profiles(dfile): return profiles +def write_torus_radial_profiles(profiles, in_filename, out_filename): + with open(out_filename, "w") as profile_file: + for key in profiles.keys(): + profile_file.write(key + "\n") + for i in range(len(profiles[key])): + profile_file.write(str(profiles[key][i]) + " ") + profile_file.write("\n") + + return out_filename + +def process_file(filename, overwrite): + + out_filename = filename[:-5] + '.profile' + + if not overwrite and os.path.exists(out_filename): + print(f"{os.path.basename(out_filename)} already exists! Skipping...") + return + + print(f"Opening file {os.path.basename(filename)}... ") + dfile = phoedf(filename) + print(f"File {os.path.basename(filename)} opened.") + + profiles = get_torus_radial_profiles(dfile) + print(f"Created profiles for file {os.path.basename(filename)}") + + out_filename = write_torus_radial_profiles(profiles, in_filename, out_filename) + print(f"Wrote profiles to file {os.path.basename(out_filename)}") + if __name__ == "__main__": parser = argparse.ArgumentParser(description="Get radial profiles from torus dump") parser.add_argument( "filenames", type=str, nargs="+", help="Files to derive radial profiles for" ) + parser.add_argument( + "--nproc", type=int, default=1, help="Number of processes to use" + ) + parser.add_argument( + "--overwrite", action="store_true", help="Whether to overwrite existing outputs" + ) args = parser.parse_args() print("Phoebus GRMHD analysis script for generating radial profiles from dumps") print(f" Number of files: {len(args.filenames)}") - for filename in args.filenames: - print(f"Opening file {os.path.basename(filename)}... ") - dfile = phoedf(filename) - print(f"File {os.path.basename(filename)} opened.") + p = Pool(processes=args.nproc) + from itertools import repeat + p.starmap( + process_file, + zip(args.filenames, repeat(args.overwrite)) + ) + + #with multiprocessing.Pool(processes = args.nproc) as pool: + # results = pool.map(process_file, args.filenames) + + + + #for filename in args.filenames: + # print(f"Opening file {os.path.basename(filename)}... ") + # dfile = phoedf(filename) + # print(f"File {os.path.basename(filename)} opened.") + + # profiles = get_torus_radial_profiles(dfile) + # print(f"Created profiles for file {os.path.basename(filename)}") - profiles = get_torus_radial_profiles(dfile) + # out_filename = write_torus_radial_profiles(profiles, filename) + # print(f"Wrote profiles to file {os.path.basename(out_name)}") diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py index 2623eb857..464d1de42 100644 --- a/scripts/python/phoedf.py +++ b/scripts/python/phoedf.py @@ -398,7 +398,7 @@ def GetGamma(self): return self.Gamma - def GetVpCon(self): + def GetVpcon(self): if self.vpcon is None: self.vpcon = np.clip( self.Get("p.velocity", flatten=False), -1.0e100, 1.0e100 @@ -407,6 +407,10 @@ def GetVpCon(self): return self.vpcon + # Backwards compatibility TODO(BRR) remove? + def GetVpCon(self): + return self.GetVpcon() + def Getucon(self): if self.ucon is None: self.ucon = np.zeros(self.FourVectorField) diff --git a/scripts/python/torus_analyzer.py b/scripts/python/torus_analyzer.py index c7ed4d61e..99aaa23e5 100644 --- a/scripts/python/torus_analyzer.py +++ b/scripts/python/torus_analyzer.py @@ -53,7 +53,7 @@ def process_file(in_filename, fig_folder, log_folder, avg_folder): and os.path.exists(log_filename) and os.path.exists(fig_filename) ): - print(" Output already exists! Skipping...") + print(f"Output for {os.path.basename(in_filename)} already exists! Skipping...") return # Load data From afac5510705462ea819a1ad313f32f5554eaab28 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Mon, 11 Dec 2023 15:09:15 -0700 Subject: [PATCH 04/11] Working --- scripts/python/get_torus_radial_profiles.py | 41 ++++++--------------- 1 file changed, 12 insertions(+), 29 deletions(-) diff --git a/scripts/python/get_torus_radial_profiles.py b/scripts/python/get_torus_radial_profiles.py index 5d0fa990c..f0c153899 100644 --- a/scripts/python/get_torus_radial_profiles.py +++ b/scripts/python/get_torus_radial_profiles.py @@ -34,9 +34,10 @@ def get_torus_radial_profiles(dfile): x1Min = min(np.array(dfile.BlockBounds)[:,0]) dx1_profile = min(dfile.Dx1) Nx1_profile = round((x1Max - x1Min) / dx1_profile) - r = np.zeros(Nx1_profile) + x1 = np.zeros(Nx1_profile) for i in range(Nx1_profile): - r[i] = x1Min + (0.5 + i) * dx1_profile + x1[i] = x1Min + (0.5 + i) * dx1_profile + r = np.exp(x1) # Simulation state rho = dfile.GetRho() @@ -113,6 +114,7 @@ def get_torus_radial_profiles(dfile): beta = Pg_sadw / Pm_sadw + profiles['x1'] = x1 profiles['r'] = r profiles['Volume'] = Volume profiles['Mass'] = Mass @@ -129,7 +131,7 @@ def get_torus_radial_profiles(dfile): return profiles -def write_torus_radial_profiles(profiles, in_filename, out_filename): +def write_torus_radial_profiles(profiles, out_filename): with open(out_filename, "w") as profile_file: for key in profiles.keys(): profile_file.write(key + "\n") @@ -137,24 +139,22 @@ def write_torus_radial_profiles(profiles, in_filename, out_filename): profile_file.write(str(profiles[key][i]) + " ") profile_file.write("\n") - return out_filename +def process_file(in_filename, overwrite): -def process_file(filename, overwrite): - - out_filename = filename[:-5] + '.profile' + out_filename = in_filename[:-5] + '.profile' if not overwrite and os.path.exists(out_filename): print(f"{os.path.basename(out_filename)} already exists! Skipping...") return - print(f"Opening file {os.path.basename(filename)}... ") - dfile = phoedf(filename) - print(f"File {os.path.basename(filename)} opened.") + print(f"Opening file {os.path.basename(in_filename)}... ") + dfile = phoedf(in_filename) + print(f"File {os.path.basename(in_filename)} opened.") profiles = get_torus_radial_profiles(dfile) - print(f"Created profiles for file {os.path.basename(filename)}") + print(f"Created profiles for file {os.path.basename(in_filename)}") - out_filename = write_torus_radial_profiles(profiles, in_filename, out_filename) + write_torus_radial_profiles(profiles, out_filename) print(f"Wrote profiles to file {os.path.basename(out_filename)}") if __name__ == "__main__": @@ -180,20 +180,3 @@ def process_file(filename, overwrite): process_file, zip(args.filenames, repeat(args.overwrite)) ) - - #with multiprocessing.Pool(processes = args.nproc) as pool: - # results = pool.map(process_file, args.filenames) - - - - #for filename in args.filenames: - # print(f"Opening file {os.path.basename(filename)}... ") - # dfile = phoedf(filename) - # print(f"File {os.path.basename(filename)} opened.") - - # profiles = get_torus_radial_profiles(dfile) - # print(f"Created profiles for file {os.path.basename(filename)}") - - # out_filename = write_torus_radial_profiles(profiles, filename) - # print(f"Wrote profiles to file {os.path.basename(out_name)}") - From c84b35535d5f241f3d4e415004c15c7b2c861823 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Tue, 12 Dec 2023 08:27:00 -0700 Subject: [PATCH 05/11] Fix v^r, add v_r --- scripts/python/get_torus_radial_profiles.py | 11 +++++++---- scripts/python/phoedf.py | 17 ++++++++++++++++- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/scripts/python/get_torus_radial_profiles.py b/scripts/python/get_torus_radial_profiles.py index f0c153899..1f887348a 100644 --- a/scripts/python/get_torus_radial_profiles.py +++ b/scripts/python/get_torus_radial_profiles.py @@ -46,6 +46,7 @@ def get_torus_radial_profiles(dfile): ucon = dfile.Getucon() ucov = dfile.Getucov() vpcon = dfile.GetVpcon() + vpcov = dfile.GetVpcov() bsq = dfile.GetPm() * 2 sigma = bsq / rho alpha = dfile.alpha @@ -54,7 +55,6 @@ def get_torus_radial_profiles(dfile): # Initialize profiles Volume = np.zeros(Nx1_profile) # For volume averages Mass = np.zeros(Nx1_profile) # For density-weighted volume averages - #Mdot = np.zeros(Nx1_profile) # Total accretion rate F_M_in = np.zeros(Nx1_profile) # Inflow mass flux F_M_out = np.zeros(Nx1_profile) # Outflow mass flux F_Eg = np.zeros(Nx1_profile) # Gas (MHD) energy flux @@ -65,7 +65,8 @@ def get_torus_radial_profiles(dfile): JetP_m = np.zeros(Nx1_profile) # MHD jet power (sigma > 1) JetP_g = np.zeros(Nx1_profile) # hydro jet power (sigma > 1) beta = np.zeros(Nx1_profile) # Plasma beta - vr = np.zeros(Nx1_profile) # SADW radial velocity + vconr = np.zeros(Nx1_profile) # SADW contravariant radial three-velocity + vcovr = np.zeros(Nx1_profile) # SADW covariant radial three-velocity # TODO(BRR) Precalculate some azimuthal averages like for Bernoulli. Store # these 2D arrays as well? @@ -110,7 +111,8 @@ def get_torus_radial_profiles(dfile): JetP_m[ip] += dA[b,k,j,i]*Tmunu_concov[1,0] JetP_g[ip] += dA[b,k,j,i]*rho[b,k,j,i]*ucon[b,1,k,j,i] - vr[ip] += np.sum(dV[b,:,:,i] * vpcon[b,0,:,:,i] / Gamma[b,:,:,i]) + vconr[ip] += np.sum(dV[b,:,:,i]*rho[b,:,:,i] * vpcon[b,0,:,:,i] / Gamma[b,:,:,i]) + vcovr[ip] += np.sum(dV[b,:,:,i]*rho[b,:,:,i] * vpcov[b,0,:,:,i] / Gamma[b,:,:,i]) beta = Pg_sadw / Pm_sadw @@ -127,7 +129,8 @@ def get_torus_radial_profiles(dfile): profiles['beta'] = beta profiles['JetP_m'] = JetP_m profiles['JetP_g'] = JetP_g - profiles['vr'] = vr + profiles['vconr'] = vconr / Mass + profiles['vcovr'] = vcovr / Mass return profiles diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py index 464d1de42..ca16d71e5 100644 --- a/scripts/python/phoedf.py +++ b/scripts/python/phoedf.py @@ -28,8 +28,11 @@ # ---------------------------------------------------------------------------- # # Phoebus-specific derived class from Parthenon's phdf datafile reader +# +# geomfile : optional geometry file to read many files more efficiently +# clip : Whether to restrict output variables to the range [-1e100, 1e100] class phoedf(phdf.phdf): - def __init__(self, filename, geomfile=None): + def __init__(self, filename, geomfile=None, clip=True): super().__init__(filename) try: @@ -173,6 +176,7 @@ def flatten_indices(mu, nu): self.Pr = None self.tau = None self.Tg = None + self.vpcov = None self.ucon = None self.ucov = None self.bcon = None @@ -398,6 +402,17 @@ def GetGamma(self): return self.Gamma + def GetVpcov(self): + if self.vpcov is None: + self.vpcov = np.zeros(self.ThreeVectorField) + + vpcon = self.GetVpCon() + for ii in range(3): + for jj in range(3): + vpcov[:,ii,:,:,:] += self.gcov[:, ii + 1,jj + 1,:,:,:]*vpcon[:,jj,:,:,:] + + return self.vpcov + def GetVpcon(self): if self.vpcon is None: self.vpcon = np.clip( From bda457b7b7315e8fc49391b24f3602fc713f050b Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Tue, 12 Dec 2023 08:31:47 -0700 Subject: [PATCH 06/11] Update clip option --- scripts/python/phoedf.py | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py index ca16d71e5..04ad4d756 100644 --- a/scripts/python/phoedf.py +++ b/scripts/python/phoedf.py @@ -1,4 +1,4 @@ -# © 2022. Triad National Security, LLC. All rights reserved. This +# © 2022-2023. Triad National Security, LLC. All rights reserved. This # program was produced under U.S. Government contract # 89233218CNA000001 for Los Alamos National Laboratory (LANL), which # is operated by Triad National Security, LLC for the U.S. Department @@ -18,14 +18,11 @@ "../../external/parthenon/scripts/python/packages/parthenon_tools/parthenon_tools/", ) import phdf - # from parthenon_tools import phdf - from phoebus_constants import * from phoebus_eos import * from phoebus_opacities import * - # ---------------------------------------------------------------------------- # # Phoebus-specific derived class from Parthenon's phdf datafile reader # @@ -35,6 +32,8 @@ class phoedf(phdf.phdf): def __init__(self, filename, geomfile=None, clip=True): super().__init__(filename) + self.clip = clip + try: self.eos_type = self.Params["eos/type"].decode() except (UnicodeDecodeError, AttributeError): @@ -305,7 +304,8 @@ def GetTau(self): self.tau[b, :, :, :] = kappaH[b, :, :, :] * dX - self.tau = np.clip(self.tau, 1.0e-100, 1.0e100) + if clip: + self.tau = np.clip(self.tau, 1.0e-100, 1.0e100) return self.tau @@ -326,7 +326,8 @@ def GetPg(self): / self.EnergyDensityCodeToCGS ) - self.Pg = np.clip(self.Pg, 1.0e-100, 1.0e100) + if self.clip: + self.Pg = np.clip(self.Pg, 1.0e-100, 1.0e100) return self.Pg @@ -348,7 +349,8 @@ def GetTg(self): / self.TemperatureCodeToCGS ) - self.Tg = np.clip(self.Tg, 1.0e-100, 1.0e100) + if self.clip: + self.Tg = np.clip(self.Tg, 1.0e-100, 1.0e100) return self.Tg @@ -370,7 +372,8 @@ def GetPm(self): ) self.Pm = bsq / 2.0 - self.Pm = np.clip(self.Pm, 1.0e-100, 1.0e100) + if self.clip: + self.Pm = np.clip(self.Pm, 1.0e-100, 1.0e100) return self.Pm @@ -382,7 +385,8 @@ def GetPr(self): for ispec in range(self.NumSpecies): self.Pr[:, :, :, :] += 1.0 / 3.0 * J[:, ispec, :, :, :] - self.Pr = np.clip(self.Pr, 1.0e-100, 1.0e100) + if self.clip: + self.Pr = np.clip(self.Pr, 1.0e-100, 1.0e100) return self.Pr @@ -415,10 +419,10 @@ def GetVpcov(self): def GetVpcon(self): if self.vpcon is None: - self.vpcon = np.clip( - self.Get("p.velocity", flatten=False), -1.0e100, 1.0e100 - ) + self.vpcon = self.Get("p.velocity", flatten=False) assert self.vpcon is not None + if self.clip: + self.vpcon = np.clip(self.vpcon, -1.0e100, 1.0e100) return self.vpcon @@ -480,7 +484,8 @@ def GetXi(self): self.xi[:, ispec, :, :, :] -= vdH * vdH self.xi = np.sqrt(self.xi) - self.xi = np.clip(self.xi, 1.0e-100, 1.0) + if self.clip: + self.xi = np.clip(self.xi, 1.0e-100, 1.0) return self.xi From d73c044481b09e0de983138b8ba26a9a0bc47208 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Tue, 12 Dec 2023 08:34:26 -0700 Subject: [PATCH 07/11] Fix bad capitalization of GetVpCon, GetVpCov --- scripts/python/get_torus_radial_profiles.py | 4 ++-- scripts/python/phoedf.py | 22 ++++++++++----------- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/scripts/python/get_torus_radial_profiles.py b/scripts/python/get_torus_radial_profiles.py index 1f887348a..d01dccc0e 100644 --- a/scripts/python/get_torus_radial_profiles.py +++ b/scripts/python/get_torus_radial_profiles.py @@ -45,8 +45,8 @@ def get_torus_radial_profiles(dfile): Pg = dfile.GetPg() ucon = dfile.Getucon() ucov = dfile.Getucov() - vpcon = dfile.GetVpcon() - vpcov = dfile.GetVpcov() + vpcon = dfile.Getvpcon() + vpcov = dfile.Getvpcov() bsq = dfile.GetPm() * 2 sigma = bsq / rho alpha = dfile.alpha diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py index 04ad4d756..b8b7dad47 100644 --- a/scripts/python/phoedf.py +++ b/scripts/python/phoedf.py @@ -234,7 +234,7 @@ def Getbcon(self): Bcon = self.GetBcon() Gamma = self.GetGamma() - vcon = self.GetVpCon() / Gamma[:, np.newaxis, :, :, :] + vcon = self.Getvpcon() / Gamma[:, np.newaxis, :, :, :] ucon = self.Getucon() gcov = self.gcov alpha = self.alpha @@ -394,7 +394,7 @@ def GetGamma(self): if self.Gamma is None: self.Gamma = np.zeros(self.ScalarField) - vpcon = self.GetVpCon() + vpcon = self.Getvpcon() for ii in range(3): for jj in range(3): self.Gamma[:, :, :, :] += ( @@ -406,18 +406,18 @@ def GetGamma(self): return self.Gamma - def GetVpcov(self): + def Getvpcov(self): if self.vpcov is None: self.vpcov = np.zeros(self.ThreeVectorField) - vpcon = self.GetVpCon() + vpcon = self.Getvpcon() for ii in range(3): for jj in range(3): vpcov[:,ii,:,:,:] += self.gcov[:, ii + 1,jj + 1,:,:,:]*vpcon[:,jj,:,:,:] return self.vpcov - def GetVpcon(self): + def Getvpcon(self): if self.vpcon is None: self.vpcon = self.Get("p.velocity", flatten=False) assert self.vpcon is not None @@ -428,13 +428,13 @@ def GetVpcon(self): # Backwards compatibility TODO(BRR) remove? def GetVpCon(self): - return self.GetVpcon() + return self.Getvpcon() def Getucon(self): if self.ucon is None: self.ucon = np.zeros(self.FourVectorField) - vpcon = self.GetVpCon() + vpcon = self.Getvpcon() Gamma = self.GetGamma() self.ucon[:, 0, :, :, :] = Gamma[:, :, :, :] / self.alpha[:, :, :, :] @@ -470,7 +470,7 @@ def GetXi(self): Hcov = self.GetHcov() / self.GetJ()[:, np.newaxis, :, :, :] Gamma = self.GetGamma() - vcon = self.GetVpCon() / Gamma[:, np.newaxis, :, :, :] + vcon = self.Getvpcon() / Gamma[:, np.newaxis, :, :, :] for ispec in range(self.NumSpecies): vdH = np.zeros(self.ScalarField) for ii in range(3): @@ -499,7 +499,7 @@ def GetE(self): ) Gamma = self.GetGamma() - vcon = self.GetVpCon() / Gamma[:, np.newaxis, :, :, :] + vcon = self.Getvpcon() / Gamma[:, np.newaxis, :, :, :] J = self.GetJ() Hcov = self.GetHcov() @@ -529,7 +529,7 @@ def GetF(self): ) Gamma = self.GetGamma() - vcon = self.GetVpCon() / Gamma[:, np.newaxis, :, :, :] + vcon = self.Getvpcon() / Gamma[:, np.newaxis, :, :, :] J = self.GetJ() Hcov = self.GetHcov() gammacon = self.gammacon @@ -572,7 +572,7 @@ def GetP(self): ) Gamma = self.GetGamma() - vcon = self.GetVpCon() / Gamma[:, np.newaxis, :, :, :] + vcon = self.Getvpcon() / Gamma[:, np.newaxis, :, :, :] J = self.GetJ() Hcov = self.GetHcov() gammacon = self.gammacon From cb1758313991bba78e655814ea6488a571b8b2b8 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Tue, 12 Dec 2023 08:35:36 -0700 Subject: [PATCH 08/11] Better banner for get_torus_radial_profiles.py run as main --- scripts/python/get_torus_radial_profiles.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/scripts/python/get_torus_radial_profiles.py b/scripts/python/get_torus_radial_profiles.py index d01dccc0e..affcd0b61 100644 --- a/scripts/python/get_torus_radial_profiles.py +++ b/scripts/python/get_torus_radial_profiles.py @@ -175,7 +175,9 @@ def process_file(in_filename, overwrite): args = parser.parse_args() print("Phoebus GRMHD analysis script for generating radial profiles from dumps") - print(f" Number of files: {len(args.filenames)}") + print(f" Number of files: {len(args.filenames)}") + print(f" Number of processors: {args.nproc}") + print(f" Overwrite? {args.overwrite}") p = Pool(processes=args.nproc) from itertools import repeat From 6ad19c9e7fb4926b88c75898ef33dbb4ff4bf3d2 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Wed, 13 Dec 2023 13:28:22 -0700 Subject: [PATCH 09/11] Add time averaging script --- scripts/python/get_torus_radial_profiles.py | 1 + scripts/python/phoedf.py | 2 +- scripts/python/tavg_torus_radial_profiles.py | 128 +++++++++++++++++++ 3 files changed, 130 insertions(+), 1 deletion(-) create mode 100644 scripts/python/tavg_torus_radial_profiles.py diff --git a/scripts/python/get_torus_radial_profiles.py b/scripts/python/get_torus_radial_profiles.py index affcd0b61..ce306ed7c 100644 --- a/scripts/python/get_torus_radial_profiles.py +++ b/scripts/python/get_torus_radial_profiles.py @@ -116,6 +116,7 @@ def get_torus_radial_profiles(dfile): beta = Pg_sadw / Pm_sadw + # TODO(BRR) Include dump time profiles['x1'] = x1 profiles['r'] = r profiles['Volume'] = Volume diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py index b8b7dad47..493a360c6 100644 --- a/scripts/python/phoedf.py +++ b/scripts/python/phoedf.py @@ -413,7 +413,7 @@ def Getvpcov(self): vpcon = self.Getvpcon() for ii in range(3): for jj in range(3): - vpcov[:,ii,:,:,:] += self.gcov[:, ii + 1,jj + 1,:,:,:]*vpcon[:,jj,:,:,:] + self.vpcov[:,ii,:,:,:] += self.gcov[:, ii + 1,jj + 1,:,:,:]*vpcon[:,jj,:,:,:] return self.vpcov diff --git a/scripts/python/tavg_torus_radial_profiles.py b/scripts/python/tavg_torus_radial_profiles.py new file mode 100644 index 000000000..7be3af438 --- /dev/null +++ b/scripts/python/tavg_torus_radial_profiles.py @@ -0,0 +1,128 @@ +#!/usr/bin/env python + +# © 2023. Triad National Security, LLC. All rights reserved. This +# program was produced under U.S. Government contract +# 89233218CNA000001 for Los Alamos National Laboratory (LANL), which +# is operated by Triad National Security, LLC for the U.S. Department +# of Energy/National Nuclear Security Administration. All rights in +# the program are reserved by Triad National Security, LLC, and the +# U.S. Department of Energy/National Nuclear Security +# Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +# license in this material to reproduce, prepare derivative works, +# distribute copies to the public, perform publicly and display +# publicly, and to permit others to do so. + +import argparse +import numpy as np +import sys +import os +from phoedf import phoedf +from multiprocessing import Pool + +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description="Time-average radial profiles from torus dumps") + parser.add_argument( + "filenames", type=str, nargs="+", help="Files to derive radial profiles for" + ) + parser.add_argument( + "--overwrite", action="store_true", help="Whether to overwrite existing outputs" + ) + parser.add_argument( + "--tmin", type=float, default=0, help="Minimum time to begin averaging" + ) + parser.add_argument( + "--dt", type=float, required=True, help="Time window over which to separate averages" + ) + parser.add_argument( + "--dt_sim", type=float, required=True, help="Dump cadence" + ) + args = parser.parse_args() + + print("Phoebus GRMHD analysis script for time-averaging radial profiles from dumps") + print(f" Number of files: {len(args.filenames)}") + print(f" Overwrite? {args.overwrite}") + + tavgs = {} + + keys_to_tavg = ['r', 'F_M', 'F_M_in', 'F_M_out', 'F_Eg', 'F_Pg', 'F_Lg', 'beta', 'vconr', 'vcovr'] + + for n, filename in enumerate(args.filenames): + base_directory = os.path.dirname(filename) + base_filename = os.path.basename(filename) + print(f"Processing file {base_filename}") + + # Overly specific way to get time + nfile = int(base_filename[11:19]) + tfile = args.dt_sim * nfile + + avg_idx = int(tfile / args.dt) + avg_key = str(avg_idx) + + if avg_key not in tavgs.keys(): + tavgs[avg_key] = {} + tavgs[avg_key]['nfiles'] = 0 + tavgs[avg_key]['tmin'] = avg_idx * args.dt + tavgs[avg_key]['dt'] = args.dt + + data_indices = {} + + with open(filename, "r") as profile: + lines = profile.readlines() + for idx, line in enumerate(lines): + line = line.strip() + for key in keys_to_tavg: + if line == key: + data_indices[key] = idx + 1 + + # Initialize data if necessary + ndata = len(np.fromstring(lines[data_indices[keys_to_tavg[0]]].strip(), dtype=float, sep=' ')) + if tavgs[avg_key]['nfiles'] == 0: + for key in keys_to_tavg: + tavgs[avg_key][key] = np.zeros(ndata) + + # Sum data + for key in keys_to_tavg: + tavgs[avg_key][key] += np.fromstring(lines[data_indices[key]].strip(), dtype=float, sep=' ') + + # Record number of files in this bin for averaging + tavgs[avg_key]['nfiles'] += 1 + + # Normalize + for avg_key in tavgs.keys(): + for key in keys_to_tavg: + tavgs[avg_key][key] /= tavgs[avg_key]['nfiles'] + + # Save dictionary as pickle + import pickle + with open(os.path.join(base_directory, 'tavgs.pickle'), 'wb') as handle: + pickle.dump(tavgs, handle, protocol=pickle.HIGHEST_PROTOCOL) + + # Save dictionary as plaintext (json) + import json + import copy + tavgs_jsonable = copy.deepcopy(tavgs) + # Convert NumPy arrays to lists in the copy + def convert_numpys_to_lists(d): + for key, value in d.items(): + if isinstance(value, np.ndarray): + d[key] = value.tolist() + elif isinstance(value, dict): + convert_numpys_to_lists(value) + convert_numpys_to_lists(tavgs_jsonable) + with open(os.path.join(base_directory, 'tavgs.json'), 'w') as handle: + json.dump(tavgs_jsonable, handle, indent=4) + + # Visualization + import matplotlib.pyplot as plt + fig, ax = plt.subplots(1,1) + for avg_key in tavgs.keys(): + ax.plot(tavgs[avg_key]['r'], tavgs[avg_key]['F_M']) + ax.set_xlabel('r') + ax.set_ylabel('FM') + ax.set_xscale('log') + plt.savefig("F_M.png") + + + From 89c23eef82bbc6bb8770046e4e4db23564d563f2 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Wed, 13 Dec 2023 14:34:00 -0700 Subject: [PATCH 10/11] Updating profile scripts now with plotting --- .../python/plot_tavg_torus_radial_profiles.py | 73 +++++++++++++++++++ scripts/python/tavg_torus_radial_profiles.py | 10 ++- 2 files changed, 80 insertions(+), 3 deletions(-) create mode 100644 scripts/python/plot_tavg_torus_radial_profiles.py diff --git a/scripts/python/plot_tavg_torus_radial_profiles.py b/scripts/python/plot_tavg_torus_radial_profiles.py new file mode 100644 index 000000000..fc589e75e --- /dev/null +++ b/scripts/python/plot_tavg_torus_radial_profiles.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python + +# © 2023. Triad National Security, LLC. All rights reserved. This +# program was produced under U.S. Government contract +# 89233218CNA000001 for Los Alamos National Laboratory (LANL), which +# is operated by Triad National Security, LLC for the U.S. Department +# of Energy/National Nuclear Security Administration. All rights in +# the program are reserved by Triad National Security, LLC, and the +# U.S. Department of Energy/National Nuclear Security +# Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +# license in this material to reproduce, prepare derivative works, +# distribute copies to the public, perform publicly and display +# publicly, and to permit others to do so. + +# This file was generated in part by OpenAI's GPT-4. + +import argparse +import numpy as np +import sys +import os +import pickle +import matplotlib.pyplot as plt +from matplotlib import cm + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description="Time-average radial profiles from torus dumps") + parser.add_argument( + "filename", type=str, help="Pickle file containing time-averaged radial profiles (from tavg_torus_radial_profiles.py)" + ) + parser.add_argument( + "key", type=str, help="Profile to plot" + ) + parser.add_argument( + "--rmin", type=float, default=None, help="Minimum radial coordinate to plot" + ) + parser.add_argument( + "--rmax", type=float, default=None, help="Maximum radial coordinate to plot" + ) + parser.add_argument( + "--ymin", type=float, default=None, help="Minimum value to enforce on data" + ) + parser.add_argument( + "--ymax", type=float, default=None, help="Maximum value to enforce on data" + ) + parser.add_argument( + "--log", action="store_true", help="Whether to use logarithmic y axis" + ) + args = parser.parse_args() + + with open(args.filename, 'rb') as handle: + tavgs = pickle.load(handle) + + fig, ax = plt.subplots(1,1) + + for n, avg_key in enumerate(tavgs.keys()): + col = cm.get_cmap("Spectral")(n / len(tavgs.keys())) + ax.plot(tavgs[avg_key]['r'], tavgs[avg_key][args.key], label=str(tavgs[avg_key]['tmin']), color=col) + if args.log: + ax.plot(tavgs[avg_key]['r'], -tavgs[avg_key][args.key], linestyle='--', color=col) + ax.legend(loc=2) + + ax.set_xlabel('r') + ax.set_ylabel(args.key) + ax.set_xscale('log') + + ax.set_xlim([args.rmin, args.rmax]) + ax.set_ylim([args.ymin, args.ymax]) + + if args.log: + ax.set_yscale('log') + + plt.savefig("tavg_profile.png") \ No newline at end of file diff --git a/scripts/python/tavg_torus_radial_profiles.py b/scripts/python/tavg_torus_radial_profiles.py index 7be3af438..c0b2aec21 100644 --- a/scripts/python/tavg_torus_radial_profiles.py +++ b/scripts/python/tavg_torus_radial_profiles.py @@ -13,12 +13,12 @@ # distribute copies to the public, perform publicly and display # publicly, and to permit others to do so. +# This file was generated in part by OpenAI's GPT-4. + import argparse import numpy as np import sys import os -from phoedf import phoedf -from multiprocessing import Pool if __name__ == "__main__": @@ -54,7 +54,11 @@ print(f"Processing file {base_filename}") # Overly specific way to get time - nfile = int(base_filename[11:19]) + try: + nfile = int(base_filename[11:19]) + except ValueError: + print(f" Can't extract dumpfile index; skipping!") + continue tfile = args.dt_sim * nfile avg_idx = int(tfile / args.dt) From bc580853f8023c6df0fbf74ed9e902e0d077d813 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Mon, 18 Dec 2023 13:30:16 -0700 Subject: [PATCH 11/11] Cusotm tavg windows, more vals --- scripts/python/get_torus_radial_profiles.py | 24 ++++ .../python/plot_tavg_torus_radial_profiles.py | 7 +- scripts/python/tavg_torus_radial_profiles.py | 120 +++++++++++++----- 3 files changed, 120 insertions(+), 31 deletions(-) diff --git a/scripts/python/get_torus_radial_profiles.py b/scripts/python/get_torus_radial_profiles.py index ce306ed7c..e60d43472 100644 --- a/scripts/python/get_torus_radial_profiles.py +++ b/scripts/python/get_torus_radial_profiles.py @@ -58,8 +58,14 @@ def get_torus_radial_profiles(dfile): F_M_in = np.zeros(Nx1_profile) # Inflow mass flux F_M_out = np.zeros(Nx1_profile) # Outflow mass flux F_Eg = np.zeros(Nx1_profile) # Gas (MHD) energy flux + F_Eg_out = np.zeros(Nx1_profile) # Outward MHD energy flux + F_Eg_in = np.zeros(Nx1_profile) # Inward MHD energy flux F_Pg = np.zeros(Nx1_profile) # Gas radial momentum flux + F_Pg_out = np.zeros(Nx1_profile) # Outward gas radial momentum flux + F_Pg_in = np.zeros(Nx1_profile) # Inward gas radial momentum flux F_Lg = np.zeros(Nx1_profile) # Gas angular momentum flux + F_Lg_out = np.zeros(Nx1_profile) # Outward gas angular momentum flux + F_Lg_in = np.zeros(Nx1_profile) # Inward gas angular momentum flux Pg_sadw = np.zeros(Nx1_profile) # SADW gas pressure Pm_sadw = np.zeros(Nx1_profile) # SADW magnetic pressure JetP_m = np.zeros(Nx1_profile) # MHD jet power (sigma > 1) @@ -67,6 +73,8 @@ def get_torus_radial_profiles(dfile): beta = np.zeros(Nx1_profile) # Plasma beta vconr = np.zeros(Nx1_profile) # SADW contravariant radial three-velocity vcovr = np.zeros(Nx1_profile) # SADW covariant radial three-velocity + vconr_out = np.zeros(Nx1_profile) # SADW v^r for outflow + vconr_in = np.zeros(Nx1_profile) # SADW v^r for inflow # TODO(BRR) Precalculate some azimuthal averages like for Bernoulli. Store # these 2D arrays as well? @@ -105,8 +113,16 @@ def get_torus_radial_profiles(dfile): F_Lg[ip] += dA[b,k,j,i] * Tmunu_concov[1, 3] if Be[k,j] > 0: F_M_out[ip] += dA[b,k,j,i]*rho[b,k,j,i]*ucon[b,1,k,j,i] + F_Eg_out[ip] += dA[b,k,j,i] * Tmunu_concov[1, 0] + F_Pg_out[ip] += dA[b,k,j,i] * Tmunu_concov[1, 1] + F_Lg_out[ip] += dA[b,k,j,i] * Tmunu_concov[1, 3] + vconr_out[ip] += dV[b,k,j,i]*rho[b,k,j,i]*vpcon[b,0,k,j,i]/Gamma[b,k,j,i] else: F_M_in[ip] += dA[b,k,j,i]*rho[b,k,j,i]*ucon[b,1,k,j,i] + F_Eg_in[ip] += dA[b,k,j,i] * Tmunu_concov[1, 0] + F_Pg_in[ip] += dA[b,k,j,i] * Tmunu_concov[1, 1] + F_Lg_in[ip] += dA[b,k,j,i] * Tmunu_concov[1, 3] + vconr_in[ip] += dV[b,k,j,i]*rho[b,k,j,i]*vpcon[b,0,k,j,i]/Gamma[b,k,j,i] if sigma[b,k,j,i] > 1.: JetP_m[ip] += dA[b,k,j,i]*Tmunu_concov[1,0] JetP_g[ip] += dA[b,k,j,i]*rho[b,k,j,i]*ucon[b,1,k,j,i] @@ -127,11 +143,19 @@ def get_torus_radial_profiles(dfile): profiles['F_Eg'] = F_Eg profiles['F_Pg'] = F_Pg profiles['F_Lg'] = F_Lg + profiles['F_Eg_out'] = F_Eg_out + profiles['F_Pg_out'] = F_Pg_out + profiles['F_Lg_out'] = F_Lg_out + profiles['F_Eg_in'] = F_Eg_in + profiles['F_Pg_in'] = F_Pg_in + profiles['F_Lg_in'] = F_Lg_in profiles['beta'] = beta profiles['JetP_m'] = JetP_m profiles['JetP_g'] = JetP_g profiles['vconr'] = vconr / Mass profiles['vcovr'] = vcovr / Mass + profiles['vconr_out'] = vconr_out / Mass + profiles['vconr_in'] = vconr_in / Mass return profiles diff --git a/scripts/python/plot_tavg_torus_radial_profiles.py b/scripts/python/plot_tavg_torus_radial_profiles.py index fc589e75e..56656cad5 100644 --- a/scripts/python/plot_tavg_torus_radial_profiles.py +++ b/scripts/python/plot_tavg_torus_radial_profiles.py @@ -46,6 +46,7 @@ parser.add_argument( "--log", action="store_true", help="Whether to use logarithmic y axis" ) + # TODO(BRR) add geometric scaling with time averaging window args = parser.parse_args() with open(args.filename, 'rb') as handle: @@ -54,8 +55,10 @@ fig, ax = plt.subplots(1,1) for n, avg_key in enumerate(tavgs.keys()): + if n == len(tavgs.keys()) - 1: + continue col = cm.get_cmap("Spectral")(n / len(tavgs.keys())) - ax.plot(tavgs[avg_key]['r'], tavgs[avg_key][args.key], label=str(tavgs[avg_key]['tmin']), color=col) + #ax.plot(tavgs[avg_key]['r'], tavgs[avg_key][args.key], label=str(tavgs[avg_key]['tmin']), color=col) if args.log: ax.plot(tavgs[avg_key]['r'], -tavgs[avg_key][args.key], linestyle='--', color=col) ax.legend(loc=2) @@ -70,4 +73,4 @@ if args.log: ax.set_yscale('log') - plt.savefig("tavg_profile.png") \ No newline at end of file + plt.savefig("tavg_profile.png") diff --git a/scripts/python/tavg_torus_radial_profiles.py b/scripts/python/tavg_torus_radial_profiles.py index c0b2aec21..9212742c9 100644 --- a/scripts/python/tavg_torus_radial_profiles.py +++ b/scripts/python/tavg_torus_radial_profiles.py @@ -22,7 +22,9 @@ if __name__ == "__main__": - parser = argparse.ArgumentParser(description="Time-average radial profiles from torus dumps") + parser = argparse.ArgumentParser( + description="Time-average radial profiles from torus dumps" + ) parser.add_argument( "filenames", type=str, nargs="+", help="Files to derive radial profiles for" ) @@ -33,20 +35,53 @@ "--tmin", type=float, default=0, help="Minimum time to begin averaging" ) parser.add_argument( - "--dt", type=float, required=True, help="Time window over which to separate averages" + "--dt", + type=float, + required=True, + help="Time window over which to separate averages", ) + parser.add_argument("--dt_sim", type=float, required=True, help="Dump cadence") parser.add_argument( - "--dt_sim", type=float, required=True, help="Dump cadence" + "--tbins", type=float, nargs="+", help="Time-averaging window bin edges to use" ) args = parser.parse_args() print("Phoebus GRMHD analysis script for time-averaging radial profiles from dumps") print(f" Number of files: {len(args.filenames)}") print(f" Overwrite? {args.overwrite}") + print(f" Time averaging bins: {args.tbins}") - tavgs = {} + using_custom_tavg = args.tbins is not None - keys_to_tavg = ['r', 'F_M', 'F_M_in', 'F_M_out', 'F_Eg', 'F_Pg', 'F_Lg', 'beta', 'vconr', 'vcovr'] + tavgs = {} + custom_tavgs = {} + + def get_custom_avg_idx(t): + for nedge in range(len(args.tbins) - 1): + if t >= args.tbins[nedge] and t < args.tbins[nedge + 1]: + return nedge + return -1 + + keys_to_tavg = [ + "r", + "F_M", + "F_M_in", + "F_M_out", + "F_Eg", + "F_Eg_out", + "F_Eg_in", + "F_Pg", + "F_Pg_out", + "F_Pg_in", + "F_Lg", + "F_Lg_out", + "F_Lg_in", + "beta", + "vconr", + "vcovr", + "vconr_out", + "vconr_in" + ] for n, filename in enumerate(args.filenames): base_directory = os.path.dirname(filename) @@ -64,11 +99,21 @@ avg_idx = int(tfile / args.dt) avg_key = str(avg_idx) + if using_custom_tavg: + custom_avg_idx = get_custom_avg_idx(tfile) + custom_avg_key = str(custom_avg_idx) + if avg_key not in tavgs.keys(): tavgs[avg_key] = {} - tavgs[avg_key]['nfiles'] = 0 - tavgs[avg_key]['tmin'] = avg_idx * args.dt - tavgs[avg_key]['dt'] = args.dt + tavgs[avg_key]["nfiles"] = 0 + tavgs[avg_key]["tmin"] = avg_idx * args.dt + tavgs[avg_key]["dt"] = args.dt + + if using_custom_tavg and custom_avg_key not in custom_tavgs.keys(): + custom_tavgs[custom_avg_key] = {} + custom_tavgs[custom_avg_key]["nfiles"] = 0 + custom_tavgs[custom_avg_key]["tmin"] = avg_idx * args.dt + custom_tavgs[custom_avg_key]["dt"] = args.dt data_indices = {} @@ -81,32 +126,55 @@ data_indices[key] = idx + 1 # Initialize data if necessary - ndata = len(np.fromstring(lines[data_indices[keys_to_tavg[0]]].strip(), dtype=float, sep=' ')) - if tavgs[avg_key]['nfiles'] == 0: + ndata = len( + np.fromstring( + lines[data_indices[keys_to_tavg[0]]].strip(), dtype=float, sep=" " + ) + ) + if tavgs[avg_key]["nfiles"] == 0: for key in keys_to_tavg: tavgs[avg_key][key] = np.zeros(ndata) + if using_custom_tavg and custom_tavgs[custom_avg_key]["nfiles"] == 0: + for key in keys_to_tavg: + custom_tavgs[custom_avg_key][key] = np.zeros(ndata) # Sum data for key in keys_to_tavg: - tavgs[avg_key][key] += np.fromstring(lines[data_indices[key]].strip(), dtype=float, sep=' ') + tavgs[avg_key][key] += np.fromstring( + lines[data_indices[key]].strip(), dtype=float, sep=" " + ) + if using_custom_tavg: + custom_tavgs[custom_avg_key][key] += np.fromstring( + lines[data_indices[key]].strip(), dtype=float, sep=" " + ) # Record number of files in this bin for averaging - tavgs[avg_key]['nfiles'] += 1 + tavgs[avg_key]["nfiles"] += 1 # Normalize for avg_key in tavgs.keys(): for key in keys_to_tavg: - tavgs[avg_key][key] /= tavgs[avg_key]['nfiles'] + tavgs[avg_key][key] /= tavgs[avg_key]["nfiles"] + if using_custom_tavg: + for custom_avg_key in custom_tavgs.keys(): + for key in keys_to_tavg: + custom_tavgs[custom_avg_key][key] /= custom_tavgs[custom_avg_key][ + "nfiles" + ] # Save dictionary as pickle import pickle - with open(os.path.join(base_directory, 'tavgs.pickle'), 'wb') as handle: + + with open(os.path.join(base_directory, "tavgs.pickle"), "wb") as handle: pickle.dump(tavgs, handle, protocol=pickle.HIGHEST_PROTOCOL) + if using_custom_tavg: + with open(os.path.join(base_directory, "custom_tavgs.pickle"), "wb") as handle: + pickle.dump(custom_tavgs, handle, protocol=pickle.HIGHEST_PROTOCOL) # Save dictionary as plaintext (json) import json import copy - tavgs_jsonable = copy.deepcopy(tavgs) + # Convert NumPy arrays to lists in the copy def convert_numpys_to_lists(d): for key, value in d.items(): @@ -114,19 +182,13 @@ def convert_numpys_to_lists(d): d[key] = value.tolist() elif isinstance(value, dict): convert_numpys_to_lists(value) + + tavgs_jsonable = copy.deepcopy(tavgs) convert_numpys_to_lists(tavgs_jsonable) - with open(os.path.join(base_directory, 'tavgs.json'), 'w') as handle: + with open(os.path.join(base_directory, "tavgs.json"), "w") as handle: json.dump(tavgs_jsonable, handle, indent=4) - - # Visualization - import matplotlib.pyplot as plt - fig, ax = plt.subplots(1,1) - for avg_key in tavgs.keys(): - ax.plot(tavgs[avg_key]['r'], tavgs[avg_key]['F_M']) - ax.set_xlabel('r') - ax.set_ylabel('FM') - ax.set_xscale('log') - plt.savefig("F_M.png") - - - + if using_custom_tavg: + custom_tavgs_jsonable = copy.deepcopy(custom_tavgs) + convert_numpys_to_lists(custom_tavgs_jsonable) + with open(os.path.join(base_directory, "custom_tavgs.json"), "w") as handle: + json.dump(custom_tavgs_jsonable, handle, indent=4)