Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
Announcements
Welcome to the Intel Community. If you get an answer you like, please mark it as an Accepted Solution to help others. Thank you!
26755 Discussions

NLS Equation using split step method in fortran

Afroja_A_
Beginner
111 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 ='
!read (*,*) length
length = 100.0d0
!print*, 'dispersion 1 for normal, -1 for anomalous'
!read (*,*) dispersion
dispersion = 1.0d0
!print*, 'Nonlinear parameter N ='
!read (*,*) N
N= 1.0d0
!print*, 'm = 0 for sech , m>0 for super-gaussian='
!read (*,*) m
m =0.0d0
 chirp =0.0d0
 
! print*, 'Enter Sampling point nt='
! read(*,*) nt
 nt =1024.0d0
! print*, 'Enter window size Tmax='
! read(*,*) 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 ='
!read (*,*) length
length = 100.0d0
!print*, 'dispersion 1 for normal, -1 for anomalous'
!read (*,*) dispersion
dispersion = 1.0d0
!print*, 'Nonlinear parameter N ='
!read (*,*) N
N= 1.0d0
!print*, 'm = 0 for sech , m>0 for super-gaussian=' 
!read (*,*) m
m =0.0d0
 chirp =0.0d0
 
! print*, 'Enter Sampling point nt='
! read(*,*) nt
 nt =1024.0d0
! print*, 'Enter window size Tmax='
! read(*,*) 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.

Thanks in advance

0 Kudos
0 Replies
Reply