Blow es un archivo principal

PROGRAM SPHEROID

USE nrtype
USE SUB_INFO
INCLUDE "/usr/local/include/fftw3.f"

INTEGER(I8B) :: plan_forward, plan_backward
INTEGER(I4B) :: i, t, int_N

REAL(DP) :: cth_i, sth_i, real_i, perturbation
REAL(DP) :: PolarEffect, dummy, x1, x2, x3

REAL(DP), DIMENSION(4096) :: dummy1, dummy2, gam, th, ph
REAL(DP), DIMENSION(4096) ::  k1, k2, k3, k4, l1, l2, l3, l4, f_in

COMPLEX(DPC), DIMENSION(2049) :: output1, output2, f_out

CHARACTER(1024) :: baseOutputFilename
CHARACTER(1024) :: outputFile, format_string

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

int_N = 4096

! File Open Section

format_string = '(I5.5)'

! Write the coodinates at t = 0

do i = 1, N
    real_i = real(i)
    gam(i) = 2d0*pi/real_N
    perturbation = 0.01d0*dsin(2d0*pi*real_i/real_N)
    ph(i) = 2d0*pi*real_i/real_N + perturbation
    th(i) = pi/3d0 + perturbation
end do

! Initialization Section for FFTW PLANS

call dfftw_plan_dft_r2c_1d(plan_forward, int_N, f_in, f_out, FFTW_ESTIMATE)
call dfftw_plan_dft_c2r_1d(plan_backward, int_N, f_out, f_in, FFTW_ESTIMATE)

! Runge-Kutta 4th Order Method Section

do t = 1, Iter_N

    call integration(th, ph, gam, k1, l1)

    do i = 1, N
        dummy1(i) = th(i) + 0.5d0*dt*k1(i)
    end do
    do i = 1, N
        dummy2(i) = ph(i) + 0.5d0*dt*l1(i)
    end do

    call integration(dummy1, dummy2, gam, k2, l2)

    do i = 1, N
        dummy1(i) = th(i) + 0.5d0*dt*k2(i)
    end do
    do i = 1, N
        dummy2(i) = ph(i) + 0.5d0*dt*l2(i)
    end do

    call integration(dummy1, dummy2, gam, k3, l3)

    do i = 1, N
        dummy1(i) = th(i) + dt*k3(i)
    end do
    do i = 1, N
        dummy2(i) = ph(i) + dt*l3(i)
    end do

    call integration(dummy1, dummy2, gam, k4, l4)

    do i = 1, N

        cth_i = dcos(th(i))
        sth_i = dsin(th(i))
        PolarEffect = (nv-sv)*dsqrt(1d0+a*sth_i**2) + (nv+sv)*cth_i
        PolarEffect = PolarEffect/(sth_i**2)
        th(i) = th(i) + dt*(k1(i) + 2d0*k2(i) + 2d0*k3(i) + k4(i))/6d0
        ph(i) = ph(i) + dt*(l1(i) + 2d0*l2(i) + 2d0*l3(i) + l4(i))/6d0
        ph(i) = ph(i) + dt*0.25d0*PolarEffect/pi

    end do

!! Fourier Filtering Section

    call dfftw_execute_dft_r2c(plan_forward, th, output1)

    do i = 1, N/2+1
       dummy = abs(output1(i))
        if (dummy.lt.threshhold) then
           output1(i) = dcmplx(0.0d0)
        end if
    end do

    call dfftw_execute_dft_c2r(plan_backward, output1, th)

    do i = 1, N
       th(i) = th(i)/real_N
    end do

    call dfftw_execute_dft_r2c(plan_forward, ph, output2)

    do i = 1, N/2+1
       dummy = abs(output2(i))
        if (dummy.lt.threshhold) then
           output2(i) = dcmplx(0.0d0)
        end if
    end do

    call dfftw_execute_dft_c2r(plan_backward, output2, ph)

    do i = 1, N
       ph(i) = ph(i)/real_N
    end do

!! Data Writing Section

    write(baseOutputFilename, format_string) t
    outputFile = "xyz" // baseOutputFilename
    open(unit=7, file=outputFile)
    outputFile = "Fsptrm" // baseOutputFilename
    open(unit=8, file=outputFile)

    do i = 1, N
        x1 = dsin(th(i))*dcos(ph(i))
        x2 = dsin(th(i))*dsin(ph(i))
        x3 = dsqrt(1d0+a)*dcos(th(i))
        write(7,*) x1, x2, x3
    end do

    do i = 1, N/2+1
        write(8,*) abs(output1(i)), abs(output2(i))
    end do

    close(7)
    close(8)

    do i = 1, N/2+1
        output1(i) = dcmplx(0.0d0)
    end do

    do i = 1, N/2+1
        output2(i) = dcmplx(0.0d0)
    end do

end do

! Destroying Process for FFTW PLANS

call dfftw_destroy_plan(plan_forward)
call dfftw_destroy_plan(plan_backward)

END PROGRAM

A continuación se muestra un archivo de subrutina para integración

! We implemented Shelly's spectrally accurate convergence method

SUBROUTINE integration(in1,in2,in3,out1,out2)

USE nrtype
USE SUB_INFO

INTEGER(I4B) :: i, j
REAL(DP) :: th_i, th_j, gi, ph_i, ph_j, gam_j, v1, v2
REAL(DP), DIMENSION(N), INTENT(INOUT) :: in1, in2, in3, out1, out2
REAL(DP) :: ui, uj, part1, part2, gj, cph, sph
REAL(DP) :: denom, numer, temp

