求大神解答，什么问题
subroutine single_step(itimestep, dt, ntotal, hsml, mass, x, vx,& u, s,rho, p, t, tdsdt, dx, dvx, du, ds, drho,itype, av)
c----------------------------------------------------------------------
c Subroutine to determine the right hand side of a differential
c equation in a single step for performing time integration
c In this routine and its subroutines the SPH algorithms are performed.
c itimestep: Current timestep number
c dt : Timestep
c ntotal :Number of particles
c hsml :Smoothing Length
c mass :Particle masses
c x :Particle position
c vx :Particle velocity
c u :Particle internal energy
c s :Particle entropy (not used here)
c rho :Density
c p :Pressure
c t :Temperature
c tdsdt :Production of viscous entropy t*ds/dt
c dx :dx = vx = dx/dt
c dvx :dvx = dvx/dt, force per unit mass
c du :du= du/dt
c ds :ds= ds/dt
c drho :drho =drho/dt
c itype :Type of particle
c av :Monaghan average velocity
implicit none
include 'param.inc'
integer itimestep, ntotal, itype(maxn)
double precision dt, hsml(maxn), mass(maxn), x(dim,maxn),
& vx(dim,maxn), u(maxn), s(maxn), rho(maxn), p(maxn),
& t(maxn), tdsdt(maxn), dx(dim,maxn), dvx(dim,maxn),
& du(maxn), ds(maxn), drho(maxn), av(dim, maxn)
integer i, d, nvirt, niac, pair_i(max_interaction),
& pair_j(max_interaction), ns(maxn)
double precision w(max_interaction), dwdx(dim,max_interaction),
& indvxdt(dim,maxn),exdvxdt(dim,maxn),ardvxdt(dim,maxn),
& avdudt(maxn), ahdudt(maxn), c(maxn), eta(maxn)
doi=1,ntotal
avdudt(i) = 0.
ahdudt(i) = 0.
c dod=1,dim
indvxdt(d,i) = 0.
ardvxdt(d,i) = 0.
exdvxdt(d,i) = 0.
c enddo
enddo
c---Positions of virtual (boundary) particles:
nvirt = 0
if (virtual_part) then
call virt_part(itimestep, ntotal,nvirt,hsml,mass,x,vx,
& rho,u,p,itype)
endif
c---Interaction parameters, calculating neighboring particles
c and optimzing smoothing length
if (nnps.eq.1) then
call direct_find(itimestep, ntotal+nvirt,hsml,x,niac,pair_i,
& pair_j,w,dwdx,ns)
else if (nnps.eq.2) then
call link_list(itimestep, ntotal+nvirt,hsml(1),x,niac,pair_i,
& pair_j,w,dwdx,ns)
else if (nnps.eq.3) then
call tree_search(itimestep, ntotal+nvirt,hsml,x,niac,pair_i,
& pair_j,w,dwdx,ns)
endif
c---Density approximation or change rate
if (summation_density) then
call sum_density(ntotal+nvirt,hsml,mass,niac,pair_i,pair_j,w,
& itype,rho)
else
call con_density(ntotal+nvirt,mass,niac,pair_i,pair_j,
& dwdx,vx, itype,x,rho, drho)
endif
c---Dynamic viscosity:
if (visc) call viscosity(ntotal+nvirt,itype,x,rho,eta)
c---Internal forces:
call int_force(itimestep,dt,ntotal+nvirt,hsml,mass,vx,niac,rho,
& eta, pair_i,pair_j,dwdx,u,itype,x,t,c,p,indvxdt,tdsdt,du)
c---Artificial viscosity:
if (visc_artificial) call art_visc(ntotal+nvirt,hsml,
& mass,x,vx,niac,rho,c,pair_i,pair_j,w,dwdx,ardvxdt,avdudt)
c---External forces:
if (ex_force) call ext_force(ntotal+nvirt,mass,x,niac,
& pair_i,pair_j,itype, hsml, exdvxdt)
c Calculating the neighboring particles and undating HSML
if (sle.ne.0) call h_upgrade(dt,ntotal, mass, vx, rho, niac,
& pair_i, pair_j, dwdx, hsml)
if (heat_artificial) call art_heat(ntotal+nvirt,hsml,
& mass,x,vx,niac,rho,u, c,pair_i,pair_j,w,dwdx,ahdudt)
c Calculating average velocity of each partile for avoiding penetration
if (average_velocity) call av_vel(ntotal,mass,niac,pair_i,
& pair_j, w, vx, rho, av)
c---Convert velocity, force, and energy to f and dfdt
do i=1,ntotal
c do d=1,dim
c dvx(d,i) = indvxdt(d,i) + exdvxdt(d,i) + ardvxdt(d,i)
c enddo
du(i) = du(i) + avdudt(i) + ahdudt(i)
enddo
if (mod(itimestep,print_step).eq.0) then
write(*,*)
write(*,*) '**** Information for particle ****',
& moni_particle
write(*,101)'internal a ','artifical a=',
& 'external a ','total a '
write(*,100)indvxdt(1,moni_particle),ardvxdt(1,moni_particle),
& exdvxdt(1,moni_particle),dvx(1,moni_particle)
endif
101 format(1x,4(2x,a12))
100 format(1x,4(2x,e12.6))
end
缺少文件 param.inc 找不到函数 tree_search
页:
[1]