update the subchandra plots to those used in the paper

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)

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,15 +76,18 @@ 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:
  nrows = 1
  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__":
 

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")