Make use of this iterative scheme for solving the neutron diffusion equation

finite differencesmatrix equationsnumerical linear algebranumerical methods

I am trying to solve the neutron diffusion equation to model neutron flux distribution in a one-dimensional two-group setting using the following iterative scheme. The governing equations of the system are:

$$-D_1\nabla^2\phi_1+(\mathcal{E}_{a_1}+\mathcal{E}_{s12})\phi_1=\frac{1}{k}[\nu_1\mathcal{E}_{f_1}\phi_1+\nu_2\mathcal{E}_{f_2}\phi_2]+\mathcal{E}_{s21}\phi_2$$

$$-D_2\nabla^2\phi_2+(\mathcal{E}_{a_2}+\mathcal{E}_{s21})\phi_2=\mathcal{E}_{s12}\phi_1$$

where $D$=Diffusion Co-efficient, $\phi$=neutron flux, $\mathcal{E}_{a}$=absorption cross-section, $\mathcal{E}_{s}$=scattering cross-section, $\nu$=neutrons per fission, $\mathcal{E}_{f}$=fission cross-section and $k$=multiplication factor. 1 and 2 refers to fast and thermal neutron groups.

I divided the geometry into 420 mesh elements and discretized the equations using Forward Difference Method (FDM). Putting the whole system into matrix form yields:

$$\Bigl[C_1\Bigr]\Bigl[\phi_1\Bigr]=\frac{1}{k}\Bigl[h\Bigr]\Bigl[\nu_1\mathcal{E}_{f_1}\Bigr]\Bigl[\phi_1\Bigr]+\frac{1}{k}\Bigl[h\Bigr]\Bigl[\nu_2\mathcal{E}_{f_2}\Bigr]\Bigl[\phi_2\Bigr]+\Bigl[h\Bigr]\Bigl[\mathcal{E}_{s21}\Bigr]\Bigl[\phi_2\Bigr]$$

$$\Bigl[C_2\Bigr]\Bigl[\phi_2\Bigr]=\Bigl[h\Bigr]\Bigl[\mathcal{E}_{s12}\Bigr]\Bigl[\phi_1\Bigr]$$

Where $\Bigl[C\Bigr]$=Co-efficient Matrix and $h$=mesh element length.

Now, for the $n$th iteration, I am supposed to use $\phi_{1}^{n-1}$, $\phi_{2}^{n-1}$ and $k^{n-1}$ to solve the two above equations and calculate $\phi_{1}^{n}$, $\phi_{2}^{n}$ and $k^{n}$ using the iterative scheme:

$$\Bigl[C_1\Bigr]\Bigl[\phi_1\Bigr]^n=\frac{1}{k}\Bigl[h\Bigr]\Bigl[\nu_1\mathcal{E}_{f_1}\Bigr]\Bigl[\phi_1\Bigr]^{n-1}+\frac{1}{k}\Bigl[h\Bigr]\Bigl[\nu_2\mathcal{E}_{f_2}\Bigr]\Bigl[\phi_2\Bigr]^{n-1}+\Bigl[h\Bigr]\Bigl[\mathcal{E}_{s21}\Bigr]\Bigl[\phi_2\Bigr]^{n-1}$$

$$\Bigl[C_2\Bigr]\Bigl[\phi_2\Bigr]^n=\Bigl[h\Bigr]\Bigl[\mathcal{E}_{s12}\Bigr]\Bigl[\phi_1\Bigr]^n$$

$$k^n=k^{n-1}\frac{\int dr(\nu_1\mathcal{E}_{f_1}\phi_1^n+\nu_2\mathcal{E}_{f_2}\phi_2^n)}{\int dr(\nu_1\mathcal{E}_{f_1}\phi_1^{n-1}+\nu_2\mathcal{E}_{f_2}\phi_2^{n-1})}$$

I am confused about how to proceed with the last equation. Do I take the differential length $dr$ as mesh element length $h$ and turn the integration into a summation over the geometry? Or do I actually try some kind of integration? I cannot imagine how an integration can be introduced here given that $\phi_1$ and $\phi_2$ are not known as functions but as matrices.

Please help! And thanks for reading through such a long question!

Best Answer

You are correct that you have to turn this into a sum over the geometry. You solved for the flux values at the gridpoints and you have the step size between the gridpoints. From this, your choice of integration scheme is limited to those that use uniformly spaced nodes. You could just do a direct sum and multiply by $h$, but in this case you can actually do a bit better by using the trapezoid rule.

Building a test case

Continuous problem with exact solution

$u_t=u_{xx}$, that is, $D=1$, with $u_x(0,t)=0$, $u(\pi,t)=1$ over $t\in[0,1]$.

Separation of variables tells us that terms $e^{-\omega^2t}\cos(\omega t)$ are simple functions that satisfy homogeneous condition automatically, using a constant for the shift gives as solution example $$u(x,t)=1+a_1e^{-t/4}\cos(x/2)+a_3e^{-9t/4}\cos(3x/2).$$

Code snippets (python, using scipy.sparse)

the test solution

a1, a3 = 2,1;
def u_exact(x,t): return 1+a1*np.exp(-0.25*t)*np.cos(0.5*x)+a3*np.exp(-2.25*t)*np.cos(1.5*x)

construction of the system matrices

I = sparse.eye( M, format="csc" )
A = sparse.diags([1, -2, 1], [-1, 0, 1], shape=(M, M), format="csc");
A[0,0]=-2.; A[0,1]=2.;
b = np.zeros(M); b[-1]=1

the time iteration

