!####################################################################### ! SUBROUTINE --- sub_calc_gravity_tree --- ### !####################################################################### ! ### ! Subroutine for Calculating the Gravitational Acceleration ### ! by Using Tree Method ### ! ### ! Output Arguments ### ! np_ave_g_tree : Average Number of Interacting Nodes ### ! ### ! 2010/08/24 ### ! 2011/06/27 ### !####################################################################### subroutine sub_calc_gravity_tree(np_ave_g_tree) !$ use omp_lib use mod_parameter use mod_variable implicit none !==== Output Arguments ================================================= real*8 np_ave_g_tree !==== Variables Used in This Subroutine ================================ integer i ! Particle Number (i-th Particle) integer j ! Particle Number (j-th Particle) integer k ! Components (x,y,z) real*8 xij(3) ! x(k,i)-x(k,j) real*8 rij2 ! Distance Squared real*8 rij,rij3 ! Dummy of Distance real*8 gc_m_rij_rev3 ! Dummy real*8 xi(3) ! Position of i-th Particle integer jj ! Dummy of Particle Number real*8 fun_m_eft_s3 ! Function of Effective Mass integer gt_tot ! Number of Interacting Node integer(8) :: gt_tot_tot ! Total Number of Interacting Node (8byte Integer) real*8 theta_g2 integer(4) :: trd real*8 dw !======================================================================= theta_g2=THETA_GRAV2 gt_tot_tot=0.0d0 trd=0 !==== Calculation of Gravitational Acceleration ======================== !$ call omp_set_num_threads(N_TRD) !$OMP parallel default(none) & !$OMP private(i,k,xi,gt_tot,j,rij2,jj, & !$OMP gc_m_rij_rev3, rij,rij3,xij,trd) & !$OMP shared(np,x,h,hg,gt_tot_tot,a_grav,theta_g2,id_nb,mc_nb,xc_nb,dis2_nb,omega,zeta) !$ trd=omp_get_thread_num() !---- Initialization ------------------ !$OMP do do i=1,np do k=1,3 a_grav(k,i)=0.0d0 enddo enddo !$OMP end do !-------------------------------------- !$OMP do reduction(+ : gt_tot_tot) do i=1,np do k=1,3 xi(k)=x(k,i) enddo ! ---- Get Interacting Nodes and Their Data --------------------- call sub_get_grav_tree_list(trd,theta_g2,xi,h(i),gt_tot) ! --------------------------------------------------------------- gt_tot_tot=gt_tot_tot+gt_tot do j=1,gt_tot rij2=dis2_nb(j,trd) if(id_nb(j,trd) .ne. -1) then jj=id_nb(j,trd) gc_m_rij_rev3 = 0.5d0*mc_nb(j,trd) & & *( fun_m_eft_s3(rij2,h(i))/h(i)**3 & & +fun_m_eft_s3(rij2,h(jj))/h(jj)**3 ) if(FLAG_hg_CONST == 1) then gc_m_rij_rev3 = gc_m_rij_rev3 & & + 0.5d0*mc_nb(j,trd) & & *( zeta(i)/omega(i)*dw(rij2,h(i)) & & +zeta(jj)/omega(jj)*dw(rij2,h(jj)) ) endif else rij=dsqrt(rij2) rij3=rij2*rij gc_m_rij_rev3=mc_nb(j,trd)/rij3 endif do k=1,3 xij(k)=xi(k)-xc_nb(k,j,trd) a_grav(k,i)=a_grav(k,i)-gc_m_rij_rev3*xij(k) enddo enddo enddo !$OMP end do !$OMP end parallel np_ave_g_tree=dble(gt_tot_tot)/dble(np) return end !####################################################################### ! SUBROUTINE --- sub_get_grav_tree_list --- ### !####################################################################### ! ### ! Subroutine for Getting the Interacting Node and Their Data ### ! ### ! Input Arguments ### ! trd : Thread Number ### ! theta_g2 : Angle**2 ### ! xi(k) : Position of i-th Particle (x,y,z) ### ! hi : Smoothing Length of i-th Particle ### ! Output Arguments ### ! gt_tot : Total Number of Interacting Nodes ### ! ### ! id_nb() : List of Particle ID ### ! mc_nb() : Mass of Interacting Node ### ! xc_nb() : COM of Interacting Node ### ! dis2_nb() : Distance Squared to Interacting Node ### ! ### ! 2010/08/24 ### ! 2011/06/27 ### !####################################################################### subroutine sub_get_grav_tree_list(trd,theta_g2,xi,hi,gt_tot) use mod_parameter use mod_variable implicit none !==== Arguments ======================================================== ! ---- Input Arguments --------------------------------------------- integer(4) :: trd real*8 theta_g2 real*8 xi(3) real*8 hi ! ---- Output Arguments -------------------------------------------- integer gt_tot !==== Variables Used in This Subroutine ================================ integer gt ! Dummy integer wait_list_head integer wait_list_tail ! Last Number of Waiting List real*8 h2i ! 2*hi integer node ! Node Number real*8 dis2_COM ! Dummy of Distance integer node_head ! First Number of Branch Node integer node_tail ! Last Number of Branch Node integer n_node ! Branch Node Number real*8 dis2,dis,dis_radius !======================================================================= gt_tot=0 gt=0 wait_list(1,trd)=1 wait_list_head=1 wait_list_tail=1 h2i=2.0d0*hi 555 continue ! ---- Case for Zero Node in Waiting List -------------------------- if(wait_list_head .gt. wait_list_tail) go to 999 ! ------------------------------------------------------------------ ! ---- Last Node in Waiting List ----------------------------------- node=wait_list(wait_list_head,trd) ! ------------------------------------------------------------------ ! ---- Calculation of Distance^2 between i-th Particle and Node ---- dis2_COM=(center_mass_node(1,node)-xi(1))**2 & +(center_mass_node(2,node)-xi(2))**2 & +(center_mass_node(3,node)-xi(3))**2 ! ------------------------------------------------------------------ ! ---- Distance between i-th Particle and COM of Node -------------- dis2=(center_node(1,node)-xi(1))**2 & +(center_node(2,node)-xi(2))**2 & +(center_node(3,node)-xi(3))**2 dis=dsqrt(dis2) ! ------------------------------------------------------------------ dis_radius=dmax1(h2i,2.0d0*h_max_node(node))+radius_node(node) ! ---- Case for One Particle in Node ------------------------------- if(flag_node(node) .ne. -1) then gt_tot=gt_tot+1 gt=gt+1 mc_nb(gt,trd)=mass_node(node) xc_nb(1,gt,trd)=center_mass_node(1,node) xc_nb(2,gt,trd)=center_mass_node(2,node) xc_nb(3,gt,trd)=center_mass_node(3,node) dis2_nb(gt,trd)=dis2_COM id_nb(gt,trd)=flag_node(node) else ! ---- Case for Particles in Node ---------------------------------- if(leng_node(node)**2 .lt. theta_g2*dis2_COM .and. & dis .gt. 1.2d0*dis_radius) then gt_tot=gt_tot+1 gt=gt+1 mc_nb(gt,trd)=mass_node(node) xc_nb(1,gt,trd)=center_mass_node(1,node) xc_nb(2,gt,trd)=center_mass_node(2,node) xc_nb(3,gt,trd)=center_mass_node(3,node) dis2_nb(gt,trd)=dis2_COM id_nb(gt,trd)=flag_node(node) else node_head=head_branch_id(node) node_tail=node_head+num_branch(node)-1 do n_node=node_head,node_tail wait_list_tail=wait_list_tail+1 wait_list(wait_list_tail,trd)=n_node enddo endif endif wait_list_head=wait_list_head+1 go to 555 999 continue return end !####################################################################### ! SUBROUTINE --- sub_calc_gravity_direct --- ### !####################################################################### ! ### ! Subroutine for Calculating the Gravitational Acceleration ### ! by Using Direct Method ### ! ### ! 2010/08/20 ### ! 2011/06/27 ### ! 2011/09/13 ### !####################################################################### subroutine sub_calc_gravity_direct !$ use omp_lib use mod_parameter use mod_variable implicit none !==== Variables Used in This Subroutine ================================ integer(4) :: i integer(4) :: j real(8) :: xij(3) real(8) :: rij ! Distance btwn i-th and j-th Particles real(8) :: rij2 ! rij**2 real(8) :: rij3 ! rij**3 real(8) :: h2i2 ! (2*h(i))**2 real(8) :: h2j2 ! (2*h(j))**2 real(8) :: h22_max ! max(h2i2,h2j2) real(8) :: gc_m_rij_rev3 real(8) :: fun_m_eft_s3 real(8) :: dw integer(4) :: trd !======================================================================= !$ call omp_set_num_threads(N_TRD) !$OMP parallel default(none) & !$OMP private(i,j,xij,rij2,h2i2,h2j2,h22_max,rij,rij3, & !$OMP gc_m_rij_rev3,trd) & !$OMP shared(np,a_grav,x,h,m,omega,zeta) !$ trd=omp_get_thread_num() !---- Initialization ------------------ !$OMP do do i=1,np a_grav(1:3,i)=0.0d0 enddo !$OMP end do !-------------------------------------- !$OMP do do i=1,np do j=1,np xij(1:3)=x(1:3,i)-x(1:3,j) rij2=xij(1)**2+xij(2)**2+xij(3)**2 h2i2=4.0d0*h(i)**2 h2j2=4.0d0*h(j)**2 h22_max=dmax1(h2i2,h2j2) if(rij2 >= h22_max) then rij=dsqrt(rij2) rij3=rij2*rij gc_m_rij_rev3=m(j)/rij3 else gc_m_rij_rev3 = 0.5d0*m(j) & & *( fun_m_eft_s3(rij2,h(i))/h(i)**3 & & +fun_m_eft_s3(rij2,h(j))/h(j)**3 ) if(FLAG_hg_CONST == 1) then gc_m_rij_rev3 = gc_m_rij_rev3 & & + 0.5d0*m(j) & & *( zeta(i)/omega(i)*dw(rij2,h(i)) & & +zeta(j)/omega(j)*dw(rij2,h(j)) ) endif endif a_grav(1:3,i)=a_grav(1:3,i)-gc_m_rij_rev3*xij(1:3) enddo enddo !$OMP end do !$OMP end parallel return end subroutine sub_calc_gravity_direct !####################################################################### ! FUNCTION --- fun_m_eft_s3 --- ### !####################################################################### ! ### ! Function for Calculating the Effective Mass of j-th Particle ### ! ### ! Input Arguments ### ! r2 : Distance Squared ### ! h : Smoothing Length ### ! ### ! 2010/08/20 ### ! 2011/06/27 ### !####################################################################### real*8 function fun_m_eft_s3(r2,h) implicit none !==== Input Arguments ================================================== real*8 r2 real*8 h !==== Variables Used in This Function ================================== real*8 s1,s2,s3,s3_rev !==== Parameters Used in This Function ================================= real(8), parameter :: C1=1.0d0/15.0d0 real(8), parameter :: C2=8.0d0/3.0d0 real(8), parameter :: C3=1.0d0/6.0d0 real(8), parameter :: C4=4.0d0/3.0d0 !======================================================================= s1=dsqrt(r2)/h if (s1 .ge. 2.0d0) then s3_rev=1.0d0/s1**3 fun_m_eft_s3=s3_rev else if (s1 .ge. 1.0d0) then s2=s1*s1 s3=s1*s2 s3_rev=1.0d0/s3 fun_m_eft_s3=-C1*s3_rev+C2-3.0d0*s1 & +1.2d0*s2-C3*s3 else s2=s1*s1 s3=s1*s2 fun_m_eft_s3=C4-1.2d0*s2+0.5d0*s3 end if return end