From 3fc2bccd4d92fb4d9eeeeeee1c0f90a2f7d6a997 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 2 Nov 2023 10:01:26 +0100 Subject: [PATCH 01/19] (shock) set units in interactive setup --- src/setup/setup_shock.F90 | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/setup/setup_shock.F90 b/src/setup/setup_shock.F90 index 943e3ef2c..a5dfe3306 100644 --- a/src/setup/setup_shock.F90 +++ b/src/setup/setup_shock.F90 @@ -406,6 +406,7 @@ subroutine choose_shock (gamma,polyk,dtg,iexist) use physcon, only:pi,Rg,au,solarm use prompting, only:prompt use units, only:udist,utime,unit_density,unit_pressure + use setunits, only:set_units_interactive real, intent(inout) :: gamma,polyk real, intent(out) :: dtg logical, intent(in) :: iexist @@ -436,6 +437,9 @@ subroutine choose_shock (gamma,polyk,dtg,iexist) yright = 0.0 zright = 0.0 const = sqrt(4.*pi) + + call set_units_interactive(gr) + ! !--list of shocks ! @@ -679,6 +683,8 @@ end function get_conserved_density subroutine write_setupfile(filename,iprint,numstates,gamma,polyk,dtg) use infile_utils, only:write_inopt use dim, only:tagline + use setunits, only:write_options_units + use part, only:gr integer, intent(in) :: iprint,numstates real, intent(in) :: gamma,polyk,dtg character(len=*), intent(in) :: filename @@ -690,6 +696,8 @@ subroutine write_setupfile(filename,iprint,numstates,gamma,polyk,dtg) write(lu,"(a)") '# '//trim(tagline) write(lu,"(a)") '# input file for Phantom shock tube setup' + call write_options_units(lu,gr) + write(lu,"(/,a)") '# shock tube' do i=1,numstates call write_inopt(leftstate(i), trim(var_label(i))//'left', trim(var_label(i))//' (left)', lu,ierr1) @@ -754,6 +762,8 @@ end subroutine write_setupfile !------------------------------------------ subroutine read_setupfile(filename,iprint,numstates,gamma,polyk,dtg,ierr) use infile_utils, only:open_db_from_file,inopts,close_db,read_inopt + use setunits, only:read_options_and_set_units + use part, only:gr character(len=*), intent(in) :: filename integer, parameter :: lu = 21 integer, intent(in) :: iprint,numstates @@ -767,6 +777,10 @@ subroutine read_setupfile(filename,iprint,numstates,gamma,polyk,dtg,ierr) write(iprint, '(1x,2a)') 'Setup_shock: Reading setup options from ',trim(filename) nerr = 0 + + ! units + call read_options_and_set_units(db,nerr,gr) + do i=1,numstates call read_inopt(leftstate(i), trim(var_label(i))//'left',db,errcount=nerr) call read_inopt(rightstate(i),trim(var_label(i))//'right',db,errcount=nerr) From 2747924e8436b3ce83d9aa56bf60b45f63207f55 Mon Sep 17 00:00:00 2001 From: Miguel Gonzalez-Bolivar Date: Thu, 16 Nov 2023 17:25:52 +1100 Subject: [PATCH 02/19] Add v_esc option for .divv files --- src/utils/analysis_common_envelope.f90 | 42 +++++++++++++++++++++++--- 1 file changed, 37 insertions(+), 5 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 6cd1c6c27..bc16bcc17 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -1400,13 +1400,14 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu) real, dimension(3) :: com_xyz,com_vxyz,xyz_a,vxyz_a real :: pC, pC2, pC2H, pC2H2, nH_tot, epsC, S real :: taustar, taugr, JstarS + real :: v_esci real, parameter :: Scrit = 2. ! Critical saturation ratio logical :: verbose = .false. allocate(quant(4,npart)) - Nquantities = 13 + Nquantities = 14 if (dump_number == 0) then - print "(13(a,/))",& + print "(14(a,/))",& '1) Total energy (kin + pot + therm)', & '2) Mach number', & '3) Opacity from MESA tables', & @@ -1419,7 +1420,8 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu) '10) Mass coordinate', & '11) Gas omega w.r.t. CoM', & '12) Gas omega w.r.t. sink 1',& - '13) JstarS' !option to calculate JstarS + '13) JstarS', & + '14) Escape velocity' quantities_to_calculate = (/1,2,4,5/) call prompt('Choose first quantity to compute ',quantities_to_calculate(1),0,Nquantities) @@ -1435,7 +1437,7 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu) com_vxyz = 0. do k=1,4 select case (quantities_to_calculate(k)) - case(0,1,2,3,6,8,9,13) ! Nothing to do + case(0,1,2,3,6,8,9,13,14) ! Nothing to do case(4,5,11,12) ! Fractional difference between gas and orbital omega if (quantities_to_calculate(k) == 4 .or. quantities_to_calculate(k) == 5) then com_xyz = (xyzmh_ptmass(1:3,1)*xyzmh_ptmass(4,1) + xyzmh_ptmass(1:3,2)*xyzmh_ptmass(4,2)) & @@ -1582,6 +1584,9 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu) case(10) ! Mass coordinate quant(k,iorder(i)) = real(i,kind=kind(time)) * particlemass + case(14) ! Escape_velocity + call calc_escape_velocities(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),xyzmh_ptmass,phii,epoti,v_esci) + quant(k,i) = v_esci case default print*,"Error: Requested quantity is invalid." stop @@ -3868,7 +3873,8 @@ subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,xyzmh_ptmass,phii,epo epoti = 2.*poten + particlemass * phii ! For individual particles, need to multiply 2 to poten to get \sum_j G*mi*mj/r ekini = particlemass * 0.5 * dot_product(vxyzu(1:3),vxyzu(1:3)) einti = particlemass * vxyzu(4) - etoti = epoti + ekini + einti + etoti = epoti + ekini + einti + end subroutine calc_gas_energies @@ -4557,4 +4563,30 @@ subroutine set_eos_options(analysis_to_perform) end subroutine set_eos_options + +!---------------------------------------------------------------- +!+ +! Calculates escape velocity for all SPH particles given the potential energy +! of the system at that time +!+ +!---------------------------------------------------------------- +subroutine calc_escape_velocities(particlemass,poten,xyzh,vxyzu,xyzmh_ptmass,phii,epoti,v_esc) + use ptmass, only:get_accel_sink_gas + use part, only:nptmass + real, intent(in) :: particlemass + real(4), intent(in) :: poten + real, dimension(4), intent(in) :: xyzh,vxyzu + real, dimension(5,nptmass), intent(in) :: xyzmh_ptmass + real :: phii,epoti + real :: fxi,fyi,fzi + real, intent(out) :: v_esc + + phii = 0.0 + call get_accel_sink_gas(nptmass,xyzh(1),xyzh(2),xyzh(3),xyzh(4),xyzmh_ptmass,fxi,fyi,fzi,phii) + + epoti = 2.*poten + particlemass * phii ! For individual particles, need to multiply 2 to poten to get \sum_j G*mi*mj/r + v_esc = sqrt(2*abs(epoti/particlemass)) + +end subroutine calc_escape_velocities + end module analysis From b17a283876eaa93f87e2450ab83c25c9159dbe88 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Mon, 20 Nov 2023 10:09:08 +0100 Subject: [PATCH 03/19] (shock) only allow setting units when using radiation --- src/setup/setup_shock.F90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/setup/setup_shock.F90 b/src/setup/setup_shock.F90 index a5dfe3306..c39e04ed5 100644 --- a/src/setup/setup_shock.F90 +++ b/src/setup/setup_shock.F90 @@ -438,7 +438,7 @@ subroutine choose_shock (gamma,polyk,dtg,iexist) zright = 0.0 const = sqrt(4.*pi) - call set_units_interactive(gr) + if (do_radiation) call set_units_interactive(gr) ! !--list of shocks @@ -682,7 +682,7 @@ end function get_conserved_density !------------------------------------------ subroutine write_setupfile(filename,iprint,numstates,gamma,polyk,dtg) use infile_utils, only:write_inopt - use dim, only:tagline + use dim, only:tagline,do_radiation use setunits, only:write_options_units use part, only:gr integer, intent(in) :: iprint,numstates @@ -696,7 +696,7 @@ subroutine write_setupfile(filename,iprint,numstates,gamma,polyk,dtg) write(lu,"(a)") '# '//trim(tagline) write(lu,"(a)") '# input file for Phantom shock tube setup' - call write_options_units(lu,gr) + if (do_radiation) call write_options_units(lu,gr) write(lu,"(/,a)") '# shock tube' do i=1,numstates @@ -764,6 +764,7 @@ subroutine read_setupfile(filename,iprint,numstates,gamma,polyk,dtg,ierr) use infile_utils, only:open_db_from_file,inopts,close_db,read_inopt use setunits, only:read_options_and_set_units use part, only:gr + use dim, only:do_radiation character(len=*), intent(in) :: filename integer, parameter :: lu = 21 integer, intent(in) :: iprint,numstates @@ -779,7 +780,7 @@ subroutine read_setupfile(filename,iprint,numstates,gamma,polyk,dtg,ierr) nerr = 0 ! units - call read_options_and_set_units(db,nerr,gr) + if (do_radiation) call read_options_and_set_units(db,nerr,gr) do i=1,numstates call read_inopt(leftstate(i), trim(var_label(i))//'left',db,errcount=nerr) From f80b5ec543c3014b2530db14dabe1d3ff6f10eb9 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 10:09:12 +1100 Subject: [PATCH 04/19] (build) use -g not -gdwarf-2 in default gfortran debugging flags --- build/Makefile_defaults_gfortran | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/build/Makefile_defaults_gfortran b/build/Makefile_defaults_gfortran index 78de58f33..02a0dfdfe 100644 --- a/build/Makefile_defaults_gfortran +++ b/build/Makefile_defaults_gfortran @@ -12,7 +12,7 @@ # endif # FC= gfortran -FFLAGS+= -O3 -Wall -Wno-unused-dummy-argument -frecord-marker=4 -gdwarf-2 \ +FFLAGS+= -O3 -Wall -Wno-unused-dummy-argument -frecord-marker=4 -g \ -finline-functions-called-once -finline-limit=1500 -funroll-loops -ftree-vectorize \ -std=f2008 -fall-intrinsics DBLFLAG= -fdefault-real-8 -fdefault-double-8 From 75c952747541096d3d25523cf6172a9e2872ac03 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 10:19:44 +1100 Subject: [PATCH 05/19] (radiation) unused variable warnings fixed --- src/main/radiation_implicit.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/main/radiation_implicit.f90 b/src/main/radiation_implicit.f90 index f93abf079..450956569 100644 --- a/src/main/radiation_implicit.f90 +++ b/src/main/radiation_implicit.f90 @@ -423,7 +423,7 @@ subroutine fill_arrays(ncompact,ncompactlocal,npart,icompactmax,dt,xyzh,vxyzu,iv real, intent(out) :: vari(:,:),EU0(6,npart),varij(2,icompactmax),varij2(4,icompactmax) integer :: n,i,j,k,icompact real :: pmi,hi,hi21,hi41,rhoi,dx,dy,dz,rij2,rij,rij1,dr,dti,& - pmj,rhoj,hj,hj21,hj41,v2i,vi,v2j,vj,dWi,dWj,rhomean,& + pmj,rhoj,hj,hj21,hj41,v2i,vi,v2j,vj,dWi,dWj,& c_code,dWidrlightrhorhom,dWjdrlightrhorhom,& xi,yi,zi,gradhi,pmjdWrijrhoi,pmjdWrunix,pmjdWruniy,pmjdWruniz,& dust_kappai,dust_cooling,heatingISRi,dust_gas @@ -432,7 +432,7 @@ subroutine fill_arrays(ncompact,ncompactlocal,npart,icompactmax,dt,xyzh,vxyzu,iv !$omp do & !$omp private(n,i,j,k,rhoi,icompact,pmi,dti) & !$omp private(dx,dy,dz,rij2,rij,rij1,dr,pmj,rhoj,hi,hj,hi21,hj21,hi41,hj41) & - !$omp private(v2i,vi,v2j,vj,dWi,dWj,rhomean) & + !$omp private(v2i,vi,v2j,vj,dWi,dWj) & !$omp private(xi,yi,zi,gradhi,dWidrlightrhorhom,pmjdWrijrhoi,dWjdrlightrhorhom) & !$omp private(pmjdWrunix,pmjdWruniy,pmjdWruniz,dust_kappai,dust_cooling,heatingISRi,dust_gas) @@ -553,12 +553,12 @@ subroutine compute_flux(ivar,ijvar,ncompact,npart,icompactmax,varij2,vari,EU0,va real, intent(out) :: varinew(3,npart) ! we use this parallel loop to set varinew to zero integer :: i,j,k,n,icompact real :: rhoi,rhoj,pmjdWrunix,pmjdWruniy,pmjdWruniz,dedx(3),dradenij,rhoiEU0 - real :: gradE1i,opacity,radRi,EU01i + real :: opacity,radRi,EU01i !$omp do schedule(runtime)& !$omp private(i,j,k,n,dedx,rhoi,rhoj,icompact)& !$omp private(pmjdWrunix,pmjdWruniy,pmjdWruniz,dradenij)& - !$omp private(gradE1i,opacity,radRi,EU01i) + !$omp private(opacity,radRi,EU01i) do n = 1,ncompact i = ivar(3,n) @@ -632,12 +632,12 @@ subroutine calc_diffusion_term(ivar,ijvar,varij,ncompact,npart,icompactmax, & real, intent(inout) :: varinew(3,npart) integer :: n,i,j,k,icompact real :: rhoi,rhoj,opacityi,opacityj,Ej,bi,bj,b1,dWdrlightrhorhom - real :: diffusion_numerator,diffusion_denominator,tempval1,tempval2 + real :: diffusion_numerator,diffusion_denominator ierr = 0 !$omp do schedule(runtime)& !$omp private(i,j,k,n,rhoi,rhoj,opacityi,opacityj,Ej,bi,bj,b1,diffusion_numerator,diffusion_denominator)& - !$omp private(dWdrlightrhorhom,tempval1,tempval2,icompact)& + !$omp private(dWdrlightrhorhom,icompact)& !$omp reduction(max:ierr) do n = 1,ncompact i = ivar(3,n) @@ -721,7 +721,7 @@ subroutine update_gas_radiation_energy(ivar,vari,npart,ncompactlocal,& logical, intent(in) :: store_drad logical, intent(out):: moresweep logical, intent(inout):: mask(npart) - integer :: i,j,n,ieqtype,ierr + integer :: i,n,ieqtype,ierr logical :: moresweep2,skip_quartic real :: dti,rhoi,diffusion_numerator,diffusion_denominator,gradEi2,gradvPi,rpdiag,rpall real :: radpresdenom,stellarradiation,gas_temp,xnH2,betaval,gammaval,tfour,betaval_d,chival @@ -729,7 +729,7 @@ subroutine update_gas_radiation_energy(ivar,vari,npart,ncompactlocal,& real :: cosmic_ray,cooling_line,photoelectric,h2form,dust_heating,dust_term,e_planetesimali real :: u4term,u1term,u0term,pcoleni,dust_cooling,heatingISRi,dust_gas real :: pres_numerator,pres_denominator,mui,U1i,E1i,Tgas,dUcomb,dEcomb - real :: residualE,residualU,xchange,maxerrU2old,Tgas4,Trad4,ck,ack + real :: residualE,residualU,maxerrU2old,Tgas4,Trad4,ck,ack real :: Ei,Ui,cvi,opacityi,eddi real :: maxerrE2i,maxerrU2i @@ -742,12 +742,12 @@ subroutine update_gas_radiation_energy(ivar,vari,npart,ncompactlocal,& !$omp end single !$omp do schedule(runtime)& - !$omp private(i,j,n,rhoi,dti,diffusion_numerator,diffusion_denominator,U1i,skip_quartic,Tgas,E1i,dUcomb,dEcomb) & + !$omp private(i,n,rhoi,dti,diffusion_numerator,diffusion_denominator,U1i,skip_quartic,Tgas,E1i,dUcomb,dEcomb) & !$omp private(gradEi2,gradvPi,rpdiag,rpall,radpresdenom,stellarradiation,dust_tempi,dust_kappai,xnH2) & !$omp private(dust_cooling,heatingISRi,dust_gas,gas_dust_val,dustgammaval,gas_dust_cooling,cosmic_ray) & !$omp private(cooling_line,photoelectric,h2form,dust_heating,dust_term,betaval,chival,gammaval,betaval_d,tfour) & !$omp private(e_planetesimali,u4term,u1term,u0term,pcoleni,pres_numerator,pres_denominator,moresweep2,mui,ierr) & - !$omp private(residualE,residualU,xchange,maxerrU2old,gas_temp,ieqtype,unit_density,Tgas4,Trad4,ck,ack) & + !$omp private(residualE,residualU,maxerrU2old,gas_temp,ieqtype,unit_density,Tgas4,Trad4,ck,ack) & !$omp private(maxerrE2i,maxerrU2i) & !$omp reduction(max:maxerrE2,maxerrU2) main_loop: do n = 1,ncompactlocal From 8b42870aa5930f8d39e0c9eb0dbed374f991c83d Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 10:50:11 +1100 Subject: [PATCH 06/19] (test_gravity) maybe-unused warning fixed --- src/tests/test_gravity.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tests/test_gravity.f90 b/src/tests/test_gravity.f90 index ec23df777..ee2bb7069 100644 --- a/src/tests/test_gravity.f90 +++ b/src/tests/test_gravity.f90 @@ -263,6 +263,7 @@ subroutine test_directsum(ntests,npass) real :: epoti,tree_acc_prev real, allocatable :: fgrav(:,:),fxyz_ptmass_gas(:,:) + maxvxyzu = size(vxyzu(:,1)) tree_acc_prev = tree_accuracy do k = 1,6 if (labeltype(k)/='bound') then @@ -282,7 +283,6 @@ subroutine test_directsum(ntests,npass) ! call init_part() np = 1000 - maxvxyzu = size(vxyzu(:,1)) totvol = 4./3.*pi*rmax**3 nx = int(np**(1./3.)) psep = totvol**(1./3.)/real(nx) From 3caee7a95b3f01851238518456979a8ded058152 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 10:50:34 +1100 Subject: [PATCH 07/19] (test_sedov; #55) remove ifdefs --- src/tests/test_sedov.F90 | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/tests/test_sedov.F90 b/src/tests/test_sedov.F90 index 268ad74e3..a70797345 100644 --- a/src/tests/test_sedov.F90 +++ b/src/tests/test_sedov.F90 @@ -30,7 +30,7 @@ module testsedov !+ !----------------------------------------------------------------------- subroutine test_sedov(ntests,npass) - use dim, only:maxp,maxvxyzu,maxalpha,use_dust,periodic,do_radiation + use dim, only:maxp,maxvxyzu,maxalpha,use_dust,periodic,do_radiation,ind_timesteps use io, only:id,master,iprint,ievfile,iverbose,real4 use boundary, only:set_boundary,xmin,xmax,ymin,ymax,zmin,zmax,dxbound,dybound,dzbound use unifdis, only:set_unifdis @@ -43,9 +43,7 @@ subroutine test_sedov(ntests,npass) use deriv, only:get_derivs_global use timestep, only:time,tmax,dtmax,C_cour,C_force,dt,tolv,bignumber use units, only:set_units -#ifndef IND_TIMESTEPS use timestep, only:dtcourant,dtforce,dtrad -#endif use testutils, only:checkval,update_test_scores use evwrite, only:init_evfile,write_evfile use energies, only:etot,totmom,angtot,mdust @@ -67,10 +65,10 @@ subroutine test_sedov(ntests,npass) real :: temp character(len=20) :: logfile,evfile,dumpfile -#ifndef PERIODIC - if (id==master) write(*,"(/,a)") '--> SKIPPING Sedov blast wave (needs -DPERIODIC)' - return -#endif + if (.not.periodic) then + if (id==master) write(*,"(/,a)") '--> SKIPPING Sedov blast wave (needs -DPERIODIC)' + return + endif #ifdef DISC_VISCOSITY if (id==master) write(*,"(/,a)") '--> SKIPPING Sedov blast wave (cannot use -DDISC_VISCOSITY)' return @@ -154,9 +152,7 @@ subroutine test_sedov(ntests,npass) ! !--now call evolve ! -#ifndef IND_TIMESTEPS - dt = min(dtcourant,dtforce,dtrad) -#endif + if (.not.ind_timesteps) dt = min(dtcourant,dtforce,dtrad) call init_step(npart,time,dtmax) iprint = 6 logfile = 'test01.log' From 45860c0892b989a16b67f6ff2d1fd7a1f10a59d2 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 10:52:08 +1100 Subject: [PATCH 08/19] (#484) fix slow tests on github actions --- .github/workflows/mpi.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/mpi.yml b/.github/workflows/mpi.yml index 9bbf3e6df..acd28d60c 100644 --- a/.github/workflows/mpi.yml +++ b/.github/workflows/mpi.yml @@ -43,7 +43,7 @@ jobs: env: OMP_STACKSIZE: 512M - OMP_NUM_THREADS: 4 + OMP_NUM_THREADS: 2 steps: From d0c1870284fc32b6511362978856f49ab274d749 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 10:58:24 +1100 Subject: [PATCH 09/19] (interp3D) cleanup --- src/utils/interpolate3D.F90 | 315 +----------------------------------- 1 file changed, 4 insertions(+), 311 deletions(-) diff --git a/src/utils/interpolate3D.F90 b/src/utils/interpolate3D.F90 index 5e1196284..00503ad41 100644 --- a/src/utils/interpolate3D.F90 +++ b/src/utils/interpolate3D.F90 @@ -16,19 +16,14 @@ module interpolations3D ! ! :Dependencies: einsteintk_utils, kernel ! - !---------------------------------------------------------------------- ! ! Module containing all of the routines required for interpolation ! from 3D data to a 3D grid (SLOW!) ! !---------------------------------------------------------------------- - use einsteintk_utils, only:exact_rendering - use kernel, only:radkern2,radkern,cnormk,wkern!,wallint ! Moved to this module - !use interpolation, only:iroll ! Moved to this module - - !use timing, only:wall_time,print_time ! Using cpu_time for now + use kernel, only:radkern2,radkern,cnormk,wkern implicit none integer, parameter :: doub_prec = kind(0.d0) real :: cnormk3D = cnormk @@ -70,7 +65,6 @@ subroutine interpolate3D(xyzh,weight,dat,itype,npart,& integer, intent(in) :: npart,npixx,npixy,npixz real, intent(in) :: xyzh(4,npart) -!real, intent(in), dimension(npart) :: x,y,z,hh ! change to xyzh() real, intent(in), dimension(npart) :: weight,dat integer, intent(in), dimension(npart) :: itype real, intent(in) :: xmin,ymin,zmin,pixwidthx,pixwidthy,pixwidthz @@ -97,7 +91,6 @@ subroutine interpolate3D(xyzh,weight,dat,itype,npart,& !logical, parameter :: exact_rendering = .true. ! use exact rendering y/n integer :: usedpart, negflag - !$ integer :: omp_get_num_threads,omp_get_thread_num integer(kind=selected_int_kind(10)) :: iprogress,j ! up to 10 digits @@ -106,16 +99,9 @@ subroutine interpolate3D(xyzh,weight,dat,itype,npart,& y(:) = xyzh(2,:) z(:) = xyzh(3,:) hh(:) = xyzh(4,:) -!print*, "smoothing length: ", hh(1:10) -! cnormk3D set the value from the kernel routine cnormk3D = cnormk radkernel = radkern radkernel2 = radkern2 -! print*, "radkern: ", radkern -! print*, "radkernel: ",radkernel -! print*, "radkern2: ", radkern2 - -! print*, "npix: ", npixx, npixy,npixz if (exact_rendering) then print "(1x,a)",'interpolating to 3D grid (exact/Petkova+2018 on subgrid) ...' @@ -160,11 +146,6 @@ subroutine interpolate3D(xyzh,weight,dat,itype,npart,& xminpix = xmin !- 0.5*pixwidthx yminpix = ymin !- 0.5*pixwidthy zminpix = zmin !- 0.5*pixwidthz -! print*, "xminpix: ", xminpix -! print*, "yminpix: ", yminpix -! print*, "zminpix: ", zminpix -! print*, "dat: ", dat(1:10) -! print*, "weights: ", weight(1:10) pixwidthmax = max(pixwidthx,pixwidthy,pixwidthz) ! !--use a minimum smoothing length on the grid to make @@ -277,15 +258,12 @@ subroutine interpolate3D(xyzh,weight,dat,itype,npart,& ! !--precalculate an array of dx2 for this particle (optimisation) ! -! Check the x position of the grid cells -!open(unit=677,file="posxgrid.txt",action='write',position='append') nxpix = 0 do ipix=ipixmin,ipixmax nxpix = nxpix + 1 ipixi = ipix if (periodicx) ipixi = iroll(ipix,npixx) xpixi = xminpix + ipix*pixwidthx - !write(677,*) ipix, xpixi !--watch out for errors with periodic wrapping... if (nxpix <= size(dx2i)) then dx2i(nxpix) = ((xpixi - xi)**2)*hi21 @@ -415,13 +393,6 @@ subroutine interpolate3D(xyzh,weight,dat,itype,npart,& call cpu_time(t_end) t_used = t_end - t_start print*, 'Interpolate3D completed in ',t_end-t_start,'s' -!if (t_used > 10.) call print_time(t_used) - -!print*, 'Number of particles in the volume: ', usedpart -! datsmooth(1,1,1) = 3.14159 -! datsmooth(32,32,32) = 3.145159 -! datsmooth(11,11,11) = 3.14159 -! datsmooth(10,10,10) = 3.145159 end subroutine interpolate3D @@ -441,7 +412,7 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& !logical, intent(in), exact_rendering real, allocatable :: datnorm(:,:,:) - integer :: i,ipix,jpix,kpix,lockindex,smoothindex + integer :: i,ipix,jpix,kpix,smoothindex integer :: iprintinterval,iprintnext integer :: ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax integer :: ipixi,jpixi,kpixi,nxpix,nwarn,threadid @@ -561,7 +532,7 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& !$omp private(ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax) & !$omp private(ipix,jpix,kpix,ipixi,jpixi,kpixi) & !$omp private(dx2i,nxpix,zpix,dz,dz2,dyz2,dy,ypix,q2,wab) & - !$omp private(pixint,wint,negflag,dfac,threadid,lockindex,smoothindex) & + !$omp private(pixint,wint,negflag,dfac,threadid,smoothindex) & !$omp firstprivate(iprintnext) & !$omp reduction(+:nwarn,usedpart) !$omp master @@ -644,15 +615,12 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& ! !--precalculate an array of dx2 for this particle (optimisation) ! - ! Check the x position of the grid cells - !open(unit=677,file="posxgrid.txt",action='write',position='append') nxpix = 0 do ipix=ipixmin,ipixmax nxpix = nxpix + 1 ipixi = ipix if (periodicx) ipixi = iroll(ipix,npixx) xpixi = xminpix + ipix*pixwidthx - !write(677,*) ipix, xpixi !--watch out for errors with periodic wrapping... if (nxpix <= size(dx2i)) then dx2i(nxpix) = ((xpixi - xi)**2)*hi21 @@ -732,24 +700,14 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& ! !--calculate data value at this pixel using the summation interpolant ! - ! Find out where this pixel sits in the lock array - ! lockindex = (k-1)*nx*ny + (j-1)*nx + i - !lockindex = (kpixi-1)*npixx*npixy + (jpixi-1)*npixx + ipixi - !!$call omp_set_lock(ilock(lockindex)) - !!$omp critical (datsmooth) do smoothindex=1, ilendat !$omp atomic datsmooth(smoothindex,ipixi,jpixi,kpixi) = datsmooth(smoothindex,ipixi,jpixi,kpixi) + term(smoothindex)*wab enddo - !!$omp end critical (datsmooth) if (normalise) then !$omp atomic - !!$omp critical (datnorm) datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab - !!$omp end critical (datnorm) - endif - - !!$call omp_unset_lock(ilock(lockindex)) + endif endif else if (q2 < radkernel2) then @@ -761,25 +719,14 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& ! !--calculate data value at this pixel using the summation interpolant ! - !!$omp atomic ! Atomic statmements only work with scalars - !!$omp set lock ! Does this work with an array? - ! Find out where this pixel sits in the lock array - ! lockindex = (k-1)*nx*ny + (j-1)*nx + i - !lockindex = (kpixi-1)*npixx*npixy + (jpixi-1)*npixx + ipixi - !!$call omp_set_lock(ilock(lockindex)) - !!$omp critical (datsmooth) do smoothindex=1,ilendat !$omp atomic datsmooth(smoothindex,ipixi,jpixi,kpixi) = datsmooth(smoothindex,ipixi,jpixi,kpixi) + term(smoothindex)*wab enddo - !!$omp end critical (datsmooth) if (normalise) then !$omp atomic - !!$omp critical (datnorm) datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab - !!$omp end critical (datnorm) endif - !!$call omp_unset_lock(ilock(lockindex)) endif endif @@ -790,11 +737,6 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& !$omp enddo !$omp end parallel -!!$ do i=1,npixx*npixy*npixz -!!$ call omp_destroy_lock(ilock(i)) -!!$ enddo -!!$ if (allocated(ilock)) deallocate(ilock) - if (nwarn > 0) then print "(a,i11,a,/,a)",' interpolate3D: WARNING: contributions truncated from ',nwarn,' particles',& ' that wrap periodic boundaries more than once' @@ -812,261 +754,12 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& endif if (allocated(datnorm)) deallocate(datnorm) - !call wall_time(t_end) call cpu_time(t_end) t_used = t_end - t_start print*, 'Interpolate3DVec completed in ',t_end-t_start,'s' - !if (t_used > 10.) call print_time(t_used) - - !print*, 'Number of particles in the volume: ', usedpart - ! datsmooth(1,1,1) = 3.14159 - ! datsmooth(32,32,32) = 3.145159 - ! datsmooth(11,11,11) = 3.14159 - ! datsmooth(10,10,10) = 3.145159 end subroutine interpolate3D_vecexact - ! subroutine interpolate3D_vec(x,y,z,hh,weight,datvec,itype,npart,& - ! xmin,ymin,zmin,datsmooth,npixx,npixy,npixz,pixwidthx,pixwidthy,pixwidthz,& - ! normalise,periodicx,periodicy,periodicz) - - ! integer, intent(in) :: npart,npixx,npixy,npixz - ! real, intent(in), dimension(npart) :: x,y,z,hh,weight - ! real, intent(in), dimension(npart,3) :: datvec - ! integer, intent(in), dimension(npart) :: itype - ! real, intent(in) :: xmin,ymin,zmin,pixwidthx,pixwidthy,pixwidthz - ! real(doub_prec), intent(out), dimension(3,npixx,npixy,npixz) :: datsmooth - ! logical, intent(in) :: normalise,periodicx,periodicy,periodicz - ! real(doub_prec), dimension(npixx,npixy,npixz) :: datnorm - - ! integer :: i,ipix,jpix,kpix - ! integer :: iprintinterval,iprintnext - ! integer :: ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax - ! integer :: ipixi,jpixi,kpixi,nxpix,nwarn - ! real :: xminpix,yminpix,zminpix - ! real, dimension(npixx) :: dx2i - ! real :: xi,yi,zi,hi,hi1,hi21,radkern,wab,q2,const,dyz2,dz2 - ! real :: termnorm,dy,dz,ypix,zpix,xpixi,ddatnorm - ! real, dimension(3) :: term - ! !real :: t_start,t_end - ! logical :: iprintprogress - ! !$ integer :: omp_get_num_threads - ! integer(kind=selected_int_kind(10)) :: iprogress ! up to 10 digits - - ! datsmooth = 0. - ! datnorm = 0. - ! if (normalise) then - ! print "(1x,a)",'interpolating to 3D grid (normalised) ...' - ! else - ! print "(1x,a)",'interpolating to 3D grid (non-normalised) ...' - ! endif - ! if (pixwidthx <= 0. .or. pixwidthy <= 0. .or. pixwidthz <= 0.) then - ! print "(1x,a)",'interpolate3D: error: pixel width <= 0' - ! return - ! endif - ! if (any(hh(1:npart) <= tiny(hh))) then - ! print*,'interpolate3D: WARNING: ignoring some or all particles with h < 0' - ! endif - - ! ! - ! !--print a progress report if it is going to take a long time - ! ! (a "long time" is, however, somewhat system dependent) - ! ! - ! iprintprogress = (npart >= 100000) .or. (npixx*npixy > 100000) - ! !$ iprintprogress = .false. - ! ! - ! !--loop over particles - ! ! - ! iprintinterval = 25 - ! if (npart >= 1e6) iprintinterval = 10 - ! iprintnext = iprintinterval - ! ! - ! !--get starting CPU time - ! ! - ! !call cpu_time(t_start) - - ! xminpix = xmin - 0.5*pixwidthx - ! yminpix = ymin - 0.5*pixwidthy - ! zminpix = zmin - 0.5*pixwidthz - - ! const = cnormk3D ! normalisation constant (3D) - ! nwarn = 0 - - ! !$omp parallel default(none) & - ! !$omp shared(hh,z,x,y,weight,datvec,itype,datsmooth,npart) & - ! !$omp shared(xmin,ymin,zmin,radkernel,radkernel2) & - ! !$omp shared(xminpix,yminpix,zminpix,pixwidthx,pixwidthy,pixwidthz) & - ! !$omp shared(npixx,npixy,npixz,const) & - ! !$omp shared(iprintprogress,iprintinterval) & - ! !$omp shared(datnorm,normalise,periodicx,periodicy,periodicz) & - ! !$omp private(hi,xi,yi,zi,radkern,hi1,hi21) & - ! !$omp private(term,termnorm,xpixi) & - ! !$omp private(iprogress,iprintnext) & - ! !$omp private(ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax) & - ! !$omp private(ipix,jpix,kpix,ipixi,jpixi,kpixi) & - ! !$omp private(dx2i,nxpix,zpix,dz,dz2,dyz2,dy,ypix,q2,wab) & - ! !$omp reduction(+:nwarn) - ! !$omp master - ! !$ print "(1x,a,i3,a)",'Using ',omp_get_num_threads(),' cpus' - ! !$omp end master - ! ! - ! !--loop over particles - ! ! - ! !$omp do schedule (guided, 2) - ! over_parts: do i=1,npart - ! ! - ! !--report on progress - ! ! - ! if (iprintprogress) then - ! iprogress = 100*i/npart - ! if (iprogress >= iprintnext) then - ! write(*,"('(',i3,'% -',i12,' particles done)')") iprogress,i - ! iprintnext = iprintnext + iprintinterval - ! endif - ! endif - ! ! - ! !--skip particles with itype < 0 - ! ! - ! if (itype(i) < 0 .or. weight(i) < tiny(0.)) cycle over_parts - - ! hi = hh(i) - ! if (hi <= 0.) cycle over_parts - - ! ! - ! !--set kernel related quantities - ! ! - ! xi = x(i) - ! yi = y(i) - ! zi = z(i) - - ! hi1 = 1./hi - ! hi21 = hi1*hi1 - ! radkern = radkernel*hi ! radius of the smoothing kernel - ! termnorm = const*weight(i) - ! term(:) = termnorm*datvec(i,:) - ! ! - ! !--for each particle work out which pixels it contributes to - ! ! - ! ipixmin = int((xi - radkern - xmin)/pixwidthx) - ! jpixmin = int((yi - radkern - ymin)/pixwidthy) - ! kpixmin = int((zi - radkern - zmin)/pixwidthz) - ! ipixmax = int((xi + radkern - xmin)/pixwidthx) + 1 - ! jpixmax = int((yi + radkern - ymin)/pixwidthy) + 1 - ! kpixmax = int((zi + radkern - zmin)/pixwidthz) + 1 - - ! if (.not.periodicx) then - ! if (ipixmin < 1) ipixmin = 1 ! make sure they only contribute - ! if (ipixmax > npixx) ipixmax = npixx ! to pixels in the image - ! endif - ! if (.not.periodicy) then - ! if (jpixmin < 1) jpixmin = 1 - ! if (jpixmax > npixy) jpixmax = npixy - ! endif - ! if (.not.periodicz) then - ! if (kpixmin < 1) kpixmin = 1 - ! if (kpixmax > npixz) kpixmax = npixz - ! endif - ! ! - ! !--precalculate an array of dx2 for this particle (optimisation) - ! ! - ! nxpix = 0 - ! do ipix=ipixmin,ipixmax - ! nxpix = nxpix + 1 - ! ipixi = ipix - ! if (periodicx) ipixi = iroll(ipix,npixx) - ! xpixi = xminpix + ipix*pixwidthx - ! !--watch out for errors with perioic wrapping... - ! if (nxpix <= size(dx2i)) then - ! dx2i(nxpix) = ((xpixi - xi)**2)*hi21 - ! endif - ! enddo - - ! !--if particle contributes to more than npixx pixels - ! ! (i.e. periodic boundaries wrap more than once) - ! ! truncate the contribution and give warning - ! if (nxpix > npixx) then - ! nwarn = nwarn + 1 - ! ipixmax = ipixmin + npixx - 1 - ! endif - ! ! - ! !--loop over pixels, adding the contribution from this particle - ! ! - ! do kpix = kpixmin,kpixmax - ! kpixi = kpix - ! if (periodicz) kpixi = iroll(kpix,npixz) - ! zpix = zminpix + kpix*pixwidthz - ! dz = zpix - zi - ! dz2 = dz*dz*hi21 - - ! do jpix = jpixmin,jpixmax - ! jpixi = jpix - ! if (periodicy) jpixi = iroll(jpix,npixy) - ! ypix = yminpix + jpix*pixwidthy - ! dy = ypix - yi - ! dyz2 = dy*dy*hi21 + dz2 - - ! nxpix = 0 - ! do ipix = ipixmin,ipixmax - ! ipixi = ipix - ! if (periodicx) ipixi = iroll(ipix,npixx) - ! nxpix = nxpix + 1 - ! q2 = dx2i(nxpix) + dyz2 ! dx2 pre-calculated; dy2 pre-multiplied by hi21 - ! ! - ! !--SPH kernel - standard cubic spline - ! ! - ! if (q2 < radkernel2) then - ! wab = wkernel(q2) - ! ! - ! !--calculate data value at this pixel using the summation interpolant - ! ! - ! !$omp atomic - ! datsmooth(1,ipixi,jpixi,kpixi) = datsmooth(1,ipixi,jpixi,kpixi) + term(1)*wab - ! !$omp atomic - ! datsmooth(2,ipixi,jpixi,kpixi) = datsmooth(2,ipixi,jpixi,kpixi) + term(2)*wab - ! !$omp atomic - ! datsmooth(3,ipixi,jpixi,kpixi) = datsmooth(3,ipixi,jpixi,kpixi) + term(3)*wab - ! if (normalise) then - ! !$omp atomic - ! datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab - ! endif - ! endif - ! enddo - ! enddo - ! enddo - ! enddo over_parts - ! !$omp enddo - ! !$omp end parallel - - ! if (nwarn > 0) then - ! print "(a,i11,a,/,a)",' interpolate3D: WARNING: contributions truncated from ',nwarn,' particles',& - ! ' that wrap periodic boundaries more than once' - ! endif - ! ! - ! !--normalise dat array - ! ! - ! if (normalise) then - ! !$omp parallel do default(none) schedule(static) & - ! !$omp shared(datsmooth,datnorm,npixz,npixy,npixx) & - ! !$omp private(kpix,jpix,ipix,ddatnorm) - ! do kpix=1,npixz - ! do jpix=1,npixy - ! do ipix=1,npixx - ! if (datnorm(ipix,jpix,kpix) > tiny(datnorm)) then - ! ddatnorm = 1./datnorm(ipix,jpix,kpix) - ! datsmooth(1,ipix,jpix,kpix) = datsmooth(1,ipix,jpix,kpix)*ddatnorm - ! datsmooth(2,ipix,jpix,kpix) = datsmooth(2,ipix,jpix,kpix)*ddatnorm - ! datsmooth(3,ipix,jpix,kpix) = datsmooth(3,ipix,jpix,kpix)*ddatnorm - ! endif - ! enddo - ! enddo - ! enddo - ! !$omp end parallel do - ! endif - - ! return - - ! end subroutine interpolate3D_vec - !------------------------------------------------------------ ! interface to kernel routine to avoid problems with openMP !----------------------------------------------------------- From b16fb767657d3551f4bdf0cbfbf6fd661bbe50b6 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 11:01:04 +1100 Subject: [PATCH 10/19] [format-bot] obsolete .gt. .lt. .ge. .le. .eq. .ne. replaced --- src/utils/analysis_radiotde.f90 | 4 ++-- src/utils/moddump_radiotde.f90 | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/utils/analysis_radiotde.f90 b/src/utils/analysis_radiotde.f90 index 4d009de92..2dd7105b0 100644 --- a/src/utils/analysis_radiotde.f90 +++ b/src/utils/analysis_radiotde.f90 @@ -194,10 +194,10 @@ subroutine tde_analysis(npart,pmass,xyzh,vxyzu) vri = dot_product(vxyz,xyz)/r vr_accum_add = vr_accum_add + vri v_accum_add = v_accum_add + v - if (r-rad_cap < drad_cap .and. (v .ge. v_min .and. v .le. v_max)) then + if (r-rad_cap < drad_cap .and. (v >= v_min .and. v <= v_max)) then thetai = atan2d(y,x) phii = atan2d(z,sqrt(x**2+y**2)) - if ((thetai .ge. theta_min .and. thetai .le. theta_max) .and. (phii .ge. phi_min .and. phii .le. phi_max)) then + if ((thetai >= theta_min .and. thetai <= theta_max) .and. (phii >= phi_min .and. phii <= phi_max)) then m_cap = m_cap + pmass n_cap = n_cap + 1 cap(i) = .true. diff --git a/src/utils/moddump_radiotde.f90 b/src/utils/moddump_radiotde.f90 index d854923a5..22784d0a0 100644 --- a/src/utils/moddump_radiotde.f90 +++ b/src/utils/moddump_radiotde.f90 @@ -132,7 +132,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) rhof_rbreak(:) = rhof_rbreak_in(1:nbreak) call calc_rhobreak() else - if (temperature .le. 0) read_temp = .true. + if (temperature <= 0) read_temp = .true. rhof => rho_tab deallocate(rhof_n,rhof_rbreak) From 0f02abfea9c0999c8b57087f821cb58a61114d8c Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 11:02:31 +1100 Subject: [PATCH 11/19] [header-bot] updated file headers --- src/main/extern_gr.F90 | 3 +-- src/main/initial.F90 | 15 ++++++++------- src/main/inject_windtunnel.f90 | 13 +++++++++---- src/main/interp_metric.F90 | 2 +- src/main/metric_et.f90 | 2 +- src/main/metric_flrw.f90 | 2 +- src/main/readwrite_dumps_fortran.F90 | 4 ++-- src/main/tmunu2grid.f90 | 2 +- src/setup/set_star.f90 | 14 +------------- src/setup/setup_asteroidwind.f90 | 4 +++- src/setup/setup_flrw.f90 | 6 +++--- src/setup/setup_flrwpspec.f90 | 6 +++--- src/setup/setup_windtunnel.f90 | 22 ++++++++++++++++++---- src/utils/analysis_dustformation.f90 | 2 +- src/utils/analysis_radiotde.f90 | 20 ++++++++++---------- src/utils/einsteintk_utils.f90 | 2 +- src/utils/einsteintk_wrapper.f90 | 2 +- src/utils/interpolate3D.F90 | 8 +------- src/utils/interpolate3Dold.F90 | 2 +- src/utils/moddump_radiotde.f90 | 26 +++++++++++++++++--------- 20 files changed, 84 insertions(+), 73 deletions(-) diff --git a/src/main/extern_gr.F90 b/src/main/extern_gr.F90 index 939d7b301..17a80ea47 100644 --- a/src/main/extern_gr.F90 +++ b/src/main/extern_gr.F90 @@ -14,8 +14,7 @@ module extern_gr ! ! :Runtime parameters: None ! -! :Dependencies: eos, fastmath, io, metric_tools, part, physcon, timestep, -! utils_gr +! :Dependencies: eos, io, metric_tools, part, physcon, timestep, utils_gr ! implicit none diff --git a/src/main/initial.F90 b/src/main/initial.F90 index 3784e7431..35906e82a 100644 --- a/src/main/initial.F90 +++ b/src/main/initial.F90 @@ -16,13 +16,14 @@ module initial ! ! :Dependencies: analysis, boundary, boundary_dyn, centreofmass, ! checkconserved, checkoptions, checksetup, cons2prim, cooling, cpuinfo, -! damping, densityforce, deriv, dim, dust, dust_formation, energies, eos, -! evwrite, extern_gr, externalforces, fastmath, fileutils, forcing, -! growth, inject, io, io_summary, krome_interface, linklist, -! metric_tools, mf_write, mpibalance, mpidomain, mpimemory, mpitree, -! mpiutils, nicil, nicil_sup, omputils, options, part, partinject, -! ptmass, radiation_utils, readwrite_dumps, readwrite_infile, timestep, -! timestep_ind, timestep_sts, timing, units, writeheader +! damping, densityforce, deriv, dim, dust, dust_formation, +! einsteintk_utils, energies, eos, evwrite, extern_gr, externalforces, +! fastmath, fileutils, forcing, growth, inject, io, io_summary, +! krome_interface, linklist, metric_tools, mf_write, mpibalance, +! mpidomain, mpimemory, mpitree, mpiutils, nicil, nicil_sup, omputils, +! options, part, partinject, ptmass, radiation_utils, readwrite_dumps, +! readwrite_infile, timestep, timestep_ind, timestep_sts, timing, +! tmunu2grid, units, writeheader ! implicit none diff --git a/src/main/inject_windtunnel.f90 b/src/main/inject_windtunnel.f90 index 7c304db07..0245f9337 100644 --- a/src/main/inject_windtunnel.f90 +++ b/src/main/inject_windtunnel.f90 @@ -9,14 +9,19 @@ module inject ! Handles injection for gas sphere in wind tunnel ! ! +! :References: None +! ! :Owner: Mike Lau ! ! :Runtime parameters: -! - lattice_type : *0: cubic distribution, 1: closepacked distribution* -! - handled_layers : *(integer) number of handled BHL wind layers* -! - v_inf : *BHL wind speed* -! - Rstar : *BHL star radius (in accretion radii)* ! - BHL_radius : *radius of the wind cylinder (in star radii)* +! - Rstar : *sphere radius (code units)* +! - handled_layers : *(integer) number of handled BHL wind layers* +! - lattice_type : *0: cubic distribution, 1: closepacked distribution* +! - nstar : *No. of particles making up sphere* +! - pres_inf : *ambient pressure (code units)* +! - rho_inf : *ambient density (code units)* +! - v_inf : *wind speed (code units)* ! - wind_injection_x : *x position of the wind injection boundary (in star radii)* ! - wind_length : *crude wind length (in star radii)* ! diff --git a/src/main/interp_metric.F90 b/src/main/interp_metric.F90 index fc4dd62bf..0d1cb7080 100644 --- a/src/main/interp_metric.F90 +++ b/src/main/interp_metric.F90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module metric_interp ! diff --git a/src/main/metric_et.f90 b/src/main/metric_et.f90 index ca792fc92..a15d185e6 100644 --- a/src/main/metric_et.f90 +++ b/src/main/metric_et.f90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module metric ! diff --git a/src/main/metric_flrw.f90 b/src/main/metric_flrw.f90 index 68152b86d..3685131b8 100644 --- a/src/main/metric_flrw.f90 +++ b/src/main/metric_flrw.f90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module metric ! diff --git a/src/main/readwrite_dumps_fortran.F90 b/src/main/readwrite_dumps_fortran.F90 index cbe7212ad..696128036 100644 --- a/src/main/readwrite_dumps_fortran.F90 +++ b/src/main/readwrite_dumps_fortran.F90 @@ -20,8 +20,8 @@ module readwrite_dumps_fortran ! ! :Dependencies: boundary, boundary_dyn, checkconserved, dim, dump_utils, ! dust, dust_formation, eos, externalforces, fileutils, io, lumin_nsdisc, -! memory, mpi, mpiutils, options, part, readwrite_dumps_common, -! setup_params, sphNGutils, timestep, units +! memory, metric_tools, mpi, mpiutils, options, part, +! readwrite_dumps_common, setup_params, sphNGutils, timestep, units ! use dump_utils, only:lenid,ndatatypes,i_int,i_int1,i_int2,i_int4,i_int8,& i_real,i_real4,i_real8,int1,int2,int1o,int2o,dump_h,lentag diff --git a/src/main/tmunu2grid.f90 b/src/main/tmunu2grid.f90 index c2ff7ab27..21b5e620b 100644 --- a/src/main/tmunu2grid.f90 +++ b/src/main/tmunu2grid.f90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module tmunu2grid ! diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index 66a5d476c..24071d154 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -15,19 +15,7 @@ module setstar ! ! :Owner: Daniel Price ! -! :Runtime parameters: -! - Mstar : *mass of star* -! - Rstar : *radius of star* -! - hsoft : *Softening length of sink particle stellar core* -! - input_profile : *Path to input profile* -! - isinkcore : *Add a sink particle stellar core* -! - isoftcore : *0=no core softening, 1=cubic core, 2=constant entropy core* -! - isofteningopt : *1=supply rcore, 2=supply mcore, 3=supply both* -! - mcore : *Mass of sink particle stellar core* -! - np : *number of particles* -! - outputfilename : *Output path for softened MESA profile* -! - rcore : *Radius of core softening* -! - ui_coef : *specific internal energy (units of GM/R)* +! :Runtime parameters: None ! ! :Dependencies: centreofmass, dim, eos, extern_densprofile, infile_utils, ! io, mpiutils, part, physcon, prompting, radiation_utils, relaxstar, diff --git a/src/setup/setup_asteroidwind.f90 b/src/setup/setup_asteroidwind.f90 index aff62f942..44f098ea0 100644 --- a/src/setup/setup_asteroidwind.f90 +++ b/src/setup/setup_asteroidwind.f90 @@ -22,13 +22,15 @@ module setup ! - ipot : *wd modelled by 0=sink or 1=externalforce* ! - m1 : *mass of white dwarf (solar mass)* ! - m2 : *mass of asteroid (ceres mass)* +! - mdot : *mass injection rate (g/s)* ! - norbits : *number of orbits* ! - npart_at_end : *number of particles injected after norbits* ! - rasteroid : *radius of asteroid (km)* ! - semia : *semi-major axis (solar radii)* ! ! :Dependencies: eos, extern_lensethirring, externalforces, infile_utils, -! io, options, part, physcon, setbinary, spherical, timestep, units +! inject, io, kernel, options, part, physcon, setbinary, spherical, +! timestep, units ! use inject, only:mdot implicit none diff --git a/src/setup/setup_flrw.f90 b/src/setup/setup_flrw.f90 index 7cdc8c868..0d577a3df 100644 --- a/src/setup/setup_flrw.f90 +++ b/src/setup/setup_flrw.f90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module setup ! @@ -23,8 +23,8 @@ module setup ! - rhozero : *initial density in code units* ! ! :Dependencies: boundary, dim, infile_utils, io, mpidomain, mpiutils, -! options, part, physcon, prompting, setup_params, stretchmap, unifdis, -! units, utils_gr +! part, physcon, prompting, setup_params, stretchmap, unifdis, units, +! utils_gr ! use dim, only:use_dust use setup_params, only:rhozero diff --git a/src/setup/setup_flrwpspec.f90 b/src/setup/setup_flrwpspec.f90 index eef12efc8..ff8db10e9 100644 --- a/src/setup/setup_flrwpspec.f90 +++ b/src/setup/setup_flrwpspec.f90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module setup ! @@ -24,8 +24,8 @@ module setup ! - rhozero : *initial density in code units* ! ! :Dependencies: boundary, dim, eos_shen, infile_utils, io, mpidomain, -! mpiutils, options, part, physcon, prompting, setup_params, stretchmap, -! unifdis, units, utils_gr +! mpiutils, part, physcon, prompting, setup_params, stretchmap, unifdis, +! units, utils_gr ! use dim, only:use_dust use setup_params, only:rhozero diff --git a/src/setup/setup_windtunnel.f90 b/src/setup/setup_windtunnel.f90 index 2c8e20ba4..29251e287 100644 --- a/src/setup/setup_windtunnel.f90 +++ b/src/setup/setup_windtunnel.f90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module setup ! @@ -10,11 +10,25 @@ module setup ! ! :References: None ! -! :Owner: Daniel Price +! :Owner: Mike Lau ! -! :Runtime parameters: None +! :Runtime parameters: +! - Mstar : *sphere mass in code units* +! - Rstar : *sphere radius in code units* +! - gamma : *adiabatic index* +! - handled_layers : *number of handled layers* +! - lattice_type : *0: cubic, 1: close-packed cubic* +! - nstar : *number of particles resolving gas sphere* +! - pres_inf : *wind pressure / dyn cm^2* +! - rho_inf : *wind density / g cm^-3* +! - v_inf : *wind speed / km s^-1* +! - wind_injection_x : *injection x in units of Rstar* +! - wind_length : *wind length in units of Rstar* +! - wind_radius : *injection radius in units of Rstar* ! -! :Dependencies: inject, part, physcon, units +! :Dependencies: dim, eos, extern_densprofile, infile_utils, inject, io, +! kernel, mpidomain, part, physcon, rho_profile, setstar_utils, setunits, +! setup_params, table_utils, timestep, unifdis, units ! use io, only:master,fatal use inject, only:init_inject,nstar,Rstar,lattice_type,handled_layers,& diff --git a/src/utils/analysis_dustformation.f90 b/src/utils/analysis_dustformation.f90 index 489c57789..353a39b1b 100644 --- a/src/utils/analysis_dustformation.f90 +++ b/src/utils/analysis_dustformation.f90 @@ -15,7 +15,7 @@ module analysis ! ! :Runtime parameters: None ! -! :Dependencies: None +! :Dependencies: dim, fileutils, part, physcon, units ! implicit none character(len=20), parameter, public :: analysistype = 'dustformation' diff --git a/src/utils/analysis_radiotde.f90 b/src/utils/analysis_radiotde.f90 index 2dd7105b0..cd17076f9 100644 --- a/src/utils/analysis_radiotde.f90 +++ b/src/utils/analysis_radiotde.f90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module analysis ! @@ -10,17 +10,17 @@ module analysis ! ! :References: None ! -! :Owner: Fitz Hu +! :Owner: fhu ! ! :Runtime parameters: -! - rad_cap : *capture shell radius* -! - drad_cap : *capture shell thickness* -! - v_max : *max velocity* -! - v_min : *min velocity* -! - theta_max : *max azimuthal angle* -! - theta_min : *min azimuthal angle* -! - phi_max : *max altitude angle* -! - phi_min : *min altitude angle* +! - drad_cap : *capture thickness (in cm) (-ve for all particles at outer radius)* +! - phi_max : *max phi (in deg)* +! - phi_min : *min phi (in deg)* +! - rad_cap : *capture inner radius (in cm)* +! - theta_max : *max theta (in deg)* +! - theta_min : *min theta (in deg)* +! - v_max : *max velocity (in c)* +! - v_min : *min velocity (in c)* ! ! :Dependencies: infile_utils, io, physcon, readwrite_dumps, units ! diff --git a/src/utils/einsteintk_utils.f90 b/src/utils/einsteintk_utils.f90 index 880ac3096..7d436fd0a 100644 --- a/src/utils/einsteintk_utils.f90 +++ b/src/utils/einsteintk_utils.f90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module einsteintk_utils ! diff --git a/src/utils/einsteintk_wrapper.f90 b/src/utils/einsteintk_wrapper.f90 index f7b5282e2..7541e5974 100644 --- a/src/utils/einsteintk_wrapper.f90 +++ b/src/utils/einsteintk_wrapper.f90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module einsteintk_wrapper ! diff --git a/src/utils/interpolate3D.F90 b/src/utils/interpolate3D.F90 index 00503ad41..e97ee8c4a 100644 --- a/src/utils/interpolate3D.F90 +++ b/src/utils/interpolate3D.F90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module interpolations3D ! @@ -16,12 +16,6 @@ module interpolations3D ! ! :Dependencies: einsteintk_utils, kernel ! -!---------------------------------------------------------------------- -! -! Module containing all of the routines required for interpolation -! from 3D data to a 3D grid (SLOW!) -! -!---------------------------------------------------------------------- use einsteintk_utils, only:exact_rendering use kernel, only:radkern2,radkern,cnormk,wkern implicit none diff --git a/src/utils/interpolate3Dold.F90 b/src/utils/interpolate3Dold.F90 index d1344fd96..c7fff7ca7 100644 --- a/src/utils/interpolate3Dold.F90 +++ b/src/utils/interpolate3Dold.F90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module interpolations3D ! diff --git a/src/utils/moddump_radiotde.f90 b/src/utils/moddump_radiotde.f90 index 22784d0a0..fa8c8ad96 100644 --- a/src/utils/moddump_radiotde.f90 +++ b/src/utils/moddump_radiotde.f90 @@ -2,7 +2,7 @@ ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! ! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module moddump ! @@ -10,17 +10,25 @@ module moddump ! ! :References: None ! -! :Owner: Fitz Hu +! :Owner: fhu ! ! :Runtime parameters: -! - temperature : *Temperature* -! - mu : *mean molecular mass* -! - ieos_in : *equation of state* -! - use_func : *use broken power law or profile date points* +! - ieos : *equation of state used* +! - ignore_radius : *tde particle inside this radius will be ignored* +! - mu : *mean molecular density of the cloud* +! - nbreak : *number of broken power laws* +! - nprof : *number of data points in the cloud profile* +! - profile_filename : *filename for the cloud profile* +! - rad_max : *outer radius of the circumnuclear gas cloud* +! - rad_min : *inner radius of the circumnuclear gas cloud* +! - remove_overlap : *remove outflow particles overlap with circum particles* +! - rhof_n_1 : *power law index of the section* +! - rhof_rho0 : *density at rad_min (in g/cm^3)* +! - temperature : *temperature of the gas cloud (-ve = read from file)* +! - use_func : *if use broken power law for density profile* ! -! :Dependencies: datafiles, eos, io, stretchmap, kernel, -! mpidomain, part, physcon, setup_params, -! spherical, timestep, units, infile_utils +! :Dependencies: eos, infile_utils, io, kernel, mpidomain, part, physcon, +! setup_params, spherical, stretchmap, timestep, units ! implicit none public :: modify_dump From 2a8b7784bece7acaf8a16b0884250b82e49df1f5 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 11:02:47 +1100 Subject: [PATCH 12/19] [space-bot] whitespace at end of lines removed --- src/main/config.F90 | 2 +- src/main/extern_gr.F90 | 8 +++--- src/main/inject_windtunnel.f90 | 10 +++---- src/main/metric_et.f90 | 2 +- src/main/readwrite_dumps_fortran.F90 | 4 +-- src/main/step_leapfrog.F90 | 2 +- src/main/tmunu2grid.f90 | 16 ++++++------ src/setup/setup_flrw.f90 | 16 ++++++------ src/setup/setup_flrwpspec.f90 | 36 +++++++++++++------------- src/setup/setup_windtunnel.f90 | 14 +++++----- src/utils/analysis_common_envelope.f90 | 2 +- src/utils/analysis_radiotde.f90 | 8 +++--- src/utils/einsteintk_wrapper.f90 | 28 ++++++++++---------- src/utils/interpolate3D.F90 | 18 ++++++------- src/utils/moddump_radiotde.f90 | 12 ++++----- 15 files changed, 89 insertions(+), 89 deletions(-) diff --git a/src/main/config.F90 b/src/main/config.F90 index 561adf30e..020e43a42 100644 --- a/src/main/config.F90 +++ b/src/main/config.F90 @@ -277,7 +277,7 @@ module dim logical, parameter :: nr = .true. #else logical, parameter :: nr = .false. -#endif +#endif !-------------------- ! Supertimestepping diff --git a/src/main/extern_gr.F90 b/src/main/extern_gr.F90 index 17a80ea47..b118c17b6 100644 --- a/src/main/extern_gr.F90 +++ b/src/main/extern_gr.F90 @@ -86,9 +86,9 @@ subroutine dt_grforce(xyzh,fext,dtf) real, intent(out) :: dtf real :: r,r2,dtf1,dtf2,f2i integer, parameter :: steps_per_orbit = 100 - + f2i = fext(1)*fext(1) + fext(2)*fext(2) + fext(3)*fext(3) - if (f2i > 0.) then + if (f2i > 0.) then dtf1 = sqrt(xyzh(4)/sqrt(f2i)) ! This is not really accurate since fi is a component of dp/dt, not da/dt else dtf1 = huge(dtf1) @@ -99,9 +99,9 @@ subroutine dt_grforce(xyzh,fext,dtf) r2 = xyzh(1)*xyzh(1) + xyzh(2)*xyzh(2) + xyzh(3)*xyzh(3) r = sqrt(r2) dtf2 = (2.*pi*sqrt(r*r2))/steps_per_orbit - case default + case default dtf2 = huge(dtf2) - end select + end select dtf = min(dtf1,dtf2) diff --git a/src/main/inject_windtunnel.f90 b/src/main/inject_windtunnel.f90 index 0245f9337..5888f288e 100644 --- a/src/main/inject_windtunnel.f90 +++ b/src/main/inject_windtunnel.f90 @@ -89,7 +89,7 @@ subroutine init_inject(ierr) if (lattice_type == 1) then psep = (sqrt(2.)*element_volume)**(1./3.) elseif (lattice_type == 0) then - psep = element_volume**(1./3.) + psep = element_volume**(1./3.) else call fatal("init_inject",'unknown lattice_type (must be 0 or 1)') endif @@ -268,16 +268,16 @@ subroutine print_summary(v_inf,cs_inf,rho_inf,pres_inf,mach,pmass,distance_betwe integer, intent(in) :: max_layers,nstar,max_particles print*, 'wind speed: ',v_inf * unit_velocity / 1e5," km s^-1" - print*, 'wind cs: ',cs_inf * unit_velocity / 1e5," km s^-1" - print*, 'wind density: ',rho_inf * unit_density," g cm^-3" - print*, 'wind pressure: ',pres_inf * unit_pressure," dyn cm^-2" + print*, 'wind cs: ',cs_inf * unit_velocity / 1e5," km s^-1" + print*, 'wind density: ',rho_inf * unit_density," g cm^-3" + print*, 'wind pressure: ',pres_inf * unit_pressure," dyn cm^-2" print*, 'wind mach number: ', mach print*, 'maximum wind layers: ', max_layers print*, 'pmass: ',pmass print*, 'nstar: ',nstar print*, 'nstar + max. wind particles: ', max_particles - print*, 'distance_between_layers: ',distance_between_layers + print*, 'distance_between_layers: ',distance_between_layers print*, 'time_between_layers: ',time_between_layers print*, 'planet crossing time: ',2*Rstar/v_inf diff --git a/src/main/metric_et.f90 b/src/main/metric_et.f90 index a15d185e6..d13454ce1 100644 --- a/src/main/metric_et.f90 +++ b/src/main/metric_et.f90 @@ -20,7 +20,7 @@ module metric character(len=*), parameter :: metric_type = 'et' integer, parameter :: imetric = 6 ! This are dummy parameters to stop the compiler complaing - ! Not used anywhere in the code - Needs a fix! + ! Not used anywhere in the code - Needs a fix! real, public :: mass1 = 1. ! mass of central object real, public :: a = 0.0 ! spin of central object contains diff --git a/src/main/readwrite_dumps_fortran.F90 b/src/main/readwrite_dumps_fortran.F90 index 696128036..03b755cc4 100644 --- a/src/main/readwrite_dumps_fortran.F90 +++ b/src/main/readwrite_dumps_fortran.F90 @@ -218,7 +218,7 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG) rad,rad_label,radprop,radprop_label,do_radiation,maxirad,maxradprop,itemp,igasP,igamma,& iorig,iX,iZ,imu,nucleation,nucleation_label,n_nucleation,tau,itau_alloc,tau_lucy,itauL_alloc,& luminosity,eta_nimhd,eta_nimhd_label - use part, only:metrics,metricderivs,tmunus + use part, only:metrics,metricderivs,tmunus use options, only:use_dustfrac,use_var_comp,icooling use dump_utils, only:tag,open_dumpfile_w,allocate_header,& free_header,write_header,write_array,write_block_header @@ -365,7 +365,7 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG) endif if (gr) then call write_array(1,pxyzu,pxyzu_label,maxvxyzu,npart,k,ipass,idump,nums,ierrs(8)) - call write_array(1,dens,'dens prim',npart,k,ipass,idump,nums,ierrs(8)) + call write_array(1,dens,'dens prim',npart,k,ipass,idump,nums,ierrs(8)) if (imetric==imet_et) then ! Output metric if imetric=iet call write_array(1,metrics(1,1,1,:), 'gtt (covariant)',npart,k,ipass,idump,nums,ierrs(8)) diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 68abc8bc9..5f8e468b8 100644 --- a/src/main/step_leapfrog.F90 +++ b/src/main/step_leapfrog.F90 @@ -659,7 +659,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) endif endif enddo iterations - + ! MPI reduce summary variables nwake = int(reduceall_mpi('+', nwake)) nvfloorp = int(reduceall_mpi('+', nvfloorp)) diff --git a/src/main/tmunu2grid.f90 b/src/main/tmunu2grid.f90 index 21b5e620b..2777e9a7d 100644 --- a/src/main/tmunu2grid.f90 +++ b/src/main/tmunu2grid.f90 @@ -104,10 +104,10 @@ subroutine get_tmunugrid_all(npart,xyzh,vxyzu,tmunus,calc_cfac) tmunugrid = 0. datsmooth = 0. - ! Vectorized tmunu calculation - + ! Vectorized tmunu calculation + ! Put tmunu into an array of form - ! tmunu(npart,16) + ! tmunu(npart,16) do k=1, 4 do j=1,4 do i=1,npart @@ -116,8 +116,8 @@ subroutine get_tmunugrid_all(npart,xyzh,vxyzu,tmunus,calc_cfac) ! print*, "Index in array is: ", (k-1)*4 + j ! print*,tmunus(k,j,1) dat(i, (k-1)*4 + j) = tmunus(k,j,i) - enddo - enddo + enddo + enddo enddo !stop ilendat = 16 @@ -139,17 +139,17 @@ subroutine get_tmunugrid_all(npart,xyzh,vxyzu,tmunus,calc_cfac) !print*, datsmooth((i-1)*4 + j, 10,10,10) enddo enddo -!stop +!stop ! do k=1,4 ! do j=1,4 ! do i=1,4 ! print*, "Lock index is: ", (k-1)*16+ (j-1)*4 + i ! enddo ! enddo -! enddo +! enddo ! tmunugrid(0,0,:,:,:) = datsmooth(1,:,:,:) - + ! TODO Unroll this loop for speed + using symmetries ! Possiblly cleanup the messy indexing ! do k=1,4 diff --git a/src/setup/setup_flrw.f90 b/src/setup/setup_flrw.f90 index 0d577a3df..e16173d2f 100644 --- a/src/setup/setup_flrw.f90 +++ b/src/setup/setup_flrw.f90 @@ -6,8 +6,8 @@ !--------------------------------------------------------------------------! module setup ! -! Setup routine for a constant density + petrubtations FLRW universe -! as described in Magnall et al. 2023 +! Setup routine for a constant density + petrubtations FLRW universe +! as described in Magnall et al. 2023 ! ! :References: None ! @@ -83,7 +83,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, procedure(rho_func), pointer :: density_func density_func => rhofunc ! desired density function - + ! !--general parameters @@ -97,7 +97,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! irrelevant for gamma = 4./3. endif - + ! ! default units ! @@ -139,7 +139,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! We assume ainit = 1, but this may not always be the case c1 = 1./(4.*pi*rhozero) !c2 = We set g(x^i) = 0 as we only want to extract the growing mode - c3 = - sqrt(1./(6.*pi*rhozero)) + c3 = - sqrt(1./(6.*pi*rhozero)) !c3 = hub/(4.d0*PI*rhozero) @@ -189,7 +189,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! general parameters ! ! time should be read in from the par file - !time = 0.08478563386065302 + !time = 0.08478563386065302 time = 0.18951066686763596 ! z~1000 lambda = perturb_wavelength*length kwave = (2.d0*pi)/lambda @@ -200,10 +200,10 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, select case(radiation_dominated) case('"yes"') - ! Set a value of rho_matter + ! Set a value of rho_matter rho_matter = 1.e-40 !rhozero = rhozero - radconst*last_scattering_temp**4 - ! Solve for temperature + ! Solve for temperature last_scattering_temp = ((rhozero-rho_matter)/radconst)**(1./4.) rhozero = rho_matter end select diff --git a/src/setup/setup_flrwpspec.f90 b/src/setup/setup_flrwpspec.f90 index ff8db10e9..f493a2766 100644 --- a/src/setup/setup_flrwpspec.f90 +++ b/src/setup/setup_flrwpspec.f90 @@ -7,8 +7,8 @@ module setup ! ! Setup routine for realistic cosmological initial conditions based -! on the Zeldovich approximation. -! Requries velocity files generated from a powerspectrum. +! on the Zeldovich approximation. +! Requries velocity files generated from a powerspectrum. ! ! :References: None ! @@ -82,7 +82,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, real :: scale_factor,gradphi(3),vxyz(3),dxgrid,gridorigin integer :: nghost, gridres, gridsize real, allocatable :: vxgrid(:,:,:),vygrid(:,:,:),vzgrid(:,:,:) - + ! !--general parameters ! @@ -117,7 +117,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, perturb = '"no"' perturb_direction = '"none"' radiation_dominated = '"no"' - ampl = 0. + ampl = 0. ! Ideally this should read the values of the box length ! and initial Hubble parameter from the par file. @@ -141,7 +141,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, isperiodic = .true. ncross = 0 - + ! Approx Temp of the CMB in Kelvins last_scattering_temp = 3000 @@ -188,8 +188,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, call set_units(dist=udist,mass=umass,G=1.) endif call set_boundary(xmini,xmaxi,ymini,ymaxi,zmini,zmaxi) - - + + allocate(vxgrid(gridsize,gridsize,gridsize)) allocate(vygrid(gridsize,gridsize,gridsize)) allocate(vzgrid(gridsize,gridsize,gridsize)) @@ -238,14 +238,14 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, pspec_filename1 = 'init_vel1_64.dat' pspec_filename2 = 'init_vel2_64.dat' pspec_filename3 = 'init_vel3_64.dat' - + ! Check if files exist otherwise skip and return flat space - if (.not. check_files(pspec_filename1,pspec_filename2,pspec_filename3)) then + if (.not. check_files(pspec_filename1,pspec_filename2,pspec_filename3)) then print*, "Velocity files not found..." print*, "Setting up flat space!" - return - endif - + return + endif + ! Read in velocities from vel file here ! Should be made into a function at some point @@ -359,7 +359,7 @@ subroutine setup_interactive(id,polyk) call prompt(' enter sound speed in code units (sets polyk)',cs0,0.) endif call bcast_mpi(cs0) - + ! ! type of lattice ! @@ -597,11 +597,11 @@ logical function check_files(file1,file2,file3) inquire(file=file1,exist=file1_exist) inquire(file=file2,exist=file2_exist) - inquire(file=file3,exist=file3_exist) - - if ((.not. file1_exist) .or. (.not. file2_exist) .or. (.not. file3_exist)) then - check_files = .false. - endif + inquire(file=file3,exist=file3_exist) + + if ((.not. file1_exist) .or. (.not. file2_exist) .or. (.not. file3_exist)) then + check_files = .false. + endif end function check_files end module setup diff --git a/src/setup/setup_windtunnel.f90 b/src/setup/setup_windtunnel.f90 index 29251e287..91e0ce7c6 100644 --- a/src/setup/setup_windtunnel.f90 +++ b/src/setup/setup_windtunnel.f90 @@ -76,7 +76,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, logical :: use_exactN,setexists character(len=30) :: lattice character(len=120) :: setupfile - + call set_units(mass=solarm,dist=solarr,G=1.) ! ! Initialise parameters, including those that will not be included in *.setup @@ -130,7 +130,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, call write_setupfile(setupfile) stop 'please check and edit .setup file and rerun phantomsetup' endif - + pmass = Mstar / real(nstar) massoftype(igas) = pmass call check_setup(pmass,ierr) @@ -160,9 +160,9 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, presi = yinterp(pres(1:npts),r(1:npts),ri) vxyzu(4,i) = presi / ( (gamma-1.) * densi) enddo - + deallocate(r,den,pres) - + print*, "udist = ", udist, "; umass = ", umass, "; utime = ", utime end subroutine setpart @@ -194,7 +194,7 @@ subroutine write_setupfile(filename) call write_inopt(nstar,'nstar','number of particles resolving gas sphere',iunit) call write_inopt(Mstar,'Mstar','sphere mass in code units',iunit) call write_inopt(Rstar,'Rstar','sphere radius in code units',iunit) - + write(iunit,"(/,a)") '# wind settings' call write_inopt(v_inf*unit_velocity/1.e5,'v_inf','wind speed / km s^-1',iunit) call write_inopt(rho_inf*unit_density,'rho_inf','wind density / g cm^-3',iunit) @@ -295,6 +295,6 @@ subroutine check_setup(pmass,ierr) endif end subroutine check_setup - + end module setup - + diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 6cd1c6c27..c64b92bbe 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -1167,7 +1167,7 @@ subroutine roche_lobe_values(time,npart,particlemass,xyzh,vxyzu) else MRL(iR1T) = MRL(iR1T) / real(nR1T) endif - + if (nFB == 0) then MRL(iFBV) = 0 else diff --git a/src/utils/analysis_radiotde.f90 b/src/utils/analysis_radiotde.f90 index cd17076f9..d4e99725c 100644 --- a/src/utils/analysis_radiotde.f90 +++ b/src/utils/analysis_radiotde.f90 @@ -151,7 +151,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) v_cap_mean, & e_accum*unit_energ, & e_cap*unit_energ - close(iunit) + close(iunit) write(*,'(I8,1X,A2,1X,I8,1X,A34)') n_cap, 'of', npart, 'particles are in the capture shell' write(*,'(I8,1X,A2,1X,I8,1X,A40)') n_accum, 'of', npart, 'particles are outside the capture radius' @@ -175,7 +175,7 @@ subroutine tde_analysis(npart,pmass,xyzh,vxyzu) vr_cap_add = 0. v_accum_add = 0. v_cap_add = 0. - + do i = 1,npart x = xyzh(1,i) y = xyzh(2,i) @@ -188,7 +188,7 @@ subroutine tde_analysis(npart,pmass,xyzh,vxyzu) r = sqrt(dot_product(xyz,xyz)) v = sqrt(dot_product(vxyz,vxyz)) if (r > rad_cap) then - m_accum = m_accum + pmass + m_accum = m_accum + pmass n_accum = n_accum + 1 e_accum = e_accum + 0.5*pmass*v**2 vri = dot_product(vxyz,xyz)/r @@ -264,7 +264,7 @@ subroutine read_tdeparams(filename,ierr) nerr = 0 ierr = 0 call open_db_from_file(db,filename,iunit,ierr) - + call read_inopt(rad_cap,'rad_cap',db,min=0.,errcount=nerr) call read_inopt(drad_cap,'drad_cap',db,errcount=nerr) diff --git a/src/utils/einsteintk_wrapper.f90 b/src/utils/einsteintk_wrapper.f90 index 7541e5974..4414c3142 100644 --- a/src/utils/einsteintk_wrapper.f90 +++ b/src/utils/einsteintk_wrapper.f90 @@ -239,7 +239,7 @@ subroutine et2phantom_tmunu() ! Correct Tmunu ! Convert to 8byte real to stop compiler warning tmunugrid = real(cfac)*tmunugrid - + end subroutine et2phantom_tmunu @@ -286,7 +286,7 @@ subroutine phantom2et_consvar() ! Correct momentum and Density ! Conversion of cfac to 8byte real to avoid - ! compiler warning + ! compiler warning rhostargrid = real(cfac)*rhostargrid pxgrid = real(cfac)*pxgrid entropygrid = real(cfac)*entropygrid @@ -426,37 +426,37 @@ subroutine et2phantom_dumphydro(time,dt_et,checkpointfile) real, intent(in) :: time, dt_et !real(kind=16) :: cfac !logical, intent(in), optional :: checkpoint - !integer, intent(in) :: checkpointno + !integer, intent(in) :: checkpointno character(*),optional, intent(in) :: checkpointfile - logical :: createcheckpoint + logical :: createcheckpoint - if (present(checkpointfile)) then + if (present(checkpointfile)) then createcheckpoint = .true. - else + else createcheckpoint = .false. - endif + endif ! Write EV_file - if (.not. createcheckpoint) then + if (.not. createcheckpoint) then call write_evfile(time,dt_et) evfilestor = getnextfilename(evfilestor) logfilestor = getnextfilename(logfilestor) dumpfilestor = getnextfilename(dumpfilestor) call write_fulldump(time,dumpfilestor) - endif + endif ! Write full dump - if (createcheckpoint) then - call write_fulldump(time,checkpointfile) - endif + if (createcheckpoint) then + call write_fulldump(time,checkpointfile) + endif ! Quick and dirty write cfac to txtfile - + ! Density check vs particles ! call check_conserved_dens(rhostargrid,cfac) ! open(unit=777, file="cfac.txt", action='write', position='append') -! print*, time, cfac +! print*, time, cfac ! write(777,*) time, cfac ! close(unit=777) diff --git a/src/utils/interpolate3D.F90 b/src/utils/interpolate3D.F90 index e97ee8c4a..80789b06f 100644 --- a/src/utils/interpolate3D.F90 +++ b/src/utils/interpolate3D.F90 @@ -695,13 +695,13 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& !--calculate data value at this pixel using the summation interpolant ! do smoothindex=1, ilendat - !$omp atomic + !$omp atomic datsmooth(smoothindex,ipixi,jpixi,kpixi) = datsmooth(smoothindex,ipixi,jpixi,kpixi) + term(smoothindex)*wab - enddo + enddo if (normalise) then !$omp atomic datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab - endif + endif endif else if (q2 < radkernel2) then @@ -714,14 +714,14 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& !--calculate data value at this pixel using the summation interpolant ! do smoothindex=1,ilendat - !$omp atomic + !$omp atomic datsmooth(smoothindex,ipixi,jpixi,kpixi) = datsmooth(smoothindex,ipixi,jpixi,kpixi) + term(smoothindex)*wab - enddo + enddo if (normalise) then !$omp atomic datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab endif - + endif endif enddo @@ -741,10 +741,10 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& if (normalise) then do i=1, ilendat where (datnorm > tiny(datnorm)) - - datsmooth(i,:,:,:) = datsmooth(i,:,:,:)/datnorm(:,:,:) + + datsmooth(i,:,:,:) = datsmooth(i,:,:,:)/datnorm(:,:,:) end where - enddo + enddo endif if (allocated(datnorm)) deallocate(datnorm) diff --git a/src/utils/moddump_radiotde.f90 b/src/utils/moddump_radiotde.f90 index fa8c8ad96..5e12a738e 100644 --- a/src/utils/moddump_radiotde.f90 +++ b/src/utils/moddump_radiotde.f90 @@ -90,7 +90,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) temperature = 10. ! Temperature in Kelvin mu = 2. ! mean molecular weight ieos_in = 2 - ignore_radius = 1.e14 ! in cm + ignore_radius = 1.e14 ! in cm use_func = .true. use_func_old = use_func remove_overlap = .true. @@ -141,7 +141,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) call calc_rhobreak() else if (temperature <= 0) read_temp = .true. - rhof => rho_tab + rhof => rho_tab deallocate(rhof_n,rhof_rbreak) allocate(dens_prof(nprof),rad_prof(nprof),masstab(nprof)) @@ -199,17 +199,17 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) call set_sphere('random',id,master,rad_min,rad_max,delta,hfact_default,npart,xyzh, & rhofunc=rhof,nptot=npart_total,exactN=.true.,np_requested=np_sphere,mask=i_belong) if (ierr /= 0) call fatal('moddump','error setting up the circumnuclear gas cloud') - + npartoftype(igas) = npart !--Set particle properties do i = npart_old+1,npart call set_particle_type(i,igas) r = dot_product(xyzh(1:3,i),xyzh(1:3,i)) - if (read_temp) temperature = get_temp_r(r,rad_prof,temp_prof) + if (read_temp) temperature = get_temp_r(r,rad_prof,temp_prof) vxyzu(4,i) = uerg(rhof(r),temperature) vxyzu(1:3,i) = 0. ! stationary for now enddo - + !--Set timesteps tmax = 10.*years/utime dtmax = tmax/1000. @@ -318,7 +318,7 @@ subroutine write_setupfile(filename) integer, parameter :: iunit = 20 integer :: i character(len=20) :: rstr,nstr - + write(*,"(a)") ' writing setup options file '//trim(filename) open(unit=iunit,file=filename,status='replace',form='formatted') write(iunit,"(a)") '# input file for setting up a circumnuclear gas cloud' From 2c73a466619a99235c8ce69ce09f32ca936a5d6b Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 11:02:47 +1100 Subject: [PATCH 13/19] [author-bot] updated AUTHORS file --- AUTHORS | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/AUTHORS b/AUTHORS index d0b207025..a647036c4 100644 --- a/AUTHORS +++ b/AUTHORS @@ -27,10 +27,10 @@ Mats Esseldeurs Stephane Michoulier Simone Ceppi MatsEsseldeurs +fhu +Spencer Magnall Caitlyn Hardiman Enrico Ragusa -Spencer Magnall -fhu Sergei Biriukov Cristiano Longarini Giovanni Dipierro @@ -38,12 +38,12 @@ Roberto Iaconi Hauke Worpel Alison Young Simone Ceppi -Amena Faruqi Stephen Neilson <36410751+s-neilson@users.noreply.github.com> +Amena Faruqi Martina Toscani Benedetta Veronesi -Sahl Rowther Simon Glover +Sahl Rowther Thomas Reichardt Jean-François Gonzalez Christopher Russell @@ -51,28 +51,28 @@ Alessia Franchini Alex Pettitt Jolien Malfait Phantom benchmark bot -Kieran Hirsh Nicole Rodrigues +Kieran Hirsh Amena Faruqi David Trevascus Farzana Meru -Chris Nixon -Megha Sharma Nicolas Cuello -Benoit Commercon -Giulia Ballabio -Joe Fisher -Maxime Lombart +Megha Sharma +Chris Nixon Megha Sharma +s-neilson <36410751+s-neilson@users.noreply.github.com> Orsola De Marco Terrence Tricco +Miguel Gonzalez-Bolivar +Benoit Commercon Zachary Pellow -s-neilson <36410751+s-neilson@users.noreply.github.com> -Alison Young -Cox, Samuel +Maxime Lombart +Joe Fisher +Giulia Ballabio Jorge Cuadra -Miguel Gonzalez-Bolivar -Nicolás Cuello -Steven Rieder Stéven Toupin +Nicolás Cuello mats esseldeurs +Alison Young +Cox, Samuel +Steven Rieder From b6f263a636a4b0e9ed4d95dfafebe58592448470 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 11:03:26 +1100 Subject: [PATCH 14/19] [indent-bot] standardised indentation --- src/main/extern_gr.F90 | 4 +- src/main/radiation_implicit.f90 | 8 +- src/main/tmunu2grid.f90 | 44 +-- src/setup/setup_flrw.f90 | 4 +- src/setup/setup_flrwpspec.f90 | 22 +- src/utils/analysis_common_envelope.f90 | 10 +- src/utils/analysis_radiotde.f90 | 18 +- src/utils/einsteintk_wrapper.f90 | 16 +- src/utils/interpolate3D.F90 | 496 ++++++++++++------------- src/utils/moddump_radiotde.f90 | 4 +- 10 files changed, 313 insertions(+), 313 deletions(-) diff --git a/src/main/extern_gr.F90 b/src/main/extern_gr.F90 index b118c17b6..17d050f42 100644 --- a/src/main/extern_gr.F90 +++ b/src/main/extern_gr.F90 @@ -89,9 +89,9 @@ subroutine dt_grforce(xyzh,fext,dtf) f2i = fext(1)*fext(1) + fext(2)*fext(2) + fext(3)*fext(3) if (f2i > 0.) then - dtf1 = sqrt(xyzh(4)/sqrt(f2i)) ! This is not really accurate since fi is a component of dp/dt, not da/dt + dtf1 = sqrt(xyzh(4)/sqrt(f2i)) ! This is not really accurate since fi is a component of dp/dt, not da/dt else - dtf1 = huge(dtf1) + dtf1 = huge(dtf1) endif select case (imetric) diff --git a/src/main/radiation_implicit.f90 b/src/main/radiation_implicit.f90 index 450956569..719111842 100644 --- a/src/main/radiation_implicit.f90 +++ b/src/main/radiation_implicit.f90 @@ -597,10 +597,10 @@ subroutine compute_flux(ivar,ijvar,ncompact,npart,icompactmax,varij2,vari,EU0,va if (dustRT) then if (dust_temp(i) < Tdust_threshold) opacity = nucleation(idkappa,i) endif - ! if (opacity < 0.) then - ! ierr = max(ierr,ierr_negative_opacity) - ! call error(label,'Negative opacity',val=opacity) - ! endif + ! if (opacity < 0.) then + ! ierr = max(ierr,ierr_negative_opacity) + ! call error(label,'Negative opacity',val=opacity) + ! endif if (limit_radiation_flux) then radRi = get_rad_R(rhoi,EU01i,dedx,opacity) diff --git a/src/main/tmunu2grid.f90 b/src/main/tmunu2grid.f90 index 2777e9a7d..bc5269940 100644 --- a/src/main/tmunu2grid.f90 +++ b/src/main/tmunu2grid.f90 @@ -109,36 +109,36 @@ subroutine get_tmunugrid_all(npart,xyzh,vxyzu,tmunus,calc_cfac) ! Put tmunu into an array of form ! tmunu(npart,16) do k=1, 4 - do j=1,4 - do i=1,npart - ! Check that this is correct!!! - ! print*,"i j is: ", k, j - ! print*, "Index in array is: ", (k-1)*4 + j - ! print*,tmunus(k,j,1) - dat(i, (k-1)*4 + j) = tmunus(k,j,i) - enddo - enddo -enddo + do j=1,4 + do i=1,npart + ! Check that this is correct!!! + ! print*,"i j is: ", k, j + ! print*, "Index in array is: ", (k-1)*4 + j + ! print*,tmunus(k,j,1) + dat(i, (k-1)*4 + j) = tmunus(k,j,i) + enddo + enddo + enddo !stop -ilendat = 16 + ilendat = 16 -call interpolate3D_vecexact(xyzh,weights,dat,ilendat,itype,npart,& + call interpolate3D_vecexact(xyzh,weights,dat,ilendat,itype,npart,& xmininterp(1),xmininterp(2),xmininterp(3), & datsmooth(:,ilower:iupper,jlower:jupper,klower:kupper),& ngrid(1),ngrid(2),ngrid(3),dxgrid(1),dxgrid(2),dxgrid(3),& normalise,periodicx,periodicy,periodicz) ! Put the smoothed array into tmunugrid -do i=1,4 - do j=1,4 - ! Check this is correct too! - !print*,"i j is: ", i, j - !print*, "Index in array is: ", (i-1)*4 + j - tmunugrid(i-1,j-1,:,:,:) = datsmooth((i-1)*4 + j, :,:,:) - !print*, "tmunugrid: ", tmunugrid(i-1,j-1,10,10,10) - !print*, datsmooth((i-1)*4 + j, 10,10,10) - enddo -enddo + do i=1,4 + do j=1,4 + ! Check this is correct too! + !print*,"i j is: ", i, j + !print*, "Index in array is: ", (i-1)*4 + j + tmunugrid(i-1,j-1,:,:,:) = datsmooth((i-1)*4 + j, :,:,:) + !print*, "tmunugrid: ", tmunugrid(i-1,j-1,10,10,10) + !print*, datsmooth((i-1)*4 + j, 10,10,10) + enddo + enddo !stop ! do k=1,4 ! do j=1,4 diff --git a/src/setup/setup_flrw.f90 b/src/setup/setup_flrw.f90 index e16173d2f..875c44de2 100644 --- a/src/setup/setup_flrw.f90 +++ b/src/setup/setup_flrw.f90 @@ -408,9 +408,9 @@ real function massfunc(x,xmin) end function massfunc real function deltaint(x) - real, intent(in) :: x + real, intent(in) :: x - deltaint = (1./kwave)*(kwave*kwave*c1 - 2)*ampl*cos(2*pi*x/lambda) + deltaint = (1./kwave)*(kwave*kwave*c1 - 2)*ampl*cos(2*pi*x/lambda) end function deltaint diff --git a/src/setup/setup_flrwpspec.f90 b/src/setup/setup_flrwpspec.f90 index f493a2766..2392255ac 100644 --- a/src/setup/setup_flrwpspec.f90 +++ b/src/setup/setup_flrwpspec.f90 @@ -241,9 +241,9 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! Check if files exist otherwise skip and return flat space if (.not. check_files(pspec_filename1,pspec_filename2,pspec_filename3)) then - print*, "Velocity files not found..." - print*, "Setting up flat space!" - return + print*, "Velocity files not found..." + print*, "Setting up flat space!" + return endif @@ -592,16 +592,16 @@ subroutine get_grid_neighbours(position,gridorigin,dx,xlower,ylower,zlower) end subroutine get_grid_neighbours logical function check_files(file1,file2,file3) - character(len=*), intent(in) :: file1,file2,file3 - logical :: file1_exist, file2_exist, file3_exist + character(len=*), intent(in) :: file1,file2,file3 + logical :: file1_exist, file2_exist, file3_exist - inquire(file=file1,exist=file1_exist) - inquire(file=file2,exist=file2_exist) - inquire(file=file3,exist=file3_exist) + inquire(file=file1,exist=file1_exist) + inquire(file=file2,exist=file2_exist) + inquire(file=file3,exist=file3_exist) - if ((.not. file1_exist) .or. (.not. file2_exist) .or. (.not. file3_exist)) then - check_files = .false. - endif + if ((.not. file1_exist) .or. (.not. file2_exist) .or. (.not. file3_exist)) then + check_files = .false. + endif end function check_files end module setup diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index c64b92bbe..02991276f 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -1163,15 +1163,15 @@ subroutine roche_lobe_values(time,npart,particlemass,xyzh,vxyzu) enddo if (nR1T == 0) then - MRL(iR1T) = 0 + MRL(iR1T) = 0 else - MRL(iR1T) = MRL(iR1T) / real(nR1T) + MRL(iR1T) = MRL(iR1T) / real(nR1T) endif if (nFB == 0) then - MRL(iFBV) = 0 + MRL(iFBV) = 0 else - MRL(iFBV) = MRL(iFBV) / real(nFB) + MRL(iFBV) = MRL(iFBV) / real(nFB) endif @@ -2549,7 +2549,7 @@ subroutine planet_profile(num,dumpfile,particlemass,xyzh,vxyzu) z(i) = dot_product(ri, vnorm) Rvec = ri - z(i)*vnorm R(i) = sqrt(dot_product(Rvec,Rvec)) - ! write(iu,"(es13.6,2x,es13.6,2x,es13.6)") R(i),z(i),rho(i) + ! write(iu,"(es13.6,2x,es13.6,2x,es13.6)") R(i),z(i),rho(i) write(iu,"(es13.6,2x,es13.6,2x,es13.6,2x,es13.6,2x,es13.6)") xyzh(1,i),xyzh(2,i),xyzh(3,i),rho(i),vxyzu(4,i) enddo diff --git a/src/utils/analysis_radiotde.f90 b/src/utils/analysis_radiotde.f90 index d4e99725c..e18da12f1 100644 --- a/src/utils/analysis_radiotde.f90 +++ b/src/utils/analysis_radiotde.f90 @@ -107,9 +107,9 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) call tde_analysis(npart,pmass,xyzh,vxyzu) if (n_cap > 0) then - open(iunit,file=output) - write(iunit,'("# ",es20.12," # TIME")') time - write(iunit,"('#',6(1x,'[',i2.2,1x,a11,']',2x))") & + open(iunit,file=output) + write(iunit,'("# ",es20.12," # TIME")') time + write(iunit,"('#',6(1x,'[',i2.2,1x,a11,']',2x))") & 1,'theta', & 2,'thetap', & 3,'phi', & @@ -117,18 +117,18 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) 5,'vtheta', & 6,'vphi' - do i = 1,npart - if (cap(i)) then - write(iunit,'(6(es18.10,1X))') & + do i = 1,npart + if (cap(i)) then + write(iunit,'(6(es18.10,1X))') & theta(i), & plot_theta(i), & phi(i), & vr(i), & vtheta(i), & vphi(i) - endif - enddo - close(iunit) + endif + enddo + close(iunit) endif deallocate(theta,plot_theta,phi,vr,vtheta,vphi,cap) diff --git a/src/utils/einsteintk_wrapper.f90 b/src/utils/einsteintk_wrapper.f90 index 4414c3142..ede060fcf 100644 --- a/src/utils/einsteintk_wrapper.f90 +++ b/src/utils/einsteintk_wrapper.f90 @@ -431,24 +431,24 @@ subroutine et2phantom_dumphydro(time,dt_et,checkpointfile) logical :: createcheckpoint if (present(checkpointfile)) then - createcheckpoint = .true. + createcheckpoint = .true. else - createcheckpoint = .false. + createcheckpoint = .false. endif ! Write EV_file if (.not. createcheckpoint) then - call write_evfile(time,dt_et) + call write_evfile(time,dt_et) - evfilestor = getnextfilename(evfilestor) - logfilestor = getnextfilename(logfilestor) - dumpfilestor = getnextfilename(dumpfilestor) - call write_fulldump(time,dumpfilestor) + evfilestor = getnextfilename(evfilestor) + logfilestor = getnextfilename(logfilestor) + dumpfilestor = getnextfilename(dumpfilestor) + call write_fulldump(time,dumpfilestor) endif ! Write full dump if (createcheckpoint) then - call write_fulldump(time,checkpointfile) + call write_fulldump(time,checkpointfile) endif ! Quick and dirty write cfac to txtfile diff --git a/src/utils/interpolate3D.F90 b/src/utils/interpolate3D.F90 index 80789b06f..ba9eac4c7 100644 --- a/src/utils/interpolate3D.F90 +++ b/src/utils/interpolate3D.F90 @@ -57,102 +57,102 @@ subroutine interpolate3D(xyzh,weight,dat,itype,npart,& xmin,ymin,zmin,datsmooth,npixx,npixy,npixz,pixwidthx,pixwidthy,pixwidthz,& normalise,periodicx,periodicy,periodicz) -integer, intent(in) :: npart,npixx,npixy,npixz -real, intent(in) :: xyzh(4,npart) -real, intent(in), dimension(npart) :: weight,dat -integer, intent(in), dimension(npart) :: itype -real, intent(in) :: xmin,ymin,zmin,pixwidthx,pixwidthy,pixwidthz -real, intent(out), dimension(npixx,npixy,npixz) :: datsmooth -logical, intent(in) :: normalise,periodicx,periodicy,periodicz + integer, intent(in) :: npart,npixx,npixy,npixz + real, intent(in) :: xyzh(4,npart) + real, intent(in), dimension(npart) :: weight,dat + integer, intent(in), dimension(npart) :: itype + real, intent(in) :: xmin,ymin,zmin,pixwidthx,pixwidthy,pixwidthz + real, intent(out), dimension(npixx,npixy,npixz) :: datsmooth + logical, intent(in) :: normalise,periodicx,periodicy,periodicz !logical, intent(in), exact_rendering -real, allocatable :: datnorm(:,:,:) - -integer :: i,ipix,jpix,kpix -integer :: iprintinterval,iprintnext -integer :: ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax -integer :: ipixi,jpixi,kpixi,nxpix,nwarn,threadid -real :: xminpix,yminpix,zminpix,hmin !,dhmin3 -real, dimension(npixx) :: dx2i -real :: xi,yi,zi,hi,hi1,hi21,wab,q2,const,dyz2,dz2 -real :: term,termnorm,dy,dz,ypix,zpix,xpixi,pixwidthmax,dfac -real :: t_start,t_end,t_used -logical :: iprintprogress -real, dimension(npart) :: x,y,z,hh -real :: radkernel, radkernel2, radkernh + real, allocatable :: datnorm(:,:,:) + + integer :: i,ipix,jpix,kpix + integer :: iprintinterval,iprintnext + integer :: ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax + integer :: ipixi,jpixi,kpixi,nxpix,nwarn,threadid + real :: xminpix,yminpix,zminpix,hmin !,dhmin3 + real, dimension(npixx) :: dx2i + real :: xi,yi,zi,hi,hi1,hi21,wab,q2,const,dyz2,dz2 + real :: term,termnorm,dy,dz,ypix,zpix,xpixi,pixwidthmax,dfac + real :: t_start,t_end,t_used + logical :: iprintprogress + real, dimension(npart) :: x,y,z,hh + real :: radkernel, radkernel2, radkernh ! Exact rendering -real :: pixint, wint + real :: pixint, wint !logical, parameter :: exact_rendering = .true. ! use exact rendering y/n -integer :: usedpart, negflag + integer :: usedpart, negflag !$ integer :: omp_get_num_threads,omp_get_thread_num -integer(kind=selected_int_kind(10)) :: iprogress,j ! up to 10 digits + integer(kind=selected_int_kind(10)) :: iprogress,j ! up to 10 digits ! Fill the particle data with xyzh -x(:) = xyzh(1,:) -y(:) = xyzh(2,:) -z(:) = xyzh(3,:) -hh(:) = xyzh(4,:) -cnormk3D = cnormk -radkernel = radkern -radkernel2 = radkern2 - -if (exact_rendering) then -print "(1x,a)",'interpolating to 3D grid (exact/Petkova+2018 on subgrid) ...' -elseif (normalise) then -print "(1x,a)",'interpolating to 3D grid (normalised) ...' -else -print "(1x,a)",'interpolating to 3D grid (non-normalised) ...' -endif -if (pixwidthx <= 0. .or. pixwidthy <= 0 .or. pixwidthz <= 0) then -print "(1x,a)",'interpolate3D: error: pixel width <= 0' -return -endif -if (any(hh(1:npart) <= tiny(hh))) then -print*,'interpolate3D: WARNING: ignoring some or all particles with h < 0' -endif + x(:) = xyzh(1,:) + y(:) = xyzh(2,:) + z(:) = xyzh(3,:) + hh(:) = xyzh(4,:) + cnormk3D = cnormk + radkernel = radkern + radkernel2 = radkern2 + + if (exact_rendering) then + print "(1x,a)",'interpolating to 3D grid (exact/Petkova+2018 on subgrid) ...' + elseif (normalise) then + print "(1x,a)",'interpolating to 3D grid (normalised) ...' + else + print "(1x,a)",'interpolating to 3D grid (non-normalised) ...' + endif + if (pixwidthx <= 0. .or. pixwidthy <= 0 .or. pixwidthz <= 0) then + print "(1x,a)",'interpolate3D: error: pixel width <= 0' + return + endif + if (any(hh(1:npart) <= tiny(hh))) then + print*,'interpolate3D: WARNING: ignoring some or all particles with h < 0' + endif !call wall_time(t_start) -datsmooth = 0. -if (normalise) then -allocate(datnorm(npixx,npixy,npixz)) -datnorm = 0. -endif + datsmooth = 0. + if (normalise) then + allocate(datnorm(npixx,npixy,npixz)) + datnorm = 0. + endif ! !--print a progress report if it is going to take a long time ! (a "long time" is, however, somewhat system dependent) ! -iprintprogress = (npart >= 100000) .or. (npixx*npixy > 100000) !.or. exact_rendering + iprintprogress = (npart >= 100000) .or. (npixx*npixy > 100000) !.or. exact_rendering ! !--loop over particles ! -iprintinterval = 25 -if (npart >= 1e6) iprintinterval = 10 -iprintnext = iprintinterval + iprintinterval = 25 + if (npart >= 1e6) iprintinterval = 10 + iprintnext = iprintinterval ! !--get starting CPU time ! -call cpu_time(t_start) + call cpu_time(t_start) -usedpart = 0 + usedpart = 0 -xminpix = xmin !- 0.5*pixwidthx -yminpix = ymin !- 0.5*pixwidthy -zminpix = zmin !- 0.5*pixwidthz -pixwidthmax = max(pixwidthx,pixwidthy,pixwidthz) + xminpix = xmin !- 0.5*pixwidthx + yminpix = ymin !- 0.5*pixwidthy + zminpix = zmin !- 0.5*pixwidthz + pixwidthmax = max(pixwidthx,pixwidthy,pixwidthz) ! !--use a minimum smoothing length on the grid to make ! sure that particles contribute to at least one pixel ! -hmin = 0.5*pixwidthmax + hmin = 0.5*pixwidthmax !dhmin3 = 1./(hmin*hmin*hmin) -const = cnormk3D ! normalisation constant (3D) + const = cnormk3D ! normalisation constant (3D) !print*, "const: ", const -nwarn = 0 -j = 0_8 -threadid = 1 + nwarn = 0 + j = 0_8 + threadid = 1 ! !--loop over particles ! @@ -177,216 +177,216 @@ subroutine interpolate3D(xyzh,weight,dat,itype,npart,& !$omp end master !$omp do schedule (guided, 2) -over_parts: do i=1,npart + over_parts: do i=1,npart ! !--report on progress ! -if (iprintprogress) then - !$omp atomic - j=j+1_8 + if (iprintprogress) then + !$omp atomic + j=j+1_8 !$ threadid = omp_get_thread_num() - iprogress = 100*j/npart - if (iprogress >= iprintnext .and. threadid==1) then - write(*,"(i3,'%.')",advance='no') iprogress - iprintnext = iprintnext + iprintinterval - endif -endif + iprogress = 100*j/npart + if (iprogress >= iprintnext .and. threadid==1) then + write(*,"(i3,'%.')",advance='no') iprogress + iprintnext = iprintnext + iprintinterval + endif + endif ! !--skip particles with itype < 0 ! -if (itype(i) < 0 .or. weight(i) < tiny(0.)) cycle over_parts - -hi = hh(i) -if (hi <= 0.) then - cycle over_parts -elseif (hi < hmin) then - ! - !--use minimum h to capture subgrid particles - ! (get better results *without* adjusting weights) - ! - termnorm = const*weight(i) !*(hi*hi*hi)*dhmin3 - if (.not.exact_rendering) hi = hmin -else - termnorm = const*weight(i) -endif + if (itype(i) < 0 .or. weight(i) < tiny(0.)) cycle over_parts + + hi = hh(i) + if (hi <= 0.) then + cycle over_parts + elseif (hi < hmin) then + ! + !--use minimum h to capture subgrid particles + ! (get better results *without* adjusting weights) + ! + termnorm = const*weight(i) !*(hi*hi*hi)*dhmin3 + if (.not.exact_rendering) hi = hmin + else + termnorm = const*weight(i) + endif ! !--set kernel related quantities ! -xi = x(i) -yi = y(i) -zi = z(i) + xi = x(i) + yi = y(i) + zi = z(i) -hi1 = 1./hi -hi21 = hi1*hi1 -radkernh = radkernel*hi ! radius of the smoothing kernel + hi1 = 1./hi + hi21 = hi1*hi1 + radkernh = radkernel*hi ! radius of the smoothing kernel !termnorm = const*weight(i) -term = termnorm*dat(i) -dfac = hi**3/(pixwidthx*pixwidthy*pixwidthz*const) + term = termnorm*dat(i) + dfac = hi**3/(pixwidthx*pixwidthy*pixwidthz*const) !dfac = hi**3/(pixwidthx*pixwidthy*const) ! !--for each particle work out which pixels it contributes to ! -ipixmin = int((xi - radkernh - xmin)/pixwidthx) -jpixmin = int((yi - radkernh - ymin)/pixwidthy) -kpixmin = int((zi - radkernh - zmin)/pixwidthz) -ipixmax = int((xi + radkernh - xmin)/pixwidthx) + 1 -jpixmax = int((yi + radkernh - ymin)/pixwidthy) + 1 -kpixmax = int((zi + radkernh - zmin)/pixwidthz) + 1 - -if (.not.periodicx) then - if (ipixmin < 1) ipixmin = 1 ! make sure they only contribute - if (ipixmax > npixx) ipixmax = npixx ! to pixels in the image -endif -if (.not.periodicy) then - if (jpixmin < 1) jpixmin = 1 - if (jpixmax > npixy) jpixmax = npixy -endif -if (.not.periodicz) then - if (kpixmin < 1) kpixmin = 1 - if (kpixmax > npixz) kpixmax = npixz -endif - -negflag = 0 + ipixmin = int((xi - radkernh - xmin)/pixwidthx) + jpixmin = int((yi - radkernh - ymin)/pixwidthy) + kpixmin = int((zi - radkernh - zmin)/pixwidthz) + ipixmax = int((xi + radkernh - xmin)/pixwidthx) + 1 + jpixmax = int((yi + radkernh - ymin)/pixwidthy) + 1 + kpixmax = int((zi + radkernh - zmin)/pixwidthz) + 1 + + if (.not.periodicx) then + if (ipixmin < 1) ipixmin = 1 ! make sure they only contribute + if (ipixmax > npixx) ipixmax = npixx ! to pixels in the image + endif + if (.not.periodicy) then + if (jpixmin < 1) jpixmin = 1 + if (jpixmax > npixy) jpixmax = npixy + endif + if (.not.periodicz) then + if (kpixmin < 1) kpixmin = 1 + if (kpixmax > npixz) kpixmax = npixz + endif + + negflag = 0 ! !--precalculate an array of dx2 for this particle (optimisation) ! -nxpix = 0 -do ipix=ipixmin,ipixmax - nxpix = nxpix + 1 - ipixi = ipix - if (periodicx) ipixi = iroll(ipix,npixx) - xpixi = xminpix + ipix*pixwidthx - !--watch out for errors with periodic wrapping... - if (nxpix <= size(dx2i)) then - dx2i(nxpix) = ((xpixi - xi)**2)*hi21 - endif -enddo + nxpix = 0 + do ipix=ipixmin,ipixmax + nxpix = nxpix + 1 + ipixi = ipix + if (periodicx) ipixi = iroll(ipix,npixx) + xpixi = xminpix + ipix*pixwidthx + !--watch out for errors with periodic wrapping... + if (nxpix <= size(dx2i)) then + dx2i(nxpix) = ((xpixi - xi)**2)*hi21 + endif + enddo !--if particle contributes to more than npixx pixels ! (i.e. periodic boundaries wrap more than once) ! truncate the contribution and give warning -if (nxpix > npixx) then - nwarn = nwarn + 1 - ipixmax = ipixmin + npixx - 1 -endif + if (nxpix > npixx) then + nwarn = nwarn + 1 + ipixmax = ipixmin + npixx - 1 + endif ! !--loop over pixels, adding the contribution from this particle ! -do kpix = kpixmin,kpixmax - kpixi = kpix - if (periodicz) kpixi = iroll(kpix,npixz) - - zpix = zminpix + kpix*pixwidthz - dz = zpix - zi - dz2 = dz*dz*hi21 - - do jpix = jpixmin,jpixmax - jpixi = jpix - if (periodicy) jpixi = iroll(jpix,npixy) - - ypix = yminpix + jpix*pixwidthy - dy = ypix - yi - dyz2 = dy*dy*hi21 + dz2 - - nxpix = 0 - do ipix = ipixmin,ipixmax - if ((kpix==kpixmin).and.(jpix==jpixmin).and.(ipix==ipixmin)) then - usedpart = usedpart + 1 - endif - - nxpix = nxpix + 1 - ipixi = ipix - if (periodicx) ipixi = iroll(ipix,npixx) - - q2 = dx2i(nxpix) + dyz2 ! dx2 pre-calculated; dy2 pre-multiplied by hi21 - - if (exact_rendering .and. ipixmax-ipixmin <= 4) then - if (q2 < radkernel2 + 3.*pixwidthmax**2*hi21) then - xpixi = xminpix + ipix*pixwidthx - - ! Contribution of the cell walls in the xy-plane - pixint = 0.0 - wint = wallint(zpix-zi+0.5*pixwidthz,xi,yi,xpixi,ypix,pixwidthx,pixwidthy,hi) - pixint = pixint + wint - - wint = wallint(zi-zpix+0.5*pixwidthz,xi,yi,xpixi,ypix,pixwidthx,pixwidthy,hi) - pixint = pixint + wint - - ! Contribution of the cell walls in the xz-plane - wint = wallint(ypix-yi+0.5*pixwidthy,xi,zi,xpixi,zpix,pixwidthx,pixwidthz,hi) - pixint = pixint + wint - - wint = wallint(yi-ypix+0.5*pixwidthy,xi,zi,xpixi,zpix,pixwidthx,pixwidthz,hi) - pixint = pixint + wint - - ! Contribution of the cell walls in the yz-plane - wint = wallint(xpixi-xi+0.5*pixwidthx,zi,yi,zpix,ypix,pixwidthz,pixwidthy,hi) - pixint = pixint + wint - - wint = wallint(xi-xpixi+0.5*pixwidthx,zi,yi,zpix,ypix,pixwidthz,pixwidthy,hi) - pixint = pixint + wint - - wab = pixint*dfac ! /(pixwidthx*pixwidthy*pixwidthz*const)*hi**3 - - if (pixint < -0.01d0) then - print*, "Error: (",ipixi,jpixi,kpixi,") -> ", pixint, term*wab - endif - - ! - !--calculate data value at this pixel using the summation interpolant - ! - !$omp atomic - datsmooth(ipixi,jpixi,kpixi) = datsmooth(ipixi,jpixi,kpixi) + term*wab - if (normalise) then - !$omp atomic - datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab - endif - endif - else - if (q2 < radkernel2) then - - ! - !--SPH kernel - standard cubic spline - ! - wab = wkernel(q2) - ! - !--calculate data value at this pixel using the summation interpolant - ! - !$omp atomic - datsmooth(ipixi,jpixi,kpixi) = datsmooth(ipixi,jpixi,kpixi) + term*wab - if (normalise) then - !$omp atomic - datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab - endif - endif - endif - enddo - enddo -enddo -enddo over_parts + do kpix = kpixmin,kpixmax + kpixi = kpix + if (periodicz) kpixi = iroll(kpix,npixz) + + zpix = zminpix + kpix*pixwidthz + dz = zpix - zi + dz2 = dz*dz*hi21 + + do jpix = jpixmin,jpixmax + jpixi = jpix + if (periodicy) jpixi = iroll(jpix,npixy) + + ypix = yminpix + jpix*pixwidthy + dy = ypix - yi + dyz2 = dy*dy*hi21 + dz2 + + nxpix = 0 + do ipix = ipixmin,ipixmax + if ((kpix==kpixmin).and.(jpix==jpixmin).and.(ipix==ipixmin)) then + usedpart = usedpart + 1 + endif + + nxpix = nxpix + 1 + ipixi = ipix + if (periodicx) ipixi = iroll(ipix,npixx) + + q2 = dx2i(nxpix) + dyz2 ! dx2 pre-calculated; dy2 pre-multiplied by hi21 + + if (exact_rendering .and. ipixmax-ipixmin <= 4) then + if (q2 < radkernel2 + 3.*pixwidthmax**2*hi21) then + xpixi = xminpix + ipix*pixwidthx + + ! Contribution of the cell walls in the xy-plane + pixint = 0.0 + wint = wallint(zpix-zi+0.5*pixwidthz,xi,yi,xpixi,ypix,pixwidthx,pixwidthy,hi) + pixint = pixint + wint + + wint = wallint(zi-zpix+0.5*pixwidthz,xi,yi,xpixi,ypix,pixwidthx,pixwidthy,hi) + pixint = pixint + wint + + ! Contribution of the cell walls in the xz-plane + wint = wallint(ypix-yi+0.5*pixwidthy,xi,zi,xpixi,zpix,pixwidthx,pixwidthz,hi) + pixint = pixint + wint + + wint = wallint(yi-ypix+0.5*pixwidthy,xi,zi,xpixi,zpix,pixwidthx,pixwidthz,hi) + pixint = pixint + wint + + ! Contribution of the cell walls in the yz-plane + wint = wallint(xpixi-xi+0.5*pixwidthx,zi,yi,zpix,ypix,pixwidthz,pixwidthy,hi) + pixint = pixint + wint + + wint = wallint(xi-xpixi+0.5*pixwidthx,zi,yi,zpix,ypix,pixwidthz,pixwidthy,hi) + pixint = pixint + wint + + wab = pixint*dfac ! /(pixwidthx*pixwidthy*pixwidthz*const)*hi**3 + + if (pixint < -0.01d0) then + print*, "Error: (",ipixi,jpixi,kpixi,") -> ", pixint, term*wab + endif + + ! + !--calculate data value at this pixel using the summation interpolant + ! + !$omp atomic + datsmooth(ipixi,jpixi,kpixi) = datsmooth(ipixi,jpixi,kpixi) + term*wab + if (normalise) then + !$omp atomic + datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab + endif + endif + else + if (q2 < radkernel2) then + + ! + !--SPH kernel - standard cubic spline + ! + wab = wkernel(q2) + ! + !--calculate data value at this pixel using the summation interpolant + ! + !$omp atomic + datsmooth(ipixi,jpixi,kpixi) = datsmooth(ipixi,jpixi,kpixi) + term*wab + if (normalise) then + !$omp atomic + datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab + endif + endif + endif + enddo + enddo + enddo + enddo over_parts !$omp enddo !$omp end parallel -if (nwarn > 0) then -print "(a,i11,a,/,a)",' interpolate3D: WARNING: contributions truncated from ',nwarn,' particles',& + if (nwarn > 0) then + print "(a,i11,a,/,a)",' interpolate3D: WARNING: contributions truncated from ',nwarn,' particles',& ' that wrap periodic boundaries more than once' -endif + endif ! !--normalise dat array ! -if (normalise) then -where (datnorm > tiny(datnorm)) - datsmooth = datsmooth/datnorm -end where -endif -if (allocated(datnorm)) deallocate(datnorm) + if (normalise) then + where (datnorm > tiny(datnorm)) + datsmooth = datsmooth/datnorm + end where + endif + if (allocated(datnorm)) deallocate(datnorm) !call wall_time(t_end) -call cpu_time(t_end) -t_used = t_end - t_start -print*, 'Interpolate3D completed in ',t_end-t_start,'s' + call cpu_time(t_end) + t_used = t_end - t_start + print*, 'Interpolate3D completed in ',t_end-t_start,'s' end subroutine interpolate3D @@ -695,8 +695,8 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& !--calculate data value at this pixel using the summation interpolant ! do smoothindex=1, ilendat - !$omp atomic - datsmooth(smoothindex,ipixi,jpixi,kpixi) = datsmooth(smoothindex,ipixi,jpixi,kpixi) + term(smoothindex)*wab + !$omp atomic + datsmooth(smoothindex,ipixi,jpixi,kpixi) = datsmooth(smoothindex,ipixi,jpixi,kpixi) + term(smoothindex)*wab enddo if (normalise) then !$omp atomic @@ -714,8 +714,8 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& !--calculate data value at this pixel using the summation interpolant ! do smoothindex=1,ilendat - !$omp atomic - datsmooth(smoothindex,ipixi,jpixi,kpixi) = datsmooth(smoothindex,ipixi,jpixi,kpixi) + term(smoothindex)*wab + !$omp atomic + datsmooth(smoothindex,ipixi,jpixi,kpixi) = datsmooth(smoothindex,ipixi,jpixi,kpixi) + term(smoothindex)*wab enddo if (normalise) then !$omp atomic @@ -739,12 +739,12 @@ subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& !--normalise dat array ! if (normalise) then - do i=1, ilendat - where (datnorm > tiny(datnorm)) + do i=1, ilendat + where (datnorm > tiny(datnorm)) - datsmooth(i,:,:,:) = datsmooth(i,:,:,:)/datnorm(:,:,:) - end where - enddo + datsmooth(i,:,:,:) = datsmooth(i,:,:,:)/datnorm(:,:,:) + end where + enddo endif if (allocated(datnorm)) deallocate(datnorm) diff --git a/src/utils/moddump_radiotde.f90 b/src/utils/moddump_radiotde.f90 index 5e12a738e..774536628 100644 --- a/src/utils/moddump_radiotde.f90 +++ b/src/utils/moddump_radiotde.f90 @@ -94,7 +94,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) use_func = .true. use_func_old = use_func remove_overlap = .true. - !--Power law default setups + !--Power law default setups rad_max = 7.1e16 ! in cm rad_min = 8.7e15 ! in cm nbreak = 1 @@ -103,7 +103,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) allocate(rhof_n(nbreak),rhof_rbreak(nbreak)) rhof_n = -1.7 rhof_rbreak = rad_min - !--Profile default setups + !--Profile default setups read_temp = .false. profile_filename = default_name nprof = 7 From 482376134f9f20ac6768fe164585fdac76fa0265 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 11:04:51 +1100 Subject: [PATCH 15/19] (interpolate) remove obsolete routine --- src/utils/interpolate3Dold.F90 | 367 --------------------------------- 1 file changed, 367 deletions(-) delete mode 100644 src/utils/interpolate3Dold.F90 diff --git a/src/utils/interpolate3Dold.F90 b/src/utils/interpolate3Dold.F90 deleted file mode 100644 index c7fff7ca7..000000000 --- a/src/utils/interpolate3Dold.F90 +++ /dev/null @@ -1,367 +0,0 @@ -!--------------------------------------------------------------------------! -! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! -! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! -! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.github.io/ ! -!--------------------------------------------------------------------------! -module interpolations3D -! -! Module containing routine for interpolation from PHANTOM data -! to 3D adaptive mesh -! -! Requires adaptivemesh.f90 module -! -! :References: None -! -! :Owner: Spencer Magnall -! -! :Runtime parameters: None -! -! :Dependencies: kernel -! - - implicit none - real, parameter, private :: dpi = 1./3.1415926536d0 - public :: interpolate3D -!$ integer(kind=8), dimension(:), private, allocatable :: ilock - -contains -!-------------------------------------------------------------------------- -! subroutine to interpolate from particle data to even grid of pixels -! -! The data is interpolated according to the formula -! -! datsmooth(pixel) = sum_b weight_b dat_b W(r-r_b, h_b) -! -! where _b is the quantity at the neighbouring particle b and -! W is the smoothing kernel, for which we use the usual cubic spline. -! -! For a standard SPH smoothing the weight function for each particle should be -! -! weight = pmass/(rho*h^3) -! -! this version is written for slices through a rectangular volume, ie. -! assumes a uniform pixel size in x,y, whilst the number of pixels -! in the z direction can be set to the number of cross-section slices. -! -! Input: particle coordinates and h : xyzh(4,npart) -! weight for each particle : weight [ same on all parts in PHANTOM ] -! scalar data to smooth : dat (npart) -! -! Output: smoothed data : datsmooth (npixx,npixy,npixz) -! -! Daniel Price, Monash University 2010 -! daniel.price@monash.edu -!-------------------------------------------------------------------------- - -subroutine interpolate3D(xyzh,weight,npart, & - xmin,datsmooth,nnodes,dxgrid,normalise,dat,ngrid,vertexcen) - use kernel, only:wkern, radkern, radkern2, cnormk - !use adaptivemesh, only:ifirstlevel,nsub,ndim,gridnodes - integer, intent(in) :: npart,nnodes,ngrid(3) - real, intent(in) :: xyzh(:,:)! ,vxyzu(:,:) - real, intent(in) :: weight !,pmass - real, intent(in) :: xmin(3),dxgrid(3) - real, intent(out) :: datsmooth(:,:,:) - logical, intent(in) :: normalise, vertexcen - real, intent(in), optional :: dat(:) - real, allocatable :: datnorm(:,:,:) -! real, dimension(nsub**ndim,nnodes) :: datnorm - integer, parameter :: ndim = 3, nsub=1 - integer :: i,ipix,jpix,kpix,isubmesh,imesh,level,icell - integer :: iprintinterval,iprintnext - integer :: ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax - integer :: ipixi,jpixi,kpixi,npixx,npixy,npixz - real :: xi,yi,zi,hi,hi1,hi21,radkernh,qq,wab,q2,const,dyz2,dz2 - real :: xorigi,yorigi,zorigi,xpix,ypix,zpix,dx,dy,dz - real :: dxcell(ndim),xminnew(ndim), dxmax(ndim) - real :: t_start,t_end - real :: termnorm - real :: term - logical :: iprintprogress -!$ integer :: omp_get_num_threads,j -#ifndef _OPENMP - integer(kind=8) :: iprogress -#endif - - print*, "size: ", size(datsmooth) - print*, "datsmooth out of bounds: ", datsmooth(35,1,1) - datsmooth = 0. - dxmax(:) = dxgrid(:) - !datnorm = 0. - if (normalise) then - print "(1x,a)",'interpolating from particles to Einstein toolkit grid (normalised) ...' - else - print "(1x,a)",'interpolating from particles to Einstein toolkit grid (non-normalised) ...' - endif -! if (any(dxmax(:) <= 0.)) then -! print "(1x,a)",'interpolate3D: error: grid size <= 0' -! return -! endif -! if (ilendat /= 0) then -! print "(1x,a)",'interpolate3D: error in interface: dat has non-zero length but is not present' -! return -! endif - if (normalise) then - allocate(datnorm(ngrid(1),ngrid(2),ngrid(3))) - datnorm = 0. - endif - -!$ allocate(ilock(0:nnodes)) -!$ do i=0,nnodes -!$ call omp_init_lock(ilock(i)) -!$ enddo - - ! - !--print a progress report if it is going to take a long time - ! (a "long time" is, however, somewhat system dependent) - ! - iprintprogress = (npart >= 100000) .or. (nnodes > 10000) - ! - !--loop over particles - ! - iprintinterval = 25 - if (npart >= 1e6) iprintinterval = 10 - iprintnext = iprintinterval - ! - !--get starting CPU time - ! - call cpu_time(t_start) - - imesh = 1 - level = 1 - dxcell(:) = dxgrid(:)/real(nsub**level) -! xminpix(:) = xmin(:) - 0.5*dxcell(:) - npixx = ngrid(1) - npixy = ngrid(2) - npixz = ngrid(3) - print "(3(a,i4))",' root grid: ',npixx,' x ',npixy,' x ',npixz - print*, "position of i cell is: ", 1*dxcell(1) + xmin(1) - print*, "npart: ", npart - - const = cnormk ! kernel normalisation constant (3D) - print*,"const: ", const - !stop - - ! - !--loop over particles - ! - !$omp parallel default(none) & - !$omp shared(npart,xyzh,dat,datsmooth,datnorm,vertexcen,const,weight) & - !$omp shared(xmin,imesh,nnodes,level) & - !$omp shared(npixx,npixy,npixz,dxmax,dxcell,normalise) & - !$omp private(i,j,hi,hi1,hi21,termnorm,term) & - !$omp private(xpix,ypix,zpix,dx,dy,dz,dz2,dyz2,qq,q2,wab,radkernh) & - !$omp private(xi,yi,zi,xorigi,yorigi,zorigi,xminnew) & - !$omp private(ipix,jpix,kpix,ipixi,jpixi,kpixi,icell,isubmesh) & - !$omp private(ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax) - !$omp master -!$ print "(1x,a,i3,a)",'Using ',omp_get_num_threads(),' cpus' - !$omp end master - !$omp do schedule(guided,10) - over_parts: do i=1,npart - ! - !--report on progress - ! - !print*, i -#ifndef _OPENMP - if (iprintprogress) then - iprogress = nint(100.*i/npart) - if (iprogress >= iprintnext) then - write(*,"('(',i3,'% -',i12,' particles done)')") iprogress,i - iprintnext = iprintnext + iprintinterval - endif - endif -#endif - ! - !--set kernel related quantities - ! - xi = xyzh(1,i); xorigi = xi - yi = xyzh(2,i); yorigi = yi - zi = xyzh(3,i); zorigi = zi - hi = xyzh(4,i) - radkernh = radkern*hi - !print*, "hi: ", hi - if (hi <= 0.) cycle over_parts - hi1 = 1./hi; hi21 = hi1*hi1 - termnorm = const*weight - ! print*, "const: ", const - ! print*, "weight: ", weight - ! print*, "termnorm: ", termnorm - - !radkern = 2.*hi ! radius of the smoothing kernel - !print*, "radkern: ", radkern - !print*, "part pos: ", xi,yi,zi - term = termnorm*dat(i) ! weight for density calculation - ! I don't understand why this doesnt involve any actual smoothing? - !dfac = hi**3/(dxcell(1)*dxcell(2)*dxcell(3)*const) - ! - !--for each particle work out which pixels it contributes to - ! - !print*, "radkern: ", radkern - ipixmin = int((xi - radkernh - xmin(1))/dxcell(1)) - jpixmin = int((yi - radkernh - xmin(2))/dxcell(2)) - kpixmin = int((zi - radkernh - xmin(3))/dxcell(3)) - - ipixmax = int((xi + radkernh - xmin(1))/dxcell(1)) + 1 - jpixmax = int((yi + radkernh - xmin(2))/dxcell(2)) + 1 - kpixmax = nint((zi + radkernh - xmin(3))/dxcell(3)) + 1 - - !if (ipixmax == 33) stop - - - !if (ipixmin == 4 .and. jpixmin == 30 .and. kpixmin == 33) print*, "particle (min): ", i - !if (ipixmax == 4 .and. jpixmax == 30 .and. kpixmax == 33) print*, "particle (max): ", i -#ifndef PERIODIC - if (ipixmin < 1) ipixmin = 1 ! make sure they only contribute - if (jpixmin < 1) jpixmin = 1 ! to pixels in the image - if (kpixmin < 1) kpixmin = 1 - if (ipixmax > npixx) ipixmax = npixx - if (jpixmax > npixy) jpixmax = npixy - if (kpixmax > npixz) kpixmax = npixz - !print*, "ipixmin: ", ipixmin - !print*, "ipixmax: ", ipixmax - !print*, "jpixmin: ", jpixmin - !print*, "jpixmax: ", jpixmax - !print*, "kpixmin: ", kpixmin - !print*, "kpixmax: ", kpixmax -#endif - !print*,' part ',i,' lims = ',ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax - ! - !--loop over pixels, adding the contribution from this particle - ! (note that we handle the periodic boundary conditions - ! entirely on the root grid) - ! - do kpix = kpixmin,kpixmax - kpixi = kpix -#ifdef PERIODIC - if (kpixi < 1) then - kpixi = kpixi + npixz - zi = zorigi !+ dxmax(3) - elseif (kpixi > npixz) then - kpixi = kpixi - npixz - zi = zorigi !- dxmax(3) - else - zi = zorigi - endif -#endif - if (vertexcen) then - zpix = xmin(3) + (kpixi-1)*dxcell(3) - else - zpix = xmin(3) + (kpixi-0.5)*dxcell(3) - endif - dz = zpix - zi - dz2 = dz*dz*hi21 - - do jpix = jpixmin,jpixmax - jpixi = jpix -#ifdef PERIODIC - if (jpixi < 1) then - jpixi = jpixi + npixy - yi = yorigi !+ dxmax(2) - elseif (jpixi > npixy) then - jpixi = jpixi - npixy - yi = yorigi !- dxmax(2) - else - yi = yorigi - endif -#endif - if (vertexcen) then - ypix = xmin(2) + (jpixi-1)*dxcell(2) - else - ypix = xmin(2) + (jpixi-0.5)*dxcell(2) - endif - dy = ypix - yi - dyz2 = dy*dy*hi21 + dz2 - - do ipix = ipixmin,ipixmax - ipixi = ipix -#ifdef PERIODIC - if (ipixi < 1) then - ipixi = ipixi + npixx - xi = xorigi !+ dxmax(1) - elseif (ipixi > npixx) then - if (ipixi == 33) then - print*,"xi old: ", xorigi - print*, "xi new: ", xorigi-dxmax(1) - print*, "ipixi new: ", ipixi - npixx - endif - ipixi = ipixi - npixx - xi = xorigi !- dxmax(1) - else - xi = xorigi - endif -#endif - icell = ((kpixi-1)*nsub + (jpixi-1))*nsub + ipixi - ! - !--particle interpolates directly onto the root grid - ! - !print*,'onto root grid ',ipixi,jpixi,kpixi - if (vertexcen) then - xpix = xmin(1) + (ipixi-1)*dxcell(1) - else - xpix = xmin(1) + (ipixi-0.5)*dxcell(1) - endif - !print*, "xpix: ", xpix - !xpix = xmin(1) + (ipixi-1)*dxcell(1) ! Since we are vertex centered from Et - dx = xpix - xi - q2 = dx*dx*hi21 + dyz2 ! dx2 pre-calculated; dy2 pre-multiplied by hi21 - ! - !--SPH kernel - standard cubic spline - ! - if (q2 < radkern2) then - ! if (q2 < 1.0) then - ! qq = sqrt(q2) - ! wab = 1.-1.5*q2 + 0.75*q2*qq - ! else - ! qq = sqrt(q2) - ! wab = 0.25*(2.-qq)**3 - ! endif - ! Call the kernel routine - qq = sqrt(q2) - wab = wkern(q2,qq) - ! - !--calculate data value at this pixel using the summation interpolant - ! - ! Change this to the access the pixel coords x,y,z - !$omp critical - datsmooth(ipixi,jpixi,kpixi) = datsmooth(ipixi,jpixi,kpixi) + term*wab - - !if (ipixi==1 .and. jpixi==1 .and. kpixi==1) print*, "x position of 1,1,1", xi,yi,zi - if (normalise) then - datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab - endif - !$omp end critical - endif - enddo - enddo - enddo - enddo over_parts - !$omp enddo - !$omp end parallel - -!$ do i=0,nnodes -!$ call omp_destroy_lock(ilock(i)) -!$ enddo -!$ if (allocated(ilock)) deallocate(ilock) - - ! - !--normalise dat array - ! - if (normalise) then - where (datnorm > tiny(datnorm)) - datsmooth = datsmooth/datnorm - end where - endif - if (allocated(datnorm)) deallocate(datnorm) - ! - !--get ending CPU time - ! - call cpu_time(t_end) - print*,'completed in ',t_end-t_start,'s' - - return - -end subroutine interpolate3D - -end module interpolations3D From 6fde5a96acd7c491738671b19d0efd97ec032cf2 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 11:59:37 +1100 Subject: [PATCH 16/19] (#484) further reduced ntasks to 2 for MPI workflow --- .github/workflows/mpi.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/mpi.yml b/.github/workflows/mpi.yml index acd28d60c..46099c9f2 100644 --- a/.github/workflows/mpi.yml +++ b/.github/workflows/mpi.yml @@ -19,7 +19,7 @@ jobs: - yes ntasks: - 1 - - 4 + - 2 input: # [SETUP, phantom_tests] - ['test', ''] - ['testkd', ''] From 3711dc7021f9172a228774e125852b4e93ee8cc1 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 12:11:23 +1100 Subject: [PATCH 17/19] updated mailmap --- .mailmap | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/.mailmap b/.mailmap index 567c60b95..b870153b5 100644 --- a/.mailmap +++ b/.mailmap @@ -21,12 +21,13 @@ Rebecca Nealon Nealon Alex Alex Pettitt Alex Pettitt -Terrence Tricco - - - - - + + + + + + +Terrence Tricco James Wurster James Wurster James Wurster jameswurster James Wurster jameswurster @@ -46,12 +47,12 @@ Stéven Toupin stoupin Guillaume Laibe glaibe Guillaume Laibe glaibe Alice Cerioli ALICE CERIOLI +Alice Cerioli Thomas Reichardt Thomas Reichardt Thomas Reichardt Thomas Reichardt Mr Thomas Reichardt Roberto Iaconi Roberto Iaconi Roberto Iaconi Roberto Iaconi -Alice Cerioli Daniel Mentiplay Daniel Mentiplay Daniel Mentiplay @@ -85,8 +86,13 @@ Fangyi (Fitz) Hu Fitz-Hu <54089891+Fitz-Hu@users.n Fangyi (Fitz) Hu root Fangyi (Fitz) Hu root Fangyi (Fitz) Hu fitzHu <54089891+Fitz-Hu@users.noreply.github.com> +Fangyi (Fitz) Hu Fitz Hu +Fangyi (Fitz) Hu fhu Megha Sharma Megha Sharma <40732335+msha0023@users.noreply.github.com> Megha Sharma megha sharma +Megha Sharma Megha Sharma +Megha Sharma Megha Sharma +Megha Sharma Megha Sharma Mike Lau Mike Lau <55525335+themikelau@users.noreply.github.com> Elisabeth Borchert emborchert <69176538+emborchert@users.noreply.github.com> Ward Homan ward @@ -103,3 +109,9 @@ Sahl Rowther Sahl Rowther sahl95 Caitlyn Hardiman caitlynhardiman <72479852+caitlynhardiman@users.noreply.github.com> Amena Faruqi <42060670+amenafaruqi@users.noreply.github.com> +Amena Faruqi Amena Faruqi +Alison Young Alison Young +Simone Ceppi Simone Ceppi +Mats Esseldeurs mats esseldeurs +Mats Esseldeurs MatsEsseldeurs +Nicolás Cuello Nicolas Cuello From 7f901e18124c20e8860d6ba4a3b2c6ce3c83323d Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 12:13:33 +1100 Subject: [PATCH 18/19] [header-bot] updated file headers --- src/utils/analysis_radiotde.f90 | 2 +- src/utils/moddump_radiotde.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utils/analysis_radiotde.f90 b/src/utils/analysis_radiotde.f90 index e18da12f1..b994822d8 100644 --- a/src/utils/analysis_radiotde.f90 +++ b/src/utils/analysis_radiotde.f90 @@ -10,7 +10,7 @@ module analysis ! ! :References: None ! -! :Owner: fhu +! :Owner: Fitz) Hu ! ! :Runtime parameters: ! - drad_cap : *capture thickness (in cm) (-ve for all particles at outer radius)* diff --git a/src/utils/moddump_radiotde.f90 b/src/utils/moddump_radiotde.f90 index 774536628..b64612288 100644 --- a/src/utils/moddump_radiotde.f90 +++ b/src/utils/moddump_radiotde.f90 @@ -10,7 +10,7 @@ module moddump ! ! :References: None ! -! :Owner: fhu +! :Owner: Fitz) Hu ! ! :Runtime parameters: ! - ieos : *equation of state used* From fee837a0c65b10899b036f69cb28dfd688784940 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 27 Nov 2023 12:13:46 +1100 Subject: [PATCH 19/19] [author-bot] updated AUTHORS file --- AUTHORS | 48 ++++++++++++++++++------------------------------ 1 file changed, 18 insertions(+), 30 deletions(-) diff --git a/AUTHORS b/AUTHORS index a647036c4..1dc027d7b 100644 --- a/AUTHORS +++ b/AUTHORS @@ -11,23 +11,19 @@ Conrad Chan James Wurster David Liptai Lionel Siess +Fangyi (Fitz) Hu Daniel Mentiplay +Megha Sharma Arnaud Vericel Mark Hutchison -Fitz Hu -Megha Sharma +Mats Esseldeurs Rebecca Nealon Elisabeth Borchert Ward Homan Christophe Pinte -Fangyi (Fitz) Hu -Megha Sharma -Terrence Tricco -Mats Esseldeurs +Terrence Tricco +Simone Ceppi Stephane Michoulier -Simone Ceppi -MatsEsseldeurs -fhu Spencer Magnall Caitlyn Hardiman Enrico Ragusa @@ -36,43 +32,35 @@ Cristiano Longarini Giovanni Dipierro Roberto Iaconi Hauke Worpel +Amena Faruqi Alison Young -Simone Ceppi Stephen Neilson <36410751+s-neilson@users.noreply.github.com> -Amena Faruqi Martina Toscani Benedetta Veronesi -Simon Glover Sahl Rowther Thomas Reichardt +Simon Glover Jean-François Gonzalez Christopher Russell +Phantom benchmark bot +Jolien Malfait Alessia Franchini Alex Pettitt -Jolien Malfait -Phantom benchmark bot Nicole Rodrigues Kieran Hirsh -Amena Faruqi -David Trevascus +Nicolás Cuello Farzana Meru -Nicolas Cuello -Megha Sharma +David Trevascus Chris Nixon -Megha Sharma -s-neilson <36410751+s-neilson@users.noreply.github.com> -Orsola De Marco -Terrence Tricco +Giulia Ballabio Miguel Gonzalez-Bolivar -Benoit Commercon -Zachary Pellow Maxime Lombart +Benoit Commercon +Orsola De Marco Joe Fisher -Giulia Ballabio -Jorge Cuadra -Stéven Toupin -Nicolás Cuello -mats esseldeurs -Alison Young +s-neilson <36410751+s-neilson@users.noreply.github.com> +Zachary Pellow Cox, Samuel Steven Rieder +Stéven Toupin +Jorge Cuadra