In most of the times i try to to run this code, it returns the error:
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0 0xffffffff
#1 0xffffffff
#2 0xffffffff
#3 0xffffffff
#4 0xffffffff
#5 0xffffffff
#6 0xffffffff
#7 0xffffffff
#8 0xffffffff
#9 0xffffffff
#10 0xffffffff
#11 0xffffffff
#12 0xffffffff
#13 0xffffffff
#14 0xffffffff
I’m still strugglin to learn the language but it’s (or it must be) a code for partciles collision using monte carlo method. I think the problem is the subroutine that change the velocities but i really have no idea. I know it is’nt the most efficient whey of writting this type of code and some of the code it’s temporary but the idea is pretty much it
program collision
use rf_mod, only: rf
implicit none
integer :: n_particles = 10, i, j, k=0, t, ts=1000
real(8), parameter :: radius=0.06, dt=0.00008, c=1.0
integer, dimension(:,:) , allocatable :: id_pairs, id_pairs_collide
real(8), dimension(:) , allocatable :: ids, d_pairs
real(8), dimension(:,:) , allocatable :: r, v, v1, v2,v1new, v2new, r1, r2
real(8), dimension(:,:,:), allocatable :: rs, vs
!GENARATING THE PARTICLES
allocate(r(2,n_particles)) !position of the particles
allocate(v(2,n_particles)) !vellocities of the particles
allocate(ids(n_particles))
allocate(rs(ts, 2, n_particles))
allocate(vs(ts, 2, n_particles))
rs = 0.0
vs = 0.0
ids = [(i, i=1, n_particles)]
call random_number(r)
!do i=1, 2 !randomly positioning each particle
!do j=1, n_particles
!r(i,j) = rf()
!end do
!end do
open(40, file = 'pos.txt', status='unknown')
do j =1, n_particles
if(r(1,j) > 0.5) then
write(40,1)r(:,j)!, 'red'
1 format(2f20.6)
else
write(40,2)r(:,j)!, 'green'
2 format(2f20.6)
end if
end do
close(40)
!call system ('gnuplot -p grafico.plt')
!initial vellocities
!if the particle is in the left corner, v=100.
!if it's in the right, v=-100
v(2,:) = 0
where(r(1,:) < 0.5)
v(1,:) = 100.
elsewhere
v(1,:) = -100.
end where
!DISTANCE BETWEEN ALL THE PAIR OF PARTICLES
!all pair of particles
id_pairs = combination(ids)
d_pairs = get_deltad_pairs(r)
do t=1, ts
do i=1, size(d_pairs)
if((d_pairs(i) <2*radius) .eqv. .true.) then
k = k+1
end if
enddo
allocate(id_pairs_collide(2,k))
j=1
do i=1, size(d_pairs)
if(d_pairs(i) < 2*radius) then
id_pairs_collide(:,j) = id_pairs(:,i)
j=j+1
end if
enddo
!Now, we need to calculate the velocities of each particle, iterate, compare distances and change the vellocities
v1 = v(:, id_pairs_collide(1,:)) !particle i
v2 = v(:, id_pairs_collide(2,:)) !particle j
r1 = r(:, id_pairs_collide(1,:))
r2 = r(:, id_pairs_collide(2,:))
call get_new_v(v1, v2, v1new, v2new)
v(:, id_pairs_collide(1,:)) = v1new
v(:, id_pairs_collide(2,:)) = v2new
where(r(1,:) > 1.0)
v(1,:) = -abs(v(1,:))
end where
where(r(1,:) < 0.0)
v(1,:) = abs(v(1,:))
end where
where(r(2,:) > 1.0)
v(2,:) = -abs(v(2,:))
end where
where(r(2,:) < 0.0)
v(2,:) = abs(v(2,:))
end where
r = r + v*dt
rs(i,:,:) = r(:,:)
vs(i,:,:) = v(:,:)
deallocate(id_pairs_collide)
deallocate(v1)
deallocate(v2)
deallocate(r1)
deallocate(r2)
!print *, r(:,1)
enddo
open(40, file = 'all_pos.txt', status='unknown')
do i = 1, ts
write(40,*)rs(i,2,:)
end do
close(40)
PRINT*, v1new
contains
function get_deltad_pairs(r)
!get the distance between all the particles
!d_ij = sqrt(x_ij ^ 2 + y_ij ^ 2)
real(8), allocatable :: get_deltad_pairs(:), r(:,:)
get_deltad_pairs = sqrt(diff(combination(r(1,:)))**2 + diff(combination(r(2,:)))**2)
end function get_deltad_pairs
function combination(arg) result(comb)
implicit none
real(8), dimension(:,:), allocatable :: comb
real(8), intent(in), dimension(:) :: arg
integer :: M, N, counter, i, j
N=size(arg)
M=int(N*(N-1)/2)
allocate(comb(2,M))
counter=1
do i=1,N
do j=i+1,N
comb(1,counter)=arg(i)
comb(2,counter)=arg(j)
counter = counter + 1
enddo
enddo
end function combination
function diff(v)
implicit none
real(8), intent(in), dimension(:,:) :: v
real(8), allocatable:: diff(:), v1(:), v2(:)
integer :: n
n = size(v)/2
allocate(v1(n))
allocate(v2(n))
v1 = v(1,:)
v2 = v(2,:)
diff = abs(v1 - v2)
end function diff
subroutine get_new_v(v1, v2, v1n, v2n)
implicit none
real(8), dimension(:,:), intent(in) :: v1, v2
real(8), dimension(:,:), allocatable :: G, gr, grnew
real(8), dimension(:,:), allocatable, intent(out) :: v1n, v2n
real(8) :: cosx, sinx, e, RF1, RF2, PI=4.D0*DATAN(1.D0)
integer :: i, j
gr = v1 - v2 !relative velocities
G = (v1 + v2)/2.0 !relative velocities o the center of mass
do i=1, 2
do j=1, size(gr)/2
call random_number(RF2)
call random_number(RF1)
cosx = 2.0*RF1 - 1.0
sinx = sqrt(1 - cosx**2)
e = 2*RF2*pi
if(i==1) then
grnew(i,j) = sqrt(sum(gr(i,:)**2))*cosx
else
grnew(i,j) = sqrt(sum(gr(i,:)**2))*sinx*sin(e)
endif
enddo
enddo
v1n = G + 0.5*grnew
v2n = G - 0.5*grnew
return
end subroutine get_new_v
end program
module rf_mod
contains
DOUBLE PRECISION FUNCTION RF()
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
data IX1,IX2,IX3,IX4,IX5,IX6/1500419,1400159,1364,1528,1,3/
RR1=1.0/FLOAT(IX1)
RR2=1.0/FLOAT(IX2)
IX5=MOD(IX5*IX3,IX1)
IX6=MOD(IX6*IX4,IX2)
RF=RR1*IX5+RR2*IX6
IF(RF.GE.1.0)RF=RF-1.0
END FUNCTION
end module