particle_ode.f90

visual fortran 源码

program odemain
real :: h,k1(6),k2(6),k3(6),k4(6),k5(6),k6(6),ubar
real :: xp(23),yp(23),zp(23),vpx(23),vpy(23),vpz(23),u(23,23),v(23,23)
real :: x(ni),y(23),xdif(23),ydif(23)
real :: puin(ni,nj,k),pvin(ni,nj,k),p3in(ni,nj,k),p4in(ni,nj,k)
real :: puout(ni,nj,k),pvout(ni,nj,k),p3out(ni,nj,k),p4out(ni,nj,k)
real :: deltx,delty
real,external::func
integer i,j,ii,jj,iii,jjj,k,num
open(unit=101,file='results.txt')
!=============================================
do i=1,23
DO J=1,23
U(I,J)=0.5
V(I,J)=0.2
XDIF(I)=0.2
yDIF(j)=0.2
END DO 
END DO 
XDIF(2)=0.1
XDIF(23)=0.1
yDIF(2)=0.1
yDIF(23)=0.1
X(1)=0
y(1)=0
DO I=2,23
X(I)=x(i-1)+xdif(i)
y(i)=y(i-1)+ydif(i)
end do
ubar=0.3
!===============================================
do k=1,num


!===================ode求解第K个颗粒：龙格-库塔法==========================
   do while(ti<tlimit)
      do j=1,6
      k1(j)=dfunc(ti,vpx(k),vpy(k),vpz(k),xp(k),yp(k),zp(k),j)
	  end do
      do j=1,6
      k2(j)=dfunc(ti+h/2,vpx(k)+k1(1)*h/2,vpy(k)+k1(2)*h/2,vpz(k)+k1(3)*h/2,xp(k)+k1(4)*h/2,yp(k)+k1(5)*h/2,zp(k)+k1(6)*h/2,j)
      end do
      do j=1,6
      k3(j)=dfunc(ti+h/2,vpx(k)+k2(1)*h/2,vpy(k)+k2(2)*h/2,vpz(k)+k2(3)*h/2,xp(k)+k2(4)*h/2,yp(k)+k2(5)*h/2,zp(k)+k2(6)*h/2,j)
      end do
      do j=1,6
      k4(j)=dfunc(ti+h/2,vpx(k)+k3(1)*h/2,vpy(k)+k3(2)*h/2,vpz(k)+k3(3)*h/2,xp(k)+k3(4)*h/2,yp(k)+k3(5)*h/2,zp(k)+k3(6)*h/2,j)
      end do
	  vpx0=vpx(k)
      vpy0=vpy(k)
      vpz0=vpz(k)
	  xp0=xp(k)
      yp0=yp(k)
      zp0=zp(k)      
	  vpx(k)=vpx(k)+h*(k1(1)+2*k2(1)+2*k3(1)+k4(1))/6
      vpy(k)=vpy(k)+h*(k1(2)+2*k2(2)+2*k3(2)+k4(2))/6
      vpz(k)=vpz(k)+h*(k1(3)+2*k2(3)+2*k3(3)+k4(3))/6
      xp(k)=xp(k)+h*(k1(4)+2*k2(4)+2*k3(4)+k4(4))/6
      yp(k)=yp(k)+h*(k1(5)+2*k2(5)+2*k3(5)+k4(5))/6
	  zp(k)=zp(k)+h*(k1(6)+2*k2(6)+2*k3(6)+k4(6))/6
	  ti=ti+h
