Discussion of fourier.f

• In order to understand what really happens inside get_diffusion_matrix and expand_temp_profile, we must learn more about PESSL subroutine fft, which has some very special requirements.  
• In order to evaluate
  
  $$\int_0^{2\pi}\int_0^{2\pi}e^{-ikx}e^{-ijy}T_0(x,y) \,\textrm{d} x\,\textrm{d} y$$
  
  
  and
  
  $$\int_0^{2\pi}\int_0^{2\pi}e^{-ikx}e^{-ijy}D(x,y) \,\textrm{d} x\,\textrm{d} y$$
  
  
  we must provide subroutine fft with values of T₀(x,y) and D(x,y) on a sufficiently dense $x\times y$ grid. We could use the same grid for this purpose as the grid $\left(x_m, y_n\right)$ which we have already defined in main in order to compute values of T(x,y, t) at those points. But those two grids don't have to be the same. Furthermore, there are special requirements imposed by subroutine fft, as to the dimensions of its own grid. And these are as follows:

  1.

  Let n_x and n_y be the number of grid points in x and y directions.

  2.

  n_x and n_y must be a multiple of a number of processors requested for the computation.

  3.

  n_x and n_y must be no greater than 37748736 and must be of the form:
  

  $$n = 2^h\times 3^i\times 5^j\times 7^k\times 11^m$$

  
  

  where h, i, j, k, and m are integers and
  

  $$\begin{eqnarray*}h &\in& [1\ldots 25] \\ i &\in& [0\ldots 2] \\ j &\in& [0, 1] \\ k &\in& [0, 1] \\ m &\in& [0, 1] \end{eqnarray*}$$

  
  

• The purpose of subroutine factor_nodes and function minpower2 is to ensure that these very special conditions are satisfied and that a suitable grid can be generated taking our original numbers n_x and n_y, which correspond to the requested numbers of harmonics in x and y, as a starting point. We shall discuss how this is done exactly towards the end of this section. For the time being let us simply assume that it's done and that we have tweaked our input n_x and n_y so that they're now suitable. In other words, for the FFT computation we generate a physical grid which has the same number of points in x and y as the requested number of harmonics in these directions.
• Subroutine fft computes the following:
  
  $$\alpha \sum_{l=0}^{n_x - 1}\sum_{m=0}^{n_y - 1} D\left(x_l, y_m\right) e^{- i k (2\pi/n_x) l} e^{- i j (2\pi/n_y) m}$$
  
  
  Now you can see why we have scaled our problem so that it would fit into $[0, 2\pi]\times[0, 2\pi]$.
• How to relate this double sum to:
  
  $$\int_0^{2\pi}\int_0^{2\pi} D(x,y)e^{-ikx}e^{-ijy} \,\textrm{d}x\,\textrm{d}y$$
  
  
  The answer is in the scaling factor $\alpha$, which must be set so that
  
  $$\alpha = \textrm{d}x\,\textrm{d}y$$
  
  
  We have n_x little segments between 0 and $2\pi$ and ditto for n_y. Therefore
  
  $$\begin{eqnarray*}\textrm{d}x &=& \frac{2\pi}{n_x}\\ \textrm{d}y &=& \frac{2\pi}{n_y} \end{eqnarray*}$$
  
  Consequently
  
  $$\alpha = \frac{4\pi^2}{n_x n_y}$$
  
  
  If we make our $\alpha = 1/\left(n_x n_y\right)$ then we'll evaluate
  
  $$\frac{1}{4\pi^2}\int_0^{2\pi}\int_0^{2\pi} D(x,y)e^{-ikx}e^{-ijy} \,\textrm{d}x\,\textrm{d}y$$
  
  
  instead.
• This is exactly what happens below:

      scale_f = 1.d0/ REAL(nx*ny, long)
      ALLOCATE(df(nx,ny))
      DO j = 1, ny
         DO i =1,  nx
            df(i,j) = diff_coef((twopi*(i-1))/nx,(twopi*(j-1))/ny)
         ENDDO
      ENDDO                 
      CALL fft(df,scale=scale_f,transpose='N')
  

  After having determined the scaling factor (our $\alpha$, theirs scale_f), we evaluate D(x_i, y_j) for $i \in [1, n_x]$ and $j \in [1, n_y]$, save the result on an auxiliary variable df (which is distributed the same way as f, i.e., (*,BLOCK), but has different dimensions, i.e., $n_x\times n_y$, whereas f has dimensions $n_x n_y\,\times\,n_x n_y$), and load the lot into fft.
• Subroutine fft will return the numbers that correspond to the Fourier transforms for $k\in[1\ldots n_x]$ and $j\in[1\ldots n_y]$, i.e.,
  
  $$\frac{1}{4\pi^2}\hat{D}_{kj} = \frac{1}{4\pi^2}\int_0^{2\pi}\int_0^{2\pi} D(x,y)e^{-ikx}e^{-ijy}\,\textrm{d}x\,\textrm{d}y,$$
  
  
  on the same matrix, df.
