DSMC on gas dynamics and looking for a bug

Hey friends, everyday i show up here again to search for help because i’m still beggining to learn the language :upside_down_face:

I don’t know if anyone of you is familiar to DSMC for gas dynamics but this is my code for now, but the thing is that it is stoping before the time it should, in the first ‘do’ and i don’t know why, if someone find what’s wrong or how to improve something please feel free to share. I think it is some parameter that is stopping the code because of some bug in the loop for changing the vellocities. N_mc is ny suspect but i’m not sure

program collisions

    implicit none

    integer               :: t, c, i, j, k, N_col=0
    integer, parameter    :: N=200, dim = 3                      !number of particles
    integer, parameter    :: Ncell=10         !number of cells
    integer               :: PPC(Ncell**3), cell_number, Nm !number of particles per cell
    integer, parameter    :: Nmft=20, Nt = 10000
    real(8), parameter    :: pi = 4.D0*DATAN(1.D0)
    real(8), parameter    :: radius = 0.06
    real(8), parameter    :: Tw = 1.0                     !temperature of the wall
    real(8), parameter    :: sigmat = pi*radius**2        !total cross section
    real(8), parameter    :: Lx = 1.0 , Ly = 1.0 , Lz = 1.0, VOLUME = Lx*Ly*Lz   !dimensions of the box
    real(8), parameter    :: n0 = N/VOLUME
    real(8), parameter    :: mfp = 1/(sqrt(2.0)*n0*sigmat**2) !mean free path
    real(8), parameter    :: Kn = mfp / Lx                    !Knudssen number
    real(8), parameter    :: v = 100.0, v_mean =  sqrt(8*Tw/pi),  v_rel_max = 200
	real(8), parameter    :: tau = mfp / v_mean, dt = 0.00001*tau         
	real(8), parameter    :: dx = Lx/Ncell, dy=Ly/Ncell, dz = Lz/Ncell                
	real(8), parameter    :: vol = dx*dy*dz, Nr = n0*VOLUME/N! !volume of each cell and number of real particles per cell
    real(8)               :: N_mc, Fn !Maximum number of model particle collisions and number of model particle per cell
    real(8), dimension(N) :: x, y, z
    real(8), dimension(N) :: vx, vy, vz
    real(8), dimension(:), allocatable :: x_c, y_c, z_c, vx_c, vy_c, vz_c
    integer, dimension(N) :: cell
    integer, allocatable, dimension(:) :: particles_in_cell
    real(8)               :: v_rel, Rf1, Rf2, Rf3, r1, r2, vx_cm, vy_cm, vz_cm, g1, g2, g3, e, cosx, sinx
    integer ::  i_prop, j_prop
    
    !setting initial positions
    call random_number(x)
    call random_number(y)
    call random_number(z)
    x = x * Lx
    y = y * Ly
    z = z * Lz

    
    !setting initial velocities
    where(x .le. Lx/2.0)
        vx = v
        else where
        vx = -v
    end where
    vy = 0
    vz = 0
    
    do t=1, Nt
        
        !drifting
        x = x + vx*dt
        y = y + vy*dt 
        z = z + vz*dt
        
        !interacting with bounderies
        
        where(x .le. 0)
            vx = abs(vx)
        endwhere
        where(x .ge. Lx)
            vx = -abs(vx)
        endwhere
        where(y .le. 0)
            vy = abs(vy)
        endwhere
        where(y .ge. Ly)
            vy = -abs(vy)
        endwhere
        where(z .le. 0)
            vz = abs(vz)
        endwhere
        where(y .ge. Lz)
            vz = -abs(vz)
        endwhere
        
        !indexing particles in cell 
        PPC = 0.0
        do c=1, N
            i = int(x(c) / dx) 
            j = int(y(c) / dy) 
            k = int(z(c) / dz) 
            cell_number = i + j * Ncell + k * Ncell**2
            cell(c) = cell_number
            PPC(cell_number) = PPC(cell_number) + 1
            !print*, c, x(c), cell_number
    
        enddo

        !passing through all cells
        do c=1, Ncell**3, 1
         
            x_c = pack(x, cell == c)
            y_c = pack(y, cell == c)
            z_c = pack(z, cell == c)
            !print*, c, PPC(c), x_c
            vx_c = pack(vx, cell == c)
            vy_c = pack(vy, cell == c)
            vz_c = pack(vz, cell == c)


            Nm = PPC(c) !number of number particles in the cell
            Fn = Nr / Nm
            N_mc = (Nm**2*Fn *sigmat * v_rel_max * dt ) / ( 2 * vol) !(Nm **2 *Fn * sigmat * v_rel_max * dt )
            !print*, N_mc
            
            do k=0, int(N_mc)
                call random_number(Rf1)
                
                !i and j prop must be integers in range of the number of particles in the cell Nm
                call random_number(r1)
                call random_number(r2)
                i_prop = 1+floor(Nm*r1) !random integer in the interval [1, Nm]
                j_prop = 1+floor(Nm*r2)

                v_rel = sqrt((vx_c(i_prop) - vx_c(j_prop))**2 + (vy_c(i_prop) - vy_c(j_prop))**2 + (vz_c(i_prop) - vz_c(j_prop))**2)
                !print*, v_rel
                !accept-reject
                if((v_rel/v_rel_max) > Rf1) then
                    
                    vx_cm = 0.5 * (vx_c(i_prop) + vy_c(j_prop))  
                    vy_cm = 0.5 * (vy_c(i_prop) + vy_c(j_prop))
                    vz_cm = 0.5 * (vz_c(i_prop) + vz_c(j_prop))

                    call random_number(Rf2)
                    call random_number(Rf3)

                    cosx = 2.0 * Rf2 - 1.0
                    sinx = sqrt(1 - cosx**2)
                    e = 2.0 * pi * Rf3

                    g1 = v_rel * cosx
                    g2 = v_rel * sinx * sin(e)
                    g3 = v_rel * sinx * cos(e)

                    vx_c(i_prop) = vx_cm + 0.5*g1
                    vy_c(i_prop) = vy_cm + 0.5*g2
                    vz_c(i_prop) = vz_cm + 0.5*g3

                    vx_c(j_prop) = vx_cm - 0.5*g1
                    vy_c(j_prop) = vy_cm - 0.5*g2
                    vz_c(j_prop) = vz_cm - 0.5*g3

                    N_col = N_col + 1
                    continue
                endif
                x = unpack(x_c, cell==c, x)
                y = unpack(y_c, cell==c, y)
                z = unpack(z_c, cell==c, z)
    
                vx = unpack(vx_c, cell==c, vx)
                vy = unpack(vy_c, cell==c, vy)
                vz = unpack(vz_c, cell==c, vz)

                deallocate(vx_c)
                deallocate(vy_c)
                deallocate(vz_c)
            enddo 
        enddo     
        write(*,*)'T=',t, 'Number of collisions = ', N_col
    enddo
