monoblock_neutronics/post_processing.py at main · RemDelaporteMathurin/monoblock_neutronics · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
import openmc
import regular_mesh_plotter as rmp
import matplotlib.pyplot as plt
import openmc_tally_unit_converter as otuc
from scipy.stats import linregress
import numpy as np


def shape_tally(tally, threshold=0):
    value = np.array(tally.mean)
    std_dev = np.array(tally.std_dev)

    voxel_volume = otuc.compute_volume_of_voxels(tally)  # cm3
    voxel_volume *= 1e-6  # m3
    value *= 1 / voxel_volume  # m-3 per neutron
    value *= source_strength  # m-3 s-1

    std_dev *= 1 / voxel_volume  # m-3 per neutron
    std_dev *= source_strength  # m-3 s-1

    value = rmp.reshape_values_to_mesh_shape(tally, value)
    std_dev = rmp.reshape_values_to_mesh_shape(tally, std_dev)

    indexes = np.where(value < threshold)

    value[indexes] = np.nan
    std_dev[indexes] = np.nan
    return value, std_dev


def plot_distribution(mean, std_dev):
    distribution = np.nanmean(mean, axis=1)
    std_dev_distrib = np.nanmean(std_dev, axis=1)
    std_dev_distrib = std_dev_distrib[np.logical_not(np.isnan(distribution))]
    distribution = distribution[np.logical_not(np.isnan(distribution))]
    x = np.linspace(0, height_mb, len(distribution), endpoint=True)

    res = linregress(x[1:], distribution[::-1][1:])

    plt.errorbar(
        x, distribution[::-1], yerr=2 * std_dev_distrib[::-1], fmt="o", alpha=0.5
    )
    (line,) = plt.plot(x, res.slope * x + res.intercept)
    x_annotation = 0.75
    y_annotation = (res.slope * x_annotation + res.intercept) * 1.1
    plt.annotate(
        "{:.1e} $x$ + {:.1e}".format(res.slope, res.intercept),
        (x_annotation, y_annotation),
        color=line.get_color(),
    )


result = openmc.StatePoint("statepoint.50.h5")

height_mb = 2.5  # cm
source_strength = otuc.find_source_strength(
    fusion_energy_per_second_or_per_pulse=500e6
)  # n/s

# He generation mesh tally
helium_generation_mesh = result.get_tally(name="(n,Xa)_on_2D_mesh_yz")

value, std_dev = shape_tally(helium_generation_mesh, threshold=4.75e18)

rmp.plot_regular_mesh_values(
    values=value,
    extent=rmp.get_tally_extent(helium_generation_mesh),
    rotate_plot=180,
    label="He generation (m$^{-3}$ s$^{-1}$)",
)
plt.gca().spines.right.set_visible(False)
plt.gca().spines.top.set_visible(False)
plt.savefig("helium_generation_in_monoblock.png")
plt.savefig("helium_generation_in_monoblock.pdf")


# Heat generation mesh tally std dev
plt.figure()

rmp.plot_regular_mesh_values(
    values=std_dev,
    extent=rmp.get_tally_extent(helium_generation_mesh),
    rotate_plot=180,
    label="He generation (m$^{-3}$ s$^{-1}$)",
)
plt.gca().spines.right.set_visible(False)
plt.gca().spines.top.set_visible(False)
plt.savefig("helium_generation_std_dev.png")
plt.savefig("helium_generation_std_dev.pdf")


# depth distribution of helium generation

plt.figure()
plot_distribution(value, std_dev)
plt.xlabel("Distance from the top surface (cm)")
plt.ylabel("He generation  (m$^{-3}$ s$^{-1}$)")
plt.ylim(0, 8e18)
plt.gca().spines.right.set_visible(False)
plt.gca().spines.top.set_visible(False)
plt.savefig("he_generation_distribution.png")
plt.savefig("he_generation_distribution.pdf")


# Heat generation mesh tally
plt.figure()
heat_generation_mesh = result.get_tally(name="heating_on_2D_mesh_yz")

value, std_dev = shape_tally(heat_generation_mesh, threshold=4.22e27)
value *= 1.602e-19 * 1e-6  # MW  m-3
std_dev *= 1.602e-19 * 1e-6  # MW  m-3

rmp.plot_regular_mesh_values(
    values=value,
    extent=rmp.get_tally_extent(heat_generation_mesh),
    rotate_plot=180,
    label="Heat generation (MW m$^{-3}$)",
)
plt.gca().get_images()[0].set_cmap("inferno")
plt.gca().spines.right.set_visible(False)
plt.gca().spines.top.set_visible(False)
plt.savefig("heat_generation.png")

# depth distribution of helium generation
plt.figure()

plot_distribution(value, std_dev)
plt.xlabel("Distance from the top surface (cm)")
plt.ylabel("Heat generation  (MW m$^{-3}$)")
plt.ylim(bottom=0)
plt.gca().spines.right.set_visible(False)
plt.gca().spines.top.set_visible(False)
plt.savefig("heat_distribution.png")