Merge pull request #52 from lindsayad/devel

Python script examples for plot generation using yt

.gitignore

@@ -885,4 +885,5 @@ HEAD
 config
 description
 hooks
-info
+info
+*kcg*

problems/constant_inlet_outlet_temp/2group.i

@@ -40,7 +40,7 @@ H=162.56
   [../]
 []
 
-[PrecursorKernel]
+[Precursors]
   [./pres]
     var_name_base = pre
     block = 'fuel'

problems/constant_inlet_outlet_temp/4group.i

@@ -40,7 +40,7 @@ H=162.56
   [../]
 []
 
-[PrecursorKernel]
+[Precursors]
   [./pres]
     var_name_base = pre
     block = 'fuel'

problems/constant_inlet_outlet_temp/axial_fluxes.py

@@ -0,0 +1,37 @@
+import yt
+import numpy as np
+
+field_label = {('all', 'group1'): r'$\phi_1\ moltres$',
+               ('all', 'group2'): r'$\phi_2\ moltres$'}
+
+ds = yt.load(
+    '/home/lindsayad/projects/moltres/problems/'
+    'constant_inlet_outlet_temp/2group.e', step=-1)
+
+plt = yt.LinePlot(ds, [('all', 'group1'), ('all', 'group2')], [0, 0, 0], [
+                  0, 162.56, 0], 1000, field_labels=field_label)
+plt.annotate_legend(('all', 'group1'))
+plt.set_x_unit('cm')
+plt.set_x_unit('inch')
+plt.set_xlabel('z (in)')
+plt.set_ylabel(r'Fluxes ($\mathrm{10^{13}cm^{-2} s^{-1}}$)')
+plt._setup_plots()
+
+msre_axial_data = np.loadtxt(
+    '/home/lindsayad/publications/figures/'
+    'msre_axial_fluxes.csv', skiprows=1, delimiter=',')
+fast_x = msre_axial_data[:, 0]
+fast_flux = msre_axial_data[:, 1]
+msre_thermal_axial_data = np.loadtxt(
+    '/home/lindsayad/publications/figures/'
+    'msre_thermal_axial_flux.csv', skiprows=1, delimiter=',')
+thermal_x = msre_thermal_axial_data[:, 0]
+thermal_flux = msre_thermal_axial_data[:, 1]
+plt.plots[('all', 'group1')].axes.plot(
+    fast_x, fast_flux, label=r'$\phi_1\ msre$')
+plt.plots[('all', 'group1')].axes.plot(
+    thermal_x, thermal_flux, label=r'$\phi_2\ msre$')
+plt.plots[('all', 'group1')].axes.legend()
+
+plt.save('/home/lindsayad/publications/figures/'
+         'combined_msre_moltres_axial.eps')

problems/constant_inlet_outlet_temp/fuel_vs_mod_temp.py

@@ -0,0 +1,42 @@
+import yt
+import numpy as np
+from yt.units import K
+
+
+def _units(field, data):
+    return data[('all', 'temp')] * K
+
+
+ds = yt.load(
+    '/home/lindsayad/projects/moltres/problems/'
+    'constant_inlet_outlet_temp/2group.e', step=-1)
+ds.add_field(('all', 'temp_K'), units='K', function=_units, take_log=False)
+lines = []
+lines.append(yt.LineBuffer(
+    ds, (0, 0, 0), (0, 162.56, 0), 1000, label='fuel moltres'))
+lines.append(yt.LineBuffer(ds, (3, 0, 0), (3, 162.56, 0),
+                           1000, label='graphite moltres'))
+plt = yt.LinePlot.from_lines(ds, [('all', 'temp_K')], lines, field_labels={
+                             ('all', 'temp_K'): ''})
+plt.annotate_legend(('all', 'temp_K'))
+plt.set_x_unit('cm')
+plt.set_x_unit('inch')
+plt.set_unit(('all', 'temp_K'), 'degF')
+plt.set_xlabel('z (in)')
+plt.set_ylabel('Temperature (F)')
+plt._setup_plots()
+
+msre_temp_axial_data = np.loadtxt(
+    '/home/lindsayad/publications/figures/'
+    'msre_axial_temps.csv', skiprows=1, delimiter=',')
+x = msre_temp_axial_data[:, 0]
+fuel_temp = msre_temp_axial_data[:, 1]
+graphite_temp = msre_temp_axial_data[:, 2]
+plt.plots[('all', 'temp_K')].axes.plot(x, fuel_temp, label=r'fuel msre')
+plt.plots[('all', 'temp_K')].axes.plot(
+    x, graphite_temp, label=r'graphite msre')
+plt.plots[('all', 'temp_K')].axes.legend()
+plt._setup_plots()
+
+plt.save('/home/lindsayad/publications/figures/'
+         'combined_msre_moltres_axial_temps.eps')

