From f1b93be753207f01c3f73b3476a7bf7903b09a73 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 27 Aug 2023 10:08:29 -0400 Subject: [PATCH] update the subchandra plots to those used in the paper (#2543) --- .../subchandra/analysis/subch_res_compare.py | 168 ++++++++++++++++++ .../subchandra/analysis/subch_sequence.py | 45 +++-- .../science/subchandra/analysis/subch_zoom.py | 52 ++++++ 3 files changed, 247 insertions(+), 18 deletions(-) create mode 100755 Exec/science/subchandra/analysis/subch_res_compare.py create mode 100755 Exec/science/subchandra/analysis/subch_zoom.py diff --git a/Exec/science/subchandra/analysis/subch_res_compare.py b/Exec/science/subchandra/analysis/subch_res_compare.py new file mode 100755 index 0000000000..e90b70f863 --- /dev/null +++ b/Exec/science/subchandra/analysis/subch_res_compare.py @@ -0,0 +1,168 @@ +#!/usr/bin/env python3 + +import argparse +import os +import sys +from functools import reduce + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from mpl_toolkits.axes_grid1 import ImageGrid + +import yt +from yt.fields.derived_field import ValidateSpatial +from yt.frontends.boxlib.api import CastroDataset +from yt.funcs import just_one +# assume that our data is in CGS +from yt.units import amu, cm + +matplotlib.use('agg') + + +clip_val = -35 +max_val = -19 + +# how much to coarsen for the contouring +blocking_factor = 8 + +def _lap_rho(field, data): + dr = just_one(data["index", "dr"]).d + r = data["index", "r"].d + rl = r - 0.5 * dr + rr = r + 0.5 * dr + + dz = just_one(data["index", "dz"]).d + dens = data["gas", "density"].d + + _lap = np.zeros_like(dens) + + lapl_field = data.ds.arr(np.zeros(dens.shape, dtype=np.float64), None) + + # r-component + _lap[1:-1, :] = 1 / (r[1:-1, :] * dr**2) * ( + - 2.0 * r[1:-1, :] * dens[1:-1:, :] + + rl[1:-1, :] * dens[:-2, :] + rr[1:-1, :] * dens[2:, :]) + + _lap[:, 1:-1] += 1 / dz**2 * (dens[:, 2:] + dens[:, :-2] - 2.0 * dens[:, 1:-1]) + lapl_field[1:-1, 1:-1] = np.log(np.abs(_lap[1:-1, 1:-1] / dens[1:-1, 1:-1])) + lapl_field[lapl_field < clip_val] = clip_val + return lapl_field + +def doit(rows, field): + + nrows = len(rows) + ncols = len(rows[0][0]) + + + fig = plt.figure() + + grid = ImageGrid(fig, 111, nrows_ncols=(nrows, ncols), + axes_pad=0.25, cbar_pad=0.05, label_mode="L", cbar_mode="single") + + + i = 0 + + for row in rows: + + plotfiles, label = row + + for irow, pf in enumerate(plotfiles): + + ds = CastroDataset(pf) + + if field == "lap_rho": + ds.force_periodicity() + ds.add_field(name=("gas", "lap_rho"), sampling_type="local", + function=_lap_rho, units="", + validators=[ValidateSpatial(1)]) + + domain_frac = 0.2 + + xmin = ds.domain_left_edge[0] + xmax = domain_frac * ds.domain_right_edge[0] + xctr = 0.5 * (xmin + xmax) + L_x = xmax - xmin + + ymin = ds.domain_left_edge[1] + ymax = ds.domain_right_edge[1] + yctr = 0.5 * (ymin + ymax) + L_y = ymax - ymin + ymin = yctr - 0.5 * domain_frac * L_y + ymax = yctr + 0.5 * domain_frac * L_y + L_y = ymax - ymin + + sp = yt.SlicePlot(ds, "theta", field, center=[xctr, yctr, 0.0*cm], width=[L_x, L_y, 0.0*cm], fontsize="14") + sp.set_buff_size((2400,2400)) + + if field == "Temp": + text_color = "white" + else: + text_color = "black" + + sp.annotate_text((0.05, 0.05), f"time = {float(ds.current_time):8.3f} s", coord_system="axis", text_args={"color": text_color, "fontsize": "12"}) + if (irow == 0): + sp.annotate_text((0.05, 0.925), f"{label}", coord_system="axis", text_args={"color": text_color, "fontsize": "14"}) + + if (irow == 0): + sp.annotate_grids(max_level=10, cmap="tab10") + + if field == "Temp": + sp.set_zlim(field, 5.e7, 4e9) + sp.set_cmap(field, "magma") + elif field == "enuc": + sp.set_log(field, True, linthresh=1.e15) + sp.set_zlim(field, -1.e22, 1.e22) + sp.set_cmap(field, "bwr") + elif field == "abar": + sp.set_zlim(field, 4, 28) + sp.set_log(field, False) + sp.set_cmap(field, "plasma_r") + elif field == "lap_rho": + sp.set_zlim(field, clip_val, max_val) + sp.set_log(field, False) + sp.set_cmap(field, "bone_r") + + sp.set_axes_unit("km") + + plot = sp.plots[field] + plot.figure = fig + plot.axes = grid[i].axes + plot.cax = grid.cbar_axes[i] + if irow < len(plotfiles)-1: + grid[i].axes.xaxis.offsetText.set_visible(False) + + sp._setup_plots() + + i += 1 + + fig.set_size_inches(10, 14) + plt.tight_layout() + plt.savefig(f"subch_{field}_res_compare.pdf") + +if __name__ == "__main__": + + + subch_40km = [("subch_sdc_40km/subch_plt02123", + "subch_sdc_40km/subch_plt04184", + "subch_sdc_40km/subch_plt08509", + "subch_sdc_40km/subch_plt17272"), "40 km"] + + subch_20km = [("subch_sdc/subch_plt02123", + "subch_sdc/subch_plt04196", + "subch_sdc/subch_plt08614", + "subch_sdc/subch_plt17412"), "20 km"] + + subch_10km = [("subch_sdc_10km_3lev/subch_plt02123", + "subch_sdc_10km_3lev/subch_plt04197", + "subch_sdc_10km_3lev/subch_plt08582", + "subch_sdc_10km_3lev/subch_plt17526"), "10 km"] + + subch_5km = [("subch_sdc_5km_4lev/subch_plt02123", + "subch_sdc_5km_4lev/subch_plt04197", + "subch_sdc_5km_4lev/subch_plt08581", + "subch_sdc_5km_4lev/subch_plt17472"), "5 km"] + + field = "Temp" + + doit([subch_40km, subch_20km, subch_10km, subch_5km], field) diff --git a/Exec/science/subchandra/analysis/subch_sequence.py b/Exec/science/subchandra/analysis/subch_sequence.py index 075c79018d..e6f6fad56c 100755 --- a/Exec/science/subchandra/analysis/subch_sequence.py +++ b/Exec/science/subchandra/analysis/subch_sequence.py @@ -1,26 +1,27 @@ #!/usr/bin/env python3 -import matplotlib -matplotlib.use('agg') - import argparse import os import sys -import yt -import matplotlib.pyplot as plt -import numpy as np from functools import reduce +import matplotlib +import matplotlib.pyplot as plt +import numpy as np from mpl_toolkits.axes_grid1 import ImageGrid -# assume that our data is in CGS -from yt.units import cm, amu +import yt +from yt.fields.derived_field import ValidateSpatial from yt.frontends.boxlib.api import CastroDataset from yt.funcs import just_one -from yt.fields.derived_field import ValidateSpatial +# assume that our data is in CGS +from yt.units import amu, cm + +matplotlib.use('agg') + #times = [0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35] -times = [0.0, 0.1, 0.2, 0.4, 0.6, 0.75] +times = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1] #times = [0.0, 0.15, 0.3, 0.45] clip_val = -35 @@ -75,7 +76,10 @@ def doit(field, add_contours, pfiles): fig = plt.figure() - if len(pfiles) > 4: + if len(pfiles) > 8: + nrows = 3 + ncols = (len(pfiles) + 1)//3 + elif len(pfiles) > 4: nrows = 2 ncols = (len(pfiles) + 1)//2 else: @@ -83,7 +87,7 @@ def doit(field, add_contours, pfiles): ncols = len(pfiles) grid = ImageGrid(fig, 111, nrows_ncols=(nrows, ncols), - axes_pad=0.75, cbar_pad=0.05, label_mode="L", cbar_mode="single") + axes_pad=0.3, cbar_pad=0.05, label_mode="L", cbar_mode="single") for i in range(nrows * ncols): @@ -102,7 +106,7 @@ def doit(field, add_contours, pfiles): function=_lap_rho, units="", validators=[ValidateSpatial(1)]) - domain_frac = 0.15 + domain_frac = 0.2 xmin = ds.domain_left_edge[0] xmax = domain_frac * ds.domain_right_edge[0] @@ -117,13 +121,18 @@ def doit(field, add_contours, pfiles): ymax = yctr + 0.5 * domain_frac * L_y L_y = ymax - ymin - sp = yt.SlicePlot(ds, "theta", field, center=[xctr, yctr, 0.0*cm], width=[L_x, L_y, 0.0*cm], fontsize="12") + sp = yt.SlicePlot(ds, "theta", field, center=[xctr, yctr, 0.0*cm], width=[L_x, L_y, 0.0*cm], fontsize="14") sp.set_buff_size((2400,2400)) - sp.annotate_text((0.05, 0.05), f"time = {float(ds.current_time):8.3f} s", coord_system="axis", text_args={"color": "black"}) + if field == "Temp": + text_color = "white" + else: + text_color = "black" + + sp.annotate_text((0.05, 0.05), f"time = {float(ds.current_time):8.3f} s", coord_system="axis", text_args={"color": text_color}) if field == "Temp": sp.set_zlim(field, 5.e7, 4e9) - sp.set_cmap(field, "magma_r") + sp.set_cmap(field, "magma") elif field == "enuc": sp.set_log(field, True, linthresh=1.e15) sp.set_zlim(field, -1.e22, 1.e22) @@ -151,9 +160,9 @@ def doit(field, add_contours, pfiles): sp._setup_plots() - fig.set_size_inches(19.2, 10.8) + fig.set_size_inches(11, 14) plt.tight_layout() - plt.savefig(f"subch_{field}_sequence.png") + plt.savefig(f"subch_{field}_sequence.pdf") if __name__ == "__main__": diff --git a/Exec/science/subchandra/analysis/subch_zoom.py b/Exec/science/subchandra/analysis/subch_zoom.py new file mode 100755 index 0000000000..38ff03cecc --- /dev/null +++ b/Exec/science/subchandra/analysis/subch_zoom.py @@ -0,0 +1,52 @@ +#!/usr/bin/env python3 + +import argparse +import os +import sys +from functools import reduce + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from mpl_toolkits.axes_grid1 import ImageGrid + +import yt +from yt.fields.derived_field import ValidateSpatial +from yt.frontends.boxlib.api import CastroDataset +from yt.funcs import just_one +# assume that our data is in CGS +from yt.units import amu, cm + +matplotlib.use('agg') + + +plotfile = "subch_plt00000" + +fig = plt.figure() + +ds = CastroDataset(plotfile) + +xmin = 0 * cm +xmax = 1.e8 * cm + +xctr = 0.5 * (xmin + xmax) +L_x = xmax - xmin + +ymin = 5.42e9 * cm +ymax = 5.58e9 * cm +yctr = 0.5 * (ymin + ymax) +L_y = ymax - ymin + +field = "Temp" + +sp = yt.SlicePlot(ds, "theta", field, center=[xctr, yctr, 0.0*cm], width=[L_x, L_y, 0.0*cm], fontsize="14") +sp.set_buff_size((2400,2400)) + +sp.set_zlim(field, 5.e7, 4e9) +sp.set_cmap(field, "magma") + +sp.annotate_contour(("gas", "density"), take_log=True, ncont=3, clim=(1.e4, 1.e6), plot_args={"colors": "white", "linestyles": ":"}) + +sp.set_axes_unit("km") + +sp.save(f"subch_{field}_zoom.pdf")