Journal Search Engine

Volume/Issue :

Year(s) : to

Search :

Title :

Author :

Keyword :

Abstract :

Figure :

Table :

Reference :

Search Advanced Search

Adode Reader(link)

View PDF Download PDF Export Citation Korean Bibliography PMC Previewer

ISSN : 1738-1894(Print)
ISSN : 2288-5471(Online)

Journal of Nuclear Fuel Cycle and Waste Technology Vol.18 No.1 pp.19-29
DOI : https://doi.org/10.7733/jnfcwt.2020.18.1.19

Efficient Computation of Radioactive Decay with Graph Algorithms

Tae-Sic Yoo*

Idaho National Laboratory, 775 MK Simpson Blvd, Idaho Falls, Idaho 83402, USA

Corresponding Author. Tae-Sic Yoo, Idaho National Laboratory, E-mail: tae-sic.yoo@inl.gov, Tel: +1-208-526-9394

Received December 4, 2019 ; Revised February 27, 2020 ; Approved March 25, 2020

Abstract

This paper gives two graph-based algorithms for radioactive decay computation. The first algorithm identifies the connected components of the graph induced from the given radioactive decay dynamics to reduce the size of the problem. The solutions are derived over the precalculated connected components, respectively and independently. The second algorithm utilizes acyclic structure of radioactive decay dynamics. The algorithm evaluates the reachable vertices of the induced system graph from the initially activated vertices and finds the minimal set of starting vertices populating the entire reachable vertices. Then, the decay calculations are performed over the reachable vertices from the identified minimal starting vertices, respectively, with the partitioned initial value over the reachable vertices. Formal arguments are given to show that the proposed graph inspired divide and conquer calculation methods perform the intended radioactive decay calculation. Empirical efforts comparing the proposed radioactive decay calculation algorithms are presented.

Key Words : Radioactive decay computation , Nuclear fuel cycle , Graph algorithms

초록

키워드 :

This article has been cited by 0 article in crossref

Cited-By

Funding:

1. Introduction

The system of linear, constant real coefficient ordinary differential equations, referred to as the linear dynamic system in this article, appears in various disciplines. The linear dynamic systems can describe reactor depletion and radioactive decay dynamics and many related literature results concerning the radioactive decay linear dynamic system are available spanning over a century [1-8].

Many nuclear engineering and science applications require numerous iterations of radioactive decay calculations. An example requiring numerous iterative radioactive decay calculations is molten salt reactors (MSR) [9,10] where the reactor discharges a part of the molten salt inventory to the used fuel processor (UFP) for chemical treatment. The UFP treats the discharged used fuel to separate the fission products and feeds the processed fuel back to the reactor along with the exogenous fuel feed to make up the discharged salt from the reactor. The performance of the MSR system requiring online used fuel processing is tightly coupled with the performance of the UFP as the daily make-up fuel feed comes from the used fuel process operation. The evaluation of the MSR system performance requires numerous iterations of reactor depletion and used fuel process simulations. Another example is a pyroprocessing facility where radioactive materials are accumulated, divided, moved, and decayed daily [11-14]. The evaluation of the pyroprocessing facility via modeling and simulation also involves tens of thousands of radioactive material decay calculations. In these applications, computationally efficient radioactive decay calculations are crucial to avoid a prolonged overall calculation time. The radioactive decay calculation involves matrix exponential that can be addressed with various numerical and mathematical techniques. A review on this topic with extensive references is available in [15]. The foci have been given to find accurate, robust, and efficient algorithms for matrix exponential.

The intent of this article is not in contributing to the prescribed aspects. Rather, specifically having the radioactive decay application in mind, the focus is on identifying the structural information of the graph induced from the radioactive decay dynamics in order to fragment the problem into the set of smaller size subproblems. The computational gain depends on the degree of the problem fragmentation and the radioactive decay application exhibits a high degree of fragmentation allowing a multiple order of magnitude computational improvement. To do so, section 2 gives relevant preliminaries regarding the linear dynamic systems and graph theories. Section 3 presents two graph-based algorithms for solving the linear dynamic system. An application to the radioactive decay calculation is detailed in Section 4. Section 5 gives concluding remarks.

2. Preliminaries

Consider the system of linear, constant real coefficient ordinary differential equations given below.

