I have this code and It does heavy calculations and I need to use open-mp library to get it done fast . This is my code :
SUBROUTINE init_random_seed()
implicit none
INTEGER :: i, n, clock
INTEGER, DIMENSION(:), ALLOCATABLE :: seed
CALL RANDOM_SEED(size = n)
ALLOCATE(seed(n))
CALL SYSTEM_CLOCK(COUNT=clock)
seed = clock + 37 * (/ (i - 1, i = 1, n) /)
CALL RANDOM_SEED(PUT = seed)
DEALLOCATE(seed)
END SUBROUTINE
!end module
Program Activ_mater
USE OMP_LIB
Implicit none
Integer,parameter :: time=1000000000, L=128, N = L**2*2
Integer,parameter:: n_thread = 8
Real(8),parameter :: pi = 3.14159265359
Real(8),parameter :: v0 = 0.50, alpha = 1.0/36.0
real(8)START,END_1 ,eta
type block_p
Integer :: partical_N
Integer :: particle_ad(10*L)
end type
Type(block_p) ,pointer,dimension(:,:) :: C
Real(8),allocatable :: x(:), y(:) , phi,angle_new(:),angle(:)
Real(8) :: sum_a, sum_b,x_in, y_in, x_out, y_out, avrage_t, r,ra,eta1(5)
Integer :: i,j,t,n_p,I_b,J_b,b_l,neighbor_i(9),neighbor_j(9),A,n_p_b,ne = 1,stateta=0,ot=0,op=175
character(len=10)::name1
call omp_set_num_threads(n_thread)
call init_random_seed()
eta1=(/2.100,2.116,2.120,2.126,2.128/) ! The value of ETA
allocate(x(2*n), y(2*n) , phi,angle_new(2*N),angle(2*N))
allocate(C(2*L,2*L))
C(:,:)%partical_N=0
do i =1,N
call random_number(ra)
x(i)=ra*L
I_b = int(x(i))+1
call random_number(ra)
y(i)=ra*L
J_b = int(y(i))+1
call random_number(ra)
angle(i)=ra*2.0*pi
C(I_b,J_b)%partical_N = C(I_b,J_b)%partical_N + 1 !Number of particle in block C(I_b,J_b)
C(I_b,J_b)%particle_ad( C(I_b,J_b)%partical_N ) = i ! The particle number in block C(I_b,J_b)
end do
! loop for eta
eta= 0.0
write(name1,'(f5.3)')eta
open(unit=10, file='Hysteresis,'//trim(name1)//'.dat')
!=====================explanation of system====================================
print*,'==========================================================================='
print*, 'eta = ', etA ,' ',' alpha = ',alpha
print*,'L=',L ,' ', 'Particle Number=', N,' ','Density=', N/L**2
print*,'==========================================================================='
!==============================================================================
START = omp_get_wtime()
do t =1,time
if (ot == 300000 )then
stateta = 0
ot = 0
op = op + 1
end if
if (stateta == 0 )then
eta = eta + ((1.0/3.0) * (10E-6))
endif
if (int(eta * 100) == op) then
stateta = 1
end if
angle_new(:)=0
!$OMP PARALLEL DEFAULT(PRIVATE) SHARED(x,y,angle,angle_new,c) firstprivate(eta)
!$OMP DO schedule(runtime)
do i =1, N
sum_a=0; sum_b=0;n_p=0
I_b = int(x(i))+1; J_b = int(y(i))+1 ! The block of particle i
! Now I should find nine neighbor of particle i-----------------------------------------------
neighbor_i=(/I_b+1, I_b, I_b-1, I_b, I_b+1, I_b-1, I_b-1, I_b+1, I_b/)
neighbor_j=(/J_b, J_b+1, J_b, J_b-1, J_b+1, J_b+1, J_b-1, J_b-1, J_b/)
do b_l = 1, 9
I_b = neighbor_i(b_l) ; J_b=neighbor_j(b_l)
if (I_b >L )I_b=1
if (I_b <1 )I_b=L
if (J_b >L )J_b=1
if (J_b <1 )J_b=L
!neighbor_i(b_l)=I_b; neighbor_j(b_l)=J_b
A = C( I_b, J_b )%partical_N ! number of particle in block C( neighbor_i(b_l), neighbor_j(b_l) )
!=============================================================================================
do n_p_b =1, A
j = C( I_b, J_b)%particle_ad(n_p_b) !particle j in the block C
if (i /= j )then
X_in = abs(max(x(i),x(j)) - min(x(i),x(j)));
Y_in = abs(max(y(i),y(j)) - min(y(i),y(j)));
X_out =L-X_in
Y_out =L-Y_in
r = sqrt(min(X_in,X_out)**2 + min(Y_in,Y_out)**2)
else
r=0.0
end if
if ( r < 1 )then
if ( j <= i )then
sum_A = sum_A + sin(angle(j));
sum_B = sum_B + cos(angle(j));
else
sum_A = sum_A + alpha*sin(angle(j));
sum_B = sum_B + alpha*cos(angle(j));
endif
n_p = n_p + 1;
endif
enddo
enddo
sum_A = sum_A/n_p; sum_B = sum_B/n_p
!if (int(sum_A*1e10) ==0 .and. int(sum_B*1e10) ==0 )print*,'zerrooo'
avrage_t=atan2(sum_A,sum_B);
if (avrage_t<0.0) then
avrage_t=avrage_t+2.0*pi;
endif
call random_number(ra)
angle_new(i)=avrage_t+eta*(ra-0.50)
if( angle_new(i)>=2*pi) angle_new(i)= angle_new(i)-2*pi
if( angle_new(i)<0) angle_new(i)= 2*pi+angle_new(i)
end do
!$OMP END DO
!$OMP END PARALLEL
angle = angle_new
C(:,:)%partical_N=0
! phi=0.0
do i=1, N
x(i) = x(i) + v0*sin(angle(i));
if (x(i)<1) x(i)=L+x(i)
if (x(i)>L) x(i)=x(i)-L
I_b = int(x(i))+1
y(i) = y(i)+ v0*cos(angle(i));
if (y(i)<1) y(i)=L+y(i)
if (y(i)>L) y(i)=y(i)-L
J_b = int(y(i))+1
C(I_b,J_b)%partical_N = C(I_b,J_b)%partical_N + 1
C(I_b,J_b)%particle_ad( C(I_b,J_b)%partical_N ) = i
end do
if (stateta == 1 )then
phi= sqrt((sum(sin(angle))**2+sum(cos(angle))**2))/N;
ot = ot + 1
end if
write(10,*)phi
if (mod(t,10)==0)then
! ave4_phi=sum(phi**4)/t;
! ave2_phi =sum(phi**2)/t;
!print* ,ave4_phi,ave2_phi
print*,'Time=',t,' ==== Eta : ',eta,"Ot : " , ot
end if
end do
END_1 = omp_get_wtime()
print*,'Run Time = ',end_1-start
End Program
As you can see , there is a parameter "L" that defines number of particles . When "L" is lower , like 64 or 32 it works fine but when I increase it , I get the error .
Also without omplib it looks like it works fine in any L size but when I use it and compile it with -g -fcheck=all , I get this error :
At line 155 of file hyster.f90 Fortran runtime error: Index '1281' of dimension 1 of array 'c%particle_ad' above upper bound of 1280
I would be grateful if you help me . Thank you all .