• Our next task is to evaluate matrix F following the formula:
  
  $$\begin{eqnarray*}F_{kjk'j'} &=& \frac{1}{4\pi^2}\big( \left(k k' + l_r^2 j j'\... ...'\right) \left(\hat{D}_{k^+j^-} - \hat{D}_{k^-j^+}\right)\big) \end{eqnarray*}$$
  
  while at the same time collapsing the double indexes into a single index according to
  
  $$\begin{eqnarray*}l &=& (j - 1) n_x + k\\ j &=& \frac{l - 1}{n_x} + 1, \> k = \hbox{mod}(l - 1, n_x) + 1 \end{eqnarray*}$$
  
  
• But instead of using
  
  $$j = \frac{l - 1}{n_x} + 1 \>\hbox{and}\> k = \hbox{mod}(l - 1, n_x) + 1$$
  
  
  we build two arrays ixi and iyi which simply remember values of j and k that correspond to a given l. This is how it's done:

      num_errors=0
      ALLOCATE(ixi(dif_nx*dif_ny), stat=istat)
      IF( istat .NE. 0 ) num_errors = num_errors + 1
      ALLOCATE(iyi(dif_nx*dif_ny), stat=istat)
      IF( istat .NE. 0 ) num_errors = num_errors + 1
      DO ix = 1, dif_nx
         DO iy = 1, dif_ny
            i = (iy-1)* dif_nx + ix
            ixi(i) = ix
            iyi(i) = iy
         ENDDO
      ENDDO
  

• Now we can use these arrays in calculating Fⁱ_j. The code below has a double-nested DO loop, whose outer index runs over j and the inner index runs over i. At the beginning of every loop we find a pair (k, l) that corresponds to a given j or i by looking up arrays iyp and ixp:

      DO j = 1, dif_nx*dif_ny
         iyp = iyi(j)
         ixp = ixi(j)
         DO i = 1, dif_nx*dif_ny
            iy = iyi(i)
            ix = ixi(i)
  

• Then we evaluate k - k' and j - j':

            ixdiff = iabs(ix-ixp) + 1
            iydiff = iabs(iy-iyp) + 1
  

  Observe that we evaluate the absolute values of these (remember the argument that F is symmetric) and then add 1. Why? That's because the fft subroutine returns results for $k\in[0, n_x-1]$ and $j\in[0, n_y-1]$, but in our case $k\in[1, n_x]$ and $j\in[1,ny]$. So, for this reason we're shifting the indexes in the matrix returned by fft by 1.
• The matrix Fⁱ_j is now calculated as follows:

            f(i,j) = ( ( ix*ixp + iy*iyp*dif_ly_ratio ) *   &
                 REAL(df(ixdiff, iydiff)    - df(ix+ixp+1,iy+iyp+1)) &
                 +     ( ix*ixp - iy*iyp*dif_ly_ratio ) * &
                 REAL(df(ix+ixp+1,iydiff )  - df(ixdiff,iy+iyp+1)))
  

• Before returning from the subroutine we free space that was allocated for auxiliary variables:

      DEALLOCATE(df)
      DEALLOCATE(ixi)
      DEALLOCATE(iyi)
      RETURN
    END SUBROUTINE get_diffusion_matrix
  

• Let us now have a look at the next subroutine, which evaluates a₀, expand_temp_profile.
• The formula for a₀ is:
  
  $$a_{kj}(0) = \frac{-1}{2\pi} \int_0^{2\pi}\int_0^{2\pi}e^{-ikx}e^{-ijy}T_0(x,y) \,\textrm{d} x\,\textrm{d} y$$
  
  

• This means that we have to scale our $\alpha$ parameter in order to get the $-1/(2\pi)$ factor. Setting $\alpha = 1/\left(n_x\times n_y\right)$ gave us $1/\left(4\pi^2\right)$, so multiplying $\alpha$ by $-2\pi$ should yield exactly what's needed:

      scale_f = -twopi / REAL( nx*ny,long)
  

• Now we create a temporary array, called atmp, which we're going to initialise to $T_0\left(x_i, y_j\right)$ on a grid in the (x, y) space with dimensions already made to suit subroutine fft:

      ALLOCATE(atmp(nx,ny), stat=istat )
      IF( istat .NE. 0 )  THEN
         WRITE(*,*) 'allocation failure in expand_temp_profile'
         STOP
      ENDIF
      DO i = 1, nx
         DO j = 1, ny
            atmp(i,j) = init_temp((twopi*(i-1))/nx, (twopi*(j-1))/ny)
         ENDDO
      ENDDO
  

• The next step is simply to call subroutine fft passing atmp and scale_f to it:

      CALL fft(atmp,scale=scale_f,transpose='N')
  

  The coefficients a_kj(0) will be returned on atmp.
• The last step is to collapse the indexes i and j onto a single index and transfer appropriate array entries from atmp to a:

      DO j =  1, dif_nx * dif_ny
         jx = MOD(j-1, dif_nx) + 2         
         jy = (j-1) / dif_nx + 2
         a(j) = REAL(atmp(jx,jy),long)
      ENDDO
  

  Observe that this time we do it differently than in subroutine get_diffusion_matrix.
• Before returning to the caller we free space used by atmp:

      DEALLOCATE(atmp)
      RETURN