end program collisions

The program is not standard Fortran. I suggest fixing this first.

Other than a few tab characters, what isn’t conforming about the code?

If you’re using gfortran compile your code with debugging options enabled. At minimum, use -g -fcheck=all -Wall -Werror -fmax-errors=1 and fix issues. You have an array indexing issue. I think you need something like

     if (cell_number < 1 .or. cell_number > ncell**3) cycle

around line 92.

2 Likes

It looks like you are deallocating some arrays, vx_c, vy_c, and vz_c, and then referencing them before they are assigned/allocated again. I’m not following all of the logic, but that is what it looks like. Also, there might be problems when the position is outside of the cell range. You allow the positions to be outside of the range 0 to Lx, etc., but you still use those values in the velocity expressions and in the cell number expression. That might be alright, but it might be causing indexing problems. You might need to clip the values at the boundaries before they are used.

What do you do with an empty cell ?
perhaps:

        do c=1, Ncell**3
            Nm = PPC(c) !number of number particles in the cell
            if ( Nm < 1 ) cycle

also
real(8), parameter :: radius = 0.06d0
this would improve accuracy as 0.06 /= 0.06d0
sqrt(2.0) /= sqrt(2.0d0)

" sinx = sqrt(1 - cosx**2)" looks to ignore the sign of sinx ?

Hi Steve

Here are the error messages from the Nag compiler.

Error: test1000.f90, line 10: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 11: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 12: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 13: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 14: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 15: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 16: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 17: KIND value (8) does not specify a valid representation method
Warning: test1000.f90, line 17: Unused PARAMETER KN
Error: test1000.f90, line 18: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 19: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 20: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 21: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 22: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 23: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 24: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 25: KIND value (8) does not specify a valid representation method
Error: test1000.f90, line 28: KIND value (8) does not specify a valid representation method

I understand that modern Fortran provides pack and unpack, but is this enabling the best coding approach ? ( especially where all unpack operations involve the same “cell==c” )

                x = unpack(x_c, cell==c, x)
                y = unpack(y_c, cell==c, y)
                z = unpack(z_c, cell==c, z)
    
                vx = unpack(vx_c, cell==c, vx)
                vy = unpack(vy_c, cell==c, vy)
                vz = unpack(vz_c, cell==c, vz)

                deallocate(vx_c)
                deallocate(vy_c)
                deallocate(vz_c)

As for pack, there is no report of the size of the resulting vector.
For me, this is a messy functionality.

What is bad about a DO loop in modern Fortran ?
Is DO the next GOTO !!

1 Like

This does not make the code non-conforming. It makes it nonportable. kind=8 works with gfortran and Intel, and IIRC, it will work with NAG with a compiler option. The only non-conforming things found with gfortran’s -std=f2018 option were 3 tab characters used where a member of the Fortran character set was expected.

1 Like

i’m more familiar with python but somethimes i forget that e clear code isn’t the same as an efficient code, coding for science is different than coding for fun and i need to remem,ber myself of that everytime

The code has numerous errors, most of which may be classified under “optimistic programming”. For example, on line 103, with c = 1, you have Nm = PPC(c) = 0, yet on this line you divide by Nm. Similarly, the code may be using zero as a subscript in some places, and may be using undefined variables when evaluating expressions.

You must fix these errors before we can discuss why the code is not yielding the kinds of results that you expect.