Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
Announcements
27660 Discussions

## NLS Equation using split step method in fortran

Beginner
192 Views

Hi Everyone,

I am trying to solve nonlinear schrodinger equation by spilt step method in fortran95. The equation is : ∂U/∂ξ=- i s/2  ∂^2 U/∂τ^ 2+ i N^ 2 |U |^ 2 U

here is my code in below:

program nls
use mkl_dfti !, FORGET=>DFTI_DOUBLE, DFTI_DOUBLE=>DFTI_DOUBLE_R

implicit none

real*8 :: length , dispersion, N , m , chirp, nt, Tmax,Input
real*8:: step_num, deltaz, dtau
real*8, parameter :: pi= 3.14159265358979323d0
integer, parameter :: j=1024, k=1024
integer :: rti, o, i
integer :: status
real*8, dimension(1024)  :: omega
real*8, dimension(1024):: tau
complex*16, dimension(1024):: disp ,  temp, f_temp,  Out,Output, F_D_I, F_D_O, u
complex*16 :: hhz, A,B=(0,1)

type(DFTI_DESCRIPTOR), POINTER :: hand

Status = DftiCreateDescriptor(hand,DFTI_double, DFTI_complex, 1, j)
Status = DftiSetValue (hand,DFTI_FORWARD_SCALE,1.0d0 )
Status = DftiSetValue (hand,DFTI_BACKWARD_SCALE, 1.0d0/dble(j))
Status = DftiCommitDescriptor(hand)

!print*, 'Enter fiber length ='
length = 100.0d0
!print*, 'dispersion 1 for normal, -1 for anomalous'
dispersion = 1.0d0
!print*, 'Nonlinear parameter N ='
N= 1.0d0
!print*, 'm = 0 for sech , m>0 for super-gaussian='
m =0.0d0
chirp =0.0d0

! print*, 'Enter Sampling point nt='
nt =1024.0d0
! print*, 'Enter window size Tmax='
Tmax= 32.0d0

step_num= 1000.0d0*length*(N**2)
deltaz= length/step_num
dtau= (2.0d0*Tmax)/nt

rti= int(nt)

!!=============

do i= 1,int(nt)
tau(i)=dble(i-nt/2-1)*dtau

if (i<= 512) then

omega(i)=dble(i-1)* pi/Tmax

else
omega(i)=dble(i-nt -1)* pi/Tmax

end if

!u(i)= 1d0/dcosh(tau(i))

end do

u = 1.0d0/dcosh(tau)

hhz= (N**2)*deltaz*B
A= hhz/2.0d0

do i = 1, 1024

temp(i)=  u(i)* exp(abs(u(i))**2 * A)

end do

do i =1, 1024

disp(i) = exp(0.50d0*B*dispersion*omega(i)**2*deltaz)

end do

do i= 1,1024
open(1, file='Input_time_domain.dat')
write(1,*) tau(i), abs(u(i))**2
end do

!!=========================================!!

!!-==========Main Program===============!!

do o = 1, int(step_num)

Status = DftiComputeBackward (hand,temp)

do i = 1,1024

f_temp(i) = temp(i)*disp(i)

end do

Status = DftiComputeForward (hand,f_temp)

do i = 1, 1024

Out(i) = f_temp(i) *exp(abs(f_temp(i))**2*hhz)

end do

end do

do i = 1, 1024

Output(i)= Out(i)*exp(-abs(f_temp(i))**2*A)

open(2, file='Output_time_domain.dat')
write(2,*) tau(i), abs(Output(i))**2

end do

!!========== FRequency domain input and Output file======!!

Status = DftiComputeBackward (hand, u)

do i = 1, 1024
F_D_I(i)  = u(i)*((nt*dtau)/dsqrt(2.0d0*pi))
end do

Status = DftiComputeBackward (hand, Output)

do i = 1,1024
F_D_O(i)= Output(i)*((nt*dtau)/sqrt(2.0d0*pi))
end do

do i =1, 1024

open(3, file='Input_Frequency_domain.dat')
write(3,*) omega(i)/2.0d0*pi, abs(F_D_I(i))**2

open(4, file='Output_Frequency_domain.dat')
write(4,*)omega(i)/(2d0*pi), abs(F_D_O(i))**2

end do

Status = DftiFreeDescriptor (hand)

end program nls

===============================

I face problem in propagation calculation portion, where I suppose to go from 1 to step_num, and get dispersion effect. but in my program

