程序收到信号 SIGSEGV：分段错误 - 大尺寸数组中的内存引用无效

我在运行此代码时遇到错误。当我使用 L=16 或 L=32 等小 L 运行代码时，我没有得到任何错误，但在 7000-8000 步后，在 L = 128 或 L=96 中，我得到以下错误：

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)

这是完整的代码：

    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)：我使用 omp lib 并行运行我的程序

P.S.(2)：我使用gfortran来编译代码

P.S.(3)：编译的代码-g -fcheck=all并给我这个错误：

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

谢谢大家

Your particle_ad数组只有空间10*L粒子，但你似乎试图存储最多N=2*L**2其中的粒子（取决于随机数如何下降）。平均而言，每个粒子都有足够的空间，但是如果太多粒子落入单个块中，您的代码将会失败，并且（如果我没记错概率）发生这种情况的可能性会增加L增加。

您可以通过替换来解决这个问题Integer :: particle_ad(10*L) with Integer :: particle_ad(N)，但这会浪费大量空间。

更好的解决方案是重新调整大小particle_ad动态数组，每次填满时都会使它们变大。为了方便起见，您可以将此行为包装在block_p class.

例如，

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

然后你可以更换线路

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)

请注意，由于每个particle_ad数组现在是allocatable，您需要先初始化每个数组，然后才能调用add_particle.