{\dot{x}}_{k} = \sum_{i = 1}^{N} a_{k, i} \cdot x_{i} (t); ​ for k = 1 to N

(1)

Denote the set of real numbers as R. With matrix and vector notations, a compact representation of the above equation set can be written as below.

\dot{X} (t) = A \cdot X (t)

(2)

where

\begin{array}{l} \dot{X} (t) = : [\begin{matrix} {\dot{x}}_{1} (t) \\ {\dot{x}}_{2} (t) \\ ⋮ \\ {\dot{x}}_{N} (t) \end{matrix}] \in R^{N \times 1}; X (t) : = [\begin{matrix} x_{1} (t) \\ x_{2} (t) \\ ⋮ \\ x_{N} (t) \end{matrix}] \in R^{N \times 1}; \\ A : = [\begin{matrix} a_{1, 1} & a_{1, 2} & \dots & a_{1, N} \\ a_{2, 1} & a_{2, 2} & \dots & a_{2, N} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ a_{N, 1} & a_{N, 2} & \dots & a_{N, N} \end{matrix}] \in R^{​ N \times N} \end{array}

(3)

The constraint is that the solution X(t) satisfies the initial state X₀, i.e., X(0) = X₀. Then, the solution is expressed with matrix multiplication of the initial state X₀ and the matrix exponential e^At as below.

X (t) = e^{A t} \cdot X_{0}

(4)

In this article, we refer matrix A as the rate matrix and e^At as the transition matrix.

A directed graph G is an ordered pair (V, E) where V is the set of vertices and E⊆V×V is the set of directed edges. A directed graph G = (V, E) is induced from rate matrix A where V is the set of vertices representing the index set of equation (1) and E is the set of edges reflecting the influence structure of rate matrix A. That is,

V = {1, 2, \dots, N} and (i, j) \in E \subseteq V \times V \Leftrightarrow a_{j, i} \neq 0

The set {V_i} is called the partition of a set V if

[V = \underset{i}{\cup} V_{i}] \land [V_{i} \cap V_{j} = \emptyset i f i \neq j]

A path of length k in a directed graph G = (V, E) is the sequence of vertex {v₀, v₁, …, v_k} such that (v_i_–1, v_i)∈E, for all i = 1, …, k. With this, we call that u∈V is reachable from v∈V if there exists a positive integer k such that {v₀, v₁, …, v_k} forms a path of length k in G = (V, E) with v = v₀, and u = v_k. Let’s define the reverse edge set of E as below.

E^{- 1} = : {(ν, u) : (u, ν) \in E}

We call that a directed graph G = (V, E) is weakly connected, for every pair of vertices u, v∈V, if u is reachable from v in G' = (V, E∪E^–1). A graph G_s = (V_s, E_s) is called a subgraph of G = (V, E) if V_s⊆V and E_s⊆E. A graph G_s = (V_s, E_s) is called a strict subgraph of G = (V, E) if V_s⊂V or E_s⊂E. A subgraph G_a = (V_a, E_a) of a directed graph G = (V, E) is called a weakly connected component of a directed graph G if G_a is weakly connected and there is no subgraph G_b of G such that G_b is weakly connected and G_a is a strict subgraph of G_b. With this, we define the weakly connected components of a directed graph G = (V, E) denoted by {G_c,i := (V_c,i, E_c,i)} such that for all i, G_c,i is a weakly connected component and {V_c,i} and {E_c,i} are the partition of V and E, respectively.

Below gives an example illustrating the definitions presented above.

Example 1. Consider the set of linear differential equations below.

{\dot{x}}_{1} = x_{2}; {\dot{x}}_{2} = x_{3}; {\dot{x}}_{3} = - x_{3}; {\dot{x}}_{4} = 0

(5)

Rate matrix A in (3) can be constructed from the linear dynamic system in (5) as below.

A = [\begin{matrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & - 1 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}]

Digraph G = (V, E) induced from rate matrix A is given in (6); Fig. 1 gives its graphical representation.

Fig. 1

Induced digraph G = (V,E) from rate matrix A.

V = {1, 2, 3, 4} a n d E = {(2, 1), (3, 2), (3, 3)}

(6)

There are two weakly connected components: G_c_,1 = ({1,2,3}, {(2,1),(3,2),(3,3)}) and G_c_,2 =({4}, Ø).

