program sph_1d cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c c c c c SPH_1D c c c c KEVIN OLSON c c USRA and NCCS, GSFC c c c c April, 1993 c c c c Smooth Particle Hydrodynamics is a gridless particle method for c c solving the Lagrangian equations of hydrodynamics. This c c implementation includes variable smoothing lengths using both the c c algorithm described in Benz (1990) and that described in Steinmetz c c and Muller (1992). Also included is a routine to compute an c c additional thermal diffusion term which is added to the energy c c equation (Monaghan 1988). I have found that in problems where c c the initial conditions are discontinuous, large deviations in the c c energy profile occur at the contact discontinuity which forms near c c the initial discontinuity. Addition of a thermal diffusion term c c removes this problem for the double shock tube problem (Woodward and c c Colella, 1984). Also included in the model are routines for c c both the molecular and bulk viscosities. c c c c This implementation of SPH has been written in Maspar Fortran to run c c on the Maspar MP-1 at GSFC. c c c c The following references where used extensively in building this c c code: c c c c Benz, W., 1990, in Numerical Modelling of Nonlinear Stellar c c Pulsations, Problems and Prospects, ed. J.R. Buchler, Kulwer c c Acdemics Publishers, Dordecht, Boston, London c c c c Hernquist, L., and Katz, N., 1989, ApJ. Suppl., Vol. 70, 419 c c c c Monaghan, J.J., 1985, Comp. Phys. Rep., Vol. 3, 71 c c c c Monaghan, J.J., 1988, "SPH Meets the Shocks of Noh", unpublished c c preprint c c c c Steinmetz, M., and Muller, E., 1992, preprint sub. to A&A c c c c Woodward, P., and Colella, P., 1984, J. Comp. Phys., Vol. 54, 115 c c c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc parameter ( nsize=400+100 ) ! sets the total no. of particles plus ! 100 for the boundries real * 4 xo(nsize),vxo(nsize),uo(nsize) real * 4 vxt(nsize) real * 4 deno(nsize) real * 4 alpha,beta real * 4 c1,c2,tend real * 4 alphab,betab,g1,g2 real * 4 hfac real * 4 ho(nsize) real * 4 sound(nsize) real * 4 maxmu(nsize) real * 4 tote,toteo real xs(128,128),ys(128,128) integer * 4 isort(nsize) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c VARIABLE LIST, THESE DEFINITIONS RETAIN THEIR MEANINGS ACROSS c c SUBROUTINE BOUNDRIES c c c c c c xx,vx,ax: -> position,velocity, and acceleration of a particle c real * 4 xx(nsize),vx(nsize),ax(nsize) c u,au: -> energy per unit mass and its rate of change, du/dt c real * 4 u(nsize),au(nsize) c den,p,mass:-> density and pressure, and mass of a particle c real * 4 den(nsize),p(nsize),mass(nsize) c hh,mass: -> the smoothing length and mass of a particle c real * 4 hh(nsize) c t,delt: -> the time and the time step real * 4 t,delt c divv: -> the divergence of v real * 4 divv(nsize) c thermdiff: -> the additional thermal diffusion term added to du/dt c real * 4 thermdiff(nsize) c gamma: -> the ratio of specific heats c real * 4 gamma c cour: -> the courant number c real * 4 cour c partno: -> the number of a particle in the list c integer * 4 partno(nsize) c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c^^^^^^^^^^^^^^^^^^^^^^BEGIN INITIALIZATIONS^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ n=nsize ! n is the total number of particles delt=1.0e-6 ! the initial value of the time step is set ! to a small value, it is modified as the ! calculation procedes c1=delt/2.0 ! c1 and c2 are constants used in the integration of c2=(delt**2)/4.0 ! the equation of motion tend=0.035 ! the ending time is set gamma=1.4 ! the ratio of specific heats is set cour=0.3 ! the courant number is set nsteps=0 ! a counter which counts the number of time steps taken xmin=0.0 ! the minimum and maximum values of x, the define the xmax=1.0 ! domain of the problem alpha=1.0 ! alpha and beta are the coefficients used in the beta=2.0 ! molecular viscosity alphab=0.5 ! alphab and betab are the coefficients used in the betab=1.5 ! bulk viscosity g1=1.0/10.0 ! g1 and g2 are the coefficients used in the g2=1.0 ! thermal diffusion term hfac=2.0 ! hfac is a factor which determines the initial ! smoothing length of particles, ! hh = hfac * (mass/den) ** (1/d), where d is the ! dimensionality c The following are switches to be set to determine if and which artificial c viscosity is to be used, as well as if the thermal diffusion term in the c energy equation is to be included. A value of 1 indicates that the c switch is on. ibulk=0 ! the bulk viscosity switch inorm=1 ! the molecular viscosity switch itherm=0 ! the thermal diffusion switch ! note: do not use thermal diffusion with entropy ientro=0 ! specifies whether or not the entropic equation ! is to be used in place of the energy equation ientro2=0 ! specifies whether entropy is used in lieu of the ! of the force equation irefl=1 ! specifies whether reflecting boundries are used c the following, twr1-twr8, define the times at which output will c be produced twr1=0.01 twr2=0.015 twr3=0.02 twr4=0.025 twr5=0.03 twr6=0.034 twr7=0.5 twr8=0.5 c^^^^^^^^^^^^^^^^^^^^END OF INITIALIZATIONS^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c^^^^^^^^^^^^^^^^^^^SET UP OF INITIAL CONDITIONS^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ c The following set of subroutines can be called individually to set up c the initial conditions for some different problems. c "double_shock_tube" sets up the initial conditions for the double shock tube c problem of Woodward and Collela c call double_shock_tube (xx,vx,p,u,den,mass,xmax,xmin,hh, c & partno,gamma,hfac,n,ientro) c "sod" sets up the initial conditions for Sod's problem call sod (xx,vx,mass,den,p,u,hh,gamma,xmax,xmin, & partno,hfac,n,ientro) c "sod_const" is used if the distance between particles is constant c call sod_const (xx,vx,mass,den,p,u,hh,gamma,xmax,xmin, c & partno,hfac,n,ientro) c "noh" sets up the initial conditions for the "shocks of noh", i.e. two c colliding streams of gas, note: turn off reflecting boundries when doing this c problem c call noh (xx,vx,mass,den,p,u,hh,gamma,xmax,xmin, c & partno,hfac,n,ientro) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c^^^^^^^^^^^^^^^^^^INITIAL SMOOTHING^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ c Initial smoothing of the pressure (energy) and velocity fields is c performed following Hernquist and Katz (1989). call update_den (xx,mass,den,hh,n) c smoothing the initial energy (entropy) field, u call smooth (xx,hh,mass,den,n,u) c smoothing the initial velocity field, vx call smooth (xx,hh,mass,den,n,vx) p=(gamma-1.0)*u*den ! computes the pressure via the equation of state if (ientro.eq.1) p=u*(den**gamma) if (ientro.eq.1) then p=u call smooth (xx,hh,mass,den,n,p)! if the entropic formulation is used ! compute the pressure correspondingly p=p*(den**gamma) end if cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c^^^^^^^^^^^^^^^^^INITIAL SORT OF PARTICLES^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ c The particles are sorted on x, this is not really necessary for one c dimensional problems since the particles do not become disordered in the c x coordinate. call r4_sort (xx,isort,partno,n) vx=vx(isort) ax=ax(isort) p=p(isort) den=den(isort) hh=hh(isort) mass=mass(isort) u=u(isort) au=au(isort) c^^^^^^^^^^^^^^^^^ALL INITIAL COMPUTATIONS NOW PERFORMED^^^^^^^^^^^^^^^^^^^^^^ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c^^^^^^^^^^^^^^^^^^^^THE MAIN TIME LOOP NOW BEGINS^^^^^^^^^^^^^^^^^^^^^^^^^^^^ do while (t.le.tend) ! do until the time (t) exceeds the ending ! time (tend) c***************************************************************************** c STORE the values of xx,vx,u and hh at the begining of the time step c to be used later in the particle pushing routine "push2" xo=xx vxo=vx uo=u ho=hh deno=den c***************************************************************************** c PUSHING THE PARTICLES 1/2 A TIME STEP AHEAD OF THE ACCELARATIONS call SM_h (xx,hh,mass,den,divv,t,delt,n) ! update values of hh ! according to ! Steinmetz and Muller if (t.eq.0.0) go to 372 ! at the first time step, ! skip the 1/2 step call push1 (xx,u,vx,ax,au,c1,c2,n) ! push the particles 1/2 a time ! step ahead c call enz_h1 (hh,divv,c1,n) ! advance hh 1/2 a time step ahead ! according to Benz's alg. c***************************************************************************** c SORTING 372 call r4_sort (xx,isort,partno,n) ! sort on x after push vx=vx(isort) ax=ax(isort) p=p(isort) den=den(isort) hh=hh(isort) mass=mass(isort) u=u(isort) au=au(isort) xo=xo(isort) vxo=vxo(isort) uo=uo(isort) ho=ho(isort) c***************************************************************************** c SETTING BOUNDRY CONDITIONS, subroutine "reflecting" creates c reflecting (solid wall) boundry conditions if (irefl.eq.1) then call reflecting (xx,vx,xo,vxo,den,mass,p,u,uo, & hh,ho,n) c where (partno.le.50.or.partno.gt.(n-50)) c vx=0.0 c ax=0.0 c end where c where (partno.le.50) c p=p(50) c u=u(50) c den=den(50) c end where c where (partno.gt.(n-50)) c p=p(n-49) c u=u(n-49) c den=den(n-49) c end where end if c***************************************************************************** c the total energy if (t.eq.0.0) then vxt=vx vxt=0.5*(vxt**2) if (ientro.eq.0) then if (irefl.eq.1) then toteo=sum(mass*u,mask=partno.gt.51.and. & partno.lt.(n-50)) else toteo=sum(mass*u) end if else if (irefl.eq.1) then toteo=sum(mass*u*(den**(gamma-1.0))/(gamma-1.0), & mask=partno.gt.51.and.partno.lt.(n-50)) else toteo=sum(mass*u*(den**(gamma-1.0))/(gamma-1.0)) end if end if if (irefl.eq.1) then toteo=toteo+sum(mass*vxt,mask=partno.gt.51.and.partno. & lt.(n-50)) else toteo=(toteo+sum(mass*vxt)) end if end if c***************************************************************************** c SOUND SPEED sound=sqrt(gamma*p/den) c***************************************************************************** c THERMAL DIFFUSION c Here, the subroutine "therm_diffusion" is used, which computes c an additional term to be added to the energy equation for each c particle as described by Monaghan (1988), unpublished preprint. if (itherm.eq.1) then ! if the switch is set, compute it call therm_diffusion (xx,vx,hh,mass,den,p,u,divv,sound,n & ,g1,g2,maxmu,gamma,thermdiff) end if c***************************************************************************** c ACCELERATIONS AND DU/DT c The subroutine "accelerations" computes the acceleration of each c particle (ax) due to pressure forces and artificial viscosity as well as c heating of the particles (au = du/dt). call accelerations (xx,vx,ax,au,hh,p,u,mass,den,divv, & sound,maxmu,n,inorm,alpha,beta,ibulk,alphab,betab, & gamma,ientro,ientro2) c if (irefl.eq.1) then c where (partno.le.50.or.partno.gt.(n-50)) ax=0.0 c end if if (ientro.eq.1) au=au*u*den*(gamma-1.0)/p c***************************************************************************** if (itherm.eq.1) then ! if the switch is set, add the thermal ! diffusion term au=au+thermdiff end if c***************************************************************************** c PUSH THE PARTICLES the rest of the way through the time step call push2 (xx,u,vx,ax,au,c1,c2,delt,xo,vxo,uo,t,n) c***************************************************************************** c UPDATE HH, following Benz (1990) c call benz_h2 (hh,delt,divv,ho,n) c***************************************************************************** c SORTING into x order after pushing particles call r4_sort (xx,isort,partno,n) vx=vx(isort) ax=ax(isort) p=p(isort) den=den(isort) hh=hh(isort) mass=mass(isort) u=u(isort) au=au(isort) c***************************************************************************** c UPDATE THE DENSITY AND PRESSURE with the new positions and new values c of the energy density. call update_den (xx,mass,den,hh,n) if (ientro.eq.1) then p=u call smooth (xx,hh,mass,den,n,p)! if the entropic formulation is used ! compute the pressure correspondingly p=p*(den**gamma) else call update_p (hh,mass,u,den,p,gamma,n,ientro) end if where (p.le.0.0) p=0.0 sound=sqrt(gamma*p/den) ! a new sound speed is computed for each ! particle c***************************************************************************** c total energy vxt=vx vxt=0.5*(vxt**2) if (ientro.eq.0) then if (irefl.eq.1) then tote=sum(mass*u,mask=partno.gt.51.and.partno.lt.(n-50)) else tote=sum(mass*u) end if else if (irefl.eq.1) then tote=sum(mass*u*(den**(gamma-1.0))/(gamma-1.0), & mask=partno.gt.51.and.partno.lt.(n-50)) else tote=sum(mass*u*(den**(gamma-1.0))/(gamma-1.0)) end if end if if (irefl.eq.1) then tote=tote+sum(mass*vxt,mask=partno.gt.51.and. & partno.lt.(n-50)) else tote=(tote+sum(mass*vxt)) end if perc=abs(toteo-tote)/toteo perc=100.0*perc c print *,' tote = ',tote,' perc = ',perc c**************************************************************************** c BOUNDRY conditions are reset since the paticles have been pushed a c second time if (irefl.eq.1) then call reflecting (xx,vx,xo,vxo,den,mass,p,u,uo, & hh,ho,n) c where (partno.le.50.or.partno.gt.(n-50)) c vx=0.0 c ax=0.0 c end where c where (partno.le.50) c p=p(50) c u=u(50) c den=den(50) c end where c where (partno.gt.(n-50)) c p=p(n-49) c u=u(n-49) c den=den(n-49) c end where c where (partno.gt.(n-50)) p=p(n-50) end if c***************************************************************************** c computing the DIVERGENCE OF THE VELOCITY FIELD following Monaghan (1988) c unpublished preprint. call div_v (xx,vx,hh,mass,den,n,divv) c***************************************************************************** c A NEW TIME STEP IS COMPUTED according to the courant condition following c Hernquist and Katz (1989), also Monaghan (1988). call new_delt (hh,divv,sound,inorm,alpha,beta, & ibulk,alphab,betab,maxmu,n,c1,c2,cour,delt) c***************************************************************************** c 62 print *,' t= ',t,' nsteps ',nsteps write (40,*) t,tote,perc c if (nt.ge.100.or.t.eq.0.0) then c xs=reshape(source=xx,shape=(/128,128/)) c ys=reshape(source=den,shape=(/128,128/)) c call SEND_VA(xs,128,128) c call SEND_VA(ys,128,128) c ys=reshape(source=p,shape=(/128,128/)) c call SEND_VA(ys,128,128) c ys=reshape(source=vx,shape=(/128,128/)) c call SEND_VA(ys,128,128) c nt=0 c end if nt=nt+1 nsteps=nsteps+1 c***************************************************************************** c OUTPUT if (t.gt.twr1.and.t.le.(twr1+delt)) then do i=1,n write (91,25) xx(i),den(i),p(i),vx(i),u(i),hh(i) end do end if if (t.gt.twr2.and.t.le.(twr2+delt)) then do i=1,n write (92,25) xx(i),den(i),p(i),vx(i),u(i),hh(i) end do end if if (t.gt.twr3.and.t.le.(twr3+delt)) then do i=1,n write (93,25) xx(i),den(i),p(i),vx(i),u(i),hh(i) end do end if if (t.gt.twr4.and.t.le.(twr4+delt)) then do i=1,n write (94,25) xx(i),den(i),p(i),vx(i),u(i),hh(i) end do end if if (t.gt.twr5.and.t.le.(twr5+delt)) then do i=1,n write (95,25) xx(i),den(i),p(i),vx(i),u(i),hh(i) end do end if if (t.gt.twr6.and.t.le.(twr6+delt)) then do i=1,n write (96,25) xx(i),den(i),p(i),vx(i),u(i),hh(i) end do end if if (t.gt.twr7.and.t.le.(twr7+delt)) then do i=1,n write (97,25) xx(i),den(i),p(i),vx(i),u(i),hh(i) end do end if if (t.gt.twr8.and.t.le.(twr8+delt)) then do i=1,n write (98,25) xx(i),den(i),p(i),vx(i),u(i),hh(i) end do end if 25 format(6(1x,f10.5)) 26 format(6(1x,e11.5)) c***************************************************************************** t=t+delt ! increment the time end do c^^^^^^^^^^^^^^^^^^^END OF TIME LOOP^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ c***************************************************************************** 2 stop end c^^^^^^^^^^^^^^^^^END OF PROGRAM SPH-1D^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc