! $Id$ ! ! This module takes care of direct N-body gravity between particles. ! !** AUTOMATIC CPARAM.INC GENERATION **************************** ! ! Declare (for generation of cparam.inc) the number of f array ! variables and auxiliary variables added by this module ! ! MPVAR CONTRIBUTION 3 ! CPARAM logical, parameter :: lparticles_nbody=.true. ! !*************************************************************** module Particles_nbody ! use Cdata use General, only: keep_compiler_quiet use Messages use Particles_cdata use Particles_map use Particles_sub ! implicit none ! include 'particles_nbody.h' ! real, dimension(nspar,mspvar) :: fsp real, dimension(nspar) :: xsp0=0.0, ysp0=0.0, zsp0=0.0 real, dimension(nspar) :: vspx0=0.0, vspy0=0.0, vspz0=0.0 real, dimension(nspar) :: pmass=0.0, r_smooth=0.0, pmass1 real, dimension(nspar) :: accrete_hills_frac=0.2, final_ramped_mass=0.0 real :: eccentricity=0.0, semimajor_axis=1.0 real :: delta_vsp0=1.0, totmass, totmass1 real :: create_jeans_constant=0.25, GNewton1 real :: GNewton=impossible, density_scale=0.001 real :: cdtpnbody=0.1 real :: hills_tempering_fraction=0.8 real, pointer :: rhs_poisson_const, tstart_selfgrav integer :: ramp_orbits=5, mspar_orig=1 integer :: iglobal_ggp=0, istar=1, iplanet=2, imass_nbody=0 integer :: maxsink=10*nspar, icreate=100 logical, dimension(nspar) :: lcylindrical_gravity_nbody=.false. logical, dimension(nspar) :: lfollow_particle=.false., laccretion=.false. logical, dimension(nspar) :: ladd_mass=.false. logical :: lcalc_orbit=.true., lbackreaction=.false., lnorm=.true. logical :: lreset_cm=.false., lnogravz_star=.false., lexclude_frozen=.false. logical :: lnoselfgrav_star=.true. logical :: lramp=.false., lcreate_sinks=.false., lcreate_gas=.true. logical :: ldt_nbody=.true., lcreate_dust=.true. logical :: linterpolate_gravity=.false., linterpolate_linear=.true. logical :: linterpolate_quadratic_spline=.false. logical :: laccrete_when_create=.true. logical :: ldust=.false. logical :: ltempering=.false. logical :: lretrograde=.false. logical :: lgas_gravity=.false.,ldust_gravity=.false. logical :: linertial_frame=.true., lvectorgravity=.false. character (len=labellen) :: initxxsp='random', initvvsp='nothing' ! namelist /particles_nbody_init_pars/ & initxxsp, initvvsp, xsp0, ysp0, zsp0, vspx0, vspy0, vspz0, delta_vsp0, & pmass, r_smooth, lcylindrical_gravity_nbody, lexclude_frozen, GNewton, & bcspx, bcspy, bcspz, ramp_orbits, lramp, final_ramped_mass, density_scale, & linterpolate_gravity, linterpolate_quadratic_spline, laccretion, & accrete_hills_frac, istar, maxsink, lcreate_sinks, icreate, lcreate_gas, & lcreate_dust, ladd_mass, laccrete_when_create, ldt_nbody, cdtpnbody, lretrograde, & linertial_frame, eccentricity, semimajor_axis, lvectorgravity ! namelist /particles_nbody_run_pars/ & dsnap_par_minor, linterp_reality_check, lcalc_orbit, lreset_cm, & lnogravz_star, lfollow_particle, lbackreaction, lexclude_frozen, & GNewton, bcspx, bcspy, bcspz, density_scale, lnoselfgrav_star, & linterpolate_quadratic_spline, laccretion, accrete_hills_frac, istar, & maxsink, lcreate_sinks, icreate, lcreate_gas, lcreate_dust, ladd_mass, & laccrete_when_create, ldt_nbody, cdtpnbody, hills_tempering_fraction, & ltempering, lgas_gravity, ldust_gravity, linertial_frame, lvectorgravity ! integer, dimension(nspar,3) :: idiag_xxspar=0,idiag_vvspar=0 integer, dimension(nspar) :: idiag_torqint=0,idiag_torqext=0 integer, dimension(nspar) :: idiag_torqext_gas=0,idiag_torqext_par=0 integer, dimension(nspar) :: idiag_torqint_gas=0,idiag_torqint_par=0 integer :: idiag_totenergy=0,idiag_totangmom=0 ! contains !*********************************************************************** subroutine register_particles_nbody() ! ! Set up indices for access to the f and fsp. ! ! 27-aug-06/wlad: adapted ! if (lroot) call svn_id( & "$Id$") ! ! Set up mass as particle index. Plus seven, since the other 6 are ! used by positions and velocities. ! imass_nbody=npvar+1 ! ! No need to solve the N-body equations for non-N-body problems. ! if (nspar < 2) then if (lroot) write(0,*) 'nspar = ', nspar call fatal_error('register_particles_nbody','the number of massive'//& ' particles is less than 2. There is no need to use the'//& ' N-body code. Consider setting gravity as a global variable'//& ' using gravity_r.f90 instead.') endif ! if (npar < nspar) then if (lroot) write(0,*) 'npar, nspar = ', npar, nspar call fatal_error('register_particles_nbody','the number of massive'//& ' particles (nspar) is less than the allocated number of particles'//& ' (npar). Increase npar to the minimum number (npar=npsar) needed'//& ' in cparam.local and recompile') endif ! ! Auxiliary variables for polar coordinates ! !if (.not.lcartesian_coords) then call append_npvar('ivpx_cart',ivpx_cart) call append_npvar('ivpy_cart',ivpy_cart) call append_npvar('ivpz_cart',ivpz_cart) !endif ! endsubroutine register_particles_nbody !*********************************************************************** subroutine initialize_particles_nbody(f) ! ! Perform any post-parameter-read initialization i.e. calculate derived ! parameters. ! ! 27-aug-06/wlad: adapted ! use FArrayManager use SharedVariables ! real, dimension (mx,my,mz,mfarray) :: f ! integer :: ks ! ! Look for initialized masses. ! mspar_orig=0 do ks=1,nspar if (pmass(ks)/=0.0) then mspar_orig=mspar_orig+1 ipar_nbody(mspar_orig)=ks endif enddo ! ! Just needs to do this starting, otherwise mspar will be overwritten after ! reading the snapshot ! if (lstart) then mspar=mspar_orig else !else read pmass from fsp pmass(1:mspar)=fsp(1:mspar,imass_nbody) endif ! ! When first called, mspar was zero, so no diagnostic index was written to ! index.pro ! call rprint_particles_nbody(.false.,LWRITE=lroot) ! ! G_Newton. Overwrite the one set by start.in if set again here, ! because I might want units in which both G and GM are 1. ! If we are solving for selfgravity, get that one instead. ! if (GNewton == impossible) then !ONLY REASSIGN THE GNEWTON !IF IT IS NOT SET IN START.IN!!!!! if (lselfgravity) then call get_shared_variable('rhs_poisson_const',rhs_poisson_const,caller='initialize_particles_nbody') GNewton=rhs_poisson_const/(4*pi) else GNewton=G_Newton endif endif ! GNewton1=1./GNewton ! ! Get the variable tstart_selfgrav from selfgrav, so we know when to create ! sinks. ! if (lselfgravity) & call get_shared_variable('tstart_selfgrav',tstart_selfgrav,caller='initialize_particles_nbody') ! ! inverse mass ! if (lstart) then if (lramp) then do ks=1,mspar if (ks/=istar) pmass(ks) = epsi enddo pmass(istar)=1-epsi*(mspar-1) endif else !read imass_nbody from the snapshot pmass(1:mspar)=fsp(1:mspar,imass_nbody) endif ! pmass1=1./max(pmass,tini) ! ! inverse total mass ! totmass=sum(pmass) totmass1=1./max(totmass, tini) ! ! Check for consistency between the cylindrical gravities switches ! from cdata and the one from nbody. ! if (((lcylindrical_gravity).and.& (.not.lcylindrical_gravity_nbody(istar))).or.& (.not.lcylindrical_gravity).and.& (lcylindrical_gravity_nbody(istar))) then call fatal_error('initialize_particles_nbody','inconsitency '//& 'between lcylindrical_gravity from cdata and the '//& 'one from nbody') endif ! if (rsmooth/=r_smooth(istar)) then print*,'rsmooth from cdata=',rsmooth print*,'r_smooth(istar)=',r_smooth(istar) call fatal_error('initialize_particles_nbody','inconsitency '//& 'between rsmooth from cdata and the '//& 'one from nbody') endif ! ! Check for consistency between the poisson and integral ! calculations of selfgravity. ! if (lbackreaction.and.lselfgravity) then print*,'self-gravity is already taken into '//& 'account. lbackreaction is a way of '//& 'integrating the self-gravity only in few '//& 'specific points of the grid' call fatal_error('initialize_particles_nbody','') endif ! ! The presence of dust particles needs to be known. ! if (npar > nspar) ldust=.true. if (linterpolate_gravity.and.(.not.ldust)) then if (lroot) print*,'interpolate gravity is just'//& ' for the dust component. No need for it if'//& ' you are not using dust particles' call fatal_error('initialize_particles_nbody','') endif if (any(laccretion).and.linterpolate_gravity) then if (lroot) print*,'interpolate gravity not yet '//& 'implemented in connection with accretion' call fatal_error('initialize_particles_nbody','') endif ! if (linterpolate_gravity) then if (lroot) print*,'initializing global array for nbody gravity' call farray_register_global('ggp',iglobal_ggp,vector=3) endif ! if (linterpolate_quadratic_spline) & linterpolate_linear=.false. ! if ((.not.linterpolate_gravity).and.& linterpolate_quadratic_spline) then if (lroot) print*,'no need for linterpolate_quadratic_spline'//& 'if linterpolate_gravity is false' call fatal_error('initialize_particles_nbody','') endif ! if ((.not.linertial_frame).and.((.not.lcartesian_coords).or.mspar>2)) & call fatal_error('initialize_particles_nbody','Fixed star only for Cartesian and mspar=2') ! call keep_compiler_quiet(f) ! endsubroutine initialize_particles_nbody !*********************************************************************** subroutine pencil_criteria_par_nbody() ! ! All pencils that the Particles_nbody module depends on are specified here. ! ! 22-sep-06/wlad: adapted ! if (ldensity) lpenc_requested(i_rho)=.true. if (ldust.and.lbackreaction) & lpenc_requested(i_rhop)=.true. ! if (lexclude_frozen) then if (lcylinder_in_a_box) then lpenc_requested(i_rcyl_mn)=.true. else lpenc_requested(i_r_mn)=.true. endif endif ! if (idiag_totenergy/=0) then lpenc_diagnos(i_u2)=.true. lpenc_diagnos(i_rho)=.true. endif ! !if (any(idiag_torqext/=0) .or. any(idiag_torqint/=0)) then ! lpenc_diagnos(i_rcyl_mn)=.true. !endif ! lpenc_diagnos(i_rcyl_mn)=.true. lpenc_diagnos(i_rhop)=.true. ! endsubroutine pencil_criteria_par_nbody !*********************************************************************** subroutine pencil_interdep_par_nbody(lpencil_in) ! ! Interdependency among pencils provided by the Particles_nbody module ! is specified here. ! ! 22-sep-06/wlad: adapted ! logical, dimension(npencils) :: lpencil_in ! call keep_compiler_quiet(lpencil_in) ! endsubroutine pencil_interdep_par_nbody !*********************************************************************** subroutine calc_pencils_par_nbody(f,p) ! ! Calculate nbody particle pencils ! ! 22-sep-06/wlad: adapted ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p ! call keep_compiler_quiet(f) call keep_compiler_quiet(p) ! endsubroutine calc_pencils_par_nbody !*********************************************************************** subroutine init_particles_nbody(f,fp) ! ! Initial positions and velocities of nbody particles. ! Overwrite the position asserted by the dust module ! ! 17-nov-05/anders+wlad: adapted ! use General, only: random_number_wrapper use Sub ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mpar_loc,mparray) :: fp real, dimension(mspar) :: kep_vel,sma real, dimension(mspar,3) :: velocity real, dimension(mspar,3) :: positions real :: tmp,parc integer :: k,ks ! intent (in) :: f intent (out) :: fp ! ! Break if there are N-body particles without mass. ! if (mspar < nspar) then call warning("init_particles_nbody", 'mspar < nspar:') print*, 'mspar = ', mspar, ', nspar = ', nspar print*, 'nspar is the number of allocated n-body particles' print*, 'mspar is the number of those particles with' print*, ' non-zero mass in this processor.' print*, '' print*, 'you should set pmass to non-zero values in start.in' endif ! ! Shortcuts ! if (mspar > 0) then positions(1:nspar,1) = xsp0 ; velocity(1:nspar,1) = vspx0 positions(1:nspar,2) = ysp0 ; velocity(1:nspar,2) = vspy0 positions(1:nspar,3) = zsp0 ; velocity(1:nspar,3) = vspz0 endif ! ! Initialize particles' positions. ! select case (initxxsp) ! case ('nothing') if (lroot) print*, 'init_particles_nbody: nothing' ! case ('origin') if (lroot) then print*, 'init_particles_nbody: All nbody particles at origin' fp(1:mspar,ixp:izp)=0.0 endif ! case ('constant') if (lroot) & print*, 'init_particles_nbody: Place nbody particles at x,y,z=', & xsp0, ysp0, zsp0 do k=1,npar_loc if (ipar(k)<=mspar) fp(k,ixp:izp)=positions(ipar(k),1:3) enddo ! case ('random') if (lroot) print*, 'init_particles_nbody: Random particle positions' do ks=1,mspar if (nxgrid/=1) call random_number_wrapper(positions(ks,ixp)) if (nygrid/=1) call random_number_wrapper(positions(ks,iyp)) if (nzgrid/=1) call random_number_wrapper(positions(ks,izp)) enddo ! if (nxgrid/=1) & positions(1:mspar,ixp)=xyz0_loc(1)+positions(1:mspar,ixp)*Lxyz_loc(1) if (nygrid/=1) & positions(1:mspar,iyp)=xyz0_loc(2)+positions(1:mspar,iyp)*Lxyz_loc(2) if (nzgrid/=1) & positions(1:mspar,izp)=xyz0_loc(3)+positions(1:mspar,izp)*Lxyz_loc(3) ! ! Loop through ipar to allocate the nbody particles. ! do k=1,npar_loc if (ipar(k) <= mspar) then if (ldebug) & print*,'initparticles_nbody. Slot for nbody particle ',& ipar(k),' was at fp position ',k,' at processor ',iproc ! fp(k,ixp:izp)=positions(ipar(k),1:3) ! ! Correct for non-existing dimensions (not really needed, I think). ! if (nxgrid==1) fp(k,ixp)=x(nghost+1) if (nygrid==1) fp(k,iyp)=y(nghost+1) if (nzgrid==1) fp(k,izp)=z(nghost+1) ! if (ldebug) & print*,'initparticles_nbody. Nbody particle ',& ipar(k),' located at ixp=',fp(k,ixp) endif enddo ! case ('fixed-cm') if (lgrav) then print*,'a gravity module is being used. Are you using '//& 'both a fixed central gravity and nbody gravity? '//& 'better stop and check' call fatal_error('init_particles_nbody','') endif ! ! Ok, I have the masses and the positions of all massive particles ! except the last, which will have a position determined to fix the ! center of mass on the center of the grid. ! if (any(ysp0/=0)) then if (lspherical_coords) then call fatal_error('init_particles_nbody','not yet generalized'//& ' for non-zero initial inclinations') else call fatal_error('init_particles_nbody','not yet generalized'//& ' for non-zero azimuthal initial position') endif endif if (any(zsp0/=0)) then if (lspherical_coords) then call fatal_error('init_particles_nbody','not yet generalized'//& ' for non-zero azimuthal initial position') else call fatal_error('init_particles_nbody','nbody code not'//& ' yet generalized to allow initial inclinations') endif endif if (lcylindrical_coords.or.lspherical_coords) then if (any(xsp0<0)) & call fatal_error('init_particles_nbody', & 'in cylindrical coordinates '//& 'all the radial positions must be positive') endif ! if (lroot) then print*,'init_particles_nbody: fixed-cm - mass and position arranged' print*,' so that the center of mass is at rest' print*,' at the center of the grid.' endif ! if (lspherical_coords) then if (lroot) print*,'put all particles in the midplane' positions(1:mspar,iyp)=pi/2 endif ! tmp = 0.;parc=0 do ks=1,mspar if (ks/=istar) then sma(ks)=abs(positions(ks,1)) tmp=tmp+pmass(ks) parc = parc - sma(ks)*pmass(ks) endif enddo ! ! Fixed-cm assumes that the total mass is always one. The mass of the ! star is adjusted to ensure this. ! pmass(istar)=1.- tmp;pmass1=1./max(pmass,tini);totmass=1.;totmass1=1. ! parc = parc*totmass1 if (tmp >= 1.0) & call fatal_error('init_particles_nbody', & 'The mass of one '//& '(or more) of the particles is too big! The masses should '//& 'never be bigger than g0. Please scale your assemble so that '//& 'the combined mass of the (n-1) particles is less than that. '//& 'The mass of the last particle in the pmass array will be '//& 'reassigned to ensure that the total mass is g0') ! do ks=1,mspar if (ks/=istar) & positions(ks,1)=sign(1.,positions(ks,1))* (sma(ks) + parc) enddo ! ! The last one (star) fixes the CM at Rcm=zero ! if (lcartesian_coords) then positions(istar,1)=parc elseif (lcylindrical_coords) then !put the star in positive coordinates, with pi for azimuth positions(istar,1)=abs(parc) positions(istar,2)=pi elseif (lspherical_coords) then positions(istar,1)=abs(parc) positions(istar,3)=pi endif ! if (ldebug) then print*,'pmass =',pmass print*,'position (x)=',positions(:,1) print*,'position (y)=',positions(:,2) print*,'position (z)=',positions(:,3) endif ! ! Loop through ipar to allocate the nbody particles ! do k=1,npar_loc if (ipar(k) <= mspar) then ! if (ldebug) & print*,'initparticles_nbody. Slot for nbody particle ',& ipar(k),' was at fp position ',k,' at processor ',iproc ! ! Here I substitute the first mspar dust particles by massive ones, ! since the first ipars are less than mspar ! fp(k,ixp:izp)=positions(ipar(k),1:3) ! ! Correct for non-existing dimensions (not really needed, I think) ! if (nxgrid==1) fp(k,ixp)=x(nghost+1) if (nygrid==1) fp(k,iyp)=y(nghost+1) if (nzgrid==1) fp(k,izp)=z(nghost+1) ! if (ldebug) & print*,'init_particles_nbody. Nbody particle ',ipar(k),& ' located at ixp=',fp(k,ixp) endif enddo ! case ('eccentric') ! ! Coded only for 2 bodies ! if (mspar /= 2) call fatal_error("init_particles_nbody",& "This initial condition is currently coded for 2 massive particles only.") ! ! Define iplanet. istar=1 and iplanet=2 is default ! if (istar == 2) iplanet=1 ! ! Radial position at barycentric coordinates. Start both at apocenter, ! ! r_i=(1+e)*a_i, where a_i = sma * m_j /(mi+mj) ! ! See, i.e., Murray & Dermott, p.45, barycentric orbits. ! positions(iplanet,1)=(1+eccentricity) * semimajor_axis * pmass( istar)/totmass positions( istar,1)=(1+eccentricity) * semimajor_axis * pmass(iplanet)/totmass ! ! Azimuthal position. Planet and star phased by pi. ! positions(iplanet,2)=0 positions( istar,2)=pi ! do k=1,npar_loc if (ipar(k) <= mspar) then fp(k,ixp:izp) = positions(ipar(k),1:3) endif enddo ! case default if (lroot) print*,'init_particles_nbody: No such such value for'//& ' initxxsp: ',trim(initxxsp) call fatal_error('init_particles_nbody','') ! endselect ! ! Redistribute particles among processors (now that positions are determined). ! call boundconds_particles(fp,ipar) ! ! Initial particle velocity. ! select case (initvvsp) ! case ('nothing') if (lroot) print*, 'init_particles_nbody: No particle velocity set' ! case ('zero') if (lroot) then print*, 'init_particles_nbody: Zero particle velocity' fp(1:mspar,ivpx:ivpz)=0. endif ! case ('constant') if (lroot) then print*, 'init_particles_nbody: Constant particle velocity' print*, 'init_particles_nbody: vspx0, vspy0, vspz0=', vspx0, vspy0, vspz0 endif do k=1,npar_loc if (ipar(k) <= mspar) & fp(k,ivpx:ivpz)=velocity(ipar(k),1:3) enddo ! case ('fixed-cm') ! ! Keplerian velocities for the planets ! parc=0. do ks=1,mspar if (ks/=istar) then kep_vel(ks)=sqrt(1./sma(ks)) !circular velocity parc = parc - kep_vel(ks)*pmass(ks) endif enddo parc = parc*totmass do ks=1,mspar if (ks/=istar) then if (lcartesian_coords) then velocity(ks,2) = sign(1.,positions(ks,1))*(kep_vel(ks) + parc) elseif (lcylindrical_coords) then !positive for the planets velocity(ks,2) = abs(kep_vel(ks) + parc) elseif (lspherical_coords) then velocity(ks,3) = abs(kep_vel(ks) + parc) endif endif enddo ! ! The last one (star) fixes the CM also with velocity zero ! if (lcartesian_coords) then velocity(istar,2)=parc elseif (lcylindrical_coords) then velocity(istar,2)=-parc elseif (lspherical_coords) then velocity(istar,3)=-parc endif ! ! Revert all velocities if retrograde ! if (lretrograde) velocity=-velocity ! ! Loop through ipar to allocate the nbody particles ! do k=1,npar_loc if (ipar(k)<=mspar) then if (ldebug) & print*,& 'init_particles_nbody. Slot for nbody particle ',ipar(k),& ' was at fp position ', k, ' at processor ', iproc ! fp(k,ivpx:ivpz) = velocity(ipar(k),1:3) ! if (ldebug) & print*,'init_particles_nbody. Nbody particle ',ipar(k),& ' has velocity y=',fp(k,ivpy) ! endif enddo ! case ('eccentric') ! ! Coded only for 2 bodies ! if (mspar /= 2) call fatal_error("init_particles_nbody",& "This initial condition is currently coded for 2 massive particles only.") ! ! Define iplanet. istar=1 and iplanet=2 is default. ! if (istar == 2) iplanet=1 velocity(iplanet,2) = sqrt((1-eccentricity)/(1+eccentricity) * GNewton/semimajor_axis) * pmass( istar)/totmass velocity( istar,2) = sqrt((1-eccentricity)/(1+eccentricity) * GNewton/semimajor_axis) * pmass(iplanet)/totmass ! ! Revert all velocities if retrograde. ! if (lretrograde) velocity=-velocity ! ! Loop through particles to allocate the velocities. ! do k=1,npar_loc if (ipar(k)<=mspar) fp(k,ivpx:ivpz) = velocity(ipar(k),1:3) enddo ! case default if (lroot) print*, 'init_particles_nbody: No such such value for'//& 'initvvsp: ',trim(initvvsp) call fatal_error('init_particles_nbody','') ! endselect ! ! Make the particles known to all processors ! call boundconds_particles(fp,ipar) call bcast_nbodyarray(fp) ! call keep_compiler_quiet(f) ! endsubroutine init_particles_nbody !*********************************************************************** subroutine dvvp_dt_nbody_pencil(f,df,fp,dfp,p,ineargrid) ! ! Add gravity from the particles to the gas and the backreaction from the ! gas onto the particles. ! ! Adding the gravity on the dust component via the grid is less accurate, ! but much faster than looping through all npar_loc particle and add the ! gravity of all massive particles to it. ! ! 07-sep-06/wlad: coded ! use Diagnostics use Sub ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp type (pencil_case) :: p integer, dimension (mpar_loc,3) :: ineargrid ! real, dimension (nx,mspar) :: rp_mn, rpcyl_mn real, dimension (mx,3) :: ggt real, dimension (nx) :: pot_energy real, dimension (3) :: xxpar,accg integer :: ks, k logical :: lintegrate, lparticle_out ! intent (in) :: f, p, fp, ineargrid intent (inout) :: df, dfp ! ! Get the total gravity field. In the case of dust, it is already ! pre-calculated ! if (linterpolate_gravity) then ! ! Interpolate the gravity to the position of the particles. ! if (npar_imn(imn)/=0) then do k=k1_imn(imn),k2_imn(imn) !loop through the pencil if (ipar(k)>mspar) then !for dust !interpolate the gravity if (linterpolate_linear) then call interpolate_linear(f,iglobal_ggp,& iglobal_ggp+2,fp(k,ixp:izp),accg,ineargrid(k,:),0,ipar(k)) else if (linterpolate_quadratic_spline) then ! ! This interpolation works also for cylindrical coordinates. ! call interpolate_quadratic_spline(f,iglobal_ggp,& iglobal_ggp+2,fp(k,ixp:izp),accg,ineargrid(k,:),0,ipar(k)) endif dfp(k,ivpx:ivpz)=dfp(k,ivpx:ivpz)+accg endif enddo endif endif ! ! Add the acceleration to the gas. ! if (lhydro) then if (linterpolate_gravity) then !already calculated, so no need to lose time !calculating again ggt=f(:,m,n,iglobal_ggp:iglobal_ggp+2) else call get_total_gravity(ggt) endif df(l1:l2,m,n,iux:iuz) = df(l1:l2,m,n,iux:iuz) + ggt(l1:l2,:) ! ! Backreaction of the gas+dust gravity onto the massive particles ! The integration is needed for these two cases: ! ! 1. We are not solving the Poisson equation (lbackreaction) ! 2. We are, but a particle is out of the box (a star, for instance) ! and therefore the potential cannot be interpolated. ! do ks=1,mspar lparticle_out=.false. if (lselfgravity) then if ((fsp(ks,ixp)< xyz0(1)).or.(fsp(ks,ixp) > xyz1(1)) .or. & (fsp(ks,iyp)< xyz0(2)).or.(fsp(ks,iyp) > xyz1(2)) .or. & (fsp(ks,izp)< xyz0(3)).or.(fsp(ks,izp) > xyz1(3))) then !particle out of box lparticle_out=.true. endif endif lintegrate=lbackreaction.or.lparticle_out ! ! Sometimes making the star feel the selfgravity of the disk leads to ! numerical troubles as the star is too close to the origin (in cylindrical ! coordinates). ! if ((ks==istar).and.lnoselfgrav_star) lintegrate=.false. ! if (lintegrate) then ! ! Get the acceleration particle ks suffers due to self-gravity. ! call get_radial_distance(rp_mn(:,ks),rpcyl_mn(:,ks),& E1_=fsp(ks,ixp),E2_=fsp(ks,iyp),E3_=fsp(ks,izp)) xxpar = fsp(ks,ixp:izp) ! if (lcylindrical_gravity_nbody(ks)) then call integrate_selfgravity(p,rpcyl_mn(:,ks),& xxpar,accg,r_smooth(ks)) else call integrate_selfgravity(p,rp_mn(:,ks),& xxpar,accg,r_smooth(ks)) endif ! ! Add it to its dfp ! do k=1,npar_loc if (ipar(k)==ks) & dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) + accg(1:3) enddo endif ! ! Diagnostic ! if (ldiagnos) then if (idiag_totenergy/=0) pot_energy=0.0 call get_radial_distance(rp_mn(:,ks),rpcyl_mn(:,ks),& E1_=fsp(ks,ixp),E2_=fsp(ks,iyp),E3_=fsp(ks,izp)) ! ! Calculate torques for output, if needed ! if ((idiag_torqext(ks)/=0).or.(idiag_torqint(ks)/=0)) & call calc_torque(p,rpcyl_mn(:,ks),ks) ! ! Total energy ! if (idiag_totenergy/=0) then !potential energy pot_energy = pot_energy - & GNewton*pmass(ks)*& (rpcyl_mn(:,ks)**2+r_smooth(ks)**2)**(-0.5) if (ks==mspar) & call sum_lim_mn_name(.5*p%rho*p%u2 + pot_energy,& idiag_totenergy,p) endif endif enddo endif !if hydro ! endsubroutine dvvp_dt_nbody_pencil !*********************************************************************** subroutine dxxp_dt_nbody(dfp) ! ! If the center of mass of the nbody particles was moved from the ! center of the grid, reset it. ! ! 22-sep-06/wlad: coded ! real, dimension (mpar_loc,mpvar) :: dfp ! if (lreset_cm) call reset_center_of_mass(dfp) ! endsubroutine dxxp_dt_nbody !*********************************************************************** subroutine dvvp_dt_nbody(f,df,fp,dfp,ineargrid) ! ! Evolution of nbody and dust particles velocities due to particle-particle ! interaction only. ! ! Coriolis and shear are already added in particles_dust. ! ! 27-aug-06/wlad: coded ! use Mpicomm, only: mpibcast_real use Sub ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp real, dimension (mpar_loc) :: x1,y1,z1,dists2 real, dimension (mpar_loc,3) :: evr_cart integer, dimension (mpar_loc,3) :: ineargrid ! real, dimension(mspar) :: sq_hills real :: rr, w2, sma2 real :: x2,y2,z2 integer :: k, ks, j, jpos, jvel logical :: lheader, lfirstcall=.true. ! intent (in) :: f, ineargrid intent (inout) :: fp, df, dfp ! ! Print out header information in first time step. ! lheader=lfirstcall .and. lroot ! ! Identify module. ! if (lheader) print*, 'dvvp_dt_nbody: Calculate dvvp_dt_nbody' ! ! Evolve massive particle positions due to the gravity of other massive ! particles. The gravity of the massive particles on the dust will be added ! inside the pencil in dvvp_dt_nbody_pencil. ! if (lramp) call get_ramped_mass ! ! Calculate Hills radius if laccretion is switched on. ! do ks=1,mspar if (laccretion(ks).and.(ks/=istar)) then if (lcartesian_coords) then rr=sqrt(fsp(ks,ixp)**2 + fsp(ks,iyp)**2 + fsp(ks,izp)**2) elseif (lcylindrical_coords) then rr=fsp(ks,ixp) if (nzgrid/=1) rr=sqrt(fsp(ks,ixp)**2+fsp(ks,izp)**2) elseif (lspherical_coords) then rr=fsp(ks,ixp) else call fatal_error('dvvp_dt_nbody','wrong coord system') rr=0.0 endif !particle velocities are non-coordinate (linear) w2 = fsp(ks,ivpx)**2 + fsp(ks,ivpy)**2 + fsp(ks,ivpz)**2 !squared semi major axis - assumes GM=1, so beware... sma2 = (rr/(2-rr*w2))**2 !squared hills radius sq_hills(ks)=sma2*(pmass(ks)*pmass1(istar)/3)**(2./3.) else sq_hills(ks)=0.0 endif enddo ! ! Add the gravity from all N-body particles. ! if (lvectorgravity) then if (lcartesian_coords) then x1(1:npar_loc) = fp(1:npar_loc,ixp) y1(1:npar_loc) = fp(1:npar_loc,iyp) z1(1:npar_loc) = fp(1:npar_loc,izp) elseif (lcylindrical_coords) then x1(1:npar_loc) = fp(1:npar_loc,ixp)*cos(fp(1:npar_loc,iyp)) y1(1:npar_loc) = fp(1:npar_loc,ixp)*sin(fp(1:npar_loc,iyp)) z1(1:npar_loc) = fp(1:npar_loc,izp) elseif (lspherical_coords) then x1(1:npar_loc) = fp(1:npar_loc,ixp)*sin(fp(1:npar_loc,iyp))*cos(fp(1:npar_loc,izp)) y1(1:npar_loc) = fp(1:npar_loc,ixp)*sin(fp(1:npar_loc,iyp))*sin(fp(1:npar_loc,izp)) z1(1:npar_loc) = fp(1:npar_loc,ixp)*cos(fp(1:npar_loc,iyp)) endif do ks=1,mspar if (lcartesian_coords) then x2 = fsp(ks,ixp) y2 = fsp(ks,iyp) z2 = fsp(ks,izp) elseif (lcylindrical_coords) then x2 = fsp(ks,ixp)*cos(fsp(ks,iyp)) y2 = fsp(ks,ixp)*sin(fsp(ks,iyp)) z2 = fsp(ks,izp) elseif (lspherical_coords) then x2 = fsp(ks,ixp)*sin(fsp(ks,iyp))*cos(fsp(ks,izp)) y2 = fsp(ks,ixp)*sin(fsp(ks,iyp))*sin(fsp(ks,izp)) z2 = fsp(ks,ixp)*cos(fsp(ks,iyp)) endif dists2(1:npar_loc) = sqrt((x1(1:npar_loc)-x2)**2.0+(y1(1:npar_loc)-y2)**2.0+(z1(1:npar_loc)-z2)**2.0+r_smooth(ks)**2.0) evr_cart(1:npar_loc,1) = x1(1:npar_loc)-x2 evr_cart(1:npar_loc,2) = y1(1:npar_loc)-y2 evr_cart(1:npar_loc,3) = z1(1:npar_loc)-z2 if (lcartesian_coords) then where(ipar(1:npar_loc) /= ks) dfp(1:npar_loc,ivpx) = & dfp(1:npar_loc,ivpy) - GNewton*pmass(ks)*dists2(1:npar_loc)**(-1.5)*evr_cart(1:npar_loc,1) where(ipar(1:npar_loc) /= ks) dfp(1:npar_loc,ivpy) = & dfp(1:npar_loc,ivpy) - GNewton*pmass(ks)*dists2(1:npar_loc)**(-1.5)*evr_cart(1:npar_loc,2) where(ipar(1:npar_loc) /= ks) dfp(1:npar_loc,ivpy) = & dfp(1:npar_loc,ivpy) - GNewton*pmass(ks)*dists2(1:npar_loc)**(-1.5)*evr_cart(1:npar_loc,3) else where(ipar(1:npar_loc) /= ks) dfp(1:npar_loc,ivpx_cart) = & dfp(1:npar_loc,ivpx_cart) - GNewton*pmass(ks)*dists2(1:npar_loc)**(-1.5)*evr_cart(1:npar_loc,1) where(ipar(1:npar_loc) /= ks) dfp(1:npar_loc,ivpy_cart) = & dfp(1:npar_loc,ivpy_cart) - GNewton*pmass(ks)*dists2(1:npar_loc)**(-1.5)*evr_cart(1:npar_loc,2) where(ipar(1:npar_loc) /= ks) dfp(1:npar_loc,ivpz_cart) = & dfp(1:npar_loc,ivpz_cart) - GNewton*pmass(ks)*dists2(1:npar_loc)**(-1.5)*evr_cart(1:npar_loc,3) endif enddo else do k=npar_loc,1,-1 ! if (linterpolate_gravity) then !only loop through massive particles if (ipar(k)<=mspar) then call loop_through_nbodies(fp,dfp,k,sq_hills,ineargrid) endif else !for all particles call loop_through_nbodies(fp,dfp,k,sq_hills,ineargrid) endif enddo endif ! ! Position and velocity diagnostics (per nbody particle). ! if (ldiagnos) then do ks=1,mspar if (lfollow_particle(ks)) then do j=1,3 jpos=j+ixp-1 ; jvel=j+ivpx-1 if (idiag_xxspar(ks,j)/=0) then call point_par_name(fsp(ks,jpos),idiag_xxspar(ks,j)) endif if (idiag_vvspar(ks,j)/=0) & call point_par_name(fsp(ks,jvel),idiag_vvspar(ks,j)) enddo endif enddo endif ! if (lheader) print*,'dvvp_dt_nbody: Finished dvvp_dt_nbody' if (lfirstcall) lfirstcall=.false. ! call keep_compiler_quiet(f,df) ! endsubroutine dvvp_dt_nbody !************************************************************ subroutine loop_through_nbodies(fp,dfp,k,sq_hills,ineargrid) ! real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp real, dimension (mspar) :: sq_hills integer, dimension (mpar_loc,3) :: ineargrid integer :: k ! if (linertial_frame) then call loop_through_nbodies_inertial(fp,dfp,k,sq_hills,ineargrid) else call loop_through_nbodies_fixstar(fp,dfp,k,sq_hills,ineargrid) endif ! endsubroutine loop_through_nbodies !************************************************************ subroutine loop_through_nbodies_inertial(fp,dfp,k,sq_hills,ineargrid) ! ! Subroutine that adds the gravity from all massive particles. ! Particle gravity has always to be added in Cartesian, for better ! conservation of the Jacobi constant. ! ! 07-mar-08/wlad:coded ! real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp integer :: k real, dimension (mspar) :: sq_hills integer, dimension (mpar_loc,3) :: ineargrid ! real :: r2_ij, rs2, invr3_ij, v_ij,tmp1,tmp2 integer :: ks ! real, dimension (3) :: evr_cart,acc_cart ! intent(inout) :: fp, dfp intent(in) :: k ! do ks=1,mspar if (ipar(k)/=ks) then ! call get_evr(fp(k,ixp:izp),fsp(ks,ixp:izp),evr_cart) ! ! Particles relative distance from each other: ! ! r_ij = sqrt(ev1**2 + ev2**2 + ev3**2) ! tmp1=sum(evr_cart**2) tmp2=r_smooth(ks)**2 r2_ij=max(tmp1,tmp2) ! ! If there is accretion, remove the accreted particles from the simulation, if any. ! if (laccretion(ks)) then rs2=(accrete_hills_frac(ks)**2)*sq_hills(ks) if (r2_ij<=rs2) then call remove_particle(fp,ipar,k,dfp,ineargrid,ks) !add mass of the removed particle to the accreting particle if (ladd_mass(ks)) pmass(ks)=pmass(ks)+mp_swarm return endif endif ! ! Shortcut: invr3_ij = r_ij**(-3) ! if (r2_ij > 0) then invr3_ij = r2_ij**(-1.5) else ! can happen during pencil_check invr3_ij = 0.0 endif ! ! Gravitational acceleration: g=-g0/|r-r0|^3 (r-r0) ! ! The acceleration is in non-coordinate basis (all have dimension of length). ! The main dxx_dt of particle_dust takes care of transforming the linear ! velocities to angular changes in position. ! acc_cart = - GNewton*pmass(ks)*invr3_ij*evr_cart(1:3) if (lcartesian_coords) then dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) + acc_cart else ! separate this N-body acceleration from other, added elsewhere in the code dfp(k,ivpx_cart:ivpz_cart) = dfp(k,ivpx_cart:ivpz_cart) + acc_cart endif ! ! Time-step constraint from N-body particles. We use both the criterion ! that the distance to the N-body particle must not change too much in ! one time-step and additionally we use the free-fall time-scale. ! if (lfirst.and.ldt.and.ldt_nbody) then v_ij=sqrt(sum((fp(k,ivpx:ivpz)-fp(ks,ivpx:ivpz))**2)) dt1_max(ineargrid(k,1)-nghost)= & max(dt1_max(ineargrid(k,1)-nghost),v_ij/sqrt(r2_ij)/cdtpnbody) dt1_max(ineargrid(k,1)-nghost)= & max(dt1_max(ineargrid(k,1)-nghost), & sqrt(GNewton*pmass(ks)*invr3_ij)/cdtpnbody) endif ! endif !if (ipar(k)/=ks) ! enddo !nbody loop ! endsubroutine loop_through_nbodies_inertial !********************************************************** subroutine loop_through_nbodies_fixstar(fp,dfp,k,sq_hills,ineargrid) ! ! Gravity acceleration for all bodies, in the reference frame of the star. ! So far, works only for two massive bodies, in Cartesian coordinates. ! ! 23-jun-14/wlad: coded ! real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp integer :: k real, dimension (mspar) :: sq_hills integer, dimension (mpar_loc,3) :: ineargrid ! real :: r2_ij, invr3_ij integer :: ks ! real, dimension (3) :: evr_cart,acc_cart ! intent(inout) :: fp, dfp intent(in) :: k ! if (ipar(k) == istar) then dfp(k,ivpx:ivpz) = 0. else acc_cart=0. do ks=1,mspar if (ipar(k)/=ks) then ! evr_cart(1)=fp(k,ixp); evr_cart(2)=fp(k,iyp); evr_cart(3)=fp(k,izp) ! r2_ij=sum(evr_cart**2) invr3_ij = r2_ij**(-1.5) ! acc_cart = - GNewton*invr3_ij*evr_cart(1:3) ! dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) + acc_cart endif enddo endif ! endsubroutine loop_through_nbodies_fixstar !********************************************************** subroutine get_evr(xxp,xxsp,evr_cart) ! ! Point-to-point vector distance, in different coordinate systems. ! Return always in Cartesian. ! ! 14-feb-14/wlad: coded ! real, dimension(3), intent(in) :: xxp,xxsp real, dimension(3), intent(out) :: evr_cart real :: x1,y1,x2,y2,z1,z2 ! if (lcartesian_coords) then x1=xxp(1) ; x2=xxsp(1) y1=xxp(2) ; y2=xxsp(2) z1=xxp(3) ; z2=xxsp(3) elseif (lcylindrical_coords) then x1=xxp(1)*cos(xxp(2)) y1=xxp(1)*sin(xxp(2)) z1=xxp(3) ! x2=xxsp(1)*cos(xxsp(2)) y2=xxsp(1)*sin(xxsp(2)) z2=xxsp(3) elseif (lspherical_coords) then x1=xxp(1)*sin(xxp(2))*cos(xxp(3)) y1=xxp(1)*sin(xxp(2))*sin(xxp(3)) z1=xxp(1)*cos(xxp(2)) ! x2=xxsp(1)*sin(xxsp(2))*cos(xxsp(3)) y2=xxsp(1)*sin(xxsp(2))*sin(xxsp(3)) z2=xxsp(1)*cos(xxsp(2)) endif ! evr_cart(1)=x1-x2 evr_cart(2)=y1-y2 evr_cart(3)=z1-z2 ! endsubroutine get_evr !********************************************************** subroutine point_par_name(a,iname) ! ! Register a, a simple scalar, as diagnostic. ! ! Works for individual particle diagnostics. ! ! 07-mar-08/wlad: adapted from sum_par_name ! real :: a integer :: iname ! if (iname/=0) then fname(iname)=a endif ! ! There is no entry in itype_name. ! endsubroutine point_par_name !*********************************************************************** subroutine read_particles_nbody_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=particles_nbody_init_pars, IOSTAT=iostat) ! endsubroutine read_particles_nbody_init_pars !*********************************************************************** subroutine write_particles_nbody_init_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=particles_nbody_init_pars) ! endsubroutine write_particles_nbody_init_pars !*********************************************************************** subroutine read_particles_nbody_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=particles_nbody_run_pars, IOSTAT=iostat) ! endsubroutine read_particles_nbody_run_pars !*********************************************************************** subroutine write_particles_nbody_run_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=particles_nbody_run_pars) ! endsubroutine write_particles_nbody_run_pars !*********************************************************************** subroutine reset_center_of_mass(dfp) ! ! If the center of mass was accelerated, reset its position ! to the center of the grid. ! ! Assumes that the total mass of the particles is one. ! ! 27-aug-06/wlad: coded ! 18-mar-08/wlad: cylindrical and spherical corrections ! real, dimension(mpar_loc,mpvar),intent(inout) :: dfp real, dimension(mspar,mspvar) :: ftmp real, dimension(3) :: vcm real :: xcm,ycm,zcm,thtcm,phicm,vxcm,vycm,vzcm integer :: k ! ftmp=fsp(1:mspar,:) ! if (lcartesian_coords) then vcm(1) = sum(ftmp(:,imass_nbody)*ftmp(:,ivpx)) vcm(2) = sum(ftmp(:,imass_nbody)*ftmp(:,ivpy)) vcm(3) = sum(ftmp(:,imass_nbody)*ftmp(:,ivpz)) else if (lcylindrical_coords) then xcm=sum(ftmp(:,imass_nbody)*(ftmp(:,ixp)*cos(ftmp(:,iyp)))) ycm=sum(ftmp(:,imass_nbody)*(ftmp(:,ixp)*sin(ftmp(:,iyp)))) phicm=atan2(ycm,xcm) ! vxcm=sum(ftmp(:,imass_nbody)*(& ftmp(:,ivpx)*cos(ftmp(:,iyp))-ftmp(:,ivpy)*sin(ftmp(:,iyp)))) vycm=sum(ftmp(:,imass_nbody)*(& ftmp(:,ivpx)*sin(ftmp(:,iyp))+ftmp(:,ivpy)*cos(ftmp(:,iyp)))) ! vcm(1)= vxcm*cos(phicm) + vycm*sin(phicm) vcm(2)=-vxcm*sin(phicm) + vycm*cos(phicm) vcm(3) = sum(ftmp(:,imass_nbody)*ftmp(:,ivpz)) ! else if (lspherical_coords) then vxcm=sum(ftmp(:,imass_nbody)*( & ftmp(:,ivpx)*sin(ftmp(:,iyp))*cos(ftmp(:,izp))& +ftmp(:,ivpy)*cos(ftmp(:,iyp))*cos(ftmp(:,izp))& -ftmp(:,ivpz)*sin(ftmp(:,izp)) )) vycm=sum(ftmp(:,imass_nbody)*( & ftmp(:,ivpx)*sin(ftmp(:,iyp))*sin(ftmp(:,izp))& +ftmp(:,ivpy)*cos(ftmp(:,iyp))*sin(ftmp(:,izp))& +ftmp(:,ivpz)*cos(ftmp(:,izp)) )) vzcm=sum(ftmp(:,imass_nbody)*(& ftmp(:,ivpx)*cos(ftmp(:,iyp))-ftmp(:,ivpy)*sin(ftmp(:,iyp)))) ! xcm=sum(ftmp(:,imass_nbody)*(ftmp(:,ixp)*sin(ftmp(:,iyp))*cos(ftmp(:,izp)))) ycm=sum(ftmp(:,imass_nbody)*(ftmp(:,ixp)*sin(ftmp(:,iyp))*sin(ftmp(:,izp)))) zcm=sum(ftmp(:,imass_nbody)*(ftmp(:,ixp)*cos(ftmp(:,iyp)))) ! thtcm=atan2(sqrt(xcm**2+ycm**2),zcm) phicm=atan2(ycm,xcm) ! vcm(1)= vxcm*sin(thtcm)*cos(phicm) + & vycm*sin(thtcm)*sin(phicm) + vzcm*cos(thtcm) vcm(2)= vxcm*cos(thtcm)*cos(phicm) + & vycm*cos(thtcm)*sin(phicm) - vzcm*sin(thtcm) vcm(3)=-vxcm*sin(phicm) + & vycm*cos(phicm) endif ! do k=1,npar_loc if (ipar(k)<=mspar) then dfp(k,ixp:izp) = dfp(k,ixp:izp) - vcm*totmass1 endif enddo ! endsubroutine reset_center_of_mass !*********************************************************************** subroutine integrate_selfgravity(p,rrp,xxpar,accg,rp0) ! ! Calculates acceleration on the point (x,y,z)=xxpar ! due to the gravity of the gas+dust. ! ! 15-sep-06/wlad : coded ! use Mpicomm ! real, dimension(nx,3) :: dist real, dimension(nx) :: rrp,selfgrav,density real, dimension(nx) :: dv,jac,dqy,tmp real :: dqx,dqz,rp0,fac real, dimension(3) :: xxpar,accg,sum_loc integer :: j type (pencil_case) :: p logical :: lfirstcall=.true. ! intent(out) :: accg ! if (lfirstcall) & print*,'Adding gas+dust gravity to the massive particles' ! ! Sanity check ! if (.not.(lgas_gravity.or.ldust_gravity)) & call fatal_error("lintegrate_selfgravity",& "No gas gravity or dust gravity to add. "//& "Switch on lgas_gravity or ldust_gravity in n-body parameters") ! if (coord_system=='cartesian') then jac=1.;dqx=dx;dqy=dy;dqz=dz dist(:,1)=x(l1:l2)-xxpar(1) dist(:,2)=y( m )-xxpar(2) dist(:,3)=z( n )-xxpar(3) elseif (coord_system=='cylindric') then jac=x(l1:l2);dqx=dx;dqy=x(l1:l2)*dy;dqz=dz dist(:,1)=x(l1:l2)-xxpar(1)*cos(y(m)-xxpar(2)) dist(:,2)= xxpar(2)*sin(y(m)-xxpar(2)) dist(:,3)=z( n )-xxpar(3) elseif (coord_system=='spherical') then call fatal_error('integrate_selfgravity', & ' not yet implemented for spherical polars') dqx=0.;dqy=0.;dqz=0. else call fatal_error('integrate_selfgravity','wrong coord_system') dqx=0.;dqy=0.;dqz=0. endif ! if (nzgrid==1) then dv=dqx*dqy else dv=dqx*dqy*dqz endif ! ! The gravity of every single cell - should exclude inner and outer radii... ! ! selfgrav = G*((rho+rhop)*dv)*mass*r*(r**2 + r0**2)**(-1.5) ! gx = selfgrav * r\hat dot x\hat ! -> gx = selfgrav * (x-x0)/r = G*((rho+rhop)*dv)*mass*(r**2+r0**2)**(-1.5) * (x-x0) ! density=0. if (lgas_gravity) density=p%rho ! ! Add the particle gravity if npar>mspar (which means dust is being used) ! if (ldust.and.ldust_gravity) density=density+p%rhop ! selfgrav = GNewton*density_scale*& density*jac*dv*(rrp**2 + rp0**2)**(-1.5) ! ! Everything inside the accretion radius of the particle should ! not exert gravity (numerical problems otherwise) ! where (rrp<=rp0) selfgrav = 0 endwhere ! ! Exclude the frozen zones ! if (lexclude_frozen) then if (lcylinder_in_a_box) then where ((p%rcyl_mn<=r_int).or.(p%rcyl_mn>=r_ext)) selfgrav = 0 endwhere else where ((p%r_mn<=r_int).or.(p%r_mn>=r_ext)) selfgrav = 0 endwhere endif endif ! ! Integrate the accelerations on this processor ! And sum over processors with mpireduce ! do j=1,3 tmp=selfgrav*dist(:,j) !take proper care of the trapezoidal rule !in the case of non-periodic boundaries fac = 1. if ((m==m1.and.lfirst_proc_y).or.(m==m2.and.llast_proc_y)) then if (.not.lperi(2)) fac = .5*fac endif ! if (lperi(1)) then sum_loc(j) = fac*sum(tmp) else sum_loc(j) = fac*(sum(tmp(2:nx-1))+.5*(tmp(1)+tmp(nx))) endif call mpireduce_sum(sum_loc(j),accg(j)) enddo ! ! Broadcast particle acceleration ! call mpibcast_real(accg,3) ! if (lfirstcall) lfirstcall=.false. ! endsubroutine integrate_selfgravity !*********************************************************************** subroutine bcast_nbodyarray(fp) ! ! Broadcast nbody particles across processors ! The non-root processors, if they find a nbody particle in their ! fp array, they: ! send it to root with a true logical ! else send a false logical ! ! The root processor receives all the logicals, and if ! they are true, then receives the value, broadcasting it ! ! 07-sep-06/wlad: coded ! use Mpicomm ! real, dimension (mpar_loc,mparray) :: fp logical, dimension(mspar) :: lnbody integer :: ks,k,tagsend,j,tagrecv ! if (lmpicomm) then ! ! Loop through the nbody particles ! do ks=1,mspar ! ! Set the logical to false initially ! lnbody(ks)=.false. ! ! Loop through the particles on this processor ! do k=1,npar_loc if (ipar_nbody(ks)==ipar(k)) then ! ! A nbody was found here. Turn the logical true and copy fp to fsp ! lnbody(ks) = .true. fsp(ks,ixp:izp) = fp(k,ixp:izp) fsp(ks,ivpx:ivpz) = fp(k,ivpx:ivpz) fsp(ks,imass_nbody) = pmass(ks) ! ! Send it to root. As there will be just mspar calls to ! mpisend, the tag can be ipar itself ! if (.not.lroot) then call mpisend_real(fsp(ks,:),mspvar,root,ks) if (ip<=6) print*,'logical for particle ',ks,& ' set to true on processor ',iproc, & ' with tag=',ks endif endif enddo ! ! Send the logicals from each processor. As all processors send mspar calls, ! the tags are now in base mspar. It assures that two logicals will not have ! the same tag. ! if (.not.lroot) then tagsend = mspar*iproc + ks call mpisend_logical(lnbody(ks),root,tagsend) else ! ! The root receives all logicals. Same tag. ! do j=1,ncpus-1 tagrecv = mspar*j + ks call mpirecv_logical(lnbody(ks),j,tagrecv) ! ! Test the received logicals ! if (lnbody(ks)) then ! ! Found a nbody particle. Get the value of fsp ! call mpirecv_real(fsp(ks,:),mspvar,j,ks) if (ip<=6) print*,'logical for particle ',ks,& ' is true on processor ',j, & ' with tag=',ks,' on root' endif enddo endif ! ! Broadcast the received fsp ! call mpibcast_real(fsp(ks,:),mspvar) ! ! Print the result in all processors ! if (ip<=8) print*,'bcast_nbodyarray: finished loop. '//& 'nbody particles in proc ',iproc,& ' are fsp(ks,:)=',fsp(ks,:) enddo else ! ! Non-parallel. Just copy fp to fsp ! do ks=1,mspar do k=1,npar_loc if (ipar_nbody(ks)==ipar(k)) then fsp(ks,ixp:izp) = fp(k,ixp:izp) fsp(ks,ivpx:ivpz) = fp(k,ivpx:ivpz) endif enddo fsp(ks,imass_nbody) = pmass(ks) enddo ! endif ! if (ldebug) print*,'bcast_nbodyarray finished' ! endsubroutine bcast_nbodyarray !*********************************************************************** subroutine particles_nbody_special ! ! Fetch fsp array to special module ! ! 01-mar-08/wlad: coded ! use Special, only: special_calc_particles_nbody ! call special_calc_particles_nbody(fsp) ! endsubroutine particles_nbody_special !*********************************************************************** subroutine get_totalmass(tmass) ! ! called from global_shear to set the initial keplerian field ! real :: tmass ! if (lnorm) then tmass=1. else tmass=sum(pmass) endif ! endsubroutine get_totalmass !*********************************************************************** subroutine get_gravity_field_nbody(grr,gg,ks) ! ! Subroutine that returns the gravity field a particle ! exterts in a pencil, respecting the coordinate system used ! ! 07-mar-08/wlad: coded ! real, dimension(mx),intent(in) :: grr real, dimension(mx,3),intent(out) :: gg real :: x0,y0,z0 integer :: ks ! x0=fsp(ks,ixp);y0=fsp(ks,iyp);z0=fsp(ks,izp) ! if (coord_system=='cartesian') then gg(:,1) = (x -x0)*grr gg(:,2) = (y(m)-y0)*grr gg(:,3) = (z(n)-z0)*grr elseif (coord_system=='cylindric') then gg(:,1) = (x -x0*cos(y(m)-y0))*grr gg(:,2) = x0*sin(y(m)-y0) *grr gg(:,3) = (z(n)-z0 )*grr elseif (coord_system=='spherical') then gg(:,1) = (x-x0*sinth(m)*sin(y0)*cos(z(n)-z0))*grr gg(:,2) = (x0*(sinth(m)*cos(y0)-& costh(m)*sin(y0)*cos(z(n)-z0)))*grr gg(:,3) = (x0*sin(y0)*sin(z(n)-z0))*grr endif ! endsubroutine get_gravity_field_nbody !*********************************************************************** subroutine calc_torque(p,dist,ks) ! ! Output torque diagnostic for nbody particle ks ! ! 05-nov-05/wlad : coded ! use Diagnostics ! type (pencil_case) :: p real, dimension(nx) :: torque_gas,torqint_gas,torqext_gas real, dimension(nx) :: torque_par,torqint_par,torqext_par real, dimension(nx) :: dist,rpre,rr_mn,tempering real :: rr,w2,smap,hills,phip,phi,pcut integer :: ks,i ! if (ks==istar) call fatal_error('calc_torque', & 'Nonsense to calculate torques for the star') ! if (lcartesian_coords) then rr = sqrt(fsp(ks,ixp)**2 + fsp(ks,iyp)**2 + fsp(ks,izp)**2) rpre = fsp(ks,ixp)*y(m) - fsp(ks,iyp)*x(l1:l2) elseif (lcylindrical_coords) then rr_mn = x(l1:l2) ; phi = y(m) rr = fsp(ks,ixp) ; phip = fsp(ks,iyp) rpre = rr_mn*rr*sin(phi-phip) elseif (lspherical_coords) then call fatal_error("calc_torque",& "not yet implemented for spherical coordinates") else call fatal_error("calc_torque",& "the world is flat and we should never gotten here") endif ! w2 = fsp(ks,ivpx)**2 + fsp(ks,ivpy)**2 + fsp(ks,ivpz)**2 smap = 1./(2./rr - w2) hills = smap*(pmass(ks)*pmass1(istar)/3.)**(1./3.) ! ! Define separate torques for gas and dust/particles ! torque_gas = GNewton*pmass(ks)*p%rho*rpre*& (dist**2 + r_smooth(ks)**2)**(-1.5) ! if (ldust) then torque_par = GNewton*pmass(ks)*p%rhop*rpre*& (dist**2 + r_smooth(ks)**2)**(-1.5) else torque_par=0. endif ! if (ldustdensity) & call fatal_error("calc_torque",& "not implemented for the dust fluid approximation") ! ! Zero torque outside r_int and r_ext in Cartesian coordinates ! if (lcartesian_coords) then do i=1,nx if (p%rcyl_mn(i) > r_ext .or. p%rcyl_mn(i) < r_int) then torque_gas(i)=0. torque_par(i)=0. endif enddo endif ! ! Exclude region inside a fraction (hills_tempering_fraction) of the Hill sphere. ! if (ltempering) then pcut=hills_tempering_fraction*hills tempering = 1./(exp(-(sqrt(dist**2)/hills - pcut)/(.1*pcut))+1.) torque_gas = torque_gas * tempering torque_par = torque_par * tempering else do i=1,nx if (dist(i)=rr) then torqext_gas(i) = torque_gas(i) torqext_par(i) = torque_par(i) else torqext_gas(i)=0.;torqext_par(i)=0. endif if (p%rcyl_mn(i)<=rr) then torqint_gas(i) = torque_gas(i) torqint_par(i) = torque_par(i) else torqint_gas(i)=0.;torqint_par(i)=0. endif enddo ! ! Sum the different torque contributions. ! if (lcartesian_coords) then call sum_lim_mn_name(torqext_gas,idiag_torqext_gas(ks),p) call sum_lim_mn_name(torqext_par,idiag_torqext_par(ks),p) call sum_lim_mn_name(torqint_gas,idiag_torqint_gas(ks),p) call sum_lim_mn_name(torqint_par,idiag_torqint_par(ks),p) ! ! Backward compatibility ! call sum_lim_mn_name(torqext_gas+torqext_par,idiag_torqext(ks),p) call sum_lim_mn_name(torqint_gas+torqint_par,idiag_torqint(ks),p) else ! ! Hack for non-cartesian coordinates. sum_lim_mn_name is lagging ! behind sum_mn_name, and whould be brought up to date. ! call integrate_mn_name(torqext_gas,idiag_torqext_gas(ks)) call integrate_mn_name(torqext_par,idiag_torqext_par(ks)) call integrate_mn_name(torqint_gas,idiag_torqint_gas(ks)) call integrate_mn_name(torqint_par,idiag_torqint_par(ks)) call integrate_mn_name(torqext_gas+torqext_par,idiag_torqext(ks)) call integrate_mn_name(torqint_gas+torqint_par,idiag_torqint(ks)) endif ! endsubroutine calc_torque !*********************************************************************** subroutine get_ramped_mass ! real :: ramping_period,tmp integer :: ks ! ramping_period=2*pi*ramp_orbits if (t <= ramping_period) then !sin ((pi/2)*(t/(ramp_orbits*2*pi)) tmp=0. !just need to do that for the original masses do ks=1,mspar_orig if (ks/=istar) then pmass(ks)= max(dble(tini),& final_ramped_mass(ks)*(sin((.5*pi)*(t/ramping_period))**2)) tmp=tmp+pmass(ks) endif enddo pmass(istar)= 1-tmp else pmass(1:mspar_orig)=final_ramped_mass(1:mspar_orig) endif ! endsubroutine get_ramped_mass !*********************************************************************** subroutine calc_nbodygravity_particles(f) ! ! For the case of interpolated gravity, add the gravity ! of all n-body particles to slots of the f-array. ! ! 08-mar-08/wlad: coded ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(mx,3) :: ggt ! if (linterpolate_gravity) then ! if (lramp) call get_ramped_mass ! ! Calculate grid - nbody particles distances ! do n=1,mz do m=1,my call get_total_gravity(ggt) f(:,m,n,iglobal_ggp:iglobal_ggp+2)=ggt enddo enddo ! ! else do nothing ! endif ! endsubroutine calc_nbodygravity_particles !*********************************************************************** subroutine get_total_gravity(ggt) ! ! Sum the gravities of all massive particles ! ! 08-mar-08/wlad: coded ! use Sub ! real, dimension (mx,mspar) :: rp_mn,rpcyl_mn real, dimension (mx,3) :: ggp,ggt real, dimension (mx) :: grav_particle,rrp integer :: ks ! intent(out) :: ggt ! ggt=0. do ks=1,mspar ! ! Spherical and cylindrical distances ! call get_radial_distance(rp_mn(:,ks),rpcyl_mn(:,ks),& e1_=fsp(ks,ixp),e2_=fsp(ks,iyp),e3_=fsp(ks,izp)) ! ! Check which particle has cylindrical gravity switched on ! if (lcylindrical_gravity_nbody(ks)) then rrp = rpcyl_mn(:,ks) else rrp = rp_mn(:,ks) endif ! ! Gravity field from the particle ks ! grav_particle =-GNewton*pmass(ks)*(rrp**2+r_smooth(ks)**2)**(-1.5) call get_gravity_field_nbody(grav_particle,ggp,ks) ! if ((ks==istar).and.lnogravz_star) & ggp(:,3) = 0. ! ! Sum up the accelerations of the massive particles ! ggt=ggt+ggp ! ! Add indirect term if the star is fixed at the center ! if (.not.linertial_frame) call add_indirect_term(ks,ggt) ! enddo ! endsubroutine get_total_gravity !*********************************************************************** subroutine add_indirect_term(ks,ggt) ! ! Add the terms due to the reference frame being away from the baricenter. ! So far, only for one perturber (two massive bodies), and in Cartesian coordinates. ! ! 23-jun-14/wlad: coded ! real, dimension(mx,3) :: ggt real :: rr1_3 integer, intent(in) :: ks ! if (ks/=istar) then rr1_3=(fsp(ks,ixp)**2 + fsp(ks,iyp)**2 + fsp(ks,izp)**2)**(-1.5) ggt(:,1) = ggt(:,1) - GNewton*pmass(ks)*fsp(ks,ixp)*rr1_3 ggt(:,2) = ggt(:,2) - GNewton*pmass(ks)*fsp(ks,iyp)*rr1_3 ggt(:,3) = ggt(:,3) - GNewton*pmass(ks)*fsp(ks,izp)*rr1_3 endif ! endsubroutine add_indirect_term !*********************************************************************** subroutine advance_particles_in_cartesian(fp,dfp) ! ! With N-body gravity, the particles should have their position advanced in ! Cartesian coordinates, for better conservation of the Jacobi constant, even ! if the grid is polar. ! ! 14-feb-14/wlad: coded ! real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp ! real, dimension (npar_loc) :: rad,phi,tht,xp,yp,zp real, dimension (npar_loc) :: vrad,vtht,vphi,vx,vy,vz real, dimension (npar_loc) :: arad,atht,aphi,ax,ay,az real, dimension (npar_loc) :: cosp,sinp,cost,sint ! if (lcylindrical_coords) then ! ! The input position, velocities, and accelerations in cylindrical coordinates. ! rad = fp(1:npar_loc,ixp) ; phi = fp(1:npar_loc,iyp) vrad = fp(1:npar_loc,ivpx) ; vphi = fp(1:npar_loc,ivpy) arad = dfp(1:npar_loc,ivpx) ; aphi = dfp(1:npar_loc,ivpy) ! ! Shortcuts. ! cosp=cos(phi) ; sinp=sin(phi) ! ! Transform the position, velocities, and accelerations to Cartesian coordinates. ! xp=rad*cosp ; yp=rad*sinp ; zp=fp(1:npar_loc,izp) ! vx=vrad*cosp - vphi*sinp vy=vrad*sinp + vphi*cosp vz=fp(1:npar_loc,izp) ! ! Add the Cartesian acceleration. ! ax=arad*cosp - aphi*sinp + dfp(1:npar_loc,ivpx_cart) ay=arad*sinp + aphi*cosp + dfp(1:npar_loc,ivpy_cart) az=dfp(1:npar_loc,ivpz) ! else if (lspherical_coords) then ! ! The input position, velocities, and accelerations in spherical coordinates. ! rad = fp(1:npar_loc,ixp) ; tht = fp(1:npar_loc,iyp) ; phi = fp(1:npar_loc,izp) vrad = fp(1:npar_loc,ivpx) ; vtht = fp(1:npar_loc,ivpy) ; vphi = fp(1:npar_loc,ivpz) arad = dfp(1:npar_loc,ivpx) ; atht = dfp(1:npar_loc,ivpy) ; aphi = dfp(1:npar_loc,ivpz) ! ! Shortcuts. ! cosp=cos(phi) ; sinp=sin(phi) cost=cos(tht) ; sint=sin(tht) ! ! Transform the position, velocities, and accelerations to Cartesian coordinates. ! xp=rad*sint*cosp ; yp=rad*sint*sinp ; zp=rad*cost ! vx=vrad*sint*cosp + vtht*cost*cosp - vphi*sinp vy=vrad*sint*sinp + vtht*cost*sinp + vphi*cosp vz=vrad*cost - vtht*sint ! ! Add the Cartesian acceleration. ! ax=arad*sint*cosp + atht*cost*cosp - aphi*sinp + dfp(1:npar_loc,ivpx_cart) ay=arad*sint*sinp + atht*cost*sinp + aphi*cosp + dfp(1:npar_loc,ivpy_cart) az=arad*cost - atht*sint + dfp(1:npar_loc,ivpz_cart) ! endif ! ! Now the time-stepping in Cartesian coordinates. ! call update_position(fp,dfp,xp,yp,zp,vx,vy,vz) call update_velocity(fp,vx,vy,vz,ax,ay,az) ! endsubroutine advance_particles_in_cartesian !*********************************************************************** subroutine update_position(fp,dfp,xp,yp,zp,vx,vy,vz) ! ! Update position if N-body is used in polar coordinates. ! ! 14-feb-14:wlad/coded ! real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp ! real, dimension (npar_loc), intent(in) :: vx,vy,vz real, dimension (npar_loc), intent(inout) :: xp,yp,zp ! ! Input vx and vy into the dfp array, for the Runge-Kutta integration. ! It is needed because of the high order of the RK integration (dfp is ! updated every subtimestep. ! dfp(1:npar_loc,ixp) = dfp(1:npar_loc,ixp) + vx dfp(1:npar_loc,iyp) = dfp(1:npar_loc,iyp) + vy dfp(1:npar_loc,izp) = dfp(1:npar_loc,izp) + vz ! ! Update. ! xp = xp + dt_beta_ts(itsub)*dfp(1:npar_loc,ixp) yp = yp + dt_beta_ts(itsub)*dfp(1:npar_loc,iyp) zp = zp + dt_beta_ts(itsub)*dfp(1:npar_loc,izp) ! ! Convert back to polar coordinates. ! if (lcylindrical_coords) then fp(1:npar_loc,ixp) = sqrt(xp**2+yp**2+zp**2) fp(1:npar_loc,iyp) = atan2(yp,xp) fp(1:npar_loc,izp) = zp else if (lspherical_coords) then fp(1:npar_loc,ixp) = sqrt(xp**2+yp**2+zp**2) fp(1:npar_loc,iyp) = acos(zp/fp(1:npar_loc,ixp)) fp(1:npar_loc,izp) = atan2(yp,xp) endif ! endsubroutine update_position !*********************************************************************** subroutine update_velocity(fp,vx,vy,vz,ax,ay,az) ! ! Update velocity if N-body is used in polar coordinates. ! ! 14-feb-14:wlad/coded ! real, dimension (mpar_loc,mparray) :: fp ! real, dimension (npar_loc), intent(in) :: ax,ay,az real, dimension (npar_loc), intent(inout) :: vx,vy,vz ! real, dimension (npar_loc) :: phi,tht real, dimension (npar_loc) :: cosp,sinp,sint,cost ! ! All accelerations in dfp are in Cartesian. ! vx = vx + dt_beta_ts(itsub)*ax !vxdot vy = vy + dt_beta_ts(itsub)*ay !vydot vz = vz + dt_beta_ts(itsub)*az !vydot ! ! Convert back to polar coordinates. ! if (lcylindrical_coords) then phi=fp(1:npar_loc,iyp) cosp=cos(phi) ; sinp=sin(phi) fp(1:npar_loc,ivpx) = vx*cosp + vy*sinp !=vrad fp(1:npar_loc,ivpy) = -vx*sinp + vy*cosp !=vphi fp(1:npar_loc,ivpz) = vz else if (lspherical_coords) then tht=fp(1:npar_loc,iyp) phi=fp(1:npar_loc,izp) cost=cos(tht) ; sint=sin(tht) cosp=cos(phi) ; sinp=sin(phi) fp(1:npar_loc,ivpx) = vx*sint*cosp + vy*sint*sinp + vz*cost !=vrad fp(1:npar_loc,ivpy) = vx*cost*cosp + vy*cost*sinp - vz*sint !=vphi fp(1:npar_loc,ivpz) = -vx*sinp + vy*cosp !=vz endif ! endsubroutine update_velocity !*********************************************************************** subroutine particles_nbody_read_snapshot(filename) ! ! Read nbody particle info ! ! 01-apr-08/wlad: dummy ! use Mpicomm, only:mpibcast_int,mpibcast_real ! character (len=*) :: filename ! if (lroot) then open(1,FILE=filename,FORM='unformatted') read(1) mspar if (mspar/=0) read(1) fsp(1:mspar,:),ipar_nbody(1:mspar) if (ip<=8) print*, 'read snapshot', filename close(1) endif ! call mpibcast_int(mspar) call mpibcast_real(fsp(1:mspar,:),(/mspar,mspvar/)) call mpibcast_int(ipar_nbody(1:mspar),mspar) ! if (ldebug) then print*,'particles_nbody_read_snapshot' print*,'mspar=',mspar print*,'fsp(1:mspar)=',fsp(1:mspar,:) print*,'ipar_nbody(1:mspar)=',ipar_nbody print*,'' endif ! endsubroutine particles_nbody_read_snapshot !*********************************************************************** subroutine particles_nbody_write_snapshot(snapbase,enum,flist) ! use General, only:safe_character_assign use Sub, only: update_snaptime, read_snaptime ! ! Input and output of information about the massive particles ! ! 01-apr-08/wlad: coded ! logical, save :: lfirst_call=.true. integer, save :: nsnap real, save :: tsnap character (len=*) :: snapbase,flist character (len=fnlen) :: snapname, filename_diag logical :: enum,lsnap character (len=intlen) :: nsnap_ch optional :: flist ! if (enum) then call safe_character_assign(filename_diag,trim(datadir)//'/tsnap.dat') if (lfirst_call) then call read_snaptime(filename_diag,tsnap,nsnap,dsnap,t) lfirst_call=.false. endif call update_snaptime(filename_diag,tsnap,nsnap,dsnap,t,lsnap,nsnap_ch) if (lsnap) then snapname=snapbase//nsnap_ch ! ! Write number of massive particles and their data ! open(lun_output,FILE=snapname,FORM='unformatted') write(lun_output) mspar if (mspar/=0) write(lun_output) fsp(1:mspar,:),ipar_nbody(1:mspar) close(lun_output) if (ip<=10 .and. lroot) & print*,'written snapshot ', snapname endif else ! ! Write snapshot without label ! snapname=snapbase open(lun_output,FILE=snapname,FORM='unformatted') write(lun_output) mspar if (mspar/=0) write(lun_output) fsp(1:mspar,:),ipar_nbody(1:mspar) close(lun_output) if (ip<=10 .and. lroot) & print*,'written snapshot ', snapname endif ! call keep_compiler_quiet(flist) ! endsubroutine particles_nbody_write_snapshot !*********************************************************************** subroutine particles_nbody_write_spdim(filename) ! ! Write nspar and mspvar to file. ! ! 01-apr-08/wlad: coded ! character (len=*) :: filename ! open(1,file=filename) write(1,'(3i9)') nspar, mspvar close(1) ! endsubroutine particles_nbody_write_spdim !*********************************************************************** subroutine rprint_particles_nbody(lreset,lwrite) ! ! Read and register print parameters relevant for nbody particles. ! ! 17-nov-05/anders+wlad: adapted ! use Diagnostics use FArrayManager, only: farray_index_append use General, only: itoa ! logical :: lreset,lwr logical, optional :: lwrite ! integer :: iname,ks,j character :: str character (len=intlen) :: sks ! ! Write information to index.pro ! lwr = .false. if (present(lwrite)) lwr=lwrite ! if (lwr) call farray_index_append('imass_nbody',imass_nbody) ! ! Reset everything in case of reset ! if (lreset) then idiag_xxspar=0;idiag_vvspar=0 idiag_torqint=0;idiag_torqext=0 idiag_totenergy=0;idiag_totangmom=0 idiag_torqext_gas=0;idiag_torqext_par=0 idiag_torqint_gas=0;idiag_torqint_par=0 endif ! ! Run through all possible names that may be listed in print.in ! if (lroot .and. ip<14) print*,'rprint_particles: run through parse list' ! ! Now check diagnostics for specific particles ! do ks=1,nspar sks=itoa(ks) do j=1,3 if (j==1) str='x';if (j==2) str='y';if (j==3) str='z' do iname=1,nname call parse_name(iname,cname(iname),cform(iname),& trim(str)//'par'//trim(sks),idiag_xxspar(ks,j)) call parse_name(iname,cname(iname),cform(iname),& 'v'//trim(str)//'par'//trim(sks),idiag_vvspar(ks,j)) enddo ! ! Run through parse list again ! if (lwr) then call farray_index_append('i_'//trim(str)//'par'//trim(sks),idiag_xxspar(ks,j)) call farray_index_append('i_v'//trim(str)//'par'//trim(sks),idiag_vvspar(ks,j)) endif ! enddo ! do iname=1,nname call parse_name(iname,cname(iname),cform(iname),& 'torqint_'//trim(sks),idiag_torqint(ks)) call parse_name(iname,cname(iname),cform(iname),& 'torqext_'//trim(sks),idiag_torqext(ks)) call parse_name(iname,cname(iname),cform(iname),& 'torqext_gas_'//trim(sks),idiag_torqext_gas(ks)) call parse_name(iname,cname(iname),cform(iname),& 'torqext_par_'//trim(sks),idiag_torqext_par(ks)) call parse_name(iname,cname(iname),cform(iname),& 'torqint_gas_'//trim(sks),idiag_torqint_gas(ks)) call parse_name(iname,cname(iname),cform(iname),& 'torqint_par_'//trim(sks),idiag_torqint_par(ks)) enddo ! if (lwr) then call farray_index_append('i_torqint_'//trim(sks),idiag_torqint(ks)) call farray_index_append('i_torqext_'//trim(sks),idiag_torqext(ks)) call farray_index_append('i_torqint_gas'//trim(sks),idiag_torqint(ks)) call farray_index_append('i_torqext_gas'//trim(sks),idiag_torqext(ks)) call farray_index_append('i_torqint_par'//trim(sks),idiag_torqint(ks)) call farray_index_append('i_torqext_par'//trim(sks),idiag_torqext(ks)) endif enddo ! ! Diagnostic related to quantities summed over all particles ! do iname=1,nname call parse_name(iname,cname(iname),cform(iname),& 'totenergy',idiag_totenergy) call parse_name(iname,cname(iname),cform(iname),& 'totangmom',idiag_totangmom) enddo ! if (lwr) then call farray_index_append('i_totenergy',idiag_totenergy) call farray_index_append('i_totangmom',idiag_totangmom) endif ! endsubroutine rprint_particles_nbody !*********************************************************************** ! !*********************************************************************** !*********************************************************************** ! ! STUFF RELATED TO EARLY IMPLEMENTATION OF PARTICLE SINKS. ! !*********************************************************************** !*********************************************************************** ! !*********************************************************************** subroutine create_particles_sink_nbody(f,fp,dfp,ineargrid) ! ! If the flow in any place of the grid has gone gravitationally ! unstable, substitute the local flow by a sink particle. By now, ! the only criterion is the local Jeans mass. The Jeans length is ! ! lambda_J = sqrt(pi*cs2/(G*rho)) ! ! We substitute lambda_J by the biggest resolution element and write ! the condition in terms of density ! ! rho_J = pi*cs2/G*Delta_x^2 ! ! The constant is a free parameter of the module. ! ! 13-mar-08/wlad: coded ! use EquationOfState, only: cs20 use FArrayManager, only: farray_use_global use Mpicomm ! real, dimension(mx,my,mz,mfarray) :: f real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp integer, dimension(mpar_loc,3) :: ineargrid ! real, dimension(maxsink,mspvar) :: fcsp,fcsp_loc real, dimension(maxsink,mspvar) :: fcsp_proc real, dimension(nx,3) :: vvpm real, dimension(nx) :: rho_jeans,rho_jeans_dust,vpm2 real, dimension(nx) :: Delta1,dvolume integer, dimension(nx,mpar_loc) :: pik integer :: i,k,kn,inx0,nc,ks,nc_loc,j integer::nc_proc integer, pointer :: iglobal_cs2 ! logical :: ltime_to_create ! real, dimension(nx,3)::puu real, dimension(nx) ::prho,prhop,pcs2 integer, dimension(nx) :: pnp,npik ! ! Just activate this routine if we want sinks to be created and self- ! gravity is used. If one does not add the t>=tstart_selfgrav flag, the ! particles initially have zero velocity dispersion and collapse ! immediately. ! if (lcreate_sinks.and.t>=tstart_selfgrav) then ! ! just do this for some specific timesteps, otherwise it takes too long! ! ltime_to_create = (mod(it-1,icreate) == 0) if (ltime_to_create.and.llast) then ! if (ldebug) print*,'Entered create_sink_particles_nbody' ! do i=1,nx if (lcartesian_coords) then Delta1(i)=max(dx_1(i),dy_1(mpoint),dz_1(npoint)) elseif (lcylindrical_coords) then Delta1(i)=max(dx_1(i),rcyl_mn1(i)*dy_1(mpoint),dz_1(npoint)) elseif (lspherical_coords) then call fatal_error('create_sink_particle_nbody', & 'not yet implemented for spherical polars') endif enddo ! ! start number of created particles ! fcsp_loc=0. nc_loc=0 ! do imn=1,ny*nz n=nn(imn) m=mm(imn) ! ! define the "pencils" ! if (ldensity_nolog) then prho=f(l1:l2,m,n,irho) else prho=exp(f(l1:l2,m,n,ilnrho)) endif ! if (ldust) then prhop=f(l1:l2,m,n,irhop) pnp =int(f(l1:l2,m,n,inp)) endif ! if (lhydro) puu=f(l1:l2,m,n,iux:iuz) ! if (llocal_iso) then call farray_use_global('cs2',iglobal_cs2) pcs2=f(l1:l2,m,n,iglobal_cs2) else pcs2=cs20 endif ! dvolume = dVol_x(l1:l2)*dVol_y(m)*dVol_z(n) ! ! Jeans analysis of the Gas ! ! The Jeans length is lambda_J = sqrt(pi*cs2/(G*rho)) ! We substitute lambda_J by the biggest resolution ! element and write the condition in terms of density ! ! rho_J = pi*cs2/G*Delta_x^2 ! ! The constant is a free parameter of the module ! if (lcreate_gas) then !test purposes ! if (ldebug) print*,"Entered create sink from gas" ! if (nzgrid/=1) then rho_jeans = create_jeans_constant*3.*pi*pcs2*GNewton1*Delta1**2/32 else rho_jeans = create_jeans_constant* pi*pcs2*GNewton1*Delta1 /8 endif ! ! start particle counter for dust particles ! do i=1,nx if (prho(i)>rho_jeans(i)) then ! ! create a new particle at the center of the grid ! nc_loc=nc_loc+1 ! ! check if we are not creating too many particles ! if (nc_loc>maxsink) then print*,'too many gas sink particles were created. '//& 'Stop and allocated more' print*,'number of created partciles (nc)=',nc_loc print*,'maximum allowed number of sinks '//& 'before merging (maxsink)=',maxsink call fatal_error('create_sink_particles_nbody','') endif ! ! store these particles in a temporary array fcsp - "F array of Created Sink Particles" ! fcsp_loc(nc_loc,ixp) = x(i+nghost) fcsp_loc(nc_loc,iyp) = y(m) fcsp_loc(nc_loc,izp) = z(n) ! ! give it the velocity of the grid cell ! fcsp_loc(nc_loc,ivpx:ivpz) = puu(i,:) ! ! the mass of the new particle is the ! collapsed mass m=[rho - rho_jeans]*dV ! fcsp_loc(nc_loc,imass_nbody) = (prho(i)-rho_jeans(i))*dvolume(i) ! ! and that amount was lost by the grid, so only rho_jeans remains ! if (ldensity_nolog) then f(i+nghost,m,n,irho) = rho_jeans(i) else f(i+nghost,m,n,ilnrho) = log(rho_jeans(i)) endif ! endif enddo endif!if lcreate_gas ! ! Jeans analysis of the dust ! if (ldust.and.lparticles_selfgravity.and.lcreate_dust) then ! if (ldebug) print*,"Entered create sink from dust" ! ! k1s,k2s and ineargrids are already defined. Substitute sound speed ! for vpm2=<(vvp-)^2>, the particle's velocity dispersion ! vvpm=0.0; vpm2=0.0 do k=k1_imn(imn),k2_imn(imn) if (.not.(any(ipar(k)==ipar_nbody))) then inx0=ineargrid(k,1)-nghost vvpm(inx0,:) = vvpm(inx0,:) + fp(k,ivpx:ivpz) endif enddo do i=1,nx if (pnp(i)>1) vvpm(i,:)=vvpm(i,:)/pnp(i) enddo ! vpm2 do k=k1_imn(imn),k2_imn(imn) if (.not.(any(ipar(k)==ipar_nbody))) then inx0=ineargrid(k,1)-nghost vpm2(inx0) = vpm2(inx0) + (fp(k,ivpx)-vvpm(inx0,1))**2 + & (fp(k,ivpy)-vvpm(inx0,2))**2 + & (fp(k,ivpz)-vvpm(inx0,3))**2 endif enddo do i=1,nx if (pnp(i)>1) vpm2(i)=vpm2(i)/pnp(i) enddo ! if (nzgrid/=1) then rho_jeans_dust = create_jeans_constant*3*pi*vpm2*GNewton1*Delta1**2/32 else rho_jeans_dust = create_jeans_constant* pi*vpm2*GNewton1*Delta1 /8 endif ! ! Now fill up an array with the identification index of the particles ! present in each grid cell ! npik=0 do k=k1_imn(imn),k2_imn(imn) if (.not.(any(ipar(k)==ipar_nbody))) then inx0=ineargrid(k,1)-nghost npik(inx0)=npik(inx0)+1 pik(inx0,npik(inx0)) = k endif enddo ! ! Now check for unstable particle concentrations within this pencil ! do i=1,nx if (prhop(i)>rho_jeans_dust(i)) then ! ! This cell is unstable. Remove all particles from it ! if (pnp(i) > 1) then !removing must always be done backwards do kn=pnp(i),1,-1 if (.not.(any(ipar(k)==ipar_nbody))) & call remove_particle(fp,ipar,pik(i,kn)) enddo ! ! create a new particle at the center of the cell ! nc_loc=nc_loc+1 ! ! check if we are not making too many particles ! if (nc_loc>maxsink) then print*,'too many dust sink particles were created. '//& 'Stop and allocated more' print*,'number of created partciles (nc)=',nc_loc print*,'maximum allowed number of sinks '//& 'before merging (maxsink)=',maxsink call fatal_error('create_sink_particle_nbody','') endif ! fcsp_loc(nc_loc,ixp) = x(i+nghost) fcsp_loc(nc_loc,iyp) = y(m) fcsp_loc(nc_loc,izp) = z(n) ! ! give it the mean velocity of the group of particles that ! was accreted ! fcsp_loc(nc_loc,ivpx:ivpz) = vvpm(i,:) ! ! the mass of the new particle is the ! collapsed mass M=m_particle*np ! fcsp_loc(nc_loc,imass_nbody) = mp_swarm*pnp(i) ! endif endif enddo ! endif !if (ldust) enddo if ((ip<=8).and.nc_loc/=0) & print*,'I, processor ',iproc,& ' created ',nc_loc,' particles' ! ! Finished creating them. Now merge them across the processors ! into a single array ! fcsp_proc=0. nc=0 if (lmpicomm) then if (.not.lroot) then if (ldebug) print*,'processor ',iproc,'sending ',nc_loc,' particles' call mpisend_int(nc_loc,root,222) if (nc_loc/=0) then print*,'sending',fcsp_loc(1:nc_loc,:) call mpisend_real(fcsp_loc(1:nc_loc,:),(/nc_loc,mspvar/),root,111) endif else !root nc_proc=nc_loc if (nc_loc/=0) & fcsp(nc+1:nc+nc_proc,:)=fcsp_loc(1:nc_proc,:) nc=nc+nc_loc !get from the other processors do j=1,ncpus-1 call mpirecv_int(nc_proc,j,222) if (ldebug) print*,'received ',nc_proc,'particles from processor,',j ! call fatal_error("","") if (nc_proc/=0) then if (ldebug) print*,'receiving ',nc_proc,' particles from processor ',j call mpirecv_real(fcsp_proc(1:nc_proc,:),(/nc_proc,mspvar/),j,111) if (ldebug) print*,'particle received=',fcsp_proc(1:nc_proc,:) fcsp(nc+1:nc+nc_proc,:)=fcsp_proc(1:nc_proc,:) nc=nc+nc_proc endif enddo endif else !serial, just copy it nc =nc_loc fcsp=fcsp_loc endif ! ! Call merge only if particles were created ! call mpibcast_int(nc) ! if (nc/=0) then call merge_and_share(fcsp,nc,fp) call bcast_nbodyarray(fp) endif ! endif ! ! print to the screen the positions and velocities of the ! newly generated sinks ! if (ldebug) then print*,'---------------------------' print*,'the massive particles present are:' do ks=1,mspar print*,'ks=',ks print*,'positions=',fsp(ks,ixp:izp) print*,'velocities=',fsp(ks,ivpx:ivpz) print*,'mass=',fsp(ks,imass_nbody) print*,'ipar_nbody=',ipar_nbody print*,'' enddo print*,'---------------------------' endif ! endif ! call keep_compiler_quiet(dfp) ! endsubroutine create_particles_sink_nbody !*********************************************************************** subroutine remove_particles_sink_nbody(f,fp,dfp,ineargrid) ! ! Subroutine for taking particles out of the simulation due to their ! proximity to a sink particle or sink point. ! ! Just a dummy routine for now. ! ! 25-sep-08/anders: coded ! real, dimension(mx,my,mz,mfarray) :: f real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp integer, dimension(mpar_loc,3) :: ineargrid ! call keep_compiler_quiet(f) call keep_compiler_quiet(fp) call keep_compiler_quiet(dfp) call keep_compiler_quiet(ineargrid) ! endsubroutine remove_particles_sink_nbody !************************************************************************** subroutine merge_and_share(fcsp,nc,fp) ! ! Subroutine to merge a gravitationally unstable ! cluster of particles into a sink particle ! ! 13-mar-08/wlad: coded ! use Mpicomm ! real, dimension (mpar_loc,mparray) :: fp real, dimension(maxsink,mspvar) :: fcsp real, dimension(nspar,mspvar) :: fleft ! integer :: nc,nf,kc,j ! ! NILS: Should fcsp_mig have dimensionality mparray????? ! real, dimension(0:ncpus-1,nspar,mpvar) :: fcsp_mig integer, dimension(0:ncpus-1) :: nsmig integer :: iz0,iy0,ipz_rec,ipy_rec,iproc_rec,npar_tot integer:: ns,i double precision :: dy1,dz1 ! ! How many particles are present in the whole grid? Needed to set ! the new ipars of the particles that will be created. ! call mpiallreduce_sum_int(npar_loc,npar_tot) ! ! The root processor is the only one that has the right fcsp (the ! array of particles that were created). Merge them with a friends of ! friends algorithm ! if (lroot) then ! call friends_of_friends(fcsp,fleft,nc,nf) ! ! Friends of friends merged the created particles into few massive ! clusters. These new clusters are truly sinks. Now we have nf new ! particles. ! pmass(mspar+1:mspar+nf)=fleft(1:nf,imass_nbody) ! ! Allocate the new ipar_nbodies ! do i=1,nf ipar_nbody(mspar+i)=mspar+i ! Activate accretion for the newly created sinks if (laccrete_when_create) then ladd_mass(mspar+i)=.true. laccretion(mspar+i)=.true. endif enddo mspar=mspar+nf ! ! But check if we did not end up with too many particles ! if (mspar>nspar) then print*,'after friends_of_friends, we still have '//& 'too many nbody particles.'//& 'Stop and allocated more' print*,'the total number of massive particles (mspar)=',mspar print*,'is bigger than the maximum number of '//& 'allowed massive particles (nspar)=',nspar call fatal_error("merge and share","") endif endif ! ! Broadcast mspar, ipar_nbody and pmass ! call mpibcast_int(mspar) call mpibcast_int(ipar_nbody(1:mspar),mspar) call mpibcast_real(pmass,mspar) call mpibcast_logical(ladd_mass,mspar) call mpibcast_logical(laccretion,mspar) ! ! Migrate the particles to their respective processors ! if (lmpicomm) then ! dy1=1/dy; dz1=1/dz !y0_mig=0.5*(y(m1)+y(m1-1)); y1_mig=0.5*(y(m2)+y(m2+1)) !z0_mig=0.5*(z(n1)+z(n1-1)); z1_mig=0.5*(z(n2)+z(n2+1)) ! ! Only the root processor knows where the particle is (fleft) ! nsmig(:)=0 ! if (lroot) then ! do kc=1,nf ! y-index ipy_rec=ipy iy0=nint((fleft(kc,iyp)-y(m1))*dy1+nygrid)-nygrid+1 do while (iy0>ny) ipy_rec=ipy_rec+1 iy0=iy0-ny enddo ! z-index ipz_rec=ipz iz0=nint((fleft(kc,izp)-z(n1))*dz1+nzgrid)-nzgrid+1 do while (iz0>nz) ipz_rec=ipz_rec+1 iz0=iz0-nz enddo ! Calculate serial index of receiving processor. iproc_rec=ipy_rec+nprocy*ipz_rec ! Check that the processor exists if (iproc_rec>=ncpus .or. iproc_rec<0) then call warning('merge_and_share','',iproc) print*, 'merge_and_share: receiving proc does not exist' print*, 'merge_and_share: iproc, iproc_rec=', & iproc, iproc_rec call fatal_error_local("","") endif ! ! Fill up the migration arrays if the particles are not at the root processor ! if (iproc_rec/=root) then nsmig(iproc_rec)=nsmig(iproc_rec)+1 fcsp_mig(iproc_rec,nsmig(iproc_rec),:)=fleft(kc,1:mpvar) else ! The particle is at the root processor, so create it here npar_loc=npar_loc+1 fp(npar_loc,1:mpvar)=fleft(kc,1:mpvar) ipar(npar_loc)=npar_tot+1 endif enddo !loop over particles ! ! Send them now ! do j=1,ncpus-1 ns=nsmig(j) call mpisend_int(ns, j, 111) if (ns/=0) then call mpisend_real(fcsp_mig(j,1:ns,:), (/ns,mpvar/), j, 222) endif enddo ! else !not lroot ! ! The other processors receive the particles ! call mpirecv_int(nsmig(iproc), root, 111) ns=nsmig(iproc) if (ns/=0) then call mpirecv_real(fp(npar_loc+1:npar_loc+ns,1:mpvar),(/ns,mpvar/),root,222) ! update the relevant quantities do kc=1,ns npar_loc=npar_loc+1 ipar(npar_loc)=npar_tot+kc enddo ! endif ! endif !if lroot ! else !serial version. Just copy it all ! do kc=1,nf npar_loc=npar_loc+1 fp(npar_loc,1:mpvar)=fleft(kc,1:mpvar) ipar(npar_loc)=npar_tot+1 enddo ! endif ! if (ldebug) then print*,'merge_and_share finished. ' print*,'We have now ',mspar,' nbody particles, ' print*,'located at ' print*,ipar_nbody(1:mspar) endif ! endsubroutine merge_and_share !*********************************************************************** subroutine friends_of_friends(fcsp,fleft,nc,nfinal) ! ! Algorithm that cluster particles together based on a proximity ! threshold. It takes a particle and find all its neighbours (=friends). ! Then loops through the neighbours finding the neighbours of the ! neighbours (friends of friends), and so on. ! ! 13-mar-08/wlad: coded ! integer, intent(in) :: nc real, dimension(maxsink,mspvar) :: fcsp real, dimension(nspar,mspvar) :: fleft integer, dimension(nc,nc) :: clusters integer, dimension(nc) :: cluster logical, dimension(nc) :: lchecked real, dimension(mspvar) :: particle ! integer :: ks,kf,nclusters,nf,i,nfinal integer :: min_number_members,kpar ! intent(out) :: nfinal,fleft ! min_number_members=3 ! ! All particles are initially unchecked ! lchecked(1:nc)=.false. ! ! No clusters yet, counter is zero ! nclusters=0 ! ! And the final number of particles is the same as the initial ! kf=0 ! ! Loop through the particles to cluster them ! do ks=1,nc if (.not.lchecked(ks)) then ! ! If a particle is unchecked, start a new cluster ! call make_cluster(ks,lchecked,cluster,nc,nf,fcsp) ! ! Cluster done, check if the number of particles is enough ! if (nf >= min_number_members) then ! ! Okay, it is a cluster, raise the counter ! nclusters=nclusters+1 clusters(1:nf,nclusters)=cluster(1:nf) ! ! Collapse the cluster into a single sink particle ! call collapse_cluster(cluster,nf,fcsp,particle) ! ! this is all a single particle ! kf=kf+1 fleft(kf,:)=particle ! else !this (or these) are isolated sinks do i=1,nf kf=kf+1 kpar=cluster(i) fleft(kf,:)=fcsp(kpar,:) enddo endif ! endif enddo ! nfinal=kf ! endsubroutine friends_of_friends !*********************************************************************** subroutine make_cluster(ks,lchecked,cluster,nc,nf,fcsp) ! ! This subroutine starts a cluster of particles, by checking the ! first one. Once a particle is checked, it will also look for its ! neighbours, checking them. In the end, all friends and friends ! of friends are checked into the cluster. ! ! 13-mar-08/wlad: coded ! integer, intent(in) :: nc real, dimension(maxsink,mspvar) :: fcsp integer, intent(inout), dimension(nc) :: cluster logical, intent(inout), dimension(nc) :: lchecked integer, intent(in) :: ks integer, intent(out) :: nf ! integer :: ic ! ! Start the cluster counter ! ic=0 ! ! Tag the particles as checked. The subroutine will also ! loop through its friends and friends of friends ! call check_particle(ks,lchecked,cluster,ic,nc,fcsp) ! ! When it ends this loop, the cluster is done. Nf is the final ! number of members in the cluster ! nf=ic ! endsubroutine make_cluster !*********************************************************************** subroutine check_particle(k,lchecked,cluster,ic,nc,fcsp) ! ! Check (tag) a particle and its neighbours ! ! 13-mar-08/wlad : coded ! integer, intent(in) :: nc real, dimension(maxsink,mspvar) :: fcsp integer, intent(inout), dimension(nc) :: cluster logical, intent(inout), dimension(nc) :: lchecked integer, intent(inout) :: ic integer, intent(in) :: k ! ! Add the particle to the cluster ! ic=ic+1 cluster(ic)=k ! ! Tag it as visited and add its friends to the cluster ! as well ! lchecked(k)=.true. call add_friends(k,lchecked,cluster,ic,nc,fcsp) ! endsubroutine check_particle !*********************************************************************** subroutine add_friends(par,lchecked,cluster,ic,nc,fcsp) ! ! Add all neighbours of a particle to the cluster. This subroutine ! is called by check_particle, but also calls it. The procedure is ! therefore recursive: the loop will be over when all friends of ! friends are added to the cluster. ! ! 13-mar-08/wlad: coded ! integer, intent(in) :: nc real, dimension(maxsink,mspvar) :: fcsp logical, intent(inout), dimension(nc) :: lchecked integer, intent(inout), dimension(nc) :: cluster integer, intent(inout) :: ic integer, intent(in) :: par ! real :: dist,link_length integer :: k ! ! linking length ! link_length=4.*dx ! ! Loop through the particle ! do k=1,nc if (.not.lchecked(k)) then call get_particles_interdistance(& fcsp(k,ixp:izp),fcsp(par,ixp:izp),& DISTANCE=dist) ! if (dist<=link_length) then ! ! if so, add it to the cluster, tag it and call its friends ! call check_particle(k,lchecked,cluster,ic,nc,fcsp) endif ! done endif enddo ! endsubroutine add_friends !*********************************************************************** subroutine collapse_cluster(cluster,nf,fcsp,particle) ! ! Collapse the cluster into a single particle ! with the total mass, momentum and energy of the ! collapsed stuff ! ! 13-mar-08/wlad: coded ! integer, intent(in) :: nf real, dimension(maxsink,mspvar) :: fcsp integer, intent(inout), dimension(nf) :: cluster real, dimension(mspvar) :: particle integer :: ic,k ! real :: mass,xcm,ycm,vxcm,vycm real, dimension(3) :: position,velocity,xxp,vvp ! intent(out):: particle ! mass=0;position=0;velocity=0 ! do ic=1,nf k=cluster(ic) mass=mass+fcsp(k,imass_nbody) if (lcartesian_coords) then xxp(1)=fcsp(k,ixp) ; vvp(1)=fcsp(k,ivpx) xxp(2)=fcsp(k,iyp) ; vvp(2)=fcsp(k,ivpy) xxp(3)=fcsp(k,izp) ; vvp(3)=fcsp(k,ivpz) else if (lcylindrical_coords) then xxp(1)=fcsp(k,ixp)*cos(fcsp(k,iyp)) xxp(2)=fcsp(k,ixp)*sin(fcsp(k,iyp)) xxp(3)=fcsp(k,izp) ! vvp(1)=fcsp(k,ivpx)*cos(fcsp(k,iyp))-fcsp(k,ivpy)*sin(fcsp(k,iyp)) vvp(2)=fcsp(k,ivpx)*sin(fcsp(k,iyp))+fcsp(k,ivpy)*cos(fcsp(k,iyp)) vvp(3)=fcsp(k,ivpz) else if (lspherical_coords) then call fatal_error("","") endif position=position+fcsp(k,imass_nbody)*xxp velocity=velocity+fcsp(k,imass_nbody)*vvp enddo ! position=position/mass if (lcylindrical_coords) then xcm=position(1);vxcm=velocity(1) ycm=position(2);vycm=velocity(2) position(1)=sqrt(xcm**2+ycm**2) position(2)=atan2(ycm,xcm) ! velocity(1)= vxcm*cos(position(2))+vycm*sin(position(2)) velocity(2)=-vxcm*sin(position(2))+vycm*cos(position(2)) endif ! ! put it onto the output array ! particle(ixp : izp)=position particle(ivpx:ivpz)=velocity particle(imass_nbody) =mass ! endsubroutine collapse_cluster !*********************************************************************** endmodule Particles_nbody