The induced digraph facilitates an intuitive examination of the influence structure of the involved variables. Typically, the Depth First Search (DFS) over the undirected graph is used for finding the weakly connected components with O(|V|+|E|) complexity. A textbook treatise of an algorithm for finding the weakly connected components is available in [16].

3. Subdynamics Fragmentation

Denote the initial state imposed on the system given in Equation (1) by

X_{0} : = {[x_{1, 0}, x_{2, 0} \dots \in x_{N, 0}]}^{T}

(7)

Consider the following index set collecting the index of the nonzero elements in Equation (7). We call this the initially activated vertices.

\hat{I} = {i : x_{i, 0} \neq 0}

The reach of i∈Î in the induced graph G = (V, E) can be defined as below.

R_{i} : = {ν \in V | ν is reachable from i \in \hat{I}}

The union of R_i gives the union of the reachable vertices from the initially activated vertices, denoted by R_Î. That is,

R_{\hat{I}} = \underset{i \in \hat{I}}{\cup} R_{i}

Then, we can construct the subgraph induced with R_Î, $G_{R_{\hat{I}}} = (R_{\hat{I}}, E_{R_{\hat{I}}})$ , where

E_{R_{\hat{I}}} = {(i, j) \in E : i, j \in R_{\hat{I}}}

The subgraph $G_{R_{\hat{I}}} = (R_{\hat{I}}, E_{R_{\hat{I}}})$ of the induced graph G = (V, E) can be constructed with an application of DFS with O(|V|+|E|) complexity.

The set of unreachable vertices from Î, denoted by $R_{\hat{I}}^{c}$ , is the complement of R_Î. That is, $R_{\hat{I}}^{c} {1, \dots, N} \ R_{\hat{I}}$ . By grouping the reachable and unreachable vertices, Equation (2) can be rewritten as below.

[\begin{matrix} {\dot{X}}_{R_{\hat{I}}} (t) \\ {\dot{X}}_{R_{\hat{I}}^{c}} (t) \end{matrix}] = [\begin{matrix} A_{R_{\hat{I}} \to R_{\hat{I}}} & A_{R_{\hat{I}}^{c} \to R_{\hat{I}}} \\ A_{R_{\hat{I}} \to R_{\hat{I}}^{c}} & A_{R_{\hat{I}}^{c} \to R_{\hat{I}}^{c}} \end{matrix}] [\begin{matrix} X_{R_{\hat{I}}} (t) \\ X_{R_{\hat{I}}^{c}} (t) \end{matrix}]

where $X_{R_{\hat{I}}} (t)$ is the collection of the states corresponding to the reach of the initially activated vertices. On the other hand, $X_{R_{\hat{I}}^{c}} (t)$ represents the remaining vertices that are not in the reach of the initially activated vertices. Note that

{\dot{X}}_{R_{\hat{I}}} (t) = A_{R_{\hat{I}} \to R_{\hat{I}}} X_{R_{\hat{I}}} (t) + A_{R_{\hat{I}}^{c} \to R_{\hat{I}}} X_{R_{\hat{I}}^{c}} (t)

(8)

and

{\dot{X}}_{R_{_{\hat{I}}}^{c}} (t) = A_{R_{\hat{I}} \to R_{\hat{I}}^{c}} X_{R_{\hat{I}}} (t) + A_{R_{\hat{I}}^{c} \to R_{\hat{I}}^{c}} X_{R_{\hat{I}}^{c}} (t)

(9)

As $R_{\hat{I}}^{c}$ is not reachable from R_Î, we have $A_{R_{\hat{I}} \to R_{\hat{I}}^{c}} = 0$ . Then, Equation (9) becomes

{\dot{X}}_{R_{_{\hat{I}}}^{c}} (t) = A_{R_{\hat{I}}^{c} \to R_{\hat{I}}^{c}} X_{R_{\hat{I}}^{c}} (t)

Thus, the unreachable part of the linear dynamic system becomes

X_{R_{_{\hat{I}}}^{c}} (t) = e^{A_{R_{\hat{I}}^{c} \to R_{\hat{I}}^{c}} t} X_{R_{\hat{I}}^{c}} (0) = e^{A_{R_{\hat{I}}^{c} \to R_{\hat{I}}^{c}} t} 0 = 0

With $X_{R_{\hat{I}}^{c}} (t) = 0$ , Equation (8) is further simplified to