!========================颗粒源项设置==>U=====================================
	  if(x(i-1)<=xp(k).and.xp(k)<=x(i))then
	     if(yv(j)<=yp(k).and.yp(k)<=yv(j+1))then
		 continue
		 elseif(yp(k)>yv(j+1))then
         delty=yp(k)-yp0+1e-30
		 fyp=(yv(j+1)-yp0)/delty
		 puin(i,j+1,k)=vpx0+fyp*(vpx(k)-vpx0)
		 puout(i,j,k)=puin(i,j+1,k)
		 j=j+1
		 else
		 delty=yp(k)-yp0-1e-30
		 fyp=(yv(j)-yp0)/delty
		 puin(i,j-1,k)=vpx0+fyp*(vpx(k)-vpx0)
		 puout(i,j,k)=puin(i,j-1,k)
		 j=j-1
		 end if
      elseif(xp(k)>x(i))then
	     deltx=xp(k)-xp0+1e-30
		 fxp=(x(i)-xp0)/deltx
		 pufxp=vpx0+fxp*(vpx(k)-vpx0)
		 puout(i,j,k)=pufxp
	     if(yv(j)<=yp(k).and.yp(k)<=yv(j+1))then
		 puin(i+1,j,k)=pufxp
		 i=i+1
		 elseif(yp(k)>yv(j+1))then
		 puin(i+1,j+1,k)=pufxp
		 i=i+1
		 j=j+1
		 else
		 puin(i+1,j-1,k)=pufxp
		 i=i+1
		 j=j-1
		 end if
	  else
	     deltx=xp(k)-xp0-1e-30
		 fxp=(x(i-1)-xp0)/deltx
		 pufxp=vpx0+fxp*(vpx(k)-vpx0)
		 puout(i,j,k)=pufxp
		 if(yv(j)<=yp(k).and.yp(k)<=yv(j+1))then
		 puin(i-1,j,k)=pufxp
		 i=i-1
		 elseif(yp(k)>yv(j+1))then
		 puin(i-1,j+1,k)=pufxp
		 i=i-1
		 j=j+1
		 else
		 puin(i-1,j-1,k)=pufxp
		 i=i-1
		 j=j-1
		 end if
	  end if
!========================颗粒源项设置==>V==========================		 
      if(y(jj-1)<=yp(k).and.yp(k)<=y(jj))then
	     if(xu(ii)<=xp(k).and.xp(k)<=xu(ii+1))then
		 continue
		 elseif(xp(k)>xu(ii+1))then
		 deltx=xp(k)-xp0+1e-30
		 fxp=(xu(ii+1)-xp0)/deltx
		 pvin(ii+1,jj,k)=vpy0+fxp*(vpy(k)-vpy0)
		 pvout(ii,jj,k)=pvin(ii+1,jj,k)
		 ii=ii+1
		 else
		 deltx=xp(k)-xp0-1e-30
		 fxp=(xu(ii)-xp0)/deltx
		 pvin(ii-1,jj,k)=vpy0+fxp*(vpy(k)-vpy0)
		 pvout(ii,jj,k)=pvin(ii-1,jj,k)
		 ii=ii-1
		 end if
	  elseif(y(jj)<yp(k))then
	     delty=yp(k)-yp0+1e-30
		 fyp=(y(jj)-yp0)/delty
		 pvfyp=vpy0+fyp*(vpy(k)-vpy0)
		 pvout(ii,jj,k)=pvfyp
		 if(xu(ii)<=xp(k).and.xp(k)<=xu(ii+1))then
		 pvin(ii,jj+1,k)=pvout(ii,jj,k)
		 jj=jj+1
		 elseif(xp(k)>xu(ii+1))then
		 pvin(ii+1,jj+1,k)=pvfyp
		 ii=ii+1
		 jj=jj+1
		 else
		 pvin(ii-1,jj+1,k)=pvfyp
		 ii=ii-1
		 jj=jj+1
		 end if
	  else
	     delty=yp(k)-yp0-1e-30
		 fyp=(y(jj-1)-yp0)/delty
		 pvfyp=vpy0+fyp*(vpy(k)-vpy0)
		 pvout(ii,jj,k)=pvfyp
		 if(xu(ii)<=xp(k).and.xp(k)<=xu(ii+1))then
		 pvin(ii,jj-1,k)=pvfyp
		 jj=jj-1
		 elseif(xp(k)>xu(ii+1))then
		 pvin(ii+1,jj-1,k)=pvfyp
		 ii=ii+1
		 jj=jj-1
		 else
		 pvin(ii-1,jj-1,k)=pvfyp
		 ii=ii-1
		 jj=jj-1
		 end if
      end if