do i = 1, N
    out1(i) = 0d0
end do
do i = 1, N
    out2(i) = 0d0
end do

do i = 1, N

    th_i = in1(i)
    ph_i = in2(i)
    ui = dcos(th_i)
    part1 = dsqrt(1d0+a)/(dsqrt(-a)*ui+dsqrt(1d0+a-a*ui*ui))
    part1 = part1**(dsqrt(-a))
    part2 = (dsqrt(1d0+a-a*ui*ui)+ui)/(dsqrt(1d0+a-a*ui*ui)-ui)
    part2 = dsqrt(part2)

    gi = dsqrt(1d0-ui*ui)*part1*part2

    do j = 1, N

        if (mod(i+j,2).eq.1) then

            th_j = in1(j)
            ph_j = in2(j)
            gam_j = in3(j)

            uj = dcos(th_j)
            part1 = dsqrt(1d0+a)/(dsqrt(-a)*uj+dsqrt(1d0+a-a*uj*uj))
            part1 = part1**(dsqrt(-a))
            part2 = (dsqrt(1d0+a-a*uj*uj)+uj)/(dsqrt(1d0+a-a*uj*uj)-uj)
            part2 = dsqrt(part2)

            gj = dsqrt(1d0-ui*ui)*part1*part2

            cph = dcos(ph_i-ph_j)
            sph = dsin(ph_i-ph_j)

            numer = dsqrt(1d0-uj*uj)*sph
            denom = (gj/gi*(1d0-ui*ui) + gi/gj*(1d0-uj*uj))*0.5d0
            denom = denom - dsqrt((1d0-ui*ui)*(1d0-uj*uj))*cph
            denom = denom + krasny_delta
            v1 = -0.25d0*gam_j*numer/denom/pi

            temp = dsqrt(1d0+(1d0-ui*ui)*a)
            numer = -(gj/gi)*(temp+ui)
            numer = numer + (gi/gj)*((1d0-uj*uj)/(1d0-ui*ui))*(temp-ui)
            numer = numer + 2d0*ui*dsqrt((1d0-uj*uj)/(1d0-ui*ui))*cph
            numer = 0.5d0*numer
            v2 = -0.25d0*gam_j*numer/denom/pi

            out1(i) = out1(i) + 2d0*v1
            out2(i) = out2(i) + 2d0*v2

        end if

    end do

end do

END

Debajo hay un archivo de módulo

module nrtype
Implicit none
!integer
INTEGER, PARAMETER :: I8B = SELECTED_INT_KIND(20)
INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4)
INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2)
!real
INTEGER, PARAMETER :: SP = KIND(1.0)
INTEGER, PARAMETER :: DP = KIND(1.0D0)
!complex
INTEGER, PARAMETER :: SPC = KIND((1.0,1.0))
INTEGER, PARAMETER :: DPC = KIND((1.0D0,1.0D0))
!defualt logical
INTEGER, PARAMETER :: LGT = KIND(.true.)
!mathematical constants
REAL(DP), PARAMETER :: pi = 3.141592653589793238462643383279502884197_dp
!derived data type s for sparse matrices,single and double precision
!User-Defined Constants
INTEGER(I4B), PARAMETER :: N = 4096, Iter_N = 20000
REAL(DP), PARAMETER :: real_N = 4096d0
REAL(DP), PARAMETER :: a = -0.1d0, dt = 0.001d0, krasny_delta = 0.01d0
REAL(DP), PARAMETER :: nv = 0d0, sv = 0d0, threshhold = 0.00000000001d0
!N : The Number of Point Vortices, Iter_N * dt = Total time, dt : Time Step
!krasny_delta : Smoothing Parameter introduced by R.Krasny
!nv : Northern Vortex Strength, sv : Southern Vortex Strength
!a : The Eccentricity in the direction of z , threshhold : Filtering Threshhold
end module nrtype

A continuación se muestra un archivo de información de subrutina

MODULE SUB_INFO

INTERFACE
  SUBROUTINE integration(in1,in2,in3,out1,out2)
  USE nrtype
  INTEGER(I4B) :: i, j
  REAL(DP) :: th_i, th_j, gi, ph_i, ph_j, gam_j, v1, v2
  REAL(DP), DIMENSION(N), INTENT(INOUT) :: in1, in2, in3, out1, out2
  REAL(DP) :: ui, uj, part1, part2, gj, cph, sph
  REAL(DP) :: denom, numer, temp
  END SUBROUTINE
END INTERFACE

END MODULE

Los compilé usando el siguiente comando

gfortran -o p0 -fbounds-check nrtype.f90 spheroid_sub_info.f90 spheroid_sub_integration.f90 spheroid_main.f90 -lfftw3 -lm -Wall -pedantic -pg

nohup ./p0 &

Tenga en cuenta que 2049 = 4096/2 + 1

Al hacer plan_backward, ¿no es correcto que usemos 2049 en lugar de 4096 ya que la dimensión de salida es 2049?

Pero cuando hago eso, explota. (Volar significa error NAN)

Si uso 4096 para hacer plan_backward, todo está bien, excepto que algunos coeficientes de Fourier son anormalmente grandes, lo que no debería suceder.

Por favor, ayúdame a usar FFTW en Fortran correctamente. Este problema me ha desanimado durante mucho tiempo.