diff --git a/AUTHORS b/AUTHORS index e2fa5818b..d2950816a 100644 --- a/AUTHORS +++ b/AUTHORS @@ -11,8 +11,8 @@ Conrad Chan James Wurster David Liptai Lionel Siess -Yann Bernard Fangyi (Fitz) Hu +Yann Bernard Megha Sharma Daniel Mentiplay Alison Young @@ -68,6 +68,7 @@ Farzana Meru Nicolás Cuello Cheryl-Lau Chris Nixon +Chunliang Mu <86601204+chunliangmu@users.noreply.github.com> Elisa Castro Martínez Miguel Gonzalez-Bolivar Sahl Rowther @@ -104,6 +105,7 @@ Shunquan Huang Steven Rieder Stéven Toupin Taj Jankovič +Yann bernard <46387267+Yrisch@users.noreply.github.com> bmcfarlane1234 eragusa <49534894+eragusa@users.noreply.github.com> fhu diff --git a/src/main/apr.f90 b/src/main/apr.f90 index 980914af6..bd16841f8 100644 --- a/src/main/apr.f90 +++ b/src/main/apr.f90 @@ -18,7 +18,7 @@ module apr ! mpiforce, neighkdtree, part, physcon, quitdump, random, relaxem, ! timestep_ind, utils_apr, vectorutils ! - use dim, only:use_apr + use dim, only:gr,use_apr use apr_region use utils_apr @@ -54,13 +54,14 @@ subroutine init_apr(apr_level,ierr) ! the resolution levels are in addition to the base resolution apr_max = apr_max_in + 1 + if (split_dir == 2) do_relax = .true. ! if we're reading in a file that already has the levels set, ! don't override these previously_set = .false. if (sum(int(apr_level(1:npart))) > npart) then previously_set = .true. - do_relax = .false. + if (split_dir /= 2) do_relax = .false. endif if (.not.previously_set) then @@ -220,7 +221,7 @@ subroutine update_apr(npart,xyzh,vxyzu,fxyzu,apr_level) ! level it does have, increment it up one if (apri > apr_current) then call splitpart(ii,npartnew) - if (do_relax .and. (apri == top_level)) then + if (do_relax .and. (gr .or. apri == top_level)) then nrelax = nrelax + 2 relaxlist(nrelax-1) = ii relaxlist(nrelax) = npartnew @@ -291,7 +292,7 @@ subroutine update_apr(npart,xyzh,vxyzu,fxyzu,apr_level) ! If we need to relax, do it here if (nrelax > 0 .and. do_relax) call relax_particles(npart,n_ref,xyzh_ref,force_ref,nrelax,relaxlist) ! Turn it off now because we only want to do this on first splits - do_relax = .false. + if (.not. gr) do_relax = .false. ! As we may have killed particles, time to do an array shuffle call shuffle_part(npart) @@ -322,18 +323,21 @@ end subroutine update_apr !+ !----------------------------------------------------------------------- subroutine splitpart(i,npartnew) - use part, only:copy_particle_all,apr_level,xyzh,vxyzu,npartoftype,igas - use part, only:set_particle_type + use part, only:copy_particle_all,apr_level,xyzh,vxyzu,npartoftype,igas,dens, & + set_particle_type,metrics,metricderivs,fext,pxyzu,eos_vars,itemp, & + igamma,igasP,aprmassoftype use physcon, only:pi use dim, only:ind_timesteps use random, only:ran2 use vectorutils, only:cross_product3D,rotatevec use utils_apr, only:apr_region_is_circle,icentre + use metric_tools, only:pack_metric,pack_metricderivs + use extern_gr, only:get_grforce integer, intent(in) :: i integer, intent(inout) :: npartnew integer :: j,npartold,next_door real :: theta,dx,dy,dz,x_add,y_add,z_add,sep,rneigh - real :: v(3),u(3),w(3),a,b,c,mag_v + real :: v(3),u(3),w(3),a,b,c,mag_v,uold,hnew,pmass integer, save :: iseed = 4 integer(kind=1) :: aprnew @@ -345,80 +349,147 @@ subroutine splitpart(i,npartnew) sep = sep_factor endif - ! Calculate the plane that the particle must be split along - ! to be tangential to the splitting region. Particles are split - ! on this plane but rotated randomly on it. - dx = xyzh(1,i) - apr_centre(1,icentre) - dy = xyzh(2,i) - apr_centre(2,icentre) - if (.not.apr_region_is_circle) then - dz = xyzh(3,i) - apr_centre(3,icentre) ! for now, let's split about the CoM - - if (directional) then - ! Calculate a vector, v, that lies on the plane - u = (/1.0,0.5,1.0/) - w = (/dx,dy,dz/) - call cross_product3D(u,w,v) - - ! rotate it around the normal to the plane by a random amount - theta = ran2(iseed)*2.*pi - call rotatevec(v,w,theta) - else - ! No directional splitting, so just create a unit vector in a random direction - a = ran2(iseed) - 0.5 - b = ran2(iseed) - 0.5 - c = ran2(iseed) - 0.5 - v = (/a, b, c/) - endif + if (gr) then + sep = sep*xyzh(4,i) + + npartold = npartnew + npartnew = npartold + 1 + npartoftype(igas) = npartoftype(igas) + 1 + apr_level(i) = apr_level(i) + int(1,kind=1) ! to prevent compiler warnings + call copy_particle_all(i,npartnew,new_part=.true.) + pmass = aprmassoftype(igas,apr_level(i)) + + uold = vxyzu(4,i) + hnew = xyzh(4,i)*(0.5**(1./3.)) + + ! new part forward + xyzh(4,npartnew) = hnew ! set new smoothing length + call integrate_geodesic_gr(pmass,xyzh(:,npartnew),vxyzu(:,npartnew),dens(npartnew),eos_vars(igasP,npartnew), & + eos_vars(igamma,npartnew),eos_vars(itemp,npartnew),pxyzu(:,npartnew),sep) + call pack_metric(xyzh(1:3,npartnew),metrics(:,:,:,npartnew)) + call pack_metricderivs(xyzh(1:3,npartnew),metricderivs(:,:,:,npartnew)) + call get_grforce(xyzh(:,npartnew),metrics(:,:,:,npartnew),metricderivs(:,:,:,npartnew), & + vxyzu(1:3,npartnew),dens(npartnew),vxyzu(4,npartnew),eos_vars(igasP,npartnew),fext(1:3,npartnew)) + if (ind_timesteps) call put_in_smallest_bin(npartnew) + + ! old part backward + ! switch direction + vxyzu(1:3,i) = -vxyzu(1:3,i) + pxyzu(1:3,i) = -pxyzu(1:3,i) + xyzh(4,i) = hnew + call integrate_geodesic_gr(pmass,xyzh(:,i),vxyzu(:,i),dens(i),eos_vars(igasP,i),eos_vars(igamma,i),eos_vars(itemp,i), & + pxyzu(:,i),sep) + ! switch direction back + vxyzu(1:3,i) = -vxyzu(1:3,i) + pxyzu(1:3,i) = -pxyzu(1:3,i) + call pack_metric(xyzh(1:3,i),metrics(:,:,:,i)) + call pack_metricderivs(xyzh(1:3,i),metricderivs(:,:,:,i)) + call get_grforce(xyzh(:,i),metrics(:,:,:,i),metricderivs(:,:,:,i), & + vxyzu(1:3,i),dens(i),vxyzu(4,i),eos_vars(igasP,i),fext(1:3,i)) + if (ind_timesteps) call put_in_smallest_bin(i) - mag_v = sqrt(dot_product(v,v)) - if (mag_v > tiny(mag_v)) then - v = v/mag_v - else - v = 0. - endif else - dz = 0. - u = 0. - w = 0. - v = 0. - theta = atan2(dy,dx) + 0.5*pi - v(1) = cos(theta) - v(2) = sin(theta) - endif - - ! Now apply it - x_add = sep*v(1)*xyzh(4,i) - y_add = sep*v(2)*xyzh(4,i) - z_add = sep*v(3)*xyzh(4,i) - - npartold = npartnew - npartnew = npartold + 1 - npartoftype(igas) = npartoftype(igas) + 1 - aprnew = apr_level(i) + int(1,kind=1) ! to prevent compiler warnings - - - !--create the new particle - do j=npartold+1,npartnew - call copy_particle_all(i,j,new_part=.true.) - xyzh(1,j) = xyzh(1,i) + x_add - xyzh(2,j) = xyzh(2,i) + y_add - xyzh(3,j) = xyzh(3,i) + z_add - vxyzu(:,j) = vxyzu(:,i) - xyzh(4,j) = xyzh(4,i)*(0.5**(1./3.)) - apr_level(j) = aprnew - if (ind_timesteps) call put_in_smallest_bin(j) - enddo - + if (split_dir == 2) then + sep = sep*xyzh(4,i) + + npartold = npartnew + npartnew = npartold + 1 + npartoftype(igas) = npartoftype(igas) + 1 + apr_level(i) = apr_level(i) + int(1,kind=1) ! to prevent compiler warnings + call copy_particle_all(i,npartnew,new_part=.true.) + pmass = aprmassoftype(igas,apr_level(i)) + + uold = vxyzu(4,i) + hnew = xyzh(4,i)*(0.5**(1./3.)) + + ! new part forward + xyzh(4,npartnew) = hnew ! set new smoothing length + call integrate_geodesic(pmass,xyzh(:,npartnew),vxyzu(:,npartnew),sep,1.) + if (ind_timesteps) call put_in_smallest_bin(npartnew) + + ! old part backward + ! switch direction + vxyzu(1:3,i) = -vxyzu(1:3,i) + xyzh(4,i) = hnew + call integrate_geodesic(pmass,xyzh(:,i),vxyzu(:,i),sep,1.) + ! switch direction back + vxyzu(1:3,i) = -vxyzu(1:3,i) + vxyzu(4,i) = uold + if (ind_timesteps) call put_in_smallest_bin(i) + else + ! Calculate the plane that the particle must be split along + ! to be tangential to the splitting region. Particles are split + ! on this plane but rotated randomly on it. + + dx = xyzh(1,i) - apr_centre(1,icentre) + dy = xyzh(2,i) - apr_centre(2,icentre) + if (.not.apr_region_is_circle) then + dz = xyzh(3,i) - apr_centre(3,icentre) ! for now, let's split about the CoM + + if (directional) then + ! Calculate a vector, v, that lies on the plane + u = (/1.0,0.5,1.0/) + w = (/dx,dy,dz/) + call cross_product3D(u,w,v) + + ! rotate it around the normal to the plane by a random amount + theta = ran2(iseed)*2.*pi + call rotatevec(v,w,theta) + else + ! No directional splitting, so just create a unit vector in a random direction + a = ran2(iseed) - 0.5 + b = ran2(iseed) - 0.5 + c = ran2(iseed) - 0.5 + v = (/a, b, c/) + endif - ! Edit the old particle that was sent in and kept - xyzh(1,i) = xyzh(1,i) - x_add - xyzh(2,i) = xyzh(2,i) - y_add - xyzh(3,i) = xyzh(3,i) - z_add - apr_level(i) = aprnew - xyzh(4,i) = xyzh(4,i)*(0.5**(1./3.)) - if (ind_timesteps) call put_in_smallest_bin(i) + mag_v = sqrt(dot_product(v,v)) + if (mag_v > tiny(mag_v)) then + v = v/mag_v + else + v = 0. + endif + else + dz = 0. + u = 0. + w = 0. + v = 0. + theta = atan2(dy,dx) + 0.5*pi + v(1) = cos(theta) + v(2) = sin(theta) + endif + ! Now apply it + x_add = sep*v(1)*xyzh(4,i) + y_add = sep*v(2)*xyzh(4,i) + z_add = sep*v(3)*xyzh(4,i) + + npartold = npartnew + npartnew = npartold + 1 + npartoftype(igas) = npartoftype(igas) + 1 + aprnew = apr_level(i) + int(1,kind=1) ! to prevent compiler warnings + + !--create the new particle + do j=npartold+1,npartnew + call copy_particle_all(i,j,new_part=.true.) + xyzh(1,j) = xyzh(1,i) + x_add + xyzh(2,j) = xyzh(2,i) + y_add + xyzh(3,j) = xyzh(3,i) + z_add + vxyzu(:,j) = vxyzu(:,i) + xyzh(4,j) = xyzh(4,i)*(0.5**(1./3.)) + apr_level(j) = aprnew + if (ind_timesteps) call put_in_smallest_bin(j) + enddo + ! Edit the old particle that was sent in and kept + xyzh(1,i) = xyzh(1,i) - x_add + xyzh(2,i) = xyzh(2,i) - y_add + xyzh(3,i) = xyzh(3,i) - z_add + apr_level(i) = aprnew + xyzh(4,i) = xyzh(4,i)*(0.5**(1./3.)) + if (ind_timesteps) call put_in_smallest_bin(i) + endif + endif end subroutine splitpart !----------------------------------------------------------------------- @@ -545,4 +616,119 @@ subroutine put_in_smallest_bin(i) end subroutine put_in_smallest_bin +!----------------------------------------------------------------------- +!+ +! Integrate particle along the geodesic +! Update vel and metric for best energy conservation +!+ +!----------------------------------------------------------------------- +subroutine integrate_geodesic_gr(pmass,xyzh,vxyzu,dens,pr,gamma,temp,pxyzu,dist) + use extern_gr, only:get_grforce + use metric_tools, only:pack_metric,pack_metricderivs + use eos, only:ieos,equationofstate + use cons2primsolver,only:conservative2primitive + use io, only:warning + use part, only:rhoh,ien_type + real, intent(inout) :: xyzh(:),vxyzu(:),pxyzu(:) + real, intent(inout) :: dens,pr,gamma,temp,pmass + real, intent(in) :: dist + real :: metrics(0:3,0:3,2),metricderivs(0:3,0:3,3),fext(3) + real :: t,tend,v,dt + real :: xyz(3),pxyz(3),eni,vxyz(1:3),uui,rho,spsoundi,pondensi + integer :: ierr + + xyz = xyzh(1:3) + pxyz = pxyzu(1:3) + eni = pxyzu(4) + vxyz = vxyzu(1:3) + uui = vxyzu(4) + rho = rhoh(xyzh(4),pmass) + + v = sqrt(dot_product(vxyz,vxyz)) + tend = dist/v + + call pack_metric(xyz,metrics(:,:,:)) + call pack_metricderivs(xyz,metricderivs(:,:,:)) + call get_grforce(xyzh(:),metrics(:,:,:),metricderivs(:,:,:),vxyz,dens,uui,pr,fext(1:3),dt) + + t = 0. + + do while (t <= tend) + dt = min(0.1,tend*0.1,dt) + t = t + dt + pxyz = pxyz + dt*fext + + call conservative2primitive(xyz,metrics(:,:,:),vxyz,dens,uui,pr,& + temp,gamma,rho,pxyz,eni,ierr,ien_type) + if (ierr > 0) call warning('cons2primsolver [in integrate_geodesic (a)]','did not converge') + + xyz = xyz + dt*vxyz + call pack_metric(xyz,metrics(:,:,:)) + call pack_metricderivs(xyz,metricderivs(:,:,:)) + + call equationofstate(ieos,pondensi,spsoundi,dens,xyzh(1),xyzh(2),xyzh(3),temp,uui) + pr = pondensi*dens + + call get_grforce(xyzh(:),metrics(:,:,:),metricderivs(:,:,:),vxyzu(1:3),dens,vxyzu(4),pr,fext(1:3),dt) + enddo + + xyzh(1:3) = xyz(1:3) + vxyzu(1:3) = vxyz(1:3) + pxyzu(1:3) = pxyz(1:3) + +end subroutine integrate_geodesic_gr + +!----------------------------------------------------------------------- +!+ +! Integrate particle along the geodesic +! Update vel and metric for best energy conservation +!+ +!----------------------------------------------------------------------- +subroutine integrate_geodesic(pmass,xyzh,vxyzu,dist,timei) + use options, only:iexternalforce + use externalforces, only:externalforce,externalforce_vdependent + use part, only:rhoh + real, intent(inout) :: xyzh(:),vxyzu(:) + real, intent(in) :: dist,timei,pmass + real :: fext(3),fextv(3) + real :: t,tend,v,dt,dens + real :: xyz(3),vxyz(1:3),poti,uui + + xyz = xyzh(1:3) + vxyz = vxyzu(1:3) + fext = 0. + uui = vxyzu(4) + dens = rhoh(xyzh(4),pmass) + + v = sqrt(dot_product(vxyz,vxyz)) + tend = dist/v + + if (iexternalforce > 0) then + call externalforce(iexternalforce,xyz(1),xyz(2),xyz(3),xyzh(4), & + timei,fext(1),fext(2),fext(3),poti,dt) + call externalforce_vdependent(iexternalforce,xyz,vxyz,fextv,poti,dens,uui) + fext = fext + fextv + endif + + t = 0. + + do while (t <= tend) + dt = min(0.1,tend*0.1,dt,0.1*dt) + t = t + dt + + vxyz = vxyz + dt*fext + xyz = xyz + dt*vxyz + + if (iexternalforce > 0) then + call externalforce(iexternalforce,xyz(1),xyz(2),xyz(3),xyzh(4), & + timei,fext(1),fext(2),fext(3),poti,dt) + call externalforce_vdependent(iexternalforce,xyz,vxyz,fextv,poti,dens,uui) + fext = fext + fextv + endif + enddo + + xyzh(1:3) = xyz(1:3) + vxyzu(1:3) = vxyz(1:3) +end subroutine integrate_geodesic + end module apr diff --git a/src/main/cons2prim.f90 b/src/main/cons2prim.f90 index 086456756..713f4b504 100644 --- a/src/main/cons2prim.f90 +++ b/src/main/cons2prim.f90 @@ -45,8 +45,10 @@ module cons2prim !+ !---------------------------------------------------------------------- subroutine prim2consall(npart,xyzh,metrics,vxyzu,pxyzu,use_dens,dens,use_sink) - use part, only:isdead_or_accreted,ien_type,eos_vars,igasP,igamma,itemp + use part, only:isdead_or_accreted,ien_type,eos_vars,igasP,igamma,itemp,igas,& + aprmassoftype,apr_level,massoftype use eos, only:gamma,ieos + use dim, only:use_apr integer, intent(in) :: npart real, intent(in) :: xyzh(:,:),metrics(:,:,:,:),vxyzu(:,:) real, intent(out) :: pxyzu(:,:) @@ -54,7 +56,7 @@ subroutine prim2consall(npart,xyzh,metrics,vxyzu,pxyzu,use_dens,dens,use_sink) logical, intent(in), optional :: use_dens, use_sink logical :: usedens integer :: i - real :: pri,tempi,xyzhi(4),vxyzui(4),densi + real :: pmassi,pri,tempi,xyzhi(4),vxyzui(4),densi ! By default, use the smoothing length to compute primitive density, and then compute the conserved variables. ! (Alternatively, use the provided primitive density to compute conserved variables. @@ -67,7 +69,8 @@ subroutine prim2consall(npart,xyzh,metrics,vxyzu,pxyzu,use_dens,dens,use_sink) !$omp parallel do default (none) & !$omp shared(xyzh,metrics,vxyzu,dens,pxyzu,npart,usedens,ien_type,eos_vars,gamma,ieos,use_sink,use_dens) & - !$omp private(i,pri,tempi,xyzhi,vxyzui,densi) + !$omp shared(massoftype,aprmassoftype,apr_level) & + !$omp private(i,pmassi,pri,tempi,xyzhi,vxyzui,densi) do i=1,npart if (present(use_sink)) then @@ -76,11 +79,17 @@ subroutine prim2consall(npart,xyzh,metrics,vxyzu,pxyzu,use_dens,dens,use_sink) vxyzui(1:3) = vxyzu(1:3,i) vxyzui(4) = 0. ! assume energy as 0. for sink densi = 1. - call prim2consi(xyzhi,metrics(:,:,:,i),vxyzui,pri,tempi,pxyzu(:,i),ien_type,& + call prim2consi(pmassi,xyzhi,metrics(:,:,:,i),vxyzui,pri,tempi,pxyzu(:,i),ien_type,& use_sink=use_sink,dens_i=densi) ! this returns temperature and pressure as 0. else if (.not.isdead_or_accreted(xyzh(4,i))) then - call prim2consi(xyzh(:,i),metrics(:,:,:,i),vxyzu(:,i),pri,tempi,pxyzu(:,i),ien_type,& + if (use_apr) then + pmassi = aprmassoftype(igas,apr_level(i)) + else + pmassi = massoftype(igas) + endif + + call prim2consi(pmassi,xyzh(:,i),metrics(:,:,:,i),vxyzu(:,i),pri,tempi,pxyzu(:,i),ien_type,& use_dens=usedens,dens_i=dens(i)) ! save eos vars for later use @@ -105,12 +114,12 @@ end subroutine prim2consall ! for a single SPH particle !+ !---------------------------------------------------------------------- -subroutine prim2consi(xyzhi,metrici,vxyzui,pri,tempi,pxyzui,ien_type,use_dens,use_sink,dens_i) +subroutine prim2consi(pmassi,xyzhi,metrici,vxyzui,pri,tempi,pxyzui,ien_type,use_dens,use_sink,dens_i) use cons2primsolver, only:primitive2conservative use utils_gr, only:h2dens use eos, only:equationofstate,ieos real, dimension(4), intent(in) :: xyzhi, vxyzui - real, intent(in) :: metrici(:,:,:) + real, intent(in) :: pmassi,metrici(:,:,:) real, intent(inout) :: pri,tempi integer, intent(in) :: ien_type real, dimension(4), intent(out) :: pxyzui @@ -140,7 +149,7 @@ subroutine prim2consi(xyzhi,metrici,vxyzui,pri,tempi,pxyzui,ien_type,use_dens,us densi = 1. ! using a value of 0. results in NaN values for the pxyzui array. pondensi = 0. else - call h2dens(densi,xyzhi,metrici,vi) ! Compute dens from h + call h2dens(densi,pmassi,xyzhi,metrici,vi) ! Compute dens from h dens_i = densi ! Feed the newly computed dens back out of the routine call equationofstate(ieos,pondensi,spsoundi,densi,xyzi(1),xyzi(2),xyzi(3),tempi,ui) endif @@ -201,7 +210,7 @@ subroutine cons2primall(npart,xyzh,metrics,pxyzu,vxyzu,dens,eos_vars) eos_vars(igamma,i) = gammai if (ierr > 0) then print*,' pmom =',pxyzu(1:3,i) - print*,' rho* =',rhoh(xyzh(4,i),massoftype(igas)) + print*,' rho* =',rhoh(xyzh(4,i),pmassi) print*,' en =',pxyzu(4,i) call fatal('cons2prim','could not solve rootfinding',i) endif diff --git a/src/main/inject_sim.f90 b/src/main/inject_sim.f90 index 3eb88b742..5a3d8bf32 100644 --- a/src/main/inject_sim.f90 +++ b/src/main/inject_sim.f90 @@ -61,7 +61,7 @@ subroutine init_inject(ierr) !--find the tde dump at the right time ! next_time = -1. - next_dump = getnextfilename(start_dump) + next_dump = start_dump call get_dump_time_npart(trim(next_dump),next_time,ierr,npart_out=npart_sim) ierr = 0 niter = 0 diff --git a/src/main/relaxem.f90 b/src/main/relaxem.f90 index 602e0bf0a..5b8b819d2 100644 --- a/src/main/relaxem.f90 +++ b/src/main/relaxem.f90 @@ -15,8 +15,9 @@ module relaxem ! :Runtime parameters: None ! ! :Dependencies: boundary, deriv, dim, eos, io, kernel, mpidomain, options, -! part +! part, utils_apr ! + use utils_apr, only:split_dir implicit none contains @@ -24,8 +25,9 @@ module relaxem ! Subroutine to relax the new set of particles to a reference particle distribution subroutine relax_particles(npart,n_ref,xyzh_ref,force_ref,nrelax,relaxlist) use deriv, only:get_derivs_global - use dim, only:mpi + use dim, only:gr,mpi use io, only:error + integer, intent(in) :: npart,n_ref,nrelax real, intent(in) :: force_ref(3,n_ref),xyzh_ref(4,n_ref) integer, intent(in) :: relaxlist(1:nrelax) @@ -34,18 +36,20 @@ subroutine relax_particles(npart,n_ref,xyzh_ref,force_ref,nrelax,relaxlist) logical :: converged integer :: ishift,nshifts - write(*,"(/,70('-'),/,/,2x,a,/,/)") 'APR: time to relax ...' + if (split_dir /= 2) write(*,"(/,70('-'),/,/,2x,a,/,/)") 'APR: time to relax ...' if (mpi) then call error('APR','relax_particles is not compatible with MPI') return endif - write(*,"(1x,1(a,i8,a,i8,a))") 'Relaxing',nrelax,' particles the heavenly way from',n_ref,' references.' + if (split_dir /= 2) write(*,"(1x,1(a,i8,a,i8,a))") & +'Relaxing',nrelax,' particles the heavenly way from',n_ref,' references.' ! Initialise for the loop converged = .false. ishift = 0 nshifts = 50 + if (split_dir /= 2) nshifts = 2 shuffle_tol = 0.05 ! a_ref stores the accelerations at the locations of the new particles as interpolated from the old ones @@ -65,7 +69,8 @@ subroutine relax_particles(npart,n_ref,xyzh_ref,force_ref,nrelax,relaxlist) if (ishift == 0) ke_init = ke - write(*,"(1x,1(a,f5.1,a,i3,a))") 'shuffle decreased to ',ke/ke_init*100.,'% of initial with',ishift,' shifts' + if (split_dir /= 2) write(*,"(1x,1(a,f5.1,a,i3,a))") & + 'shuffle decreased to ',ke/ke_init*100.,'% of initial with',ishift,' shifts' ! Todo: cut-off criteria if (ishift >= nshifts .or. (ke/ke_init < shuffle_tol)) converged = .true. @@ -76,7 +81,7 @@ subroutine relax_particles(npart,n_ref,xyzh_ref,force_ref,nrelax,relaxlist) ! Tidy up deallocate(a_ref) - write(*,"(/,/,2x,a,/,/,70('-'))") 'APR: relaxing finished.' + if (split_dir /= 2) write(*,"(/,/,2x,a,/,/,70('-'))") 'APR: relaxing finished.' end subroutine relax_particles @@ -157,9 +162,9 @@ end subroutine get_reference_accelerations !---------------------------------------------------------------- subroutine shift_particles(npart,a_ref,nrelax,relaxlist,ke,maxshift) - use dim, only:periodic + use dim, only:periodic,gr use part, only:xyzh,vxyzu,fxyzu,igas,aprmassoftype,rhoh, & - apr_level + apr_level,fext use eos, only:get_spsound use options, only:ieos use boundary, only:cross_boundary @@ -180,7 +185,7 @@ subroutine shift_particles(npart,a_ref,nrelax,relaxlist,ke,maxshift) !$omp parallel do schedule(guided) default(none) & !$omp shared(npart,xyzh,vxyzu,fxyzu,ieos,a_ref,maxshift) & - !$omp shared(apr_level,aprmassoftype) & + !$omp shared(split_dir,fext,apr_level,aprmassoftype) & !$omp shared(isperiodic,ncross,relaxlist,nrelax) & !$omp private(i,dx,dti,cs,rhoi,hi,vi,err,pmassi) & !$omp reduction(+:nlargeshift,ke) @@ -192,11 +197,16 @@ subroutine shift_particles(npart,a_ref,nrelax,relaxlist,ke,maxshift) rhoi = rhoh(hi,pmassi) cs = get_spsound(ieos,xyzh(:,i),rhoi,vxyzu(:,i)) dti = 0.3*hi/cs ! h/cs + if (split_dir == 2) dti = 0.1*hi/cs ! h/cs - dx = 0.5*dti**2*(fxyzu(1:3,i) - a_ref(1:3,i)) + dx = 0.5*dti**2*(fxyzu(1:3,i) + fext(1:3,i) - a_ref(1:3,i)) if (sqrt(dot_product(dx,dx)) > maxshift) maxshift = sqrt(dot_product(dx,dx)) - if (dot_product(dx,dx) > hi**2) then - + if (split_dir == 2) then + do while (dot_product(dx,dx) > hi**2) + dti = 0.1*dti + dx = 0.5*dti**2*(fxyzu(1:3,i) + fext(1:3,i) - a_ref(1:3,i)) + enddo + elseif (dot_product(dx,dx) > hi**2) then dx = dx / sqrt(dot_product(dx,dx)) * hi ! Avoid large shift in particle position !check with what James has done nlargeshift = nlargeshift + 1 endif diff --git a/src/main/substepping.F90 b/src/main/substepping.F90 index c3c9bc652..eb162e6a4 100644 --- a/src/main/substepping.F90 +++ b/src/main/substepping.F90 @@ -1222,7 +1222,7 @@ subroutine kickdrift_gr(dt,npart,nptmass,ntypes,xyzh,vxyzu,pxyzu,dens,metrics,me pri = eos_vars(igasP,i) gammai = eos_vars(igamma,i) tempi = eos_vars(itemp,i) - rhoi = rhoh(hi,massoftype(igas)) + rhoi = rhoh(hi,pmassi) ! since fext includes both the sink-gas interaction and the external force, ! we need to work out the "previous" force from the metric derivatives in order ! to perform the pmom_iterations diff --git a/src/main/utils_apr.f90 b/src/main/utils_apr.f90 index c5c6b8146..1e34d423a 100644 --- a/src/main/utils_apr.f90 +++ b/src/main/utils_apr.f90 @@ -19,6 +19,7 @@ module utils_apr ! - apr_type : *1: static, 2: sink, 3: clumps, 4: sequential sinks, 5: com, 6: vertical* ! - ref_dir : *increase (1) or decrease (-1) resolution* ! - rho_crit_cgs : *density above which apr zones are created (g/cm^3)* +! - split_dir : *split particle (1) tangential to split boundary, (2) along trajectory* ! - track_part : *number of sink to track* ! ! :Dependencies: infile_utils, io, part, ptmass @@ -34,6 +35,7 @@ module utils_apr integer :: apr_max_in = 3, ref_dir = 1, apr_type = 1, apr_max = 4 integer :: top_level = 1, ntrack = 1, ntrack_max = 10 integer :: icentre = 1, track_part_in = 1 + integer :: split_dir = 1 integer, allocatable :: npart_regions(:), track_part(:) real :: apr_rad = 1.0, apr_drad = 0.1, apr_centre_in(3) = 0. real, allocatable :: apr_regions(:), apr_centre(:,:) @@ -140,6 +142,9 @@ subroutine read_options_apr(name,valstring,imatch,igotall,ierr) case('ref_dir') read(valstring,*,iostat=ierr) ref_dir ngot = ngot + 1 + case('split_dir') + read(valstring,*,iostat=ierr) split_dir + ngot = ngot + 1 case('apr_type') read(valstring,*,iostat=ierr) apr_type ngot = ngot + 1 @@ -266,6 +271,7 @@ subroutine write_options_apr(iunit) write(iunit,"(/,a)") '# options for adaptive particle refinement' call write_inopt(apr_max_in,'apr_max','number of additional refinement levels (3 -> 2x resolution)',iunit) call write_inopt(ref_dir,'ref_dir','increase (1) or decrease (-1) resolution',iunit) + call write_inopt(split_dir,'split_dir','split particle (1) tangential to split boundary, (2) along trajectory',iunit) call write_inopt(apr_type,'apr_type','1: static, 2: sink, 3: clumps, 4: sequential sinks, 5: com, 6: vertical',iunit) select case (apr_type) diff --git a/src/main/utils_gr.f90 b/src/main/utils_gr.f90 index 76bf593a4..fe3eb05d3 100644 --- a/src/main/utils_gr.f90 +++ b/src/main/utils_gr.f90 @@ -90,15 +90,15 @@ end subroutine get_bigv ! the conserved density rho* (stored as the smoothing length) !+ !---------------------------------------------------------------- -subroutine h2dens(dens,xyzh,metrici,v) - use part, only: rhoh,massoftype,igas - real, intent(in) :: xyzh(1:4),metrici(:,:,:),v(1:3) +subroutine h2dens(dens,pmass,xyzh,metrici,v) + use part, only: rhoh + real, intent(in) :: pmass,xyzh(1:4),metrici(:,:,:),v(1:3) real, intent(out):: dens real :: rho, h, xyz(1:3) xyz = xyzh(1:3) h = xyzh(4) - rho = rhoh(h,massoftype(igas)) + rho = rhoh(h,pmass) call rho2dens(dens,rho,xyz,metrici,v) end subroutine h2dens @@ -283,4 +283,5 @@ subroutine perturb_metric(phi,gcovper,gcov) end subroutine perturb_metric + end module utils_gr diff --git a/src/utils/analysis_radiotde.f90 b/src/utils/analysis_radiotde.f90 index 06ee1f9fb..c9b9d46c8 100644 --- a/src/utils/analysis_radiotde.f90 +++ b/src/utils/analysis_radiotde.f90 @@ -41,9 +41,9 @@ module analysis real :: v_accum_mean, v_cap_mean real :: e_accum, e_cap integer :: n_accum, n_cap - real :: shock_v, rad_min, rad_max, shock_e, shock_m!, shock_rho - real :: shock_v_tde, rad_min_tde, rad_max_tde, shock_e_tde, shock_m_tde!, shock_rho - real :: shock_v_cnm, rad_min_cnm, rad_max_cnm, shock_e_cnm, shock_m_cnm!, shock_rho + real :: shock_v, rad_min, rad_max, shock_e, shock_m, shock_u, shock_g!, shock_rho + real :: shock_v_tde, rad_min_tde, rad_max_tde, shock_e_tde, shock_m_tde, shock_u_tde!, shock_rho + real :: shock_v_cnm, rad_min_cnm, rad_max_cnm, shock_e_cnm, shock_m_cnm, shock_u_cnm!, shock_rho !---- These can be changed in the params file real :: rad_cap = 1.e16 ! radius where the outflow in captured (in cm) @@ -202,7 +202,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) case ('shock') write(*,'(a)') ' Analysing the shock ...' - call shock_analysis(npart,pmass,rad_all,vr_all,pxyzu(4,:)) + call shock_analysis(npart,pmass,rad_all,vr_all,vxyzu(4,:),pxyzu(4,:)) deallocate(rad_all,vr_all,v_all) @@ -211,16 +211,20 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) open(iunit,file='shock',status='old',position='append') else open(iunit,file='shock',status='new') - write(iunit,'(17(A,1x))') '#', 'time', 'rad_min[cm]', 'rad_max[cm]', 'velocity[c]', 'mass[Msun]', 'energy[erg]', & !'density[g/cm-3]' - 'rad_min_tde[cm]', 'rad_max_tde[cm]', 'vel_tde[c]', 'mass_tde[Msun]', 'ene_tde[erg]', & - 'rad_min_cnm[cm]', 'rad_max_cnm[cm]', 'vel_cnm[c]', 'mass_cnm[Msun]', 'ene_cnm[erg]' + write(iunit,'(21(A,1x))') '#', 'time', 'rad_min[cm]', 'rad_max[cm]', 'velocity[c]', 'mass[Msun]', 'ki_ene[erg]', 'u[erg]', & !'density[g/cm-3]' + 'grav[erg]', 'rad_min_tde[cm]', 'rad_max_tde[cm]', 'vel_tde[c]', 'mass_tde[Msun]', & + 'ki_ene_tde[erg]', 'u_tde[erg]', & + 'rad_min_cnm[cm]', 'rad_max_cnm[cm]', 'vel_cnm[c]', 'mass_cnm[Msun]', & + 'ki_ene_cnm[erg]', 'u_cnm[erg]' endif if (rad_max > 0.) then - write(iunit,'(16(es18.10,1x))') & + write(iunit,'(20(es18.10,1x))') & time*todays, & - rad_min*udist, rad_max*udist, shock_v, shock_m*umass/solarm, shock_e*unit_energ, & - rad_min_tde*udist, rad_max_tde*udist, shock_v_tde, shock_m_tde*umass/solarm, shock_e_tde*unit_energ, & - rad_min_cnm*udist, rad_max_cnm*udist, shock_v_cnm, shock_m_cnm*umass/solarm, shock_e_cnm*unit_energ !shock_rho*unit_density + rad_min*udist, rad_max*udist, shock_v, shock_m*umass/solarm, shock_e*unit_energ, shock_u*unit_energ, & + shock_g*unit_energ, rad_min_tde*udist, rad_max_tde*udist, shock_v_tde, & + shock_m_tde*umass/solarm, shock_e_tde*unit_energ, shock_u_tde*unit_energ, & + rad_min_cnm*udist, rad_max_cnm*udist, shock_v_cnm, & + shock_m_cnm*umass/solarm, shock_e_cnm*unit_energ, shock_u_cnm*unit_energ !shock_rho*unit_density endif close(iunit) @@ -328,19 +332,19 @@ subroutine record_background(ent,npart_old,npart_new,ent_bg) print*, 'Record background entropy of ', npart_new, ' particles' do i=1,npart_new - ent_bg(npart_old+i) = ent(npart_old+i)*1.1 ! give some range for self evolution - !(is there a reasonable choice instead of arbitrary?) + ent_bg(npart_old+i) = ent(npart_old+i)*1.3 ! give some range for self evolution + !(is there a reasonable choice instead of arbitrary?) enddo end subroutine record_background -subroutine shock_analysis(npart,pmass,rad_all,vr_all,ent) +subroutine shock_analysis(npart,pmass,rad_all,vr_all,u,ent) use units, only: udist use physcon, only: au,pi integer, intent(in) :: npart - real, intent(in) :: pmass,rad_all(:),vr_all(:),ent(:) + real, intent(in) :: pmass,rad_all(:),vr_all(:),ent(:),u(:) integer :: i,n,n_cnm,n_tde - real :: ri,half_m,ei,vi + real :: ri,half_m,ei,vi,ui ! !------Determine the shock ! @@ -350,6 +354,9 @@ subroutine shock_analysis(npart,pmass,rad_all,vr_all,ent) shock_e = 0. shock_e_cnm = 0. shock_e_tde = 0. + shock_u = 0. + shock_u_cnm = 0. + shock_u_tde = 0. shock_v = 0. ! take max vel shock_v_cnm = 0. shock_v_tde = 0. @@ -366,11 +373,14 @@ subroutine shock_analysis(npart,pmass,rad_all,vr_all,ent) ri = rad_all(i) vi = vr_all(i) ei = half_m*vi**2 + ui = u(i)*pmass n = n + 1 if (vi > shock_v) shock_v = vi if (ri < rad_min) rad_min = ri if (ri > rad_max) rad_max = ri shock_e = shock_e + ei + shock_u = shock_u + ui + shock_g = pmass/ri if (i > npart_cnm) then ! tde outflow @@ -379,6 +389,7 @@ subroutine shock_analysis(npart,pmass,rad_all,vr_all,ent) if (ri < rad_min_tde) rad_min_tde = ri if (ri > rad_max_tde) rad_max_tde = ri shock_e_tde = shock_e_tde + ei + shock_u_tde = shock_u_tde + ui else ! cnm n_cnm = n_cnm + 1 @@ -386,6 +397,7 @@ subroutine shock_analysis(npart,pmass,rad_all,vr_all,ent) if (ri < rad_min_cnm) rad_min_cnm = ri if (ri > rad_max_cnm) rad_max_cnm = ri shock_e_cnm = shock_e_cnm + ei + shock_u_cnm = shock_u_cnm + ui endif endif enddo diff --git a/src/utils/moddump_radiotde.f90 b/src/utils/moddump_radiotde.f90 index 01e3b61ef..710fdea3c 100644 --- a/src/utils/moddump_radiotde.f90 +++ b/src/utils/moddump_radiotde.f90 @@ -64,7 +64,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) eos_vars,itemp,igamma,igasP use io, only:fatal,master,id use units, only:umass,udist,utime,set_units,unit_density - use timestep, only:dtmax,tmax + use timestep, only:dtmax,tmax,idtmax_frac,dtmax_ifactor,idtmax_n use eos, only:ieos,gmw use kernel, only:hfact_default use stretchmap, only:get_mass_r,rho_func @@ -253,6 +253,9 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) !--Set timesteps tmax = 3.*years/utime dtmax = tmax/1000. + dtmax_ifactor = 0 + idtmax_frac = 0 ! so don't write to .restart + idtmax_n = 1 end subroutine modify_dump diff --git a/src/utils/phantomanalysis.f90 b/src/utils/phantomanalysis.f90 index abebb9bee..47babb2f6 100644 --- a/src/utils/phantomanalysis.f90 +++ b/src/utils/phantomanalysis.f90 @@ -14,11 +14,11 @@ program phantomanalysis ! ! :Usage: phantomanalysis dumpfile(s) ! -! :Dependencies: analysis, dim, eos, externalforces, fileutils, +! :Dependencies: analysis, apr, dim, eos, externalforces, fileutils, ! infile_utils, io, kernel, part, readwrite_dumps ! - use dim, only:tagline,do_nucleation,inucleation - use part, only:xyzh,hfact,massoftype,vxyzu,npart !,npartoftype + use dim, only:tagline,do_nucleation,inucleation,use_apr + use part, only:xyzh,hfact,massoftype,vxyzu,npart,apr_level !,npartoftype use io, only:set_io_unit_numbers,iprint,idisk1,ievfile,ianalysis use readwrite_dumps, only:read_dump,read_smalldump,is_small_dump use infile_utils, only:open_db_from_file,inopts,read_inopt,close_db @@ -27,6 +27,7 @@ program phantomanalysis use eos, only:ieos use kernel, only:hfact_default use externalforces, only:mass1,accradius1 + use apr, only:init_apr implicit none integer :: nargs,iloc,ierr,iarg,i,idust_opacity real :: time @@ -48,6 +49,9 @@ program phantomanalysis stop endif + ! initialise apr if it is being used + if (use_apr) call init_apr(apr_level,ierr) + print "(/,a,/)",' Phantom analysis ('//trim(analysistype)//'): You data, we analyse' over_args: do iarg=1,nargs