{\dot{X}}_{R_{_{\hat{I}}}} (t) = A_{R_{\hat{I}} \to R_{\hat{I}}} X_{R_{\hat{I}}} (t)

(10)

The solution of Equation (10) is

X_{R_{_{\hat{I}}}} (t) = e^{A_{R_{\hat{I}} \to R_{\hat{I}}} t} \cdot X_{R_{\hat{I}}} (0)

(11)

Given the subgraph $G_{R_{\hat{I}}} = (R_{\hat{I}}, E_{R_{\hat{I}}})$ , building the rate matrix $A_{R_{\hat{I}} \to R_{\hat{I}}}$ takes O(|R_Î|²). The matrix exponential $e^{A_{R_{\hat{I}} \to R_{\hat{I}}} t}$ can be calculated with O(|R_Î|³) [15]. The multiplication, $e^{A_{R_{\hat{I}} \to R_{\hat{I}}} t} X_{R_{\hat{I}}} (0)$ , can be done in O(|R_Î|^2.3737) with the algorithm given in [17]. Thus, overall complexity is governed by the matrix exponential computation step and O(|R_Î|³).

3.1 Algorithm with Weakly Connected Component (WCC)

Evaluating the activated part of the system with graph reachability analysis before deriving the solution provides an opportunity to reduce the size of the problem as in Equation (10). In this subsection, for the problem formulated in Equation (10), we attempt to divide the problem further into smaller subproblems with the application of the weakly connected components graph analysis described in Section 2, which can be performed with $O (| R_{\hat{I}} | + | E_{R_{\hat{I}}} |)$ complexity. In a rough sense, two vertices of a directed graph, G, are weakly connected if a path can be found from one to the other in the undirected graph created by ignoring the direction of the edges in the directed graph, G.

First, denote the i^th weakly connected component of $G_{R_{\hat{I}}} (R_{\hat{I}}, E_{R_{\hat{I}}})$ by G_c,i := (V_c,i, E_i) where V_c,i⊆R_Î and E_c,i is the edges defined over V_c,i. That is,

E_{c, i} : = {(u, ν) \in E : u, ν \in V_{c, i}}

Suppose that there are M weakly connected components of G_Î. By rearranging the index, one can get the following block diagonal rate matrix where rate matrix A_c,i is constructed from the connected component G_c,i := (V_c,i, E_i). We note that the value of off-diagonals elements representing the connection between the weakly connected components is zero. That is,

A_{R_{\hat{I}}} = [\begin{matrix} A_{c, 1} & 0 & \dots & 0 \\ 0 & A_{c, 2} & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & A_{c, M} \end{matrix}]

With the above rate matrix, the corresponding linear dynamic system can be described as below.

[\begin{matrix} {\dot{X}}_{c, 1} \\ {\dot{X}}_{c, 2} \\ ⋮ \\ {\dot{X}}_{c, M} \end{matrix}] = [\begin{matrix} A_{c, 1} & 0 & \dots & 0 \\ 0 & A_{c, 2} & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & A_{c, M} \end{matrix}] [\begin{matrix} X_{c, 1} \\ X_{c, 2} \\ ⋮ \\ X_{c, M} \end{matrix}]

Then, for i = 1, …, M,

X_{c, i} (t) = e^{A_{c, i} t} X_{c, i} (0)

(12)

Rather than solving the linear system with Equations (4) or (11) in a monolithic manner, the solution can be derived for each diagonal block induced from the weakly connected components.

Note that synthesizing rate matrix A_c,i takes O(|G_i|²) and building all diagonal block matrices takes $\sum_{i}^{M} O ({| G_{i} |}^{2})$ . As $| R_{\hat{I}} | = \sum_{i}^{M} | G_{i} |$ , one can see that

\sum_{i}^{M} O ({| G_{i} |}^{2}) \leq O ({| R_{\hat{I}} |}^{2})

(13)

Equation (13) implies that the overall complexity of fragmenting the whole dynamics into the weakly connected component based subdynamics is not worse than building the rate matrix $A_{R_{\hat{I}} \to R_{\hat{I}}}$ from $G_{R_{\hat{I}}} = (R_{\hat{I}}, E_{R_{\hat{I}}})$ Also, the overall complexity of building matrix exponentials of all subdynamics is $\sum_{i}^{M} O ({| G_{i} |}^{3})$ , which is not worse than O(|R_Î|³), the complexity of computing matrix exponential $e^{A_{R_{\hat{I}} \to R_{\hat{I}}} t}$ . That is,

