求大神解答，什么问题

wx_yxZLpy1y · 发表于 2020-10-16 10:59:51

[Fortran] 纯文本查看 复制代码

001      subroutine single_step(itimestep, dt, ntotal, hsml, mass, x, vx, 
002     &           u, s,rho, p, t, tdsdt, dx, dvx, du, ds, drho,itype, av)
003 
004c----------------------------------------------------------------------
005c   Subroutine to determine the right hand side of a differential 
006c   equation in a single step for performing time integration 
007 
008c   In this routine and its subroutines the SPH algorithms are performed.
009c     itimestep: Current timestep number                            [in]
010c     dt       : Timestep                                           [in]
011c     ntotal   :  Number of particles                               [in]
012c     hsml     :  Smoothing Length                                  [in]
013c     mass     :  Particle masses                                   [in]
014c     x        :  Particle position                                 [in]
015c     vx       :  Particle velocity                                 [in]
016c     u        :  Particle internal energy                          [in]
017c     s        :  Particle entropy (not used here)                  [in]
018c     rho      :  Density                                       [in/out]
019c     p        :  Pressure                                         [out]
020c     t        :  Temperature                                   [in/out]
021c     tdsdt    :  Production of viscous entropy t*ds/dt            [out]
022c     dx       :  dx = vx = dx/dt                                  [out]
023c     dvx      :  dvx = dvx/dt, force per unit mass                [out]
024c     du       :  du  = du/dt                                      [out]
025c     ds       :  ds  = ds/dt                                      [out]    
026c     drho     :  drho =  drho/dt                                  [out]
027c     itype    :  Type of particle                                 [in]
028c     av       :  Monaghan average velocity                        [out]
029 
030      implicit none
031      include 'param.inc'
032 
033      integer itimestep, ntotal, itype(maxn)   
034      double precision dt, hsml(maxn), mass(maxn), x(dim,maxn),
035     &       vx(dim,maxn), u(maxn), s(maxn), rho(maxn), p(maxn),
036     &       t(maxn), tdsdt(maxn), dx(dim,maxn), dvx(dim,maxn),
037     &       du(maxn), ds(maxn), drho(maxn), av(dim, maxn)             
038      integer i, d, nvirt, niac, pair_i(max_interaction),
039     &        pair_j(max_interaction), ns(maxn)
040      double precision w(max_interaction), dwdx(dim,max_interaction), 
041     &       indvxdt(dim,maxn),exdvxdt(dim,maxn),ardvxdt(dim,maxn), 
042     &       avdudt(maxn), ahdudt(maxn), c(maxn), eta(maxn)            
043 
044      do  i=1,ntotal
045        avdudt(i) = 0.
046        ahdudt(i) = 0.
047c       do  d=1,dim
048          indvxdt(d,i) = 0.
049          ardvxdt(d,i) = 0.
050          exdvxdt(d,i) = 0.
051c        enddo
052      enddo 
053  
054c---  Positions of virtual (boundary) particles:
055 
056      nvirt = 0
057      if (virtual_part) then
058        call virt_part(itimestep, ntotal,nvirt,hsml,mass,x,vx,
059     &       rho,u,p,itype)
060      endif
061      
062c---  Interaction parameters, calculating neighboring particles
063c     and optimzing smoothing length
064 
065      if (nnps.eq.1) then
066        call direct_find(itimestep, ntotal+nvirt,hsml,x,niac,pair_i,
067     &       pair_j,w,dwdx,ns)
068      else if (nnps.eq.2) then
069        call link_list(itimestep, ntotal+nvirt,hsml(1),x,niac,pair_i,
070     &       pair_j,w,dwdx,ns)
071      else if (nnps.eq.3) then
072        call tree_search(itimestep, ntotal+nvirt,hsml,x,niac,pair_i,
073     &       pair_j,w,dwdx,ns)
074      endif        
075                         
076c---  Density approximation or change rate
077       
078      if (summation_density) then     
079        call sum_density(ntotal+nvirt,hsml,mass,niac,pair_i,pair_j,w,
080     &       itype,rho)         
081      else            
082        call con_density(ntotal+nvirt,mass,niac,pair_i,pair_j,
083     &       dwdx,vx, itype,x,rho, drho)        
084      endif
085 
086c---  Dynamic viscosity:
087 
088      if (visc) call viscosity(ntotal+nvirt,itype,x,rho,eta)
089        
090c---  Internal forces:
091  
092      call int_force(itimestep,dt,ntotal+nvirt,hsml,mass,vx,niac,rho,
093     &     eta, pair_i,pair_j,dwdx,u,itype,x,t,c,p,indvxdt,tdsdt,du)
094                   
095c---  Artificial viscosity:
096 
097      if (visc_artificial) call art_visc(ntotal+nvirt,hsml,
098     &      mass,x,vx,niac,rho,c,pair_i,pair_j,w,dwdx,ardvxdt,avdudt)
099       
100c---  External forces:
101 
102      if (ex_force) call ext_force(ntotal+nvirt,mass,x,niac,
103     &                   pair_i,pair_j,itype, hsml, exdvxdt)
104 
105 
106c     Calculating the neighboring particles and undating HSML
107       
108      if (sle.ne.0) call h_upgrade(dt,ntotal, mass, vx, rho, niac,
109     &                   pair_i, pair_j, dwdx, hsml)
110 
111      if (heat_artificial) call art_heat(ntotal+nvirt,hsml,
112     &         mass,x,vx,niac,rho,u, c,pair_i,pair_j,w,dwdx,ahdudt)
113      
114c     Calculating average velocity of each partile for avoiding penetration
115 
116      if (average_velocity) call av_vel(ntotal,mass,niac,pair_i,
117     &                           pair_j, w, vx, rho, av)
118 
119c---  Convert velocity, force, and energy to f and dfdt  
120 
121      do i=1,ntotal
122c        do d=1,dim
123c          dvx(d,i) = indvxdt(d,i) + exdvxdt(d,i) + ardvxdt(d,i)
124c        enddo
125        du(i) = du(i) + avdudt(i) + ahdudt(i)
126      enddo
127 
128      if (mod(itimestep,print_step).eq.0) then     
129          write(*,*)
130          write(*,*) '**** Information for particle ****',
131     &            moni_particle         
132          write(*,101)'internal a ','artifical a=',
133     &            'external a ','total a '  
134          write(*,100)indvxdt(1,moni_particle),ardvxdt(1,moni_particle),
135     &                exdvxdt(1,moni_particle),dvx(1,moni_particle)    
136      endif
137101   format(1x,4(2x,a12))     
138100   format(1x,4(2x,e12.6))     
139 
140      end