problems/constant_inlet_outlet_temp/radial_fluxes.py

@@ -0,0 +1,38 @@
+import yt
+import numpy as np
+
+field_label = {('all', 'group1'): r'$\phi_1\ moltres$',
+               ('all', 'group2'): r'$\phi_2\ moltres$'}
+
+ds = yt.load(
+    '/home/lindsayad/projects/moltres/problems/'
+    'constant_inlet_outlet_temp/2group.e', step=-1)
+
+plt = yt.LinePlot(ds, [('all', 'group1'), ('all', 'group2')], [0, 81.28, 0], [
+                  73, 81.28, 0], 1000, field_labels=field_label)
+plt.annotate_legend(('all', 'group1'))
+plt.set_x_unit('cm')
+plt.set_x_unit('inch')
+plt.set_xlabel('r (in)')
+plt.set_ylabel(r'Fluxes ($\mathrm{10^{13}cm^{-2} s^{-1}}$)')
+plt.save("moltres_radial_fluxes.eps")
+
+msre_fast_radial_data = np.loadtxt(
+    '/home/lindsayad/publications/figures/'
+    'msre_fast_radial_flux.csv', skiprows=1, delimiter=',')
+fast_x = msre_fast_radial_data[:, 0]
+fast_flux = msre_fast_radial_data[:, 1]
+msre_thermal_radial_data = np.loadtxt(
+    '/home/lindsayad/publications/figures/'
+    'msre_thermal_radial_flux.csv', skiprows=1, delimiter=',')
+thermal_x = msre_thermal_radial_data[:, 0]
+thermal_flux = msre_thermal_radial_data[:, 1]
+plt.plots[('all', 'group1')].axes.plot(
+    fast_x, fast_flux, label=r'$\phi_1\ msre$')
+plt.plots[('all', 'group1')].axes.plot(
+    thermal_x, thermal_flux, label=r'$\phi_2\ msre$')
+plt.plots[('all', 'group1')].axes.legend()
+plt._setup_plots()
+
+plt.save('/home/lindsayad/publications/figures/'
+         'combined_msre_moltres_radial.eps')

python/calc_fuel_pitch_fraction.py

@@ -14,16 +14,16 @@ def calc_spacing(R, n):
     x = sp.symbols('x')
     Af = 0
     Am = 0
-    for m in range(2*n - 1):
+    for m in range(2 * n - 1):
         if m % 2 == 0:
-            Af += sp.pi * (R/n * (m/2 + x))**2
+            Af += sp.pi * (R / n * (m / 2 + x))**2
         if m % 2 == 1:
-            Af -= sp.pi * (R/n * (m+1)/2)**2
-    for m in range(2*n):
+            Af -= sp.pi * (R / n * (m + 1) / 2)**2
+    for m in range(2 * n):
         if m % 2 == 0:
-            Am -= sp.pi * (R/n * (m/2 + x))**2
+            Am -= sp.pi * (R / n * (m / 2 + x))**2
         if m % 2 == 1:
-            Am += sp.pi * (R/n * (m+1)/2)**2
+            Am += sp.pi * (R / n * (m + 1) / 2)**2
     return sp.solve(Af / (Af + Am) - .225, x)
 
 