\sum_{i}^{M} O ({| G_{i} |}^{3}) \leq O ({| R_{\hat{I}} |}^{3})

(14)

Other miscellaneous computation steps such as matrix multiplications and combining subdynamics solution can be performed with complexities less than O(|R_Î|³). Thus, the proposed solution strategy with the weakly connected component algorithm takes O(|R_Î|³). For brevity, we shall refer to the algorithm described in this section as WCC.

3.2 Algorithm with Minimum Starting Vertices (MSV)

The section assumes that $G_{R_{\hat{I}}} = (R_{\hat{I}}, E_{R_{\hat{I}}})$ is acyclic. The minimum starting vertex set of $G_{R_{\hat{I}}}$ is the collection of vertices that do not have incoming edges. That is,

S (G_{R_{\hat{I}}}) = {ν \in R_{\hat{I}} : \forall u \in V \ {ν}, \exists (u, ν) \in E_{R_{\hat{I}}}}

As $G_{R_{\hat{I}}}$ is acyclic, all vertices in R_Î should be reachable from the minimum starting vertex set. The reachable vertices from the initially activated vertices Î can be expressed with $S (G_{R_{\hat{I}}})$ as below.

R_{\hat{I}} = \underset{i \in S (G_{R_{\hat{I}}})}{\cup} R_{i}

For $i \in S (G_{R_{\hat{I}}})$ , one can construct the disjoint set of the reachable vertices as below.

R_{i}^{-} = R_{i} \ \underset{j < i}{\cup} R_{j}

Then, ${R_{i}^{-} : i \in S (G_{R_{\hat{I}}})}$ becomes a partition of R_Î. That is,

[R_{\hat{I}} = \underset{i \in S (G_{R_{\hat{I}}})}{\cup} R_{i}^{-}] \land [R_{i}^{-} \cap R_{j}^{-} = \emptyset if i \neq j]

For A⊆{1, …, N}, the projection map P_A : R^N → R^N, is defined as:

P_{A} ({[x_{1} x_{2} \dots x_{N}]}^{T}) : = {[z_{1} z_{2} \dots z_{N}]}^{T}

where