as in:  do o = 1, int(step_num)

Status = DftiComputeBackward (hand,temp)

do i = 1,1024

f_temp(i) = temp(i)*disp(i)

end do

Status = DftiComputeForward (hand,f_temp)

do i = 1, 1024

Out(i) = f_temp(i) *exp(abs(f_temp(i))**2*hhz)

end do

end do

==== this do loop is i calculate from 1 to 1 in do loop as highlited

```program nls
use mkl_dfti !, FORGET=>DFTI_DOUBLE, DFTI_DOUBLE=>DFTI_DOUBLE_R

implicit none

real*8 :: length , dispersion, N , m , chirp, nt, Tmax,Input
real*8:: step_num, deltaz, dtau
real*8, parameter :: pi= 3.14159265358979323d0
integer, parameter :: j=1024, k=1024
integer :: rti, o, i
integer :: status
real*8, dimension(1024)  :: omega
real*8, dimension(1024):: tau
complex*16, dimension(1024):: disp ,  temp, f_temp,  Out,Output, F_D_I, F_D_O, u
complex*16 :: hhz, A,B=(0,1)

type(DFTI_DESCRIPTOR), POINTER :: hand

Status = DftiCreateDescriptor(hand,DFTI_double, DFTI_complex, 1, j)
Status = DftiSetValue (hand,DFTI_FORWARD_SCALE,1.0d0 )
Status = DftiSetValue (hand,DFTI_BACKWARD_SCALE, 1.0d0/dble(j))
Status = DftiCommitDescriptor(hand)

!print*, 'Enter fiber length ='
length = 100.0d0
!print*, 'dispersion 1 for normal, -1 for anomalous'
dispersion = 1.0d0
!print*, 'Nonlinear parameter N ='
N= 1.0d0
!print*, 'm = 0 for sech , m>0 for super-gaussian='
m =0.0d0
chirp =0.0d0

! print*, 'Enter Sampling point nt='
nt =1024.0d0
! print*, 'Enter window size Tmax='
Tmax= 32.0d0

step_num= 1000.0d0*length*(N**2)
deltaz= length/step_num
dtau= (2.0d0*Tmax)/nt

rti= int(nt)

!!=============

do i= 1,int(nt)
tau(i)=dble(i-nt/2-1)*dtau

if (i<= 512) then

omega(i)=dble(i-1)* pi/Tmax

else
omega(i)=dble(i-nt -1)* pi/Tmax

end if

!u(i)= 1d0/dcosh(tau(i))

end do

u = 1.0d0/dcosh(tau)

hhz= (N**2)*deltaz*B
A= hhz/2.0d0

do i = 1, 1024

temp(i)=  u(i)* exp(abs(u(i))**2 * A)

end do

do i =1, 1024

disp(i) = exp(0.50d0*B*dispersion*omega(i)**2*deltaz)

end do

do i= 1,1024
open(1, file='Input_time_domain.dat')
write(1,*) tau(i), abs(u(i))**2
end do

!!=========================================!!

!!-==========Main Program===============!!

do o = 1, int(step_num)

Status = DftiComputeBackward (hand,temp)

do i = 1,1024

f_temp(i) = temp(i)*disp(i)

end do

Status = DftiComputeForward (hand,f_temp)

do i = 1, 1024

Out(i) = f_temp(i) *exp(abs(f_temp(i))**2*hhz)

end do

end do

do i = 1, 1024

Output(i)= Out(i)*exp(-abs(f_temp(i))**2*A)

open(2, file='Output_time_domain.dat')
write(2,*) tau(i), abs(Output(i))**2

end do

!!========== FRequency domain input and Output file======!!

Status = DftiComputeBackward (hand, u)

do i = 1, 1024
F_D_I(i)  = u(i)*((nt*dtau)/dsqrt(2.0d0*pi))
end do

Status = DftiComputeBackward (hand, Output)

do i = 1,1024
F_D_O(i)= Output(i)*((nt*dtau)/sqrt(2.0d0*pi))
end do

do i =1, 1024

open(3, file='Input_Frequency_domain.dat')
write(3,*) omega(i)/2.0d0*pi, abs(F_D_I(i))**2

open(4, file='Output_Frequency_domain.dat')
write(4,*)omega(i)/(2d0*pi), abs(F_D_O(i))**2

end do

Status = DftiFreeDescriptor (hand)

end program nls
```

, then data is ok, but if it increase from 1 to 2 and above it does not give any data. plesae help.