! $Id: pointmasses.f90,v 1.1 2019/02/02 03:54:41 wlyra Exp $ ! ! This module takes care of direct N-body gravity between point masses. ! !** AUTOMATIC CPARAM.INC GENERATION **************************** ! ! Declare (for generation of cparam.inc) the number of f array ! variables and auxiliary variables added by this module ! ! MQVAR CONTRIBUTION 7 ! MQAUX CONTRIBUTION 0 ! MPVAR CONTRIBUTION 3 ! CPARAM logical, parameter :: lpointmasses=.true. ! !*************************************************************** module PointMasses ! use Cdata use General, only: keep_compiler_quiet use Messages use Cparam ! implicit none ! include 'pointmasses.h' ! character (len=labellen), dimension(mqarray) :: qvarname real, dimension(nqpar,mqarray) :: fq=0.0 real, dimension(nqpar,mqvar) :: dfq=0.0 real, dimension(nqpar,3) :: dfq_cart=0.0 real, dimension(nqpar) :: xq0=0.0, yq0=0.0, zq0=0.0 real, dimension(nqpar) :: vxq0=0.0, vyq0=0.0, vzq0=0.0 real, dimension(nqpar) :: pmass=0.0, r_smooth=impossible, pmass1 real, dimension(nqpar) :: r1_smooth real, dimension(nqpar) :: frac_smooth=0.4 real, dimension(nqpar) :: accrete_hills_frac=0.2, final_ramped_mass=0.0 real, dimension(nqpar) :: StokesNumber=1. real, pointer :: rhs_poisson_const, tstart_selfgrav real :: eccentricity=0.0, semimajor_axis=1.0 real :: totmass, totmass1 real :: GNewton1, GNewton=impossible, density_scale=0.001 real :: cdtq=0.1 real :: hills_tempering_fraction=0.8 real :: ugas=0.0,Omega_coriolis=0.0 real :: tau_accretion=1.0 ! integer :: ramp_orbits=5 integer :: iprimary=1, isecondary=2 integer :: ivpx_cart=0,ivpy_cart=0,ivpz_cart=0 integer :: imass=0, ixq=0, iyq=0, izq=0, ivxq=0, ivyq=0, ivzq=0 integer :: nqvar=0 ! logical, dimension(nqpar) :: lcylindrical_gravity_nbody=.false. logical, dimension(nqpar) :: lfollow_particle=.false., laccretion=.false. logical :: llive_secondary=.false., lnorm=.true. logical :: lreset_cm=.true., lnogravz_star=.false., lexclude_frozen=.false. logical :: lramp=.false. logical :: ldt_pointmasses=.true. logical :: ldust=.false. logical :: lretrograde=.false. logical :: lnoselfgrav_primary=.true. logical :: lgas_gravity=.true.,ldust_gravity=.false. logical :: lcorrect_gasgravity_lstart=.false. logical :: lexclude_hills=.false. logical :: lgas_removal=.false.,lmomentum_removal=.false. logical :: ladd_dragforce=.false.,lquadratic_drag=.false.,llinear_drag=.true. logical :: lcoriolis_force=.false. logical :: l2D,l3D ! character (len=labellen) :: initxxq='random', initvvq='nothing' character (len=labellen), dimension (nqpar) :: ipotential_pointmass='newton' character (len=2*bclen+1) :: bcqx='p', bcqy='p', bcqz='p' ! type IndexDustParticles integer :: ixw=0,iyw=0,izw=0 integer :: ivxw=0,ivyw=0,ivzw=0 endtype IndexDustParticles type (IndexDustParticles), save :: index ! namelist /pointmasses_init_pars/ & initxxq, initvvq, xq0, yq0, zq0, vxq0, vyq0, vzq0, & pmass, r_smooth, lcylindrical_gravity_nbody, lexclude_frozen, GNewton, & bcqx, bcqy, bcqz, ramp_orbits, lramp, final_ramped_mass, & laccretion, & accrete_hills_frac, iprimary, & ldt_pointmasses, cdtq, lretrograde, & eccentricity, semimajor_axis, & ipotential_pointmass, density_scale,& lgas_gravity,ldust_gravity,lcorrect_gasgravity_lstart,& frac_smooth ! namelist /pointmasses_run_pars/ & lreset_cm, & lnogravz_star, lfollow_particle, llive_secondary, lexclude_frozen, & GNewton, bcqx, bcqy, bcqz, & laccretion, accrete_hills_frac, iprimary, & ldt_pointmasses, cdtq, hills_tempering_fraction, & ipotential_pointmass, density_scale,& lgas_gravity,ldust_gravity,& ladd_dragforce,ugas,StokesNumber,& lquadratic_drag,llinear_drag,lcoriolis_force,Omega_coriolis,& frac_smooth,lexclude_hills,tau_accretion,lgas_removal,lmomentum_removal ! integer, dimension(nqpar,3) :: idiag_xxq=0,idiag_vvq=0 integer, dimension(nqpar) :: idiag_torqint=0,idiag_torqext=0 integer, dimension(nqpar) :: idiag_torqext_gas=0,idiag_torqext_par=0 integer, dimension(nqpar) :: idiag_torqint_gas=0,idiag_torqint_par=0 integer, dimension(nqpar) :: idiag_period=0,idiag_torque=0 integer :: idiag_totenergy=0,idiag_mdot_pt=0 ! contains !*********************************************************************** subroutine register_pointmasses() ! ! Set up indices for access to the f and fq. ! ! 27-aug-06/wlad: adapted ! use Particles_main, only: append_particle_index ! integer :: iqvar ! if (lroot) call svn_id( & "$Id: pointmasses.f90,v 1.1 2019/02/02 03:54:41 wlyra Exp $") ! ! No need to solve the N-body equations for non-N-body problems. ! if (nqpar < 2) then if (lroot) write(0,*) 'nqpar = ', nqpar call fatal_error('register_pointmasses','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 ! ! Auxiliary variables for polar coordinates ! ixq = nqvar+1 qvarname(ixq)='ixq' iyq = nqvar+2 qvarname(iyq)='iyq' izq = nqvar+3 qvarname(izq)='izq' nqvar=nqvar+3 ! ivxq = nqvar+1 qvarname(ivxq)='ivxq' ivyq = nqvar+2 qvarname(ivyq)='ivyq' ivzq = nqvar+3 qvarname(ivzq)='ivzq' nqvar=nqvar+3 ! ! Set up mass as particle index. Plus seven, since the other 6 are ! used by positions and velocities. ! imass=nqvar+1 qvarname(imass)='imass' nqvar=nqvar+1 ! ! Check that the fq and dfq arrays are big enough. ! if (nqvar > mqvar) then if (lroot) write(0,*) 'nqvar = ', nqvar, ', mqvar = ', mqvar call fatal_error('register_pointmasses','nqvar > mqvar') endif ! if (lparticles) then call append_particle_index('ivpx_cart',ivpx_cart) call append_particle_index('ivpy_cart',ivpy_cart) call append_particle_index('ivpz_cart',ivpz_cart) endif ! if (lroot) then open(3,file=trim(datadir)//'/qvarname.dat',status='replace') do iqvar=1,mqarray write(3,"(i4,2x,a)") iqvar, qvarname(iqvar) enddo close(3) endif ! endsubroutine register_pointmasses !*********************************************************************** subroutine initialize_pointmasses(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 logical :: lxpresent,lypresent,lzpresent,l2Dcyl,l2Dsph ! integer :: ks ! ! Look for initialized masses. ! do ks=1,nqpar if (pmass(ks)/=0.0) then fq(ks,imass)=pmass(ks) else call fatal_error("initialize_pointmasses",& "one of the bodies has zero mass") endif enddo ! ! Just needs to do this starting, otherwise mqar will be overwritten after ! reading the snapshot ! if (.not.lstart) pmass(1:nqpar)=fq(1:nqpar,imass) ! ! When first called, nqpar was zero, so no diagnostic index was written to ! index.pro ! call rprint_pointmasses(.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_pointmasses') 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_pointmasses') ! ! inverse mass ! if (lstart) then if (lramp) then do ks=1,nqpar if (ks/=iprimary) pmass(ks) = epsi enddo pmass(iprimary)=1-epsi*(nqpar-1) endif else !read imass from the snapshot pmass(1:nqpar)=fq(1:nqpar,imass) endif ! pmass1=1./max(pmass,tini) ! ! inverse total mass ! totmass=sum(pmass) if (totmass/=1) & call warning("initialize_pointmasses","The masses do not sum up to one!") totmass1=1./max(totmass, tini) ! ! Check for consistency between the cylindrical gravities switches ! from cdata and the one from point mass. ! if (((lcylindrical_gravity).and.& (.not.lcylindrical_gravity_nbody(iprimary))).or.& (.not.lcylindrical_gravity).and.& (lcylindrical_gravity_nbody(iprimary))) then call fatal_error('initialize_pointmasses','inconsistency '//& 'between lcylindrical_gravity from cdata and the '//& 'one from point mass') endif ! ! Smoothing radius ! if (any(r_smooth == impossible)) then do ks=1,nqpar if (ks/=iprimary) then r_smooth(ks) = frac_smooth(ks) * xq0(ks) * (pmass(ks)/3.)**(1./3) else r_smooth(ks) = rsmooth endif enddo endif where (r_smooth/=0.) r1_smooth=1./r_smooth elsewhere r1_smooth=0. !MR: something better here? endwhere ! if (rsmooth/=r_smooth(iprimary)) then print*,'rsmooth from cdata=',rsmooth print*,'r_smooth(iprimary)=',r_smooth(iprimary) call fatal_error('initialize_pointmasses','inconsistency '//& 'between rsmooth from cdata and the '//& 'one from point mass') endif ! ! The presence of dust particles needs to be known. ! if (npar > nqpar) ldust=.true. ! ! Check if the run is 2D or 3D to use cylindrical gravity or not ! lxpresent=(nxgrid/=1) lypresent=(nygrid/=1) lzpresent=(nzgrid/=1) ! l3D = lxpresent.and. lypresent .and. lzpresent l2Dcyl= lcylindrical_coords.and.(lxpresent.and. lypresent .and.(.not.lzpresent)) l2Dsph= lspherical_coords .and.(lxpresent.and.(.not.lypresent).and. lzpresent ) l2D=l2Dcyl.or.l2Dsph ! if (lcylindrical_coords.and.bcqy=='p2pi'.and.lselfgravity) & call fatal_error("initialize_pointmasses",& "p2pi assumes the range is -pi/pi; selfgravity_logspirals" // & "assumes the range is 0/2pi. Change bcqy to 'p'") ! call keep_compiler_quiet(f) ! endsubroutine initialize_pointmasses !*********************************************************************** subroutine pencil_criteria_pointmasses() ! ! All pencils that the pointmasses module depends on are specified here. ! ! TODO: (08/2022) Change the criteria so that p%rho is called only when needed, ! instead of every time the density module is used. ! ! 22-sep-06/wlad: adapted ! if (ldensity) lpenc_requested(i_rho)=.true. if (ldust.and.llive_secondary) & lpenc_requested(i_rhop)=.true. ! if (llive_secondary.and.(.not.lselfgravity)) then if (l2D.or.(l3D.and.lcylindrical_gravity)) then lpenc_requested(i_rcyl_mn)=.true. else lpenc_requested(i_r_mn)=.true. endif endif ! if (lmomentum_removal) lpenc_requested(i_uu)=.true. ! if (idiag_totenergy/=0) then lpenc_diagnos(i_u2)=.true. lpenc_diagnos(i_rho)=.true. endif ! if (l2D.or.(l3D.and.lcylindrical_gravity)) then lpenc_diagnos(i_rcyl_mn)=.true. else lpenc_diagnos(i_r_mn)=.true. endif if (ldust) lpenc_diagnos(i_rhop)=.true. ! endsubroutine pencil_criteria_pointmasses !*********************************************************************** subroutine pencil_interdep_pointmasses(lpencil_in) ! ! Interdependency among pencils provided by the pointmasses module ! is specified here. ! ! 22-sep-06/wlad: adapted ! logical, dimension(npencils) :: lpencil_in ! call keep_compiler_quiet(lpencil_in) ! endsubroutine pencil_interdep_pointmasses !*********************************************************************** subroutine calc_pencils_pointmasses(f,p) ! ! Calculate point mass 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_pointmasses !*********************************************************************** subroutine init_pointmasses(f)!,fp) ! ! Initial positions and velocities of point mass particles. ! Overwrite the position asserted by the dust module ! ! 17-nov-05/anders+wlad: adapted ! use General, only: random_number_wrapper use Sub use Mpicomm, only: mpibcast_real ! real, dimension (mx,my,mz,mfarray) :: f !real, dimension (mpar_loc,mparray) :: fp real, dimension(nqpar) :: kep_vel,sma real, dimension(nqpar,3) :: velocity real, dimension(nqpar,3) :: positions real :: tmp,parc real :: absolute_offset_star,baricenter_secondaries real :: velocity_baricenter_secondaries,mass_secondaries integer :: k,ks ! intent (in) :: f !intent (out) :: fp ! ! Shortcuts ! positions(1:nqpar,1) = xq0 ; velocity(1:nqpar,1) = vxq0 positions(1:nqpar,2) = yq0 ; velocity(1:nqpar,2) = vyq0 positions(1:nqpar,3) = zq0 ; velocity(1:nqpar,3) = vzq0 ! ! Initialize particles' positions. ! select case (initxxq) ! case ('nothing') if (lroot) print*, 'init_pointmasses: nothing' ! case ('origin') if (lroot) then print*, 'init_pointmasses: All point mass particles at origin' fq(1:nqpar,ixq:izq)=0.0 endif ! case ('constant') if (lroot) & print*, 'init_pointmasses: Place point mass particles at x,y,z=', & xq0, yq0, zq0 do k=1,nqpar fq(k,ixq:izq)=positions(k,1:3) enddo ! case ('random') if (lroot) print*, 'init_pointmasses: Random particle positions' do ks=1,nqpar if (nxgrid/=1) call random_number_wrapper(positions(ks,ixq)) if (nygrid/=1) call random_number_wrapper(positions(ks,iyq)) if (nzgrid/=1) call random_number_wrapper(positions(ks,izq)) enddo ! if (nxgrid/=1) & positions(1:nqpar,ixq)=xyz0_loc(1)+positions(1:nqpar,ixq)*Lxyz_loc(1) if (nygrid/=1) & positions(1:nqpar,iyq)=xyz0_loc(2)+positions(1:nqpar,iyq)*Lxyz_loc(2) if (nzgrid/=1) & positions(1:nqpar,izq)=xyz0_loc(3)+positions(1:nqpar,izq)*Lxyz_loc(3) ! do k=1,nqpar ! fq(k,ixq:izq)=positions(k,1:3) ! ! Correct for non-existing dimensions (not really needed, I think). ! if (nxgrid==1) fq(k,ixq)=x(nghost+1) if (nygrid==1) fq(k,iyq)=y(nghost+1) if (nzgrid==1) fq(k,izq)=z(nghost+1) ! enddo ! case ('fixed-cm') if (lgrav) then print*,'a gravity module is being used. Are you using '//& 'both a fixed central gravity and point mass gravity? '//& 'better stop and check' call fatal_error('init_pointmasses','') 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(yq0/=0)) then if (lspherical_coords) then call fatal_error('init_pointmasses','not yet generalized'//& ' for non-zero initial inclinations') else call fatal_error('init_pointmasses','not yet generalized'//& ' for non-zero azimuthal initial position') endif endif if (any(zq0/=0)) then if (lspherical_coords) then call fatal_error('init_pointmasses','not yet generalized'//& ' for non-zero azimuthal initial position') else call fatal_error('init_pointmasses','point mass code not'//& ' yet generalized to allow initial inclinations') endif endif if (lcylindrical_coords.or.lspherical_coords) then if (any(xq0<0)) & call fatal_error('init_pointmasses', & 'in cylindrical coordinates '//& 'all the radial positions must be positive') endif ! if (lroot) then print*,'init_pointmasses: 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:nqpar,iyq)=pi/2 endif ! mass_secondaries = 0. baricenter_secondaries=0. do ks=1,nqpar if (ks/=iprimary) then sma(ks) = abs(positions(ks,1)) mass_secondaries = mass_secondaries + pmass(ks) baricenter_secondaries = baricenter_secondaries + positions(ks,1)*pmass(ks) endif enddo absolute_offset_star = abs(baricenter_secondaries) ! ! Fixed-cm assumes that the total mass is always one. The mass of the ! star is adjusted to ensure this. ! pmass(iprimary)=1.- mass_secondaries pmass1=1./max(pmass,tini) totmass=1.;totmass1=1. ! absolute_offset_star = absolute_offset_star*totmass1 if (mass_secondaries >= 1.0) & call fatal_error('init_pointmasses', & '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') ! ! Correct the semimajor of the secondaries by the offset they generate. ! do ks=1,nqpar ! sign(A,B) returns the value of A with the sign of B if (ks/=iprimary) & positions(ks,1)=sign(1.,positions(ks,1))* (sma(ks) - absolute_offset_star) enddo ! ! The last one (star) fixes the CM at Rcm=zero ! if (lcartesian_coords) then !put the star opposite to the baricenter of planets positions(iprimary,1)=-sign(1.,baricenter_secondaries)*absolute_offset_star elseif (lcylindrical_coords) then !put the star in positive coordinates, with pi for azimuth positions(iprimary,1)=absolute_offset_star positions(iprimary,2)=pi elseif (lspherical_coords) then positions(iprimary,1)=absolute_offset_star positions(iprimary,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 ! do k=1,nqpar ! ! Here we substitute the first nqpar dust particles by massive ones, ! since the first ipars are less than nqpar ! fq(k,ixq:izq)=positions(k,1:3) ! ! Correct for non-existing dimensions (not really needed, I think) ! if (nxgrid==1) fq(k,ixq)=x(nghost+1) if (nygrid==1) fq(k,iyq)=y(nghost+1) if (nzgrid==1) fq(k,izq)=z(nghost+1) ! enddo ! case ('eccentric') ! ! Coded only for 2 bodies ! if (nqpar /= 2) call fatal_error("init_pointmasses",& "This initial condition is currently coded for 2 massive particles only.") ! ! Define isecondary. iprimary=1 and isecondary=2 is default ! if (iprimary == 2) isecondary=1 ! ! Reassign total mass of the star so that totmass=1 ! pmass(iprimary)=1-pmass(isecondary) pmass1=1./max(pmass,tini) totmass=1.;totmass1=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(isecondary,1)=(1+eccentricity) * semimajor_axis * pmass( iprimary)/totmass positions( iprimary,1)=(1+eccentricity) * semimajor_axis * pmass(isecondary)/totmass ! ! Azimuthal position. Planet and star phased by pi. ! positions(isecondary,2)=0 positions( iprimary,2)=pi ! do k=1,nqpar !if (ipar(k) <= nqpar) then fq(k,ixq:izq) = positions(k,1:3) !endif enddo ! case default if (lroot) print*,'init_pointmasses: No such such value for'//& ' initxxq: ',trim(initxxq) call fatal_error('init_pointmasses','') ! endselect ! ! Redistribute particles among processors (now that positions are determined). ! call boundconds_pointmasses if (lroot) then print*,'init_pointmasses: particle positions' do ks=1,nqpar print*,'particle',ks, 'position x,y,z:',fq(ks,ixq:izq) enddo endif ! ! Initial particle velocity. ! select case (initvvq) ! case ('nothing') if (lroot) print*, 'init_pointmasses: No particle velocity set' ! case ('zero') if (lroot) then print*, 'init_pointmasses: Zero particle velocity' fq(1:nqpar,ivxq:ivzq)=0. endif ! case ('constant') if (lroot) then print*, 'init_pointmasses: Constant particle velocity' print*, 'init_pointmasses: vxq0, vyq0, vzq0=', vxq0, vyq0, vzq0 endif do k=1,nqpar !if (ipar(k) <= nqpar) & fq(k,ivxq:ivzq)=velocity(k,1:3) enddo ! case ('fixed-cm') ! ! Keplerian velocities for the planets ! velocity_baricenter_secondaries=0. do ks=1,nqpar if (ks/=iprimary) then kep_vel(ks)=sqrt(1./sma(ks)) !circular velocity velocity_baricenter_secondaries = velocity_baricenter_secondaries + kep_vel(ks)*pmass(ks) endif enddo velocity_baricenter_secondaries = velocity_baricenter_secondaries*totmass do ks=1,nqpar if (ks/=iprimary) then if (lcartesian_coords) then velocity(ks,2) = sign(1.,positions(ks,1))*(kep_vel(ks) - velocity_baricenter_secondaries) elseif (lcylindrical_coords) then !positive for the planets velocity(ks,2) = kep_vel(ks) - velocity_baricenter_secondaries elseif (lspherical_coords) then velocity(ks,3) = kep_vel(ks) - velocity_baricenter_secondaries endif endif enddo ! if (lcorrect_gasgravity_lstart) call initcond_correct_gasgravity(f,velocity,isecondary) ! ! The last one (star) fixes the CM also with velocity zero ! if (lcartesian_coords) then velocity(iprimary,2)= -sign(1.,baricenter_secondaries)*velocity_baricenter_secondaries elseif (lcylindrical_coords) then velocity(iprimary,2)= velocity_baricenter_secondaries elseif (lspherical_coords) then velocity(iprimary,3)= velocity_baricenter_secondaries endif ! ! Revert all velocities if retrograde ! if (lretrograde) velocity=-velocity ! ! Allocate the point mass particles ! do k=1,nqpar fq(k,ivxq:ivzq) = velocity(k,1:3) enddo ! case ('eccentric') ! ! Coded only for 2 bodies ! if (nqpar /= 2) call fatal_error("init_pointmasses",& "This initial condition is currently coded for 2 massive particles only.") ! ! Define isecondary. iprimary=1 and isecondary=2 is default. ! if (iprimary == 2) isecondary=1 velocity(isecondary,2) = sqrt((1-eccentricity)/(1+eccentricity) * GNewton/semimajor_axis) * pmass( iprimary)/totmass ! ! Correct secondary by gas gravity ! if (lcorrect_gasgravity_lstart) call initcond_correct_gasgravity(f,velocity,isecondary) ! velocity( iprimary,2) = velocity(isecondary,2) * pmass(isecondary)/pmass(iprimary) ! ! Revert all velocities if retrograde. ! ! if (lretrograde) velocity=-velocity ! ! Loop through particles to allocate the velocities. ! do k=1,nqpar fq(k,ivxq:ivzq) = velocity(k,1:3) enddo ! case default if (lroot) print*, 'init_pointmasses: No such such value for'//& 'initvvq: ',trim(initvvq) call fatal_error('init_pointmasses','') ! endselect ! ! Make the particles known to all processors ! call boundconds_pointmasses()!fq)!,ipar) call mpibcast_real(fq,(/nqpar,mqarray/)) ! if (lroot) then print*,'init_pointmasses: particle velocities' do ks=1,nqpar print*,'particle',ks, 'velocities x,y,z:',fq(ks,ivxq:ivzq) enddo endif ! call keep_compiler_quiet(f) ! endsubroutine init_pointmasses !*********************************************************************** subroutine pointmasses_pde_pencil(f,df,p) ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p ! !call dxxpq_dt_pointmasses_pencil(f,df,fp,dfp,p,ineargrid) call dvvq_dt_pointmasses_pencil(f,df,p) ! endsubroutine pointmasses_pde_pencil !*********************************************************************** subroutine dvvq_dt_pointmasses_pencil(f,df,p) ! ! 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 use Mpicomm, only: mpireduce_sum ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p ! real, dimension (nx,nqpar) :: rp_mn, rpcyl_mn real, dimension (mx,3) :: ggt real, dimension (3) :: xxq,rpsecondary,accg real, dimension (nx) :: pot_energy,torque real :: accg_local integer :: ks,j,ju logical :: lintegrate, lparticle_out ! intent (in) :: f, p intent (inout) :: df ! ! Get the pre-calculated gravity field. ! lhydroif: if (lhydro) then call get_total_gravity(ggt) df(l1:l2,m,n,iux:iuz) = df(l1:l2,m,n,iux:iuz) + ggt(l1:l2,:) ! if (lgas_removal.or.lmomentum_removal) call gas_accretion_by_pointmass(df,p) ! ! Add the gas gravity to the pointmasses by integration if llive_secondary ! is true. ! ! The integration is needed for these two cases: ! ! 1. We are not solving the Poisson equation ! 2. We are, but a particle is out of the box (a star, for instance) ! and therefore the potential cannot be interpolated. ! if (llive_secondary) then pointmasses1: do ks=1,nqpar lparticle_out=.false. ! ! Check if the particle is out of the box. ! if ((fq(ks,ixq)< xyz0(1)).or.(fq(ks,ixq) > xyz1(1)) .or. & (fq(ks,iyq)< xyz0(2)).or.(fq(ks,iyq) > xyz1(2)) .or. & (fq(ks,izq)< xyz0(3)).or.(fq(ks,izq) > xyz1(3))) then !particle out of box lparticle_out=.true. endif ! ! A live particle needs to be integrated if selfgravity is not being used (poisson ! equation not being solved) or if it is not in the grid (in which case the selfgravity ! cannot be interpolated. ! lintegrate=(.not.lselfgravity).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==iprimary).and.lnoselfgrav_primary) lintegrate=.false. ! integrategas: if (lintegrate) then ! ! Get the acceleration particle ks suffers due to self-gravity. ! call get_radial_distance(rp_mn(:,ks),rpcyl_mn(:,ks),& E1_=fq(ks,ixq),E2_=fq(ks,iyq),E3_=fq(ks,izq)) xxq = fq(ks,ixq:izq) ! if (lcylindrical_gravity_nbody(ks)) then call integrate_gasgravity(p,rpcyl_mn(:,ks),& xxq,accg,r_smooth(ks)) else call integrate_gasgravity(p,rp_mn(:,ks),& xxq,accg,r_smooth(ks)) endif ! ! Add it to its dfp ! dfq(ks,ivxq:ivzq) = dfq(ks,ivxq:ivzq) + accg(1:3) ! ! Calculate torques for output, if needed ! if (ldiagnos.and.idiag_torque(ks)/=0) then rpsecondary=(/fq(ks,ixq)*cos(fq(ks,iyq)-y(m)),& fq(ks,ixq)*sin(fq(ks,iyq)-y(m)),& fq(ks,izq) & /) ! ! Integrate will add the cell volume, so we first remove it from the ! calculation. ! !torque=(rpsecondary(1)*accg_mn(:,2)-rpsecondary(2)*accg_mn(:,1))*& ! dVol1_x(l1:l2)*dVol1_y(m)*dVol1_z(n) call cross(rpsecondary,accg,torque) call integrate_mn_name(pmass(ks)*torque,idiag_torque(ks)) endif ! endif integrategas enddo pointmasses1 endif ! ! Diagnostic ! diagnos: if (ldiagnos) then pointmasses2: do ks=1,nqpar ! if (idiag_totenergy/=0.or.& idiag_torqext(ks)/=0.or.& idiag_torqint(ks)/=0) & call get_radial_distance(rp_mn(:,ks),rpcyl_mn(:,ks),& E1_=fq(ks,ixq),E2_=fq(ks,iyq),E3_=fq(ks,izq)) ! ! Total energy ! if (idiag_totenergy/=0) then pot_energy=0.0 !potential energy pot_energy = pot_energy - & GNewton*pmass(ks)*(rpcyl_mn(:,ks)**2+r_smooth(ks)**2)**(-0.5) if (ks==nqpar) then if (lcartesian_coords) then call sum_lim_mn_name(.5*p%rho*p%u2 + pot_energy,idiag_totenergy,p) else call integrate_mn_name(.5*p%rho*p%u2 + pot_energy,idiag_totenergy) endif endif endif ! ! Calculate torques splitting inner and outer, for backward compatibility ! if ((idiag_torqext(ks)/=0).or.(idiag_torqint(ks)/=0)) & call calc_torque_split_int_ext(p,rpcyl_mn(:,ks),ks) ! enddo pointmasses2 endif diagnos endif lhydroif ! call keep_compiler_quiet(f) ! endsubroutine dvvq_dt_pointmasses_pencil !*********************************************************************** subroutine gas_accretion_by_pointmass(df,p) ! ! Compute gas accretion rate, and remove it from the gas. ! ! Add it accordingly to Eq 6 of Duffell et al. 2020 (ApJ, 901, 25) ! ! 01-aug-22/wlad: coded ! use Sub, only: get_radial_distance use Diagnostics ! real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p real, dimension (nx,3) :: momentum_removal_rate real, dimension (nx) :: prefactor,mass_removal_rate,argexp real, dimension (nx) :: rp_sph,rp_cyl integer :: j,ks ! intent (in) :: p intent (inout) :: df ! do ks=1,nqpar call get_radial_distance(rp_sph,rp_cyl,& e1_=fq(ks,ixq),e2_=fq(ks,iyq),e3_=fq(ks,izq)) argexp = min(rp_sph*r1_smooth(ks),3.0) prefactor = tau_accretion*exp(-argexp**4.0) ! if (lgas_removal) then mass_removal_rate = - prefactor * p%rho df(l1:l2,m,n,irho) = df(l1:l2,m,n,irho) + mass_removal_rate if (ldiagnos.and.idiag_mdot_pt/=0) & call sum_lim_mn_name(mass_removal_rate,idiag_mdot_pt,p) endif ! if (lmomentum_removal) then do j=1,3 momentum_removal_rate(:,j) = - prefactor * p%uu(:,j) enddo df(l1:l2,m,n,iux:iuz) = df(l1:l2,m,n,iux:iuz) + momentum_removal_rate endif ! enddo ! endsubroutine gas_accretion_by_pointmass !*********************************************************************** subroutine pointmasses_pde(f,df) ! use Mpicomm, only: mpibcast_real ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df ! real, dimension(nqpar) :: hill_radius_square logical :: ldiagnostic_only ! ! Add gas selfgravity if present ! if (lselfgravity) then ldiagnostic_only=ldiagnos.and.idiag_torque(isecondary)/=0 if (llive_secondary.or.ldiagnostic_only) & call selfgravity_gas_on_pointmass endif ! ! Add gravity. Both dvvq_dt_pointmasses and dvvp_dt_dustparticles call ! gravity_pointmasses, which call the Courant condition. ! if (lroot) then call calc_hill_radius(hill_radius_square) call dxxq_dt_pointmasses call dvvq_dt_pointmasses(hill_radius_square) endif ! if (lparticles) then call mpibcast_real(hill_radius_square,nqpar) call dvvp_dt_dustparticles(hill_radius_square) endif ! call keep_compiler_quiet(f,df) ! endsubroutine pointmasses_pde !*********************************************************************** subroutine calc_hill_radius(hill_radius_square) real, dimension(nqpar), intent(out) :: hill_radius_square integer :: ks real :: rr,w2,sma2 ! do ks=1,nqpar if (laccretion(ks).and.(ks/=iprimary)) then if (lcartesian_coords) then rr=sqrt(fq(ks,ixq)**2 + fq(ks,iyq)**2 + fq(ks,izq)**2) elseif (lcylindrical_coords) then rr=fq(ks,ixq) if (nzgrid/=1) rr=sqrt(fq(ks,ixq)**2+fq(ks,izq)**2) elseif (lspherical_coords) then rr=fq(ks,ixq) else call fatal_error('calc_hill_radius','wrong coord system') rr=0.0 endif ! ! particle velocities are non-coordinate (linear) ! w2 = fq(ks,ivxq)**2 + fq(ks,ivyq)**2 + fq(ks,ivzq)**2 ! ! squared semi major axis - assumes GM=1, so beware... ! if (totmass/=1) call fatal_error('calc_hill_radius',& 'can only calculate semimajor axis for normalized total mass') sma2 = (rr/(2-rr*w2))**2 ! ! squared hills radius ! hill_radius_square(ks)=sma2*(pmass(ks)*pmass1(iprimary)/3)**(2./3.) else hill_radius_square(ks)=0.0 endif enddo endsubroutine calc_hill_radius !*********************************************************************** subroutine dxxq_dt_pointmasses ! ! If the center of mass of the point masses was moved from the ! center of the grid, reset it. ! ! 22-sep-06/wlad: coded ! if (lreset_cm) call reset_center_of_mass ! endsubroutine dxxq_dt_pointmasses !*********************************************************************** subroutine dvvq_dt_pointmasses(hill_radius_square) ! ! Evolution of point masses 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 Sub ! real, dimension(nqpar) :: hill_radius_square integer :: k, ks, j, jpos, jvel logical :: lheader, lfirstcall=.true. ! ! Print out header information in first time step. ! lheader=lfirstcall .and. lroot ! ! Identify module. ! if (lheader) print*, 'dvvq_dt_pointmasses: Calculate dvvq_dt_pointmasses' ! ! 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_pointmasses_pencil. ! if (lramp) call get_ramped_mass ! ! To the N^2 problem of adding gravity of each massive particle to another. ! do k=1,nqpar call gravity_pointmasses(k,hill_radius_square,.true.) enddo ! ! Position and velocity diagnostics (per massive particle). ! if (ldiagnos) then do ks=1,nqpar if (lfollow_particle(ks)) then do j=1,3 jpos=j+ixq-1 ; jvel=j+ivxq-1 if (idiag_xxq(ks,j)/=0) then call point_par_name(fq(ks,jpos),idiag_xxq(ks,j)) endif if (idiag_vvq(ks,j)/=0) & call point_par_name(fq(ks,jvel),idiag_vvq(ks,j)) enddo if (idiag_period(ks)/=0) call point_par_name(2*pi*fq(ks,ixq)/fq(ks,ivyq),idiag_period(ks)) endif enddo endif ! if (lheader) print*,'dvvp_dt_pointmasses: Finished dvvp_dt_pointmasses' if (lfirstcall) lfirstcall=.false. ! endsubroutine dvvq_dt_pointmasses !************************************************************ subroutine dvvp_dt_dustparticles(hill_radius_square) use Particles_main, only: fetch_nparloc,fetch_fp_array,return_fp_array real, dimension (mpar_loc,mparray) :: fp_aux real, dimension (mpar_loc,mpvar) :: dfp_aux integer :: np_aux,k logical, dimension(mpar_loc) :: flag real, dimension(nqpar) :: hill_radius_square ! call fetch_nparloc(np_aux) if (np_aux/=0) then call fetch_fp_array(fp_aux,dfp_aux,& index%ixw,index%iyw,index%izw,& index%ivxw,index%ivyw,index%ivzw) do k=1,np_aux flag(k)=.false. call gravity_pointmasses(k,hill_radius_square,.false.,& fp_aux(k,:),dfp_aux(k,:),flag(k)) enddo call return_fp_array(fp_aux,dfp_aux,flag) endif ! endsubroutine dvvp_dt_dustparticles !************************************************************ subroutine gravity_pointmasses(k,hill_radius_square,& lcallpointmass,fp_pt,dfp_pt,flag_pt) ! ! 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 ! integer, intent(in) :: k real, dimension (nqpar) :: hill_radius_square ! real :: rr2, r2_ij, rs2, rsmooth2, Omega2_pm, rhill1, invr3_ij real :: dt1_nbody, r_ij, v_ij, a_ij integer :: ks ! real, dimension (3) :: evr_cart,evr,positions ! logical, intent(in) :: lcallpointmass ! real, dimension (mparray), optional :: fp_pt real, dimension (mpvar), optional :: dfp_pt logical, optional :: flag_pt ! if (lcallpointmass) then positions = fq(k,ixq:izq) else positions = fp_pt(index%ixw:index%izw) endif ! do ks=1,nqpar if (.not.(lcallpointmass.and.k==ks)) then !prevent self-acceleration ! ! r_ij = sqrt(ev1**2 + ev2**2 + ev3**2) ! call get_evr(positions,fq(ks,ixq:izq),evr_cart)!,.true.) rr2=sum(evr_cart**2) rsmooth2=r_smooth(ks)**2 ! ! Particles relative distance from each other: ! select case (ipotential_pointmass(ks)) ! case ('newton-hill','newton','newtonian') ! ! Newtonian potential. Should be used for a particle outside the grid (that does not need smoothing). ! r2_ij=max(rr2,rsmooth2) if (r2_ij > 0) then invr3_ij = r2_ij**(-1.5) else ! can happen during pencil_check invr3_ij = 0.0 endif Omega2_pm = GNewton*pmass(ks)*invr3_ij ! case ('plummer') ! ! Potential of a Plummer sphere. It's better to use Boley's (below) but here for testing purposes. ! r2_ij=rr2+rsmooth2 Omega2_pm = GNewton*pmass(ks)*r2_ij**(-1.5) ! case ('boley') ! ! Potential that Aaron Boley uses (should ask him for reference). Exactly matches Newtonian outside Hill radius. ! r2_ij=rr2 if (r2_ij .gt. hill_radius_square(ks)) then Omega2_pm = GNewton*pmass(ks)*r2_ij**(-1.5) else rhill1=1./sqrt(hill_radius_square(ks)) Omega2_pm = -GNewton*pmass(ks)*(3*sqrt(r2_ij)*rhill1 - 4)*rhill1**3 endif ! case default if (lroot) print*, 'gravity_pointmasses: '//& 'No such value for ipotential_pointmass: ', trim(ipotential_pointmass(ks)) call fatal_error("","") endselect ! ! If there is accretion, remove the accreted particles from the simulation, if any. ! if (.not.(lcallpointmass).and.laccretion(ks)) then rs2=(accrete_hills_frac(ks)**2)*hill_radius_square(ks) if (r2_ij<=rs2) then !flag particle for removal flag_pt=.true. return endif 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. ! if (lcallpointmass) then dfq_cart(k,:) = dfq_cart(k,:) - Omega2_pm*evr_cart(1:3) else dfp_pt(ivpx_cart:ivpz_cart) = & dfp_pt(ivpx_cart:ivpz_cart) - Omega2_pm*evr_cart(1:3) endif ! if (ladd_dragforce) call dragforce_pointmasses(k) ! if (lcoriolis_force) then dfq(k,ivxq) = dfq(k,ivxq) + 2*Omega_Coriolis*fq(k,ivyq) dfq(k,ivyq) = dfq(k,ivyq) - 2*Omega_Coriolis*fq(k,ivxq) 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_pointmasses) then if (.not.lcallpointmass) then v_ij = sqrt(sum((fp_pt(ivxq:ivzq)-fq(ks,ivxq:ivzq))**2)) else v_ij = sqrt(sum((fq(k,ivxq:ivzq)-fq(ks,ivxq:ivzq))**2)) endif a_ij = sqrt(sum(Omega2_pm*evr_cart**2)) r_ij = sqrt(r2_ij) dt1_nbody = max(v_ij/r_ij,sqrt(a_ij/r_ij)) dt1_max(l1-nghost) = max(dt1_max(l1-nghost),dt1_nbody/cdtq) endif ! endif !if (ipar(k)/=ks) ! enddo !nbody loop ! endsubroutine gravity_pointmasses !********************************************************** subroutine dragforce_pointmasses(k) ! ! Adds dragforce on massive particles. ! 06-feb-18/wlad: coded ! real, dimension (3) :: uup integer, intent(in) :: k ! ! Supports only Cartesian with ugas=uy so far. ! uup=(/0.,ugas,0./) if (llinear_drag) then dfq(k,ivxq:ivzq) = dfq(k,ivxq:ivzq) - (fq(k,ivxq:ivzq)-uup)/StokesNumber(k) else if (lquadratic_drag) then dfq(k,ivxq:ivzq) = dfq(k,ivxq:ivzq) - & abs(fq(k,ivxq:ivzq)-uup)*(fq(k,ivxq:ivzq)-uup)/StokesNumber(k) else call fatal_error("drag should be linear or quadratic","") endif ! endsubroutine dragforce_pointmasses !********************************************************** subroutine selfgravity_gas_on_pointmass() ! ! Adds selfgravity on massive particles. This subroutine ! fetchs the (global) acceleration array calculated by the ! logspirals method and interpolates the acceleration to ! the position of the pointmass. TODO: other poisson methods ! should get the potential, take the gradient and interpolate. ! ! 14-may-20/wlad: coded ! use Mpicomm use Sub, only: get_radial_distance,cross use Poisson, only:inverse_laplacian,get_acceleration ! integer :: ks real, dimension (nx,ny,nz,3) :: acceleration real, dimension (3) :: accg,torque real :: a_ij,r_ij ! call get_acceleration(acceleration) ! do ks=1,nqpar if (ks/=iprimary) then call bilinear_interpolate(acceleration,fq(ks,ixq:izq),accg) if (lroot.and.llive_secondary) & dfq(ks,ivxq:ivzq) = dfq(ks,ivxq:ivzq) + accg ! if (ldiagnos) call calc_torque(ks,accg) ! if (lfirst.and.ldt.and.ldt_pointmasses) then a_ij = sqrt(sum(accg**2)) r_ij = r_ref dt1_max(l1-nghost) = max(dt1_max(l1-nghost),sqrt(a_ij/r_ij)/cdtq) endif ! endif enddo ! endsubroutine selfgravity_gas_on_pointmass !********************************************************** subroutine get_evr(xxp,xxq,evr_output) ! ! Point-to-point vector distance, in different coordinate systems. ! Return always in Cartesian. ! ! 14-feb-14/wlad: coded ! real, dimension(3), intent(in) :: xxp,xxq real, dimension(3), intent(out) :: evr_output real :: x1,y1,x2,y2,z1,z2 real :: e1,e2,e3,e10,e20,e30 ! e1=xxp(1);e10=xxq(1) e2=xxp(2);e20=xxq(2) e3=xxp(3);e30=xxq(3) ! if (lcartesian_coords) then x1=e1 ; x2=e10 y1=e2 ; y2=e20 z1=e3 ; z2=e30 elseif (lcylindrical_coords) then ! x1=e1*cos(e2) y1=e1*sin(e2) z1=e3 ! x2=e10*cos(e20) y2=e10*sin(e20) z2=e30 ! elseif (lspherical_coords) then x1=e1*sin(e2)*cos(e3) y1=e1*sin(e2)*sin(e3) z1=e1*cos(e2) ! x2=e10*sin(e20)*cos(e30) y2=e10*sin(e20)*sin(e30) z2=e10*cos(e20) endif ! evr_output(1)=x1-x2 evr_output(2)=y1-y2 evr_output(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_pointmasses_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=pointmasses_init_pars, IOSTAT=iostat) ! endsubroutine read_pointmasses_init_pars !*********************************************************************** subroutine write_pointmasses_init_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=pointmasses_init_pars) ! endsubroutine write_pointmasses_init_pars !*********************************************************************** subroutine read_pointmasses_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=pointmasses_run_pars, IOSTAT=iostat) ! endsubroutine read_pointmasses_run_pars !*********************************************************************** subroutine write_pointmasses_run_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=pointmasses_run_pars) ! endsubroutine write_pointmasses_run_pars !*********************************************************************** subroutine reset_center_of_mass() ! ! 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(nqpar,mqarray) :: ftmp real, dimension(3) :: vcm real :: xcm,ycm,zcm,thtcm,phicm,vxcm,vycm,vzcm integer :: ks ! ftmp=fq(1:nqpar,:) ! if (lcartesian_coords) then vcm(1) = sum(ftmp(:,imass)*ftmp(:,ivxq)) vcm(2) = sum(ftmp(:,imass)*ftmp(:,ivyq)) vcm(3) = sum(ftmp(:,imass)*ftmp(:,ivzq)) else if (lcylindrical_coords) then xcm=sum(ftmp(:,imass)*(ftmp(:,ixq)*cos(ftmp(:,iyq)))) ycm=sum(ftmp(:,imass)*(ftmp(:,ixq)*sin(ftmp(:,iyq)))) phicm=atan2(ycm,xcm) ! vxcm=sum(ftmp(:,imass)*(& ftmp(:,ivxq)*cos(ftmp(:,iyq))-ftmp(:,ivyq)*sin(ftmp(:,iyq)))) vycm=sum(ftmp(:,imass)*(& ftmp(:,ivxq)*sin(ftmp(:,iyq))+ftmp(:,ivyq)*cos(ftmp(:,iyq)))) ! vcm(1)= vxcm*cos(phicm) + vycm*sin(phicm) vcm(2)=-vxcm*sin(phicm) + vycm*cos(phicm) vcm(3) = sum(ftmp(:,imass)*ftmp(:,ivzq)) ! else if (lspherical_coords) then vxcm=sum(ftmp(:,imass)*( & ftmp(:,ivxq)*sin(ftmp(:,iyq))*cos(ftmp(:,izq))& +ftmp(:,ivyq)*cos(ftmp(:,iyq))*cos(ftmp(:,izq))& -ftmp(:,ivzq)*sin(ftmp(:,izq)) )) vycm=sum(ftmp(:,imass)*( & ftmp(:,ivxq)*sin(ftmp(:,iyq))*sin(ftmp(:,izq))& +ftmp(:,ivyq)*cos(ftmp(:,iyq))*sin(ftmp(:,izq))& +ftmp(:,ivzq)*cos(ftmp(:,izq)) )) vzcm=sum(ftmp(:,imass)*(& ftmp(:,ivxq)*cos(ftmp(:,iyq))-ftmp(:,ivyq)*sin(ftmp(:,iyq)))) ! xcm=sum(ftmp(:,imass)*(ftmp(:,ixq)*sin(ftmp(:,iyq))*cos(ftmp(:,izq)))) ycm=sum(ftmp(:,imass)*(ftmp(:,ixq)*sin(ftmp(:,iyq))*sin(ftmp(:,izq)))) zcm=sum(ftmp(:,imass)*(ftmp(:,ixq)*cos(ftmp(:,iyq)))) ! 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 ks=1,nqpar dfq(ks,ixq:izq) = dfq(ks,ixq:izq) - vcm*totmass1 enddo ! endsubroutine reset_center_of_mass !*********************************************************************** 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_pointmasses(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=fq(ks,ixq);y0=fq(ks,iyq);z0=fq(ks,izq) ! 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_pointmasses !*********************************************************************** subroutine calc_torque(ks,accg) ! use Sub, only: cross ! integer, intent(in) :: ks real, dimension(3), intent(in) :: accg real, dimension(3) :: torque ! if (lroot) then if (idiag_torque(ks)/=0) then call cross((/fq(ks,ixq),0.,0./),accg,torque) call point_par_name(pmass(ks)*torque(3),idiag_torque(ks)) endif endif ! endsubroutine calc_torque !*********************************************************************** subroutine calc_torque_split_int_ext(p,dist,ks) ! ! Output torque diagnostic for nbody particle ks ! (maintained for backward compatibility). As of May 2020 ! the torque diagnostic should be computed after integrate_gasgravity, ! from the same acceleration used on the planet, instead of integrating ! separately. ! ! 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==iprimary) call fatal_error('calc_torque', & ! 'Nonsense to calculate torques for the star') ! if (lcartesian_coords) then rr = sqrt(fq(ks,ixq)**2 + fq(ks,iyq)**2 + fq(ks,izq)**2) rpre = fq(ks,ixq)*y(m) - fq(ks,iyq)*x(l1:l2) elseif (lcylindrical_coords) then rr_mn = x(l1:l2) ; phi = y(m) rr = fq(ks,ixq) ; phip = fq(ks,iyq) 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 = fq(ks,ivxq)**2 + fq(ks,ivyq)**2 + fq(ks,ivzq)**2 smap = 1./(2./rr - w2) hills = smap*(pmass(ks)*pmass1(iprimary)/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 (lexclude_hills) then pcut=hills_tempering_fraction*hills tempering = 1./(exp(-(dist/hills - pcut)/(.1*pcut))+1.) torque_gas = torque_gas * tempering torque_par = torque_par * tempering endif ! ! Separate internal and external torques ! do i=1,nx if (p%rcyl_mn(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_split_int_ext !*********************************************************************** 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,nqpar !_orig if (ks/=iprimary) then pmass(ks)= max(dble(tini),& final_ramped_mass(ks)*(sin((.5*pi)*(t/ramping_period))**2)) tmp=tmp+pmass(ks) endif enddo pmass(iprimary)= 1-tmp else pmass(1:nqpar)=final_ramped_mass(1:nqpar) endif ! endsubroutine get_ramped_mass !*********************************************************************** subroutine get_total_gravity(ggt) ! ! Sum the gravities of all massive particles ! ! 08-mar-08/wlad: coded ! use Sub ! real, dimension (mx,nqpar) :: rp_mn,rpcyl_mn real, dimension (mx,3) :: ggp,ggt ! real, dimension (mx) :: Omega2_pm,rrp real :: rr,rp1,rhill,rhill1 integer :: ks,i !real, dimension (mx) :: grav_particle,rrp !integer :: ks ! intent(out) :: ggt ! ggt=0. do ks=1,nqpar ! ! Hill radius ! if (lcartesian_coords) then rp1 = sqrt(fq(ks,ixq)**2+fq(ks,iyq)**2+fq(ks,izq)**2) elseif (lcylindrical_coords) then rp1 = sqrt(fq(ks,ixq)**2+ fq(ks,izq)**2) elseif (lspherical_coords) then rp1 = fq(ks,ixq) endif ! rhill = rp1*(GNewton*pmass(ks)/3.)**(1./3) rhill1 = 1./rhill ! ! Spherical and cylindrical distances ! call get_radial_distance(rp_mn(:,ks),rpcyl_mn(:,ks),& e1_=fq(ks,ixq),e2_=fq(ks,iyq),e3_=fq(ks,izq)) ! ! 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 ! select case (ipotential_pointmass(ks)) ! case ('plummer') ! ! Potential of a Plummer sphere ! if (ks==iprimary .and. headtt) call warning("get_total_gravity",& "The primary is not newtonian, make sure you know what you are doing.") Omega2_pm =-GNewton*pmass(ks)*(rrp**2+r_smooth(ks)**2)**(-1.5) ! case ('boley') ! ! Correct potential outside Hill sphere ! if (ks==iprimary .and. headtt) call warning("get_total_gravity",& "The primary is not newtonian, make sure you know what you are doing.") do i=1,mx if (rrp(i) .gt. rhill) then Omega2_pm(i) = -GNewton*pmass(ks)*rrp(i)**(-3) else Omega2_pm(i) = GNewton*pmass(ks)*(3*rrp(i)*rhill1 - 4)*rhill1**3 endif enddo ! case ('newton-hill','newton','newtonian') ! ! Newtonian potential; same as boley but constant inside rsmooth ! if (ks==iprimary.and.r_smooth(ks)/=0) call fatal_error("get_total_gravity",& "Use r_smooth=0 for the primary's potential") ! do i=1,mx rr=max(rrp(i),r_smooth(ks)) if (rr > 0) then Omega2_pm(i) = -GNewton*pmass(ks)*rr**(-3) else ! can happen during pencil_check Omega2_pm(i) = 0. endif enddo ! case default ! ! Catch unknown values ! if (lroot) print*, 'get_total_gravity: '//& 'No such value for ipotential_pointmass: ', trim(ipotential_pointmass(ks)) call fatal_error("","") endselect ! call get_gravity_field_pointmasses(Omega2_pm,ggp,ks) ! if ((ks==iprimary).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 ! enddo ! endsubroutine get_total_gravity !*********************************************************************** subroutine initcond_correct_gasgravity(f,velocity,k) ! ! Calculates acceleration on the point (x,y,z)=xxpar ! due to the gravity of gas. ! ! 14-may-20/wlad : coded ! real, dimension (mx,my,mz,mfarray) :: f real, dimension(nqpar,3), intent(inout) :: velocity real, dimension(3) :: accg real :: vphi,phidot2,OO,rr integer, intent(in) :: k ! if (lselfgravity) then call correct_gasgravity_selfgravity(f,k,accg) else call correct_gasgravity_integrate(f,k,accg) endif ! ! Correct original velocity by this acceleration ! ! phidot**2 = OmegaK**2 + 1/r (d/dr PHI_sg) ! phidot**2 = OmegaK**2 - 1/r gr ! ! Current azimuthal velocity and azimuthal frequency ! if (lcylindrical_coords) then vphi = velocity(k,2) rr = fq(k,ixq) elseif (lspherical_coords) then vphi = velocity(k,3) rr = fq(k,ixq)*sin(fq(k,iyq)) endif OO = vphi/rr ! ! Update angular velocity by selfgravitational acceleration ! This line assumes axisymmetry ! phidot2 = OO**2 - accg(1)/rr ! ! Corrected velocity ! vphi = sqrt(phidot2) * rr ! if (lcylindrical_coords) then velocity(k,2) = vphi elseif (lspherical_coords) then velocity(k,3) = vphi endif ! endsubroutine initcond_correct_gasgravity !*********************************************************************** subroutine correct_gasgravity_selfgravity(f,k,accg) ! ! Calculates acceleration on the point (x,y,z)=xxpar ! due to the gravity of gas, by interpolating from the ! gravitational acceleration from the selfgravity module. ! ! 14-may-20/wlad : coded ! use Mpicomm use Sub, only: get_radial_distance use Poisson, only:inverse_laplacian,get_acceleration ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (nx,ny,nz,3) :: acceleration real, dimension (nx,ny,nz) :: rho real, dimension (3) :: accg integer, intent (in) :: k ! if (ldensity_nolog) then rho=f(l1:l2,m1:m2,n1:n2,irho) else rho=exp(f(l1:l2,m1:m2,n1:n2,ilnrho)) endif call inverse_laplacian(rho) call get_acceleration(acceleration) call bilinear_interpolate(acceleration,fq(k,ixq:izq),accg) ! endsubroutine correct_gasgravity_selfgravity !*********************************************************************** subroutine correct_gasgravity_integrate(f,k,accg) ! ! Calculates acceleration on the point (x,y,z)=xxpar ! due to the gravity of gas, by integrating the whole grid. ! ! 29-aug-18/wlad : coded ! use Mpicomm use Sub, only: get_radial_distance ! implicit none real, dimension (mx,my,mz,mfarray) :: f real, dimension(3) :: xxpar,sum_loc,accg real, dimension(nx,3) :: dist real, dimension(nx) :: rrp,rp_mn,rpcyl_mn,gasgravity,density,cellmass real, dimension(nx) :: dv,tmp real :: vphi,phidot2,OO,rr,rp0,fac integer :: j,k ! ! Sanity check ! !if (.not.(lgas_gravity.or.ldust_gravity)) & ! call fatal_error("initcond_correct_selfgravity",& ! "No gas gravity or dust gravity to add. "//& ! "Switch on lgas_gravity or ldust_gravity in n-body parameters") ! xxpar = fq(k,ixq:izq) rp0=r_smooth(k) ! sum_loc=0. ! mloop: do m=m1,m2 nloop: do n=n1,n2 if (coord_system=='cartesian') then dist(:,1)=x(l1:l2)-xxpar(1) dist(:,2)=y( m )-xxpar(2) dist(:,3)=z( n )-xxpar(3) elseif (coord_system=='cylindric') then 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('correct_gasgravity_integrate', & ' not yet implemented for spherical polars') else call fatal_error('correct_gasgravity_integrate','wrong coord_system') endif ! dv=1 if (nxgrid/=1) dv=dv/dline_1(:,1) if (nygrid/=1) dv=dv/dline_1(:,2) if (nzgrid/=1) dv=dv/dline_1(:,3) ! ! 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 (ldensity_nolog) then density=f(l1:l2,m,n,irho) else density=exp(f(l1:l2,m,n,ilnrho)) endif ! ! Add the particle gravity if npar>mspar (which means dust is being used) ! !if (ldust.and.ldust_gravity) density=density+f(l1:l2,m,n,irhop) ! call get_radial_distance(rp_mn,rpcyl_mn,E1_=xxpar(1),E2_=xxpar(2),E3_=xxpar(3)) if (lcylindrical_gravity_nbody(k)) then rrp=rpcyl_mn else rrp=rp_mn endif ! cellmass = density_scale*density*dv gasgravity = GNewton*cellmass*(rrp**2 + rp0**2)**(-1.5) ! ! Exclude the frozen zones ! if (lcartesian_coords) call fatal_error("","") ! ! Integrate the accelerations on this processor ! And sum over processors with mpireduce ! do j=1,3 tmp=gasgravity*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) = sum_loc(j)+fac*sum(tmp) else sum_loc(j) = sum_loc(j)+fac*(sum(tmp(2:nx-1))+.5*(tmp(1)+tmp(nx))) endif enddo enddo nloop enddo mloop ! do j=1,3 call mpireduce_sum(sum_loc(j),accg(j)) enddo ! ! Broadcast particle acceleration ! call mpibcast_real(accg,3) ! endsubroutine correct_gasgravity_integrate !*********************************************************************** subroutine bilinear_interpolate(v,q,vp) ! ! Interpolate the gravity vector with bilinear interpolation to the position ! of the pointmass(es). ! use Mpicomm, only: mpisend_real,mpirecv_real,mpibcast_real ! real, dimension (nx,ny,nz,3) :: v real, dimension (3) :: q real, dimension (3) :: vp integer :: j real :: rdx1,rdy1,rdz1 real :: rdx2,rdy2,rdz2 integer :: ix0_global,iy0_global,iz0_global integer :: ix1_global,iy1_global,iz1_global integer :: ipx0,ipy0,ipz0 integer :: ipx1,ipy1,ipz1 integer :: ix0,iy0,iz0 integer :: ix1,iy1,iz1 ! integer :: iproc_q11,iproc_q21,iproc_q12,iproc_q22 real, dimension(3) :: Q11,Q21,Q12,Q22 ! intent(in) :: v,q intent(out) :: vp ! call calc_indices_and_processors(q(1),xgrid,dx1grid,nprocx,nx,& lperi(1),ix0,ix1,ipx0,ipx1,rdx1,rdx2) call calc_indices_and_processors(q(2),ygrid,dy1grid,nprocy,ny,& lperi(2),iy0,iy1,ipy0,ipy1,rdy1,rdy2) ! if (nz==1) then iz0=1 else call fatal_error("trilinear_interpolate","logpsirals works only for 2D") !call calc_indices_and_processors(q(3),zgrid,dz1grid,nprocz,nz,& !lperi(3),iz0,iz1,ipz0,ipz1,rdz1,rdz2) endif ! iproc_q11 = ipx0 + nprocx*ipy0 + nprocx*nprocy*ipz iproc_q21 = ipx1 + nprocx*ipy0 + nprocx*nprocy*ipz iproc_q12 = ipx0 + nprocx*ipy1 + nprocx*nprocy*ipz iproc_q22 = ipx1 + nprocx*ipy1 + nprocx*nprocy*ipz ! if (iproc_q11 /= root) then if (iproc==iproc_q11) call mpisend_real(v(ix0,iy0,iz0,:),3,root ,111) if (lroot) call mpirecv_real(Q11 ,3,iproc_q11,111) else Q11 = v(ix0,iy0,iz0,:) endif ! if (iproc_q21 /= root) then if (iproc==iproc_q21) call mpisend_real(v(ix1,iy0,iz0,:),3,root ,121) if (lroot) call mpirecv_real(Q21 ,3,iproc_q21,121) else Q21 = v(ix1,iy0,iz0,:) endif ! if (iproc_q12 /= root) then if (iproc==iproc_q12) call mpisend_real(v(ix0,iy1,iz0,:),3,root ,112) if (lroot) call mpirecv_real(Q12 ,3,iproc_q12,112) else Q12 = v(ix0,iy1,iz0,:) endif ! if (iproc_q22 /= root) then if (iproc==iproc_q22) call mpisend_real(v(ix1,iy1,iz0,:),3,root ,122) if (lroot) call mpirecv_real(Q22 ,3,iproc_q22,122) else Q22 = v(ix1,iy1,iz0,:) endif ! if (lroot) then do j=1,3 ! ! Interpolation formula ! vp(j) & = Q11(j) * rdx2* rdy2 & ! Q11 + Q21(j) * rdx1* rdy2 & ! Q21 + Q12(j) * rdx2* rdy1 & ! Q12 + Q22(j) * rdx1* rdy1 ! Q22 enddo endif ! call mpibcast_real(vp,3) ! endsubroutine bilinear_interpolate !*********************************************************************** subroutine calc_indices_and_processors(q,qgrid,dq1grid,nprocq,nq,lperiodic_direction,& iq0,iq1,ipq0,ipq1,rdq1,rdq2) ! use Sub, only: find_index_by_bisection ! real, dimension(:) :: qgrid real, dimension(:) :: dq1grid real :: q,rdq1,rdq2,dq1,dq integer :: iq0_global,iq1_global,ipq0,ipq1,iq0,iq1,nprocq,nq logical :: lperiodic_direction integer :: nqgrid ! intent(in) :: q,qgrid,lperiodic_direction,nprocq,nq intent(out) :: iq0,iq1,ipq0,ipq1,rdq1,rdq2 lperiodic: if (.not.lperiodic_direction) then ! ! Non-periodic grid ! call find_index_by_bisection(q,qgrid,iq0_global) if (qgrid(iq0_global)>q) iq0_global = iq0_global-1 iq1_global=iq0_global+1 ! dq1=1./(qgrid(iq1_global)-qgrid(iq0_global)) rdq1 = (q-qgrid(iq0_global))*dq1 rdq2 = (qgrid(iq1_global)-q)*dq1 ! if (nprocq/=1) then ipq0 = (iq0_global-1)/nq ipq1 = (iq1_global-1)/nq ! iq0 = iq0_global-ipq0*nq iq1 = iq1_global-ipq1*nq else ipq0=0 ipq1=0 ! iq0=iq0_global iq1=iq1_global endif ! else ! ! Periodic grid ! call find_index_by_bisection(q,qgrid,iq0_global) dq1=dq1grid(iq0_global) dq=1./dq1 if (qgrid(iq0_global)>q) iq0_global = iq0_global-1 iq1_global=iq0_global+1 ! nqgrid=size(qgrid) if (iq0_global < 1) then iq0_global = iq0_global + nqgrid rdq1 = (q-(qgrid(1)-dq))*dq1 else rdq1 = (q-qgrid(iq0_global))*dq1 endif ! if (iq1_global > nqgrid) then iq1_global = iq1_global - nqgrid rdq2 = (qgrid(nqgrid)+dq-q)*dq1 else rdq2 = (qgrid(iq1_global)-q)*dq1 endif ! if (nprocq/=1) then ipq0=(iq0_global-1)/nq ipq1=(iq1_global-1)/nq ! iq0=iq0_global-ipq0*nq iq1=iq1_global-ipq1*nq else ipq0=0 ipq1=0 ! iq0=iq0_global iq1=iq1_global endif ! endif lperiodic ! endsubroutine calc_indices_and_processors !*********************************************************************** subroutine integrate_gasgravity(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,rr,gasgravity,density,cellmass real, dimension(nx) :: dv,tmp real :: rp0,fac real, dimension(3) :: xxpar,accg,sum_loc integer :: j type (pencil_case) :: p logical :: lfirstcall=.true. ! intent(out) :: accg ! if (lfirstcall.and.lroot) & print*,'Adding gas+dust gravity to the massive particles' ! ! Sanity check ! if (.not.(lgas_gravity.or.ldust_gravity)) & call fatal_error("lintegrate_gasgravity",& "No gas gravity or dust gravity to add. "//& "Switch on lgas_gravity or ldust_gravity in n-body parameters") ! if (coord_system=='cartesian') then dist(:,1)=x(l1:l2)-xxpar(1) dist(:,2)=y( m )-xxpar(2) dist(:,3)=z( n )-xxpar(3) elseif (coord_system=='cylindric') then !re = r-rp !rp = rp*cos(phip-phi) \hat{r} + rp*sin(phi-phi) \hat{phi} dist(:,1)=x(l1:l2)-xxpar(1)*cos(y(m)-xxpar(2)) dist(:,2)= xxpar(1)*sin(y(m)-xxpar(2)) dist(:,3)=z( n )-xxpar(3) elseif (coord_system=='spherical') then call fatal_error('integrate_gasgravity', & ' not yet implemented for spherical polars') else call fatal_error('integrate_gasgravity','wrong coord_system') endif ! dv=1 if (nxgrid/=1) dv=dv/dline_1(:,1) if (nygrid/=1) dv=dv/dline_1(:,2) if (nzgrid/=1) dv=dv/dline_1(:,3) ! ! 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 ! cellmass = density_scale*density*dv gasgravity = GNewton*cellmass*(rrp**2 + rp0**2)**(-1.5) ! ! 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)) gasgravity = 0 endwhere else if (l2D.or.(l3D.and.lcylindrical_gravity)) then rr=p%rcyl_mn else rr=p%r_mn endif where ((rr<=r_int).or.(rr>=r_ext)) gasgravity = 0 endwhere endif endif ! ! Integrate the accelerations on this processor ! And sum over processors with mpireduce ! do j=1,3 tmp=gasgravity*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_gasgravity !*********************************************************************** subroutine pointmasses_read_snapshot(filename) ! ! Read nbody particle info ! ! 01-apr-08/wlad: dummy ! use IO, only: input_pointmass ! character (len=*), intent(in) :: filename ! call input_pointmass(filename, qvarname, fq, nqpar, mqarray) ! if (ldebug) then print*,'pointmasses_read_snapshot' print*,'nqpar=',nqpar print*,'fq =',fq print*,'' endif ! endsubroutine pointmasses_read_snapshot !*********************************************************************** subroutine pointmasses_write_snapshot(file,enum,flist) ! use General, only:safe_character_assign use Sub, only: update_snaptime, read_snaptime use IO, only: log_filename_to_file, output_pointmass ! ! Input and output of information about the massive particles ! ! 01-apr-08/wlad: coded ! character (len=*), intent(in) :: file logical, intent(in) :: enum character (len=*), optional, intent(in) :: flist ! logical, save :: lfirst_call=.true. integer, save :: nsnap real, save :: tsnap character (len=fnlen) :: filename, filename_diag logical :: lsnap character (len=intlen) :: nsnap_ch ! ! Input is either file=qvar if enum=F or file=QVAR if enum=T. ! If the latter, append a number to QVAR. ! filename = trim(file) if (enum) then ! ! Prepare to read tsnap.dat ! call safe_character_assign(filename_diag,trim(datadir)//'/tsnap.dat') ! ! Read the tsnap.dat to get N for QVARN (N=nsnap). ! if (lfirst_call) then call read_snaptime(filename_diag,tsnap,nsnap,dsnap,t) lfirst_call=.false. endif ! ! Update snaptime to set nsnap_ch=N and append N to filename QVARN. ! call update_snaptime(filename_diag,tsnap,nsnap,dsnap,t,lsnap,nsnap_ch,nowrite=.true.) call safe_character_assign(filename,trim(filename)//trim(nsnap_ch)) endif ! if (lsnap .or. .not. enum ) then call output_pointmass (filename, qvarname, fq, nqpar, mqarray) if (lroot .and. present(flist)) call log_filename_to_file(filename,flist) endif ! endsubroutine pointmasses_write_snapshot !*********************************************************************** subroutine pointmasses_write_qdim(filename) ! ! Write nqpar and mqvar to file. ! ! 01-apr-08/wlad: coded ! character (len=*) :: filename ! open(1,file=filename) write(1,'(3i9)') nqpar, mqvar close(1) ! endsubroutine pointmasses_write_qdim !*********************************************************************** subroutine rprint_pointmasses(lreset,lwrite) ! ! Read and register print parameters relevant for nbody particles. ! ! 17-nov-05/anders+wlad: adapted ! use Diagnostics use HDF5_IO, only: pointmass_index_append use General, only: itoa ! logical :: lreset logical, optional :: lwrite ! integer :: iname,ks,j character :: str character (len=intlen) :: sks logical :: lwr ! lwr = .false. if (present(lwrite)) lwr=lwrite ! ! Reset everything in case of reset ! if (lreset) then idiag_xxq=0;idiag_vvq=0 idiag_torqint=0;idiag_torqext=0 idiag_totenergy=0 idiag_torqext_gas=0;idiag_torqext_par=0 idiag_torqint_gas=0;idiag_torqint_par=0 idiag_period=0 idiag_torque=0 idiag_mdot_pt=0 endif ! ! Run through all possible names that may be listed in print.in ! if (lroot .and. ip<14) print*,'rprint_pointmasses: run through parse list' ! ! Now check diagnostics for specific particles ! do ks=1,nqpar 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)//'q'//trim(sks),idiag_xxq(ks,j)) call parse_name(iname,cname(iname),cform(iname),& 'v'//trim(str)//'q'//trim(sks),idiag_vvq(ks,j)) enddo ! ! Run through parse list again ! if (lwr) then call pointmass_index_append('i_'//trim(str)//'q'//trim(sks),idiag_xxq(ks,j)) call pointmass_index_append('i_v'//trim(str)//'q'//trim(sks),idiag_vvq(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)) call parse_name(iname,cname(iname),cform(iname),& 'period'//trim(sks),idiag_period(ks)) call parse_name(iname,cname(iname),cform(iname),& 'torque_'//trim(sks),idiag_torque(ks)) enddo ! if (lwr) then call pointmass_index_append('i_torqint_'//trim(sks),idiag_torqint(ks)) call pointmass_index_append('i_torqext_'//trim(sks),idiag_torqext(ks)) call pointmass_index_append('i_torqint_gas'//trim(sks),idiag_torqint(ks)) call pointmass_index_append('i_torqext_gas'//trim(sks),idiag_torqext(ks)) call pointmass_index_append('i_torqint_par'//trim(sks),idiag_torqint(ks)) call pointmass_index_append('i_torqext_par'//trim(sks),idiag_torqext(ks)) call pointmass_index_append('i_period'//trim(sks),idiag_period(ks)) call pointmass_index_append('i_torque_'//trim(sks),idiag_torque(ks)) endif enddo ! ! Diagnostic related to quantities summed over all point masses ! do iname=1,nname call parse_name(iname,cname(iname),cform(iname),'totenergy',idiag_totenergy) call parse_name(iname,cname(iname),cform(iname),'mdot_pt',idiag_mdot_pt) enddo ! if (lwr) then call pointmass_index_append('i_totenergy',idiag_totenergy) call pointmass_index_append('i_mdot_pt',idiag_mdot_pt) endif ! endsubroutine rprint_pointmasses !*********************************************************************** subroutine boundconds_pointmasses ! ! Global boundary conditions for particles. ! ! 30-dec-04/anders: coded ! use Mpicomm use General, only: random_number_wrapper use SharedVariables, only: get_shared_variable ! integer :: k character (len=2*bclen+1) :: boundx, boundy, boundz ! do k=1,nqpar ! boundx=bcqx; boundy=bcqy; boundz=bcqz ! ! Cartesian boundaries: Boundary condition in the x-direction. The physical ! domain is in the interval ! ! x \in [x0,x1[ ! y \in [y0,y1[ ! z \in [z0,z1[ ! if (nxgrid/=1) then if (boundx=='p') then ! xp < x0 if (fq(k,ixq)< xyz0(1)) then fq(k,ixq)=fq(k,ixq)+Lxyz(1) ! ! Particle position must never need more than one addition of Lx to get back ! in the box. Often a NaN or Inf in the particle position will show up as a ! problem here. if (fq(k,ixq)< xyz0(1)) then print*, 'boundconds_pointmasses: ERROR - particle ', k, & ' was further than Lx outside the simulation box!' print*, 'This must never happen.' print*, 'xxq=', fq(k,ixq:izq) call fatal_error_local('boundconds_pointmasses','') endif endif ! xp > x1 if (fq(k,ixq)>=xyz1(1)) then fq(k,ixq)=fq(k,ixq)-Lxyz(1) ! if (fq(k,ixq)>=xyz1(1)) then print*, 'boundconds_pointmasses: ERROR - particle ', k, & ' was further than Lx outside the simulation box!' print*, 'This must never happen.' print*, 'xxq=', fq(k,ixq:izq) call fatal_error_local('boundconds_pointmasses','') endif endif elseif (boundx=='out') then ! ! Do nothing. A massive particle, can be out of the box. A star, for example, ! in a cylindrical simulation ! else print*, 'boundconds_pointmasses: No such boundary condition =', boundx call stop_it('boundconds_pointmasses') endif endif ! ! Boundary condition in the y-direction. ! if (nygrid/=1) then if (boundy=='p') then ! yp < y0 if (fq(k,iyq)< xyz0(2)) then fq(k,iyq)=fq(k,iyq)+Lxyz(2) if (fq(k,iyq)< xyz0(2)) then print*, 'boundconds_pointmasses: ERROR - particle ', k, & ' was further than Ly outside the simulation box!' print*, 'This must never happen.' print*, 'xxq=', fq(k,ixq:izq) call fatal_error_local('boundconds_pointmasses','') endif endif ! yp > y1 if (fq(k,iyq)>=xyz1(2)) then fq(k,iyq)=fq(k,iyq)-Lxyz(2) if (fq(k,iyq)>=xyz1(2)) then print*, 'boundconds_pointmasses: ERROR - particle ', k, & ' was further than Ly outside the simulation box!' print*, 'This must never happen.' print*, 'xxq=', fq(k,ixq:izq) call fatal_error_local('boundconds_pointmasses','') endif endif elseif (boundy=='p2pi') then ! yp < y0 if (fq(k,iyq)< -pi) then fq(k,iyq)=fq(k,iyq)+2*pi if (fq(k,iyq)<-pi) then print*, 'boundconds_pointmasses: ERROR - particle ', k, & ' was further than Ly outside the simulation box!' print*, 'This must never happen.' print*, 'xxq=', fq(k,ixq:izq) call fatal_error_local('boundconds_pointmasses','') endif endif ! yp > y1 if (fq(k,iyq)>=pi) then fq(k,iyq)=fq(k,iyq)-2*pi if (fq(k,iyq)>=pi) then print*, 'boundconds_pointmasses: ERROR - particle ', k, & ' was further than Ly outside the simulation box!' print*, 'This must never happen.' print*, 'xxq=', fq(k,ixq:izq) call fatal_error_local('boundconds_pointmasses','') endif endif elseif (boundy=='out') then ! massive particles can be out of the box ! the star, for example, in a cylindrical simulation else print*, 'boundconds_pointmasses: No such boundary condition =', boundy call stop_it('boundconds_pointmasses') endif endif ! ! Boundary condition in the z-direction. ! if (nzgrid/=1) then if (boundz=='p') then ! zp < z0 if (fq(k,izq)< xyz0(3)) then fq(k,izq)=fq(k,izq)+Lxyz(3) if (fq(k,izq)< xyz0(3)) then print*, 'boundconds_pointmasses: ERROR - particle ', k, & ' was further than Lz outside the simulation box!' print*, 'This must never happen.' print*, 'xxq=', fq(k,ixq:izq) call fatal_error_local('boundconds_pointmasses','') endif endif ! zp > z1 if (fq(k,izq)>=xyz1(3)) then fq(k,izq)=fq(k,izq)-Lxyz(3) if (fq(k,izq)>=xyz1(3)) then print*, 'boundconds_pointmasses: ERROR - particle ', k, & ' was further than Lz outside the simulation box!' print*, 'This must never happen.' print*, 'xxq=', fq(k,ixq:izq) call fatal_error_local('boundconds_pointmasses','') endif endif elseif (boundz=='out') then ! massive particles can be out of the box ! the star, for example, in a cylindrical simulation else print*, 'boundconds_pointmasses: No such boundary condition=', boundz call stop_it('boundconds_pointmasses') endif endif enddo ! endsubroutine boundconds_pointmasses !*********************************************************************** subroutine pointmasses_timestep_first(f) ! use Particles_main, only: fetch_nparloc ! real, dimension (mx,my,mz,mfarray) :: f integer :: np_aux ! if (lroot) then if (itsub==1) then dfq(1:nqpar,:) = 0.0 dfq_cart(1:nqpar,:)=0.0 else dfq(1:nqpar,:)=alpha_ts(itsub)*dfq(1:nqpar,:) dfq_cart(1:nqpar,:)=alpha_ts(itsub)*dfq_cart(1:nqpar,:) endif endif ! call keep_compiler_quiet(f) ! endsubroutine pointmasses_timestep_first !*********************************************************************** subroutine pointmasses_timestep_second(f) ! ! Time evolution of particle variables. ! ! 07-jan-05/anders: coded ! use Mpicomm, only: mpibcast_real use Particles_main, only: fetch_nparloc,fetch_fp_array,return_fp_array ! real, dimension (mx,my,mz,mfarray) :: f real, dimension(3) :: xx_polar,vv_polar,xxdot_polar,vvdot_polar,aa_nbody_cart integer :: ks ! real, dimension (mpar_loc,mparray) :: fp_aux real, dimension (mpar_loc,mpvar) :: dfp_aux integer :: np_aux,k ! if (lroot) then do ks=1,nqpar xx_polar = fq(ks,ixq:izq) vv_polar = fq(ks,ivxq:ivzq) xxdot_polar = dfq(ks,ixq:izq) vvdot_polar = dfq(ks,ivxq:ivzq) ! aa_nbody_cart=dfq_cart(ks,:) call advance_particles_in_cartesian(xx_polar,vv_polar,xxdot_polar,vvdot_polar,aa_nbody_cart) fq(ks,ixq:izq)=xx_polar fq(ks,ivxq:ivzq)=vv_polar dfq(ks,ixq:izq)=xxdot_polar enddo endif call mpibcast_real(fq,(/nqpar,mqarray/)) ! if (lparticles) then call fetch_nparloc(np_aux) if (np_aux/=0) then call fetch_fp_array(fp_aux,dfp_aux,& index%ixw,index%iyw,index%izw,& index%ivxw,index%ivyw,index%ivzw) do k=1,np_aux xx_polar = fp_aux(k,index%ixw:index%izw) vv_polar = fp_aux(k,index%ivxw:index%ivzw) xxdot_polar = dfp_aux(k,index%ixw:index%izw) vvdot_polar = dfp_aux(k,index%ivxw:index%ivzw) aa_nbody_cart=dfp_aux(k,ivpx_cart:ivpz_cart) ! call advance_particles_in_cartesian(xx_polar,vv_polar,xxdot_polar,vvdot_polar,aa_nbody_cart) ! fp_aux(k,index%ixw:index%izw)=xx_polar fp_aux(k,index%ivxw:index%ivzw)=vv_polar dfp_aux(k,index%ixw:index%izw)=xxdot_polar enddo call return_fp_array(fp_aux,dfp_aux) endif ! endif ! call keep_compiler_quiet(f) ! endsubroutine pointmasses_timestep_second !*********************************************************************** subroutine advance_particles_in_cartesian(xx_polar,vv_polar,xxdot_polar,vvdot_polar,aa_nbody_cart) ! ! 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 :: rad,phi,tht,zed,xp,yp,zp real :: vrad,vtht,vphi,vzed,vx,vy,vz real :: raddot,thtdot,phidot,zeddot,xdot,ydot,zdot real :: vraddot,vthtdot,vphidot,vzeddot,vxdot,vydot,vzdot real :: cosp,sinp,cost,sint ! real, dimension(3), intent(inout) :: xx_polar,vv_polar real, dimension(3) :: aa_nbody_cart,vv_cart,dvv_cart,xx_cart,dxx_cart,xxdot_polar,vvdot_polar ! if (lcartesian_coords) then ! xp = xx_polar(1); yp = xx_polar(2); zp = xx_polar(3) vx = vv_polar(1); vy = vv_polar(2); vz = vv_polar(3) xdot = xxdot_polar(1) + vx ydot = xxdot_polar(2) + vy zdot = xxdot_polar(3) + vz ! vxdot = vvdot_polar(1) + aa_nbody_cart(1) ! ivqx_cart vydot = vvdot_polar(2) + aa_nbody_cart(2) ! ivqy_cart vzdot = vvdot_polar(3) + aa_nbody_cart(3) ! else if (lcylindrical_coords) then ! ! The input position, velocities, and accelerations in cylindrical coordinates. ! rad = xx_polar(1) ; phi = xx_polar(2) ; zed = xx_polar(3) vrad = vv_polar(1) ; vphi = vv_polar(2) ; vzed = vv_polar(3) raddot = xxdot_polar(1) ; phidot = xxdot_polar(2) ; zeddot = xxdot_polar(3) vraddot = vvdot_polar(1) ; vphidot = vvdot_polar(2) ; vzeddot = vvdot_polar(3) ! ! Shortcuts. ! cosp = cos(phi) ; sinp = sin(phi) ! ! Transform the positions and velocities to Cartesian coordinates. ! xp = rad*cosp ; yp = rad*sinp ; zp = zed ! vx = vrad*cosp - vphi*sinp vy = vrad*sinp + vphi*cosp vz = vzed ! ! Transform the position and velocity derivatives to Cartesian coordinates. ! Add the cartesian velocity to position and cartesian acceleration to velocity. ! xdot = raddot*cosp - phidot*sinp + vx ydot = raddot*sinp + phidot*cosp + vy zdot = zeddot + vz ! vxdot = vraddot*cosp - vphidot*sinp + aa_nbody_cart(1) ! ivqx_cart vydot = vraddot*sinp + vphidot*cosp + aa_nbody_cart(2) ! ivqy_cart vzdot = aa_nbody_cart(3) ! ivqz_cart ! else if (lspherical_coords) then ! ! The input position, velocities, and accelerations in spherical coordinates. ! rad = xx_polar(1) ; tht = xx_polar(2) ; phi = xx_polar(3) raddot = xxdot_polar(1) ; thtdot = xxdot_polar(2) ; phidot = xxdot_polar(3) ! vrad = vv_polar(1) ; vtht = vv_polar(2) ; vphi = vv_polar(3) vraddot = vvdot_polar(1) ; vthtdot = vvdot_polar(2) ; vphidot = vvdot_polar(3) ! ! Shortcuts. ! cosp = cos(phi) ; sinp = sin(phi) cost = cos(tht) ; sint = sin(tht) ! ! Transform the positions and velocities 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 ! ! Transform the position and velocity derivatives to Cartesian coordinates. ! Add the cartesian velocity to position and cartesian acceleration to velocity. ! xdot = raddot*sint*cosp + thtdot*cost*cosp - phidot*sinp + vx ydot = raddot*sint*sinp + thtdot*cost*sinp + phidot*cosp + vy zdot = raddot*cost - thtdot*sint + vz ! vxdot = vraddot*sint*cosp + vthtdot*cost*cosp - vphidot*sinp + aa_nbody_cart(1) !ivxq_cart) vydot = vraddot*sint*sinp + vthtdot*cost*sinp + vphidot*cosp + aa_nbody_cart(2) !ivyq_cart) vzdot = vraddot*cost - vthtdot*sint + aa_nbody_cart(3) !ivzq_cart) ! endif ! ! Now the time-stepping in Cartesian coordinates. ! xx_cart(1) = xp ; xx_cart(2) = yp ; xx_cart(3) = zp dxx_cart(1) = xdot ; dxx_cart(2) = ydot ; dxx_cart(3) = zdot vv_cart(1) = vx ; vv_cart(2) = vy ; vv_cart(3) = vz dvv_cart(1) = vxdot ; dvv_cart(2) = vydot ; dvv_cart(3) = vzdot call update_position(xx_cart,dxx_cart,xx_polar,xxdot_polar) call update_velocity(vv_cart,dvv_cart,vv_polar,xx_polar) ! endsubroutine advance_particles_in_cartesian !*********************************************************************** subroutine update_position(xx_cart,dxx_cart,xx_polar,xxdot_polar) ! ! Update position if N-body is used in polar coordinates. ! ! 14-feb-14:wlad/coded ! real, dimension(3), intent(inout) :: xx_polar,xx_cart,dxx_cart real, dimension(3), intent(out) :: xxdot_polar real :: xp,yp,zp,xdot,ydot,zdot,cosp,sinp,sint,cost,phi,tht ! ! Update. ! xx_cart = xx_cart + dt_beta_ts(itsub)*dxx_cart ! ! Convert back to polar coordinates. ! xp=xx_cart(1); yp=xx_cart(2); zp=xx_cart(3) if (lcartesian_coords) then xx_polar(1) = xp xx_polar(2) = yp xx_polar(3) = zp else if (lcylindrical_coords) then xx_polar(1) = sqrt(xp**2+yp**2+zp**2) xx_polar(2) = atan2(yp,xp) xx_polar(3) = zp else if (lspherical_coords) then xx_polar(1) = sqrt(xp**2+yp**2+zp**2) xx_polar(2) = acos(zp/xx_polar(1)) xx_polar(3) = atan2(yp,xp) endif ! xdot=dxx_cart(1); ydot=dxx_cart(2); zdot=dxx_cart(3) if (lcartesian_coords) then xxdot_polar(1) = xdot xxdot_polar(2) = ydot xxdot_polar(3) = zdot elseif (lcylindrical_coords) then phi=xx_polar(2) cosp=cos(phi) ; sinp=sin(phi) xxdot_polar(1) = xdot*cosp + ydot*sinp !=vrad xxdot_polar(2) = -xdot*sinp + ydot*cosp !=vphi xxdot_polar(3) = zdot else if (lspherical_coords) then tht=xx_polar(2) phi=xx_polar(3) cost=cos(tht) ; sint=sin(tht) cosp=cos(phi) ; sinp=sin(phi) xxdot_polar(1) = xdot*sint*cosp + ydot*sint*sinp + zdot*cost !=vrad xxdot_polar(2) = xdot*cost*cosp + ydot*cost*sinp - zdot*sint !=vphi xxdot_polar(3) = -xdot*sinp + ydot*cosp !=vz endif endsubroutine update_position !*********************************************************************** subroutine update_velocity(vv_cart,dvv_cart,vv_polar,xx_polar) ! ! Update velocity if N-body is used in polar coordinates. ! ! 14-feb-14:wlad/coded ! !real, intent(in) :: ax,ay,az real :: vx,vy,vz real, dimension(3), intent(inout) :: vv_polar,vv_cart real, dimension(3), intent(in) :: dvv_cart,xx_polar ! real :: phi,tht real :: cosp,sinp,sint,cost ! ! dvv_cart contain the accelerations transformed to Cartesian, plus ! the n-body acceleration, also in Cartesian. Do the RK timestepping. ! vv_cart = vv_cart + dt_beta_ts(itsub)*dvv_cart ! ! Convert back to polar coordinates. ! vx=vv_cart(1); vy=vv_cart(2); vz=vv_cart(3) if (lcartesian_coords) then vv_polar(1) = vx vv_polar(2) = vy vv_polar(3) = vz else if (lcylindrical_coords) then phi=xx_polar(2) cosp=cos(phi) ; sinp=sin(phi) vv_polar(1) = vx*cosp + vy*sinp !=vrad vv_polar(2) = -vx*sinp + vy*cosp !=vphi vv_polar(3) = vz else if (lspherical_coords) then tht=xx_polar(2) phi=xx_polar(3) cost=cos(tht) ; sint=sin(tht) cosp=cos(phi) ; sinp=sin(phi) vv_polar(1) = vx*sint*cosp + vy*sint*sinp + vz*cost !=vrad vv_polar(2) = vx*cost*cosp + vy*cost*sinp - vz*sint !=vphi vv_polar(3) = -vx*sinp + vy*cosp !=vz endif ! endsubroutine update_velocity !*********************************************************************** endmodule PointMasses