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

NLS Equation using split step method in fortran

Afroja_A_
Beginner
‎05-15-2017 11:21 PM
363 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

Link Copied

0 Replies
Reply

Community support is provided during standard business hours (Monday to Friday 7AM - 5PM PST). Other contact methods are available here.

Intel does not verify all solutions, including but not limited to any file transfers that may appear in this community. Accordingly, Intel disclaims all express and implied warranties, including without limitation, the implied warranties of merchantability, fitness for a particular purpose, and non-infringement, as well as any warranty arising from course of performance, course of dealing, or usage in trade.

For more complete information about compiler optimizations, see our Optimization Notice.