3

I am getting an error while I run this code. When I run the code with small L's like L=16 or L=32 I get no error but in L = 128 or L=96 after 7000-8000 steps I get following error :

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0x7FBA5CAC3E08
#1  0x7FBA5CAC2F90
#2  0x7FBA5C1E84AF
#3  0x402769 in MAIN__ at newhys.f90:?
Segmentation fault (core dumped)

This is the full 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

P.S.(1) : I use omp lib to run my program parallel

P.S.(2) : I use gfortran to compile the code

P.S.(3) : Code compiled with -g -fcheck=all and gives me this error :

At line 155 of file z.f90 Fortran runtime error: Index '1281' of dimension 1 of array 'c%particle_ad' above upper bound of 1280

Thanks to you all

Mohammad Jafari
  • 37
  • 1
  • 5
  • 1
    How are you compiling your code? If you use the [`-g` flag](https://gcc.gnu.org/onlinedocs/gcc/Debugging-Options.html) then the backtrace will contain useful information rather than just pointers. – veryreverie Feb 06 '22 at 15:44
  • 2
    When debugging with gfortran, always compile your code with error checking. Try `gfortran -g -Wall -fcheck=all`. An alternative option can be the address sanitizations `gfortran -g -Wall -fsanitize=address,undefined`. Please try it and tell us what happened. – Vladimir F Героям слава Feb 06 '22 at 16:39
  • 3
    If you compile with `-g -fcheck=all` and let it run you eventually get `At line 155 of file z.f90 Fortran runtime error: Index '1281' of dimension 1 of array 'c%particle_ad' above upper bound of 1280` . Note `Real( 8 )` is not portable, not guaranteed to compile, and may not do what you expect it to do - https://stackoverflow.com/questions/838310/fortran-90-kind-parameter – Ian Bush Feb 06 '22 at 16:43
  • 1
    Also, always first try without OpenMP and then with OpenMP. – Vladimir F Героям слава Feb 06 '22 at 17:09
  • @VladimirF I did and I got Index above upper bound error . Also it looks like it works fine without OpenMP. – Mohammad Jafari Feb 19 '22 at 07:13
  • @MohammadJafari Use [edit] to shoe us the exact error message. Even if it is the same as Ian Bush shows in his comment, put it in your question. Also, think why it gets the value of 1280 and whether you should increase the array size. Perhaps 1280 is too small. Why do you use 1280? – Vladimir F Героям слава Feb 19 '22 at 07:34

1 Answers1

0

Your particle_ad arrays only have space for 10*L particles, but you appear to be trying to store up to N=2*L**2 particles in them (depending on how the random numbers fall). On average, each will have enough space, but your code will fail if too many particles fall into a single block, and (if I'm remembering the probabilities right) the chances of this happening increase as L increases.

You could solve this problem by replacing Integer :: particle_ad(10*L) with Integer :: particle_ad(N), but this will waste an awful lot of space.

A better solution would be to re-size the particle_ad arrays on the fly, making them bigger every time they get full. For convenience, you could wrap this behaviour in a method of the block_p class.

For example,

type block_p
  Integer :: partical_N
  Integer, allocatable :: particle_ad(:)
contains
  subroutine add_particle(this, index)
    class(block_p), intent(inout) :: this
    integer, intent(in) :: index
    
    integer, allocatable :: temp
    
    ! Update partical_N.
    this%partical_N = this%partical_N + 1
    
    ! Resize particle_ad if it is full.
    if (size(this%particle_ad)<this%partical_N) then
      temp = this%particle_ad
      deallocate(this%particle_ad)
      allocate(this%particle_ad(2*this%partical_N))
      this%particle_ad(:size(temp)) = temp
    endif
    
    ! Add the new index to particle_ad.
    this%particle_ad(this%partical_N) = index
  end subroutine
end type

You could then replace the lines

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

with

call C(I_b,J_b)%add_particle(i)

Note that as each particle_ad array is now allocatable, you will need to initialise each array before you can call add_particle.

veryreverie
  • 2,871
  • 2
  • 13
  • 26