Friday, October 10, 2008

Method Sn - Discrete Ordinates with finite Differences method for solving Radiative Transfer Equation

Things are getting clear now.

There are a few methods that can be used to solve RTEs. The Sn Method is supposed to provide a fair balance between computational requirements and accuracy.

From the mathematical perspective, I do not feel confident enough to provide a very detailed explanation. In this case, I will tell you what I know from the computational perspective.

This method uses two special math operations in order to solve the equation. As one may see in the original equation, there are two particularly hard components to implement on a computer. One derivative and one integral.
The first component is not very hard. We can replace the derivative by an approximation. The basic idea of this derivative is finding the intensity of photons that could represent the intensity of photons in a particular region/depth of the water. We can find this information by splitting the studied space into an arbitrary number of regions, we can assume that the derivative value in each region will be equal to the intensity of photons right on the middle of each region. A fair enough approximation would be calculated by a simple average of the photon intensity in each border/interface of the desired region.

Calculating this for each of the regions would give us the derivative for each region.

The integral is much more complicated. On the integral case, we use a numeric quadrature. We replace the single equation by a set of attached differential equations. In this case it is used the Gaussian Quadrature.

The reason why we have several equations instead of just one is because the original equation accumulated the intensity of all photons coming from all angles into a single value. With the quadrature approach, there is a weight representation for each angle. Each new equation represents the photon intensity for each evaluated angle. It is possible to use quadratures with lower or higher degrees. As high the degree of the quadrature is, as high the accuracy is going to be. The degree of the quadrature represents the number of different angles being evaluated for the equation.

After simplifying the general equation by removing the derivative and the integral components, one will get to these two versions of the same equation:

also considering :

After getting to these equations, the general algorithm for solving this problem is:

----- Variables -------
precision <- Desired precision
I <- hold RTE results for each angle in each interface - N (angles) lines by (# interfaces) columns
H <- previous iteration of I

---- Initialization -------
set first column of I with the external source
set last column of I with the Internal Source
set second column of I with 0.5 (*)

------ Body ----------
do {

Copy all values from I to H
for ( each layer interface )
for ( each angle on gaussian quadrature table )
if ( angle cosine > 0)

} while (I-H > precision)

At this point, even if you do not have a clue about Math and can understand algorithms, you should be asking yourself: "OK! On equation (2), I'm calculating the value for [k] interface based on interface [k+1]. How is that possible if interface [k+1] has never been calculated?"

This is a fair question. That is why we have to get back to the initialization of the second column (*) with 0.5 values. This method uses several iterations in order to reach the desired results. For the very first iteration, calculating the very first column, you are right! There is no value on column number 2. We solve this by guessing a number for the second position, any value between the interval ]0 1] (any number between 0 and 1, but 0). This way, the algorithm will be able to run in order to populate the entire matrix. After this initial guessing, the matrix will be entirely populated and the next iteration will use values from previous iterations to keep going.

The stop condition is found when the values from the current iteration, compared to the values from the previous iteration matches the desired precision.

At this point, I'm comfortable with the algorithm itself. This convergence to a single result after many iterations is the part that is still bugging me. On next class, I will spend a lot of time with the professor trying to understand this behavior because for me it still looks like magic. I could not find a good reason yet on why the values will converge to the result we seek.

I will post again when I have more information.