About misleading hints

Hint 9 above is actually somewhat misleading. At some stage I got confused myself by those sizing manipulations, which looked quite puzzling to me (and for a good reason, as it turned out: they have lose ends).

First, the Fourier tranform of function D(x,y)is going to be real, because D(x,y)is even about $\pi$. You can see it as follows:

$$\begin{eqnarray*}\int_0^{2\pi}D(x,y) e^{-ikx}\,\textrm{d}x &=& \int_0^{2\pi}D(x,y) \left(\cos{kx} - i \sin{kx}\right)\,\textrm{d}x \\ \end{eqnarray*}$$

If D(x,y) is even about $\pi$ and $\sin{kx}$, obviously, is odd about $\pi$, then the imaginary part of the integral must vanish. In both functions, expand_temp_profile and in get_diffusion_matrix, subroutine fft is called in the same way, i.e.,

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

This form of calling assumes that both the input matrix to fft and the output, which will be written on the same matrix are going to be complex. In other words, we do not make any use of our knowledge that Fourier Transforms of both the initial temperature distribution and of the diffusion coefficient are going to be real. Instead we simply use the complex-to-complex form of fft and then extract the real part from the result.

The reasons for oversizing the transform have to do with negative wave numbers.

Still the first part of the hint, referring you to the discussion of subroutine fft in the PESSL manual was fine, and if you read that discussion carefully, you would have noticed that both invocations to fft in our program are done in the same, complex-to-complex, way. Consequently the sizing differences associated with special calling conventions for a complex-to-real case wouldn't apply here anyway.

It is instructive to have a look at the Fourier transform of the diffusion coefficient for the following input:

&input
ly_ratio = 1.0d0,
delta_t = 0.2,
numx = 5,
numy = 5,
nx = 7,
ny = 7,
numt = 10,
init_f = 8,
diff_f = 2,
/

Diffusion coefficient model 2 is

$$D(x,y) = \frac{(1 + x)(\pi - x + 1)(y + \pi)}{3 \pi}$$

The values of matrix df before calling fft look as follows:

    df( 1, 1) = (0.165E+02, 0.000E+00)
    df( 1, 2) = (0.175E+02, 0.000E+00)
    df( 1, 3) = (0.184E+02, 0.000E+00)
    df( 1, 4) = (0.194E+02, 0.000E+00)
    df( 1, 5) = (0.204E+02, 0.000E+00)
    df( 1, 6) = (0.214E+02, 0.000E+00)
    df( 1, 7) = (0.223E+02, 0.000E+00)
    df( 1, 8) = (0.233E+02, 0.000E+00)
...
    df(31,25) = (0.223E+02, 0.000E+00)
    df(31,26) = (0.214E+02, 0.000E+00)
    df(31,27) = (0.204E+02, 0.000E+00)
    df(31,28) = (0.194E+02, 0.000E+00)
    df(31,29) = (0.184E+02, 0.000E+00)
    df(31,30) = (0.175E+02, 0.000E+00)
    df(31,31) = (0.165E+02, 0.000E+00)

On exit from fft the non-zero coefficients look as follows:

    df( 2, 1) = (0.339E+00, 0.158E-16)
    df( 2, 3) = (0.386E-01, -.640E-16)
    df( 2, 5) = (0.146E-01, 0.166E-17)
    df( 2, 7) = (0.809E-02, -.433E-19)
    df( 2, 9) = (0.545E-02, -.206E-16)
    df( 2,11) = (0.418E-02, 0.673E-16)
    df( 2,13) = (0.355E-02, -.405E-16)
    df( 2,15) = (0.329E-02, -.555E-16)
    df( 2,17) = (0.329E-02, -.665E-16)
    df( 2,19) = (0.355E-02, -.736E-16)
    df( 2,21) = (0.418E-02, 0.347E-16)
    df( 2,23) = (0.545E-02, -.394E-16)
    df( 2,25) = (0.809E-02, 0.717E-17)
    df( 2,27) = (0.146E-01, 0.396E-16)
    df( 2,29) = (0.386E-01, -.342E-16)
    df( 2,31) = (0.339E+00, 0.231E-16)
    df( 4, 1) = (0.880E-01, -.301E-16)
    df( 4, 3) = (0.100E-01, 0.105E-15)
    df( 4, 5) = (0.381E-02, -.190E-16)
    df( 4, 7) = (0.210E-02, 0.445E-16)
    df( 4, 9) = (0.142E-02, -.327E-16)
    df( 4,11) = (0.109E-02, 0.409E-16)
    df( 4,13) = (0.924E-03, 0.615E-16)
    df( 4,15) = (0.854E-03, -.151E-16)
    df( 4,17) = (0.854E-03, -.167E-16)
    df( 4,19) = (0.924E-03, 0.642E-16)
    df( 4,21) = (0.109E-02, 0.461E-16)
    df( 4,23) = (0.142E-02, -.322E-16)
    df( 4,25) = (0.210E-02, 0.436E-16)
    df( 4,27) = (0.381E-02, -.181E-16)
    df( 4,29) = (0.100E-01, 0.104E-15)
    df( 4,31) = (0.880E-01, -.316E-16)
...

There are quite a few of them, and they are all indeed real (with ``zero'' evaluating to something like $1.0\times10^{-16}$). The input function here is not a simple trigonometric function, so the Fourier spectrum gets scattered all over the place. But we can make things cleaner by adding the following model for the diffusion coefficient:

    CASE (4)
       diff_coef = COS(x1) * COS(y1)

and changing the input namelist to

&input
ly_ratio = 1.0d0,
delta_t = 0.2,
numx = 5,
numy = 5,
nx = 7,
ny = 7,
numt = 10,
init_f = 8,
diff_f = 4,
/

This is a very strange kind of a diffusion coefficient, because it is negative in places. But the Fourier transform of the diffusion coefficient now has terms different from zero only for:

    df( 1, 1) = (0.250E+00, -.573E-17)
    df( 1,31) = (0.250E+00, 0.187E-17)
    df(31, 1) = (0.250E+00, -.632E-19)
    df(31,31) = (0.250E+00, 0.526E-17)

The scale factor here is different than it was for the initial temperature. There we had:

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

whereas here it is

    scale_f = 1.d0/ REAL(nx*ny, long)

which means that the result should be $2\pi$ smaller, i.e.,

$$\frac{\pi}{2} \frac{1}{2 \pi} = \frac{1}{4} = 0.25$$

Last, you may still wonder, why the Fourier coefficients obtained both for the initial temperature distribution and for the diffusion coefficient are positive in this case, even though the scale factor for the initial temperature distribution is negative. The reason for this is that T₀(x,y) is odd about $\pi$, whereas D(x,y) is even about $\pi$. This means that the integral

$$\int_0^{2\pi}T_0(x,y)e^{-ikx}\,\textrm{d}x = \int_0^{2\pi}T_0(x,y)\left( \cos{kx} - i\sin{kx}\right)\,\textrm{d}x$$

will select the imaginary component this time. But, the same will happen in the y direction so the two -i factors will multiply yielding -1, which will produce a positive result on multiplication by $-2\pi/\left(n_x n_y\right)$.