- #1
Telemachus
- 835
- 30
Hi. I have written a code to solve the one dimensional one group (a fixed velocity is considered for the particles) time independent transport equation. The code uses the ##S_N## discrete ordinates method, a Gauss-Legendre quadrature in the angular directions, and a Diamond Difference formula for the spatial discretization.
The one dimensional transport equation is:
##\mu \frac{d\Psi(x,\mu)}{dx}+\sigma_T \Psi(x,\mu)=\sigma_s \int \Psi(x,\mu') d\mu' + q(x,\mu)##.I am interested in solving it in the domain ##[0,10]## with vacuum boundary conditions:
##\Psi(0,\mu)=0## for all ##\mu>0##
##\Psi(10,\mu)=0## for all ##\mu<0##
In the ##S_N## discrete ordinates method the equation is solved in the discrete directions ##\mu_m##, which for my discretization are given by the Gauss-Legendre abscissas.
The scalar flux is then ##\phi(x)=\int \Psi(x,\mu') d\mu' \sim \sum_n \Psi_n(x) w_n##
Where ##\Psi_n(x)=\Psi(x,\mu_n)## and ##w_n## are the weights in the Gauss-Legendre quadrature.
The Diamond discretization subdivided the domain in meshes, and there is an equation that connects the values in the mesh boundaries, and the value for the mesh centers is the average value of the value at the boundaries.
The complete discretized equation reads:
##\frac{\mu_n}{\Delta x} (\Psi_n,{j+1/2}-\Psi_n,{j-1/2})+\sigma_T \Psi_n,j=\sigma_s \sum_m \Psi_m w_m + q_n,j##.
With the auxiliary relation that connects the boundaries with the mesh centers:
##\Psi_n,j=\frac{1}{2}(\Psi_n,{j+1/2}+\Psi_n,{j-1/2})##.
Replacing this last equation in the one above it, one obtain a discrete relation which together with the boundary conditions provides all the relations needed to solve the discretized equation. I am using an iteration scheme that starts with an initial value of zero for the function inside the integral, and iterates until convergence.
Finally, the relation reads for the iteration step l+1:
##\Psi^{l+1}_n,{j+1/2}=\frac{(2\mu_n-\sigma_T\Delta x)\Psi_n,{j-1/2}+2 \Delta x(\sigma_s \Phi^l_j+q_j)}{2\mu_n+\sigma_t \Delta x}## for ##\mu_n>0##
And:
##\Psi^{l+1}_n,{j-1/2}=\frac{(-2\mu_n-\sigma_T\Delta x)\Psi_n,{j-1/2}+2 \Delta x(\sigma_s \Phi^l_j+q_j)}{-2\mu_n+\sigma_t \Delta x}## for ##\mu_n<0##.
I have written a code in Fortran77 to solve this problem using the steps described above. I've been trying to reproduce the results in this paper (page 312 and 313), which uses ##q=0.01, \sigma_T=100,\sigma_s=100##: http://www.sciencedirect.com/science/article/pii/0021999187901707
However, my code doesn't converge to the solution. I have tried with manufactured solutions, and the code converges. But I don't know why it doesn't converge in this particular case. I know my solution is wrong, because in the first place I don't obtain what is in this paper, and in the second place because if I try to recover the source using finite difference in my solution, I don't obtain the source I have put in there.
The one dimensional transport equation is:
##\mu \frac{d\Psi(x,\mu)}{dx}+\sigma_T \Psi(x,\mu)=\sigma_s \int \Psi(x,\mu') d\mu' + q(x,\mu)##.I am interested in solving it in the domain ##[0,10]## with vacuum boundary conditions:
##\Psi(0,\mu)=0## for all ##\mu>0##
##\Psi(10,\mu)=0## for all ##\mu<0##
In the ##S_N## discrete ordinates method the equation is solved in the discrete directions ##\mu_m##, which for my discretization are given by the Gauss-Legendre abscissas.
The scalar flux is then ##\phi(x)=\int \Psi(x,\mu') d\mu' \sim \sum_n \Psi_n(x) w_n##
Where ##\Psi_n(x)=\Psi(x,\mu_n)## and ##w_n## are the weights in the Gauss-Legendre quadrature.
The Diamond discretization subdivided the domain in meshes, and there is an equation that connects the values in the mesh boundaries, and the value for the mesh centers is the average value of the value at the boundaries.
The complete discretized equation reads:
##\frac{\mu_n}{\Delta x} (\Psi_n,{j+1/2}-\Psi_n,{j-1/2})+\sigma_T \Psi_n,j=\sigma_s \sum_m \Psi_m w_m + q_n,j##.
With the auxiliary relation that connects the boundaries with the mesh centers:
##\Psi_n,j=\frac{1}{2}(\Psi_n,{j+1/2}+\Psi_n,{j-1/2})##.
Replacing this last equation in the one above it, one obtain a discrete relation which together with the boundary conditions provides all the relations needed to solve the discretized equation. I am using an iteration scheme that starts with an initial value of zero for the function inside the integral, and iterates until convergence.
Finally, the relation reads for the iteration step l+1:
##\Psi^{l+1}_n,{j+1/2}=\frac{(2\mu_n-\sigma_T\Delta x)\Psi_n,{j-1/2}+2 \Delta x(\sigma_s \Phi^l_j+q_j)}{2\mu_n+\sigma_t \Delta x}## for ##\mu_n>0##
And:
##\Psi^{l+1}_n,{j-1/2}=\frac{(-2\mu_n-\sigma_T\Delta x)\Psi_n,{j-1/2}+2 \Delta x(\sigma_s \Phi^l_j+q_j)}{-2\mu_n+\sigma_t \Delta x}## for ##\mu_n<0##.
I have written a code in Fortran77 to solve this problem using the steps described above. I've been trying to reproduce the results in this paper (page 312 and 313), which uses ##q=0.01, \sigma_T=100,\sigma_s=100##: http://www.sciencedirect.com/science/article/pii/0021999187901707
However, my code doesn't converge to the solution. I have tried with manufactured solutions, and the code converges. But I don't know why it doesn't converge in this particular case. I know my solution is wrong, because in the first place I don't obtain what is in this paper, and in the second place because if I try to recover the source using finite difference in my solution, I don't obtain the source I have put in there.