x = np.linspace(0,np.pi,M+1); dx = x[1]-x[0];
u1 = u_exact(x,0)[:-1];
t = np.linspace(0,T,N+1); dt = t[1]-t[0];
p = dt/dx**2;
for n in range(N):
    u1 = sparse.linalg.spsolve(I-0.5*p*A, (I+0.5*p*A).dot(u1)+p*b)

Numerical results

Using the left side BC treatment like in the question, I get the following table of results. Horizontally $M$ varies over the size of the space subdivision, while vertically $N$ gives the number of time steps.

In each cell the upper value is the error against a run with half the time step, which is an estimate of the error against exact time integration. The lower value is the error against the exact solution. The errors displayed is an approximation of the $L^2$ norm, that is, the Euclidean norm divided by $\sqrt M$. One can see that where $\frac{Δt}{Δx^2}$ gets small, the error against the exact solution becomes constant.

\begin{array}{|c|l|l|l|l|l|}\hline & ~M=10& ~M=50& ~M=200& ~M=1000& ~M=5000 \\ \hline N=10 & \begin{split} 5.5166e-04 \\ 2.6801e-03 \end{split} & \begin{split} 5.3780e-04 \\ 5.9059e-04 \end{split} & \begin{split} 5.3395e-04 \\ 7.0361e-04 \end{split} & \begin{split} 5.3286e-04 \\ 7.0966e-04 \end{split} & \begin{split} 5.3264e-04 \\ 7.0967e-04 \end{split} \\ \hline N=20 & \begin{split} 1.3755e-04 \\ 3.2229e-03 \end{split} & \begin{split} 1.3409e-04 \\ 5.6237e-05 \end{split} & \begin{split} 1.3313e-04 \\ 1.6967e-04 \end{split} & \begin{split} 1.3286e-04 \\ 1.7680e-04 \end{split} & \begin{split} 1.3280e-04 \\ 1.7703e-04 \end{split} \\ \hline N=40 & \begin{split} 3.4365e-05 \\ 3.3587e-03 \end{split} & \begin{split} 3.3500e-05 \\ 8.4139e-05 \end{split} & \begin{split} 3.3260e-05 \\ 3.6555e-05 \end{split} & \begin{split} 3.3193e-05 \\ 4.3943e-05 \end{split} & \begin{split} 3.3179e-05 \\ 4.4224e-05 \end{split} \\ \hline N=80 & \begin{split} 8.5898e-06 \\ 3.3927e-03 \end{split} & \begin{split} 8.3736e-06 \\ 1.1693e-04 \end{split} & \begin{split} 8.3137e-06 \\ 3.5075e-06 \end{split} & \begin{split} 8.2968e-06 \\ 1.0751e-05 \end{split} & \begin{split} 8.2933e-06 \\ 1.1045e-05 \end{split} \\ \hline N=160 & \begin{split} 2.1474e-06 \\ 3.4011e-03 \end{split} & \begin{split} 2.0933e-06 \\ 1.2518e-04 \end{split} & \begin{split} 2.0783e-06 \\ 5.1960e-06 \end{split} & \begin{split} 2.0741e-06 \\ 2.4544e-06 \end{split} & \begin{split} 2.0733e-06 \\ 2.7519e-06 \end{split} \\ \hline N=320 & \begin{split} 5.3683e-07 \\ 3.4033e-03 \end{split} & \begin{split} 5.2333e-07 \\ 1.2725e-04 \end{split} & \begin{split} 5.1958e-07 \\ 7.2298e-06 \end{split} & \begin{split} 5.1852e-07 \\ 3.8307e-07 \end{split} & \begin{split} 5.1832e-07 \\ 6.7867e-07 \end{split} \\ \hline \end{array}

Speculation

Without knowing details of the code used, there are two typical places that come to mind that might introduce a first order error:

a wrong coefficient somewhere, reducing the method to first order. This is unlikely as the matrix and function in general has a simple structure. If there were an error one would expect more grave consequences.
an error in the implementation of the time iteration. Such as
- iterating N times but setting dt=T/(N-1), or
- having T and dt as input and iterating while t<T without caring that the final t can be anywhere between T and T+dt.
Such error results in a time difference between exact and numerical solution of order O(dt) which gives a difference in the function values of also O(dt).

Convergence conditions for an iterative scheme

A little hint.

As I said in comments, consider the eigenvalue basis. The basis vectors are orthogonal and may be scaled to form an orthonormal basis: $$ A \phi_m = \lambda_m \phi_m, \quad m = 1, \dots, m\\ (\phi_m, \phi_{m'}) = \delta_{mm'}. $$

Expanding error vectors over the basis $e_k = \sum_{m=1}^n c_{k,m} \phi_m$ allows to rewrite convergence condition using expansion coefficients. Using Parseval's identity $$ \|e_k\|_2^2 = \sum_{m} c_{k,m}^2 $$ we obtain that $e_k \to 0$ happens only if for all $m$ each coefficient converges to zero, that is $$ \lim_{k \to \infty} c_{k,m} = 0, \quad m = 1,\dots,n. $$

Acting with $(I - \alpha A)^k$ on $e_0$ acts on each eigenvalue separately: $$ e_k = (I - \alpha A)^k e_0 = (I - \alpha A)^k \sum_{m=1}^n c_{0,m} \phi_m = \\ = \sum_{m=1}^n c_{0,m} (I - \alpha A)^k \phi_m = \sum_{m=1}^n c_{0,m} (1 - \alpha \lambda_m)^k \phi_m. $$

Comparing the right hand side with $\sum_{m=1}^n c_{k,m} \phi_m$ we immediately get the relation $$ c_{k,m} = (1 - \alpha \lambda_m)^k c_{0,m}. $$

Now it's up to you to find the conditions when $\lim_{k \to \infty} c_{k,m} = 0$ for every $m = 1,\dots,n$.