程序目的：一维波动过程仿真

程序语言：Fortran90程序

编译方式：ifort -qopenmp XX.f90 -o xx 编译，（Gcc icc）

运行方式：openmp程序，网格多时计算时间较长

program test

include omp_lib.h

!参数定义

!!定义浮点变量，包括第一π值、速度区间、

real(8) PI,v_max

!!定义格点特征变量,动能E_f,动量P_f，速度分布W_f,振子速度u_t

real(8) ,allocatable ::E_f(:,:),P_f(:,:),W_f(:,:),u_t(:)

!!定义速度分布，v_f(xyz,N_total,n_g),n为第n个格点，x,y,x为格点粒子速度分量,m为第每个粒子

real(8) , allocatable ::v_f(:,:,:);

!!

!!定义整数变量，包括碰撞循环次数，每次格点间碰撞次数（定为粒子数），网格数，统计网格数

integer(4) n_r,N_total,n_g,n_t

!!定义临时变量

integer(4) l_l,m_l,n_l,i,j,k,l,m,n,z_m

real(8) vx1,vx2,vy1,vy2,v_rela;

real(8) v_para1,v_para2,v_vert1,v_vert2,angle_collide;

!!定义随机数

real(8) , allocatable ::rand1(:,:),rand2(:,:),rand3(:,:),rand4(:,:);

!参数初始化

!!初始化π值

PI=2.0d0*dasin(1.0d0)

!!初始化碰撞循环次数

n_r=512

!!初始化统计格点数

n_t=101

!!初始化速度区间，速度区间为0-v_max

v_max=20000.0d0

!!初始化格点数目

n_g=100

!!初始化粒子速度

!!! open file 初始化粒子m速度

open(11,file=/data/lizi/pingheng/init.txt)

read(11,*)N_total

print*,N_total

!!!初始化矩阵

allocate (v_f(3,N_total,n_g))

allocate (E_f(3,n_g))

allocate (P_f(3,n_g))

allocate (W_f(101,n_g))

v_f=0.0d0

do i=1,N_total

read(11,*)v_f(1,i,1),v_f(2,i,1),v_f(3,i,1)

end do

close(11)

!!!初始化其他格点粒子速度

do i=2,n_g

v_f(:,:,i)=v_f(:,:,1)

end do

print*,end load data

!!数据统计

E_f=0.0d0

P_f=0.0d0

W_f=0.0d0

do z_n=1,n_g

do i=1,N_total

E_f(1,z_n)=E_f(1,z_n)+v_f(1,i,z_n)**2

E_f(2,z_n)=E_f(2,z_n)+v_f(2,i,z_n)**2

P_f(1,z_n)=P_f(1,z_n)+v_f(1,i,z_n)

P_f(2,z_n)=P_f(2,z_n)+v_f(2,i,z_n)

l_l=nint(v_f(1,i,z_n)/1000.0d0)+51

W_f(l_l,z_n)=W_f(l_l,z_n)+1

end do

end do

print*,数据统计结束

!数据输出

!输出粒子分布

open(12,file=out.txt)

write(12,*)初始数据输出

do i=1,n_g

write(12,*)网格,i,Ex、Ey=,E_f(1,i),E_f(2,i)

end do

do i=1,n_g

write(12,*)网格,i,Px、Py=,P_f(1,i),P_f(2,i)

end do

! do i=1,n_g

! write(12,*)网格,i,速度分布为

! do j=1,n_t

! write(12,*)W_f(j,i)

! end do

! end do

write(12,*)-------------------end output------------------------

print*,数据输出结束

!定义振子速度

allocate (u_t(n_r))

do i=1,n_r

u_t(i)=1000.0d0*dsin(2*PI*i/2/128)

write(12,*)振子速度为,u_t(i)

end do

!初始化随机数

allocate (rand1(n_g+1,N_total))

allocate (rand2(n_g+1,N_total))

allocate (rand3(n_g+1,N_total))

allocate (rand4(n_g+1,N_total))

!粒子碰撞,循环次数为z_l

do z_l=1,n_r

!!与振子碰撞

!!!产生随机数

call random_seed()

call random_number(rand1(:,:))