python/moltres_functions.py

@@ -0,0 +1,7 @@
+def actual_center_and_widths(ds):
+    actual_domain_widths = (ds.domain_right_edge.to_ndarray(
+    ) - ds.domain_left_edge.to_ndarray()) / 1.2
+    actual_le = ds.domain_left_edge.to_ndarray() + actual_domain_widths * .1
+    actual_re = ds.domain_right_edge.to_ndarray() - actual_domain_widths * .1
+    actual_center = (actual_le + actual_re) / 2
+    return actual_domain_widths, actual_center

python/process_serpent_data_build_least_squares_fits.py

@@ -3,28 +3,30 @@
 import numpy as np
 import subprocess
 
-xsecs = ["FLUX", "REMXS", "FISSXS", "NUBAR", "NSF", "FISSE", "DIFFCOEF", \
+xsecs = ["FLUX", "REMXS", "FISSXS", "NUBAR", "NSF", "FISSE", "DIFFCOEF",
          "RECIPVEL", "CHI", "GTRANSFXS", "BETA_EFF", "DECAY_CONSTANT"]
 num_groups = 2
 num_precursor_groups = 8
-xsec_dict = {"FLUX" : num_groups, "REMXS" : num_groups, "FISSXS" : num_groups, \
-             "NUBAR" : num_groups, "NSF" : num_groups, "FISSE" : num_groups, \
-             "DIFFCOEF" : num_groups, "RECIPVEL" : num_groups, "CHI" : num_groups, \
-             "GTRANSFXS" : num_groups * num_groups, "BETA_EFF" : num_precursor_groups, \
-             "DECAY_CONSTANT" : num_precursor_groups}
+xsec_dict = {"FLUX": num_groups, "REMXS": num_groups, "FISSXS": num_groups,
+             "NUBAR": num_groups, "NSF": num_groups, "FISSE": num_groups,
+             "DIFFCOEF": num_groups, "RECIPVEL": num_groups, "CHI": num_groups,
+             "GTRANSFXS": num_groups * num_groups,
+             "BETA_EFF": num_precursor_groups,
+             "DECAY_CONSTANT": num_precursor_groups}
 materials = ["fuel", "mod"]
 
 for material in materials:
-    write_file_name =  subprocess.call(["touch", "write_" + material])
-    file_base = "msr2g_Th_U_two_mat_homogenization_dens_func_" + material + "_data_func_of_" + material + "_temp"
+    write_file_name = subprocess.call(["touch", "write_" + material])
+    file_base = "msr2g_Th_U_two_mat_homogenization_dens_func_" + \
+        material + "_data_func_of_" + material + "_temp"
     write_file = open("write_" + material, 'r+')
     for xsec in xsecs:
         file_name = file_base + "_" + xsec + ".txt"
         data = np.loadtxt(file_name)
         temperature = data[:, 0]
         A = np.vstack([temperature, np.ones(len(temperature))]).T
         for k in range(xsec_dict[xsec]):
-            dep_var = data[:, k+1]
+            dep_var = data[:, k + 1]
             m, c = np.linalg.lstsq(A, dep_var)[0]
             write_file.write(str(m) + " " + str(c))
             write_file.write('\n')

python/slice_plot_from_2d_sim_ex.py

@@ -0,0 +1,50 @@
+import yt
+from moltres_functions import actual_center_and_widths
+
+field_label = {
+    'group1': r'$\phi_1\cdot$10$^{13}$ cm$^{-2}$ s$^{-1}$',
+    'group2': r'$\phi_2\cdot$10$^{13}$ cm$^{-2}$ s$^{-1}$',
+    'temp': 'T (K)',
+    'pre1': r'C$_1$ cm$^{-3}$',
+    'pre2': r'C$_2$ cm$^{-3}$',
+    'pre3': r'C$_3$ cm$^{-3}$',
+    'pre4': r'C$_4$ cm$^{-3}$',
+    'pre5': r'C$_5$ cm$^{-3}$',
+    'pre6': r'C$_6$ cm$^{-3}$'}
+names = [name for name in field_label.keys() if 'pre' in name]
+for name in names:
+    field_label[name + '_scaled'] = field_label[name]
+
+
+def _scale(field, data):
+    new_field = field.name
+    ftype, fname = new_field
+    original_fname = fname.split('_')[0]
+    return 1e13 * data[('all', original_fname)]
+
+
+ds = yt.load(
+    '/home/lindsayad/projects/moltres/tests/'
+    'twod_axi_coupled/auto_diff_rho.e', step=-1)
+for i in range(1, 7):
+    fname = 'pre%d_scaled' % i
+    ds.add_field(('all', fname), function=_scale, take_log=False)
+actual_domain_widths, actual_center = actual_center_and_widths(ds)
+
+for field in ds.field_info.keys():
+    field_type, field_name = field
+    if field_type != 'all':
+        continue
+    if 'scaled' not in field_name:
+        continue
+    slc = yt.SlicePlot(ds, 'z', field, origin='native', center=actual_center)
+    slc.set_log(field, False)
+    slc.set_width((actual_domain_widths[0], actual_domain_widths[1]))
+    slc.set_xlabel("r (cm)")
+    slc.set_ylabel("z (cm)")
+    if field_name == 'temp':
+        slc.set_zlim(field, 922, ds.all_data()[('all', 'temp')].max())
+    slc.set_colorbar_label(field, field_label[field_name])
+    slc.set_figure_size(3)
+    slc.save('/home/lindsayad/publications/figures/2d_gamma_heating_' +
+             field_name + '.eps')

python/slice_plot_from_3d_sim_ex.py

@@ -0,0 +1,36 @@
+import yt
+from moltres_functions import actual_center_and_widths
+
+field_label = {'group1': r'$\phi_1\cdot$10$^{13}$ cm$^{-2}$ s$^{-1}$',
+               'group2': r'$\phi_2\cdot$10$^{13}$ cm$^{-2}$ s$^{-1}$',
+               'temp': 'T (K)',
+               'pre1': r'C$_1$ cm$^{-3}$',
+               'pre2': r'C$_2$ cm$^{-3}$',
+               'pre3': r'C$_3$ cm$^{-3}$',
+               'pre4': r'C$_4$ cm$^{-3}$',
+               'pre5': r'C$_5$ cm$^{-3}$',
+               'pre6': r'C$_6$ cm$^{-3}$'}
+
+ds = yt.load(
+    '/home/lindsayad/Dropbox/MSR_data/converged_gamma_heating_3d_steady_state/'
+    'gamma_heating_newton.e', step=-1)
+actual_domain_widths, actual_center = actual_center_and_widths(ds)
+
+for field in ds.field_list:
+    field_type, field_name = field
+    import pdb
+    pdb.set_trace()
+    if field_type != 'all' or field_name != 'temp':
+        continue
+    slc = yt.SlicePlot(ds, 'z', field, origin='native', center=(
+        actual_center[0], actual_center[1], actual_domain_widths[2]))
+    slc.set_log(field, False)
+    slc.set_width((actual_domain_widths[1], actual_domain_widths[2]))
+    slc.set_xlabel("x (cm)")
+    slc.set_ylabel("y (cm)")
+    if field_name == 'temp':
+        slc.set_zlim(field, 922, ds.all_data()[('all', 'temp')].max())
+    slc.set_colorbar_label(field, field_label[field_name])
+    slc.set_figure_size(3)
+    slc.save('/home/lindsayad/publications/figures/3d_gamma_heating_z_slice_' +
+             field_name + '.eps')

tests/laminar_flow/tests

@@ -48,5 +48,6 @@
     type = 'Exodiff'
     input = 'boussinesq_square.i'
     exodiff = 'boussinesq_square_out.e'
+    abs_zero = 5e-5
   [../]
 []