!==================其他颗粒源项设置============================
      if(xu(iii)<=xp(k).and.xp(k)<=xu(iii+1)then
	     if(yv(jjj)<=yp(k).and.yp(k)<=yv(jjj+1))then
		 continue
		 elseif(yp(k)>yv(jjj+1))then
		 delty=yp(k)-yp0+1e-30
		 fyp=(yv(jjj+1)-yp0)/delty
		 !==========待选择变量================
		 p3in(iii,jjj+1,k)=vpz0+fyp*(vpz(k)-vpz0)
		 p4in(iii,jjj+1,k)=zp0+fyp*(zp(k)-zp0)
		 p3out(iii,jjj,k)=p3in(iii,jjj+1,k)
		 p4out(iii,jjj,k)=p4in(iii,jjj+1,k)
		 !====================================
		 jjj=jjj+1
		 else
		 delty=yp(k)-yp0-1e-30
		 fyp=(yv(jjj)-yp0)/delty
		 !==========待选择变量================
		 p3in(iii,jjj-1,k)=vpz0+fyp*(vpz(k)-vpz0)
		 p4in(iii,jjj-1,k)=zp0+fyp*(zp(k)-zp0)
		 p3out(iii,jjj,k)=p3in(iii,jjj-1,k)
		 p4out(iii,jjj,k)=p4in(iii,jjj-1,k)
		 !====================================
		 jjj=jjj-1
		 end if
	  elseif(xp(k)>xu(iii+1))then
		 if(yv(jjj)<=yp(k).and.yp(k)<=yv(jjj+1))then
		 deltx=xp(k)-xp0+1e-30
		 fxp=(xu(iii+1)-xp0)/deltx
		 p3in(iii+1,jjj,k)=vpz0+fxp*(vpz(k)-vpz0)
         p3out(iii,jjj,k)=p3in(iii+1,jjj,k)
         p4in(iii+1,jjj,k)=zp0+fxp*(zp(k)-zp0)
         p4out(iii,jjj,k)=p4in(iii+1,jjj,k)
		 iii=iii+1
		 elseif(yp(k)>yv(j+1))then
		 p3in(iii+1,jjj+1,k)=(vpz0+vpz(k))/2
         p3out(iii,jjj,k)=p3in(iii+1,jjj+1,k)
		 p4in(iii+1,jjj+1,k)=((zp0+zp(k))/2
         p4out(iii,jjj,k)=p4in(iii+1,jjj+1,k)
		 iii=iii+1
		 jjj=jjj+1
		 else
		 p3in(iii+1,jjj-1,k)=(vpz0+vpz(k))/2
         p3out(iii,jjj,k)=p3in(iii+1,jjj-1,k)
         p4in(iii+1,jjj-1,k)=((zp0+zp(k))/2
         p4out(iii,jjj,k)=p4in(iii+1,jjj-1,k)
		 iii=iii+1
		 jjj=jjj-1
		 end if 
	  else
		 if(yv(jjj)<=yp(k).and.yp(k)<=yv(jjj+1))then
	     deltx=xp(k)-xp0-1e-30
		 fxp=(xu(iii)-xp0)/deltx
		 p3in(iii-1,jjj,k)=vpz0+fxp*(vpz(k)-vpz0)
		 p3out(iii,jjj,k)=p3in(iii-1,jjj,k)
         p4in(iii-1,jjj,k)=zp0+fxp*(zp(k)-zp0)
		 p4out(iii,jjj,k)=p4in(iii-1,jjj,k)
		 iii=iii-1
		 elseif(yp(k)>yv(j+1))then
		 p3in(iii-1,jjj+1,k)=(vpz0+vpz(k))/2
         p3out(iii,jjj,k)=p3in(iii-1,jjj+1,k)
		 p4in(iii-1,jjj+1,k)=((zp0+zp(k))/2
         p4out(iii,jjj,k)=p4in(iii-1,jjj+1,k)
		 iii=iii-1
		 jjj=jjj+1
		 else
		 p3in(iii-1,jjj-1,k)=(vpz0+vpz(k))/2
         p3out(iii,jjj,k)=p3in(iii-1,jjj-1,k)
         p4in(iii-1,jjj-1,k)=((zp0+zp(k))/2
         p4out(iii,jjj,k)=p4in(iii-1,jjj-1,k)
		 iii=iii-1
		 jjj=jjj-1
		 end if
	  end if
   end do
!============================================================================================          
   

end program

real function dfunc(x1,x2,x3,x4,x5,x6,x7,k)
real x1,x2,x3,x4,x5,x6,x7
integer k
select case(k)
case(1) 
dfunc=1
case(2)
dfunc=1
case(3)
dfunc=0
case(4)
dfunc=1 
case(5)
dfunc=1
case(6)
dfunc=0
case default
write(*,*) "data error!"
end select
return
end function dfunc