z_{i} = {\begin{array}{l} x_{i} & i f i \in A \\ 0 & o . w . \end{array}

Denote the complement of R_Î by $R_{\hat{I}}^{c} : = {1, \dots, N} \ R_{\hat{I}}$ . As the non-zero initial condition index Î and $R_{\hat{I}}^{c}$ are disjoint, we have

P_{R_{\hat{I}}^{c}} ({[x_{1, 0}, x_{2, 0} \dots x_{N, 0}]}^{T}) = 0

Note that, with any partition {V₁, …, V_M} of R_Î, ${V_{1}, \dots, V_{M}, R_{\hat{I}}^{c}}$ becomes the partition of the index set {1, …, N}. The initial condition can then be expressed with the summation of the projected initial values with the partition ${V_{1}, \dots, V_{M}, R_{\hat{I}}^{c}}$ of R_Î. That is,

X (0) = P_{R_{\hat{I}}^{c}} (X (0)) + \sum_{i = 1}^{M} P_{V_{i}} (X (0)) = \sum_{i = 1}^{M} P_{V_{i}} (X (0))

Note that X(t) can be expressed as the sum of respective state evolution with the projected initial state X₀ induced from any partition {V_i} of the index set V. That is,

X (t) = e^{A t} \sum_{i} P_{V_{i}} (X (0)) = \sum_{i} e^{A t} P_{V_{i}} (X (0))

As ${R_{i}^{-} : i \in S (G_{R_{\hat{I}}})}$ is a partition of R_Î, ${R_{i}^{-} : i \in S (G_{R_{\hat{I}}})}$ and $R_{\hat{I}}^{c}$ become the partition of the index set V. Then, we have

X (0) = P_{R_{\hat{I}}^{c}} (X (0)) + \sum_{i \in S (G_{R_{\hat{I}}})} P_{R_{i}^{-}} (X (0))

As $P_{R_{\hat{I}}^{c}} (X (0)) = 0$ , we have

X (t) = \sum_{i \in S (G_{R_{\hat{I}}})} e^{A t} P_{R_{i}^{-}} (X (0))

Let us focus on the following dynamics below.

{\dot{X}}_{i} (t) = A X_{i} (t); X_{i} (0) = P_{R_{i}^{-}} (X (0))

(15)

The solution of the dynamics given in Equation (15) is

X_{i} (t) = e^{A t} P_{R_{i}^{-}} (X (0))

The set of unreachable vertices of G from $i \in S (G_{R_{\hat{I}}})$ , denoted by $R_{\hat{I}}^{c}$ , is the complement of R_Î, that is, $R_{\hat{I}}^{c} {1, \dots, N} \ R_{\hat{I}}$ . Note that $R_{i}^{-} \subseteq R_{\hat{I}}$ . Thus, we have $R_{i}^{-} \cap R_{\hat{I}}^{c} = \emptyset$ and all components of X(0) belonging to $R_{\hat{I}}^{c}$ become zero with the projection map $P_{R_{\hat{I}}^{-}}$ . Grouping the reachable and unreachable vertices and rewriting Equation (15) result in

[\begin{matrix} {\dot{X}}_{R_{i}} (t) \\ {\dot{X}}_{R_{i}^{c}} (t) \end{matrix}] = [\begin{matrix} A_{R_{i} \to R_{i}} & A_{R_{i}^{c} \to R_{i}} \\ A_{R_{i} \to R_{i}^{c}} & A_{R_{i}^{c} \to R_{i}^{c}} \end{matrix}] [\begin{matrix} X_{R_{i}} (t) \\ X_{R_{i}^{c}} (t) \end{matrix}]; [\begin{matrix} X_{R_{i}} (0) \\ X_{R_{i}^{c}} (0) \end{matrix}] = [\begin{matrix} X_{R_{i}} (0) \\ 0 \end{matrix}]

With $A_{R_{\hat{I}} \to R_{\hat{I}}^{c}} = 0$ , one can see that

\begin{matrix} {\dot{X}}_{R_{i}^{c}} (t) = A_{R_{i} \to R_{i}^{c}} X_{R_{i}} (t) + A_{R_{i}^{c} \to R_{i}^{c}} X_{R_{i}^{c}} (t) \\ = A_{R_{i}^{c} \to R_{i}^{c}} X_{R_{i}^{c}} (t) \end{matrix}

Then, the solution of the unreachable part of the dynamics becomes for all t ≥ 0,

X_{R_{i}^{c}} (t) = e^{A_{R_{i}^{c} \to R_{i}^{c}} t} X_{R_{i}^{c}} (0) + e^{A_{R_{i}^{c} \to R_{i}^{c}} t} 0 = 0

The dynamics of the reachable part is

\begin{matrix} {\dot{X}}_{R_{i}} (t) = A_{R_{i} \to R_{i}} X_{R_{i}} (t) + A_{R_{i}^{c} \to R_{i}} X_{R_{i}^{c}} (t) \\ = A_{R_{i} \to R_{i}} X_{R_{i}} (t) \end{matrix}

The solution of the reachable part dynamics becomes

X_{R_{i}} (t) = e^{A_{R_{i} \to R_{i}} t} X_{R_{i}} (0)

(16)

The above derivation implies that the state evolution activated with the initial state X₀ can be expressed by combining subdynamics solutions described in Equation (16) that are governed by the reachable vertices from a given starting vertex and an appropriately partitioned initial value. For brevity, we shall refer to the algorithm described in this section as MSV. Below gives an example that illustrates the partition concepts described in this section.

Example 2. Consider the digraph G = (V, E) in Fig. 2. The shaded vertices in the figure are the starting vertices of the graph G. The reachable vertices from each starting vertex are R₁ and R₅. The partition of the vertex set V is {{1,2,3,4,6},{5}}. The corresponding initial value partition are $P_{{\bar{R}}_{1}} (X_{0})$ and $P_{{\bar{R}}_{5}} (X_{0})$ .

Fig. 2

Illustrative example describing a value partition induced from the starting vertices.

4. Application to Decay Chain Graph

The set of differential equations representing radioactive decay with no external neutron flux is shown in Equation (17). For k = 1, …, N

{\dot{x}}_{k} = - λ_{k} x_{k} + \sum_{m \in (1, \dots, N) \ {k}} f_{k, m} \cdot λ_{m} \cdot x_{m}

(17)

In Equation (17), x_k(t) represents the amount of the k^th nuclide at decayed to time t and ẋ_k is its time derivative. λ_k is the decay constants for the k^th nuclide and f_k,m is the moles of the k^th nuclide created from the decay of one mole of the m^th nuclide, i.e., branching fraction. Then, the rate matrix can be constructed as below:

A = [\begin{matrix} - λ_{1} & f_{1, 2} λ_{2} & \dots & f_{1, N} λ_{N} \\ f_{2, 1} λ_{1} & - λ_{2} & \dots & f_{2, N} λ_{N} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ f_{N, 1} λ_{1} & f_{N, 2} λ_{2} & \dots & - λ_{N} \end{matrix}] \in R^{N \times N}

Rate matrix A is constructed with Idaho National Laboratory Mass Tracking System Database built from ENDF/ B-VII database [18]. Decay digraph G = (V, E) that is induced from rate matrix A has the following properties:

|V| = N = 3631 vertices
|E| = 2367 transitions
2063 weakly connected components
Largest weakly connected component: 1299 vertices (mostly actinides)
Smallest component size: 1
Average component size: 1.76
Acyclic

Fig. 3 below shows decay digraph G = (V, E) clustering weakly connected components. Fig. 4 gives the number of reachable vertices from each vertex. Contrary to the size of the largest weakly connected component (1299 vertices), the size of the largest reachable vertex set is substantially small (36 vertices). The transition structure of the graph allures that an algorithm building the subdynamics with the minimum starting vertices is likely to be efficient. Fig. 5 shows the largest decay chain directed graph with 36 vertices (isotopes). The diagram includes neutron as a vertex as well.

Fig. 3

Isotope decay chain graph (3631 isotopes).

Fig. 4

Number of reachable isotopes from each isotope.

Fig. 5

Largest decay chain (36 isotopes counting neutron as isotope).

MATLAB calculations were conducted to compare the performance of the algorithms described in this article.

The monolithic calculation is performed with expm MATLAB function where a scaling and squaring algorithm with a Pade approximation [15,19,20] is used. WCC and MSV algorithms also use the same expm function for their subsystem transition matrix computation.

For the above numerical calculations, the following resources are used:

Processor: Intel® Core™ i7-6600U CPU @2.60GHz
RAM: 16.0GB
MATLAB Version 9.3 (R2017b)

Table 1 and Fig. 6 give the average computational time of the calculation methods described in this article for 100 trials. Given the number of nonzero initial isotopes to be simulated, selection of isotopes was randomized and given the value of 1. The values in the table are relative to the monolithic calculation time with 10 random nonzero initial isotopes. One can observe that the monolithic calculation ignoring the reachable vertices shows a consistent computation time. Once the reachable part of graph is considered, a significant reduction of computational time is observed for all algorithms. All three methods accounting the reachable part of the graph show a gradual increase in computational time as the number of reachable vertices increases. MSV algorithm shows the least steep increase of computational time followed by WCC algorithm and expm with reachability consideration. The superiority of expm function with reachability consideration over WCC and MSV algorithms for the smaller size initially activated vertexes may come from the implementation inefficiency of WCC and MSV algorithms where multiple graph searches were done due to the adaptation of MATLAB generic functions, which can be consolidated with a customized graph search algorithm.

Table 1

Relative average computational time with varying number of nonzero initial variables

Fig. 6

Computation time comparison with respect to the number of activated initial variables.

5. Concluding Remarks

An attempt to view the system of linear, constant real coefficient ordinary differential equations via a graph structure allows various graph algorithm treatments to enable efficient computation of the solution. Connected component analysis, reachability analysis, and the synthesis of the minimal starting vertex set were used to break the problem into the smaller size subproblems. From a theoretical perspective, the worst-case computational complexity does not exhibit an improvement. However, empirical results demonstrated with radioactive decay application show a significant saving of computational efforts. Furthermore, we note that the subproblems are independent and allow parallel computations.

6. Acknowledgements

Rate matrix A in Section 4 is adapted from the MTG decay rate matrix implemented by DeeEarl Vaden.

This submitted manuscript was authored by a contractor of the U.S. Government under DOE Contract No. DE-AC07-05ID14517. Accordingly, the U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for U.S. Government purposes.

This information was prepared as an account of work sponsored by an agency of the U.S. Government. Neither the U.S. Government nor any agency thereof, nor any of their employees, makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness, of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. References herein to any specific commercial product, process, or service by trade name, trade mark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agency thereof.

Figures

Tables

References

H. Bateman, “The solution of a system of differential equations occurring in the theory of radio-active transformations”, Proc. Cambridge Phil. Soc., 1908 15, 423- 427 (1908).
D.R. Vondy, Development of a General Method of Explicit Solution to the Nuclide Chain Equations for Digital Machine Calculations (Thesis), No. ORNL-TM-361, Oak Ridge National Lab., Tenn., 1962.
M.J. Bell, ORIGEN: the ORNL isotope generation and depletion code. No. ORNL--4628. Oak Ridge National Lab., (1973).
J. Cetnar. “General solution of Bateman equations for nuclear transmutations”, Annals of Nuclear Energy, 33 (7), 640-645 (2006).
M. Pusa and J. Leppänen, “Computing the matrix exponential in burnup calculations”, Nuclear science and engineering, 164(2), 140-150 (2010).
A.E. Isotalo and P.A. Aarnio, “Comparison of depletion algorithms for large systems of nuclides”, Annals of Nuclear Energy, 38(2-3), 261-268 (2011).
R. Dreher, “Modified Bateman solution for identical eigenvalues”, Annals of Nuclear Energy, 53, 427-438 (2013).
D. Vaden and T. Yoo, “Radioactive Decay Computation with Dynamic Source Terms”, Nuclear Science and Engineering, 193(5), 549-553 (2019).
J. Serp, M. Allibert, O. Beneš, S. Delpech, O. Feynberg, V. Ghetta, D. Heuer, D. Holcomb, V. Ignatiev, J.L. Kloosterman, and L. Luzzi, “The molten salt reactor (MSR) in generation IV: overview and perspectives”, Progress in Nuclear Energy, 77, 308-319 (2014).
T.J. Dolan ed., Molten Salt Reactors and Thorium Energy, Woodhead Publishing (2017).
H. Lee, G. Park, K. Kang, J. Hur, J. Kim, D. Ahn, Y. Cho, and E. Kim, “Pyroprocessing technology development at KAERI”, Nuclear engineering and technology, 43(4), 317-328 (2011).
T. Inoue, T. Koyama, and Y. Arai, “State of the art of pyroprocessing technology in Japan”, Energy Procedia, 7, 405-413 (2011).
T.Y. Karlsson, G.L. Fredrickson, T. Yoo, D. Vaden, M.N. Patterson, and V. Utgikar, “Thermal analysis of projected molten salt compositions during FFTF and EBR-II used nuclear fuel processing”, Journal of Nuclear Materials, 520, 87-95 (2019).
T. Yoo and D. Vaden, “A new inventory tracking method for Mark-V electrorefiner”, Annals of Nuclear Energy, 128, 406-413 (2019).
C. Moler and C. Van Loan, “Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later”, SIAM review, 45(1), 3-49 (2003).
T.H. Cormen, C.E. Leiserson, R.L. Rivest, and C. Stein, Introduction to algorithms, MIT press (2009).
D. Coppersmith and S. Winograd, “Matrix multiplication via arithmetic progressions”, Journal of symbolic computation, 9(3), 251-280 (1990).
M. B. Chadwick, M. Herman, P. Obložinský, M.E. Dunn, Y. Danon, A.C. Kahler, D.L. Smith, B. Pritychenko, G. Arbanas, R. Arcilla, and R. Brewer, “ENDF/B-VII. 1 nuclear data for science and technology: cross sections, covariances, fission product yields and decay data”, Nuclear data sheets, 112(12), 2887- 2996 (2011).
H. Padé, Sur la représentation approchée d'une fonction par des fractions rationnelles, No. 740, Gauthier- Villars et fils, 1892.
R.C. Ward, “Numerical computation of the matrix exponential with accuracy estimate”, SIAM Journal on Numerical Analysis, 14(4), 600-610 (1977).

JNFCWT

Online Submission

Korean Radioactive
Waste Society (KRS)

Editorial Office
Contact Information

- Tel: +82-42-861-5851, 866-4157
- Fax: +82-42-861-5852
- E-mail: krs@krs.or.kr

SITE MAP