call random_number(rand2(:,:))

call random_number(rand3(:,:))

call random_number(rand4(:,:))

!omp

!$omp parallel

!private(z_m,n_m1,n_m2,vx1,vy1,vx2,vy2,v_para1,v_vert1,v_para2,v_vert2,v_rela,velocity_mc,angle_collide)

!$omp do

do i=1,N_total

z_m=0

n_m1=INT(N_total*rand1(z_m+1,i))+1

vx1=v_f(1,n_m1,1)

v_rela=u_t(z_l)-vx1;

velocity_mc=1.65d0*v_max*rand2(z_m+1,i);

if(velocity_mc

v_f(1,n_m1,1)=2.0d0*u_t(z_l)-vx1

end if

!!1--(n_g-1)间波动传递

do z_m=1,n_g-1

!!!产生随机数

!!!产生随机两个格点粒子

n_m1=INT(N_total*rand1(z_m+1,i))+1

n_m2=INT(N_total*rand2(z_m+1,i))+1

angle_collide=PI*(rand3(z_m+1,i)-0.5d0)

vx1=v_f(1,n_m1,z_m)

vy1=v_f(2,n_m1,z_m)

vx2=v_f(1,n_m2,z_m+1)

vy2=v_f(2,n_m2,z_m+1)

v_para1=vx1*dcos(angle_collide)+vy1*dsin(angle_collide);

v_vert1=-vx1*dsin(angle_collide)+vy1*dcos(angle_collide)

v_para2=vx2*dcos(angle_collide)+vy2*dsin(angle_collide);

v_vert2=-vx2*dsin(angle_collide)+vy2*dcos(angle_collide)

v_rela=v_para1-v_para2;

velocity_mc=v_max*rand4(z_m+1,i);

if(velocity_mc

vx1=v_para2*dcos(angle_collide)-v_vert1*dsin(angle_collide);

vy1=v_para2*dsin(angle_collide)+v_vert1*dcos(angle_collide);

vx2=v_para1*dcos(angle_collide)-v_vert2*dsin(angle_collide);

vy2=v_para1*dsin(angle_collide)+v_vert2*dcos(angle_collide);

v_f(1,n_m1,z_m)=vx1

v_f(2,n_m1,z_m)=vy1

v_f(1,n_m2,z_m+1)=vx2

v_f(2,n_m2,z_m+1)=vy2

end if

end do

!!边界n_g格点速度分布

!!!产生随机数

!!!产生随机两个格点粒子

z_m=n_g

n_m1=INT(N_total*rand1(z_m+1,i))+1

vx1=v_f(1,n_m1,z_m)

v_rela=vx1;

velocity_mc=1.65d0*v_max*rand2(z_m+1,i);

if(velocity_mc

v_f(1,n_m1,z_m)=-vx1

end if

end do

!end omp

!$omp end do

!$OMP END PARALLEL

print*,第,z_l,次循环结束

!!数据统计

E_f=0.0d0

P_f=0.0d0

W_f=0.0d0

do z_n=1,n_g

do i=1,N_total

E_f(1,z_n)=E_f(1,z_n)+v_f(1,i,z_n)**2

E_f(2,z_n)=E_f(2,z_n)+v_f(2,i,z_n)**2

P_f(1,z_n)=P_f(1,z_n)+v_f(1,i,z_n)

P_f(2,z_n)=P_f(2,z_n)+v_f(2,i,z_n)

l_l=nint(v_f(1,i,z_n)/1000.0d0)+51

W_f(l_l,z_n)=W_f(l_l,z_n)+1

end do

end do

!!数据输出

!输出粒子分布

write(12,*)第,z_l,次数据输出

do i=1,n_g

write(12,*)网格,i,Ex、Ey=,E_f(1,i),E_f(2,i)

end do

do i=1,n_g

write(12,*)网格,i,Px、Py=,P_f(1,i),P_f(2,i)

end do

! do i=1,n_g

! write(12,*)网格,i,速度分布为

! do j=1,n_t

! write(12,*)W_f(j,i)

! end do

! end do

write(12,*)-------------------end output------------------------

!结束循环

end do

close(12);

end program