From d39f1d6cc8218800187ccd104573ce26142f898e Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Fri, 18 Oct 2024 12:17:01 +0200 Subject: [PATCH] Fix build warnings, reduce duplicate code, use 64 bit indices --- src/main/evolve.F90 | 2 +- src/utils/libphantom-amuse.F90 | 1481 ++++++++++++++++---------------- 2 files changed, 746 insertions(+), 737 deletions(-) diff --git a/src/main/evolve.F90 b/src/main/evolve.F90 index b7a3801f7..bc5fb5a78 100644 --- a/src/main/evolve.F90 +++ b/src/main/evolve.F90 @@ -139,11 +139,11 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag) integer :: j,nskip,nskipped,nevwrite_threshold,nskipped_sink,nsinkwrite_threshold character(len=120) :: dumpfile_orig + tzero = time if (.not. initialized) then tprint = 0. nsteps = 0 nsteplast = 0 - tzero = time dtlast = 0. dtinject = huge(dtinject) dtrad = huge(dtrad) diff --git a/src/utils/libphantom-amuse.F90 b/src/utils/libphantom-amuse.F90 index be8f8167a..bccb50c6d 100644 --- a/src/utils/libphantom-amuse.F90 +++ b/src/utils/libphantom-amuse.F90 @@ -3,10 +3,9 @@ module AmusePhantom - use iso_fortran_env integer(kind=8), allocatable :: amuse_id_lookup(:) - integer :: new_particles_since_last_update = 0 - integer :: particles_added_by_amuse = 0 + integer(kind=8) :: new_particles_since_last_update = 0 + integer(kind=8) :: particles_added_by_amuse = 0 contains subroutine construct_id_lookup() @@ -15,7 +14,7 @@ subroutine construct_id_lookup() use dim, only: maxp, maxp_hard use part, only: iorig, norig implicit none - integer :: i, j + integer (kind=8) :: i, j write(*,*) "Rebuilding lookup table" do i=1, maxp amuse_id_lookup(i) = 0 @@ -30,7 +29,6 @@ subroutine construct_id_lookup() enddo write(*,*) size(amuse_id_lookup), "?=", norig, maxp, maxp_hard write(*,*) "Lookup table rebuilt" - flush(output_unit) end subroutine construct_id_lookup subroutine amuse_initialize_code() @@ -268,10 +266,10 @@ end subroutine amuse_initialize_wind subroutine amuse_commit_parameters() use memory, only:allocate_memory use dim, only:maxp_hard - use initial, only:initialise use io, only:set_io_unit_numbers - use units, only:set_units,utime,umass,udist,set_units_extra - use physcon, only:gram,seconds,solarm,pc,au + use initial, only:initialise + use units, only:set_units,set_units_extra + use physcon, only:solarm,au call allocate_memory(int(maxp_hard,kind=8)) call set_io_unit_numbers() @@ -281,18 +279,24 @@ subroutine amuse_commit_parameters() ! Hard coded defaults, should be set in a different way... call set_units(dist=1 * au,mass=1.*solarm,G=1.) call set_units_extra() - ! write(*,*) "UNITS: ", udist, utime, umass end subroutine amuse_commit_parameters subroutine amuse_commit_particles() use part, only: norig use initial, only:startrun - use energies, only: mtot + use timestep, only:nmax,nsteps + use evolve, only:evol implicit none character(len=120) :: infile,logfile,evfile,dumpfile + integer :: nsteps_orig call startrun(infile,logfile,evfile,dumpfile, .true.) - call amuse_evol(.true.) + ! Make sure the evol call only initialises and doesn't do an evolve step + nsteps_orig = nsteps + nsteps = nmax + call evol(infile,logfile,evfile,dumpfile) + nsteps = nsteps_orig + call construct_id_lookup() new_particles_since_last_update = norig - particles_added_by_amuse end subroutine amuse_commit_particles @@ -304,7 +308,6 @@ subroutine amuse_recommit_particles() eos_vars,metrics,pxyzu,rad use timestep, only:time,dtmax implicit none - integer :: ierr integer :: nactive double precision :: dt, dtnew_first dtnew_first = dtmax @@ -327,7 +330,7 @@ end subroutine amuse_cleanup_code subroutine amuse_get_norig(norig_out) use part, only:norig implicit none - integer, intent(out) :: norig_out + integer(kind=8), intent(out) :: norig_out norig_out = norig end subroutine amuse_get_norig @@ -354,7 +357,7 @@ end subroutine amuse_get_npart ! Set default parameters subroutine amuse_set_defaults() - use options, only: set_default_options,iexternalforce,write_files + use options, only: set_default_options,write_files implicit none call set_default_options() @@ -363,635 +366,635 @@ subroutine amuse_set_defaults() end subroutine amuse_set_defaults ! This initialises things. This really should only be called once, before the first step. -subroutine amuse_evol(amuse_initialise) - use io, only:iprint,iwritein,id,master,iverbose,& - flush_warnings,nprocs,fatal,warning - use timestep, only:time,tmax,dt,dtmax,nmax,nout,nsteps,dtextforce,rhomaxnow,& - dtmax_ifactor,dtmax_ifactorWT,dtmax_dratio,check_dtmax_for_decrease,& - idtmax_n,idtmax_frac,idtmax_n_next,idtmax_frac_next,dtlast - use evwrite, only:write_evfile,write_evlog - use energies, only:etot,totmom,angtot,mdust,np_cs_eq_0,np_e_eq_0,hdivBonB_ave,hdivBonB_max,compute_energies - use checkconserved, only:etot_in,angtot_in,totmom_in,mdust_in,& - init_conservation_checks,check_conservation_error,& - check_magnetic_stability - use dim, only:maxvxyzu,mhd,periodic,idumpfile - use fileutils, only:getnextfilename - use options, only:nfulldump,twallmax,nmaxdumps,rhofinal1,iexternalforce,rkill - use readwrite_infile, only:write_infile - use readwrite_dumps, only:write_smalldump,write_fulldump - use step_lf_global, only:step - use timing, only:get_timings,print_time,timer,reset_timer,increment_timer,& - setup_timers,timers,reduce_timers,ntimers,& - itimer_fromstart,itimer_lastdump,itimer_step,itimer_io,itimer_ev - use mpiutils, only:reduce_mpi,reduceall_mpi,barrier_mpi,bcast_mpi -#ifdef IND_TIMESTEPS - use part, only:ibin,iphase - use timestep_ind, only:istepfrac,nbinmax,set_active_particles,update_time_per_bin,& - write_binsummary,change_nbinmax,nactive,nactivetot,maxbins,& - print_dtlog_ind,get_newbin,print_dtind_efficiency - use timestep, only:dtdiff - use timestep_sts, only:sts_get_dtau_next,sts_init_step - use step_lf_global, only:init_step -#else - use timestep, only:dtforce,dtcourant,dterr,print_dtlog -#endif - use timestep_sts, only:use_sts - use supertimestep, only:step_sts -#ifdef DRIVING - use forcing, only:write_forcingdump -#endif -#ifdef CORRECT_BULK_MOTION - use centreofmass, only:correct_bulk_motion -#endif - use part, only:ideadhead,shuffle_part -#ifdef INJECT_PARTICLES - use inject, only:inject_particles - use part, only:npartoftype - use partinject, only:update_injected_particles -#endif - use dim, only:do_radiation - use options, only:exchange_radiation_energy,implicit_radiation - use part, only:rad,radprop - use radiation_utils, only:update_radenergy - use timestep, only:dtrad -#ifdef LIVE_ANALYSIS - use analysis, only:do_analysis - use part, only:igas - use fileutils, only:numfromfile - use io, only:ianalysis -#endif - use part, only:npart,nptmass,xyzh,vxyzu,fxyzu,fext,divcurlv,massoftype, & - xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dptmass,gravity,iboundary, & - fxyz_ptmass_sinksink,ntot,poten,ndustsmall,accrete_particles_outside_sphere - use quitdump, only:quit - use ptmass, only:icreate_sinks,ptmass_create,ipart_rhomax,pt_write_sinkev,calculate_mdot, & - set_integration_precision - use io_summary, only:iosum_nreal,summary_counter,summary_printout,summary_printnow - use externalforces, only:iext_spiral - use boundary_dyn, only:dynamic_bdy,update_boundaries -#ifdef MFLOW - use mf_write, only:mflow_write -#endif -#ifdef VMFLOW - use mf_write, only:vmflow_write -#endif -#ifdef BINPOS - use mf_write, only:binpos_write -#endif - implicit none - logical, intent(in) :: amuse_initialise - logical :: initialized - - integer :: flag - character(len=256) :: infile - character(len=256) :: logfile,evfile,dumpfile - integer :: i,noutput,noutput_dtmax,nsteplast,ncount_fulldumps - real :: dtnew,timecheck,rhomaxold,dtmax_log_dratio - real :: tprint,tzero,dtmaxold,dtinject - real(kind=4) :: t1,t2,tcpu1,tcpu2,tstart,tcpustart - real(kind=4) :: twalllast,tcpulast,twallperdump,twallused -#ifdef IND_TIMESTEPS - integer :: nalive,inbin - integer(kind=1) :: nbinmaxprev - integer(kind=8) :: nmovedtot,nalivetot - real :: tlast,tcheck,dtau - real(kind=4) :: tall - real(kind=4) :: timeperbin(0:maxbins) - logical :: dt_changed -#else - real :: tlast - real :: dtprint - integer :: nactive,istepfrac - integer(kind=1) :: nbinmax - logical, parameter :: dt_changed = .false. -#endif -#ifdef INJECT_PARTICLES - integer :: npart_old -#endif - logical :: fulldump,abortrun,abortrun_bdy,at_dump_time,writedump - logical :: should_conserve_energy,should_conserve_momentum,should_conserve_angmom - logical :: should_conserve_dustmass - logical :: use_global_dt - integer :: j,nskip,nskipped,nevwrite_threshold,nskipped_sink,nsinkwrite_threshold - character(len=120) :: dumpfile_orig - - initialized = .not. amuse_initialise - tzero = time - tlast = time - write(*,*) "AMUSE_EVOL called, amuse_initialise=", amuse_initialise - write(*,*) "dtlast", dtlast - write(*,*) "tzero, tmax: ", tzero, dtmax, tmax - if (.not. initialized) then - tprint = 0. - nsteps = 0 - nsteplast = 0 - dtlast = 0. - dtinject = huge(dtinject) - dtrad = huge(dtrad) - np_cs_eq_0 = 0 - np_e_eq_0 = 0 - abortrun_bdy = .false. - - call init_conservation_checks(should_conserve_energy,should_conserve_momentum,& - should_conserve_angmom,should_conserve_dustmass) - - noutput = 1 - noutput_dtmax = 1 - ncount_fulldumps = 0 - tprint = tzero + dtmax - rhomaxold = rhomaxnow - if (dtmax_dratio > 0.) then - dtmax_log_dratio = log10(dtmax_dratio) - else - dtmax_log_dratio = 0.0 - endif - - ! - ! Set substepping integration precision depending on the system (default is FSI) - ! - call set_integration_precision - -#ifdef IND_TIMESTEPS - write(*,*) "Using individual timesteps" - use_global_dt = .false. - istepfrac = 0 - tlast = time - write(*,*) "***********tlast: ", tlast - dt = dtmax/2.**nbinmax ! use 2.0 here to allow for step too small - nmovedtot = 0 - tall = 0. - tcheck = time - timeperbin(:) = 0. - dt_changed = .false. - call init_step(npart,time,dtmax) - if (use_sts) then - call sts_get_dtau_next(dtau,dt,dtmax,dtdiff,nbinmax) - call sts_init_step(npart,time,dtmax,dtau) ! overwrite twas for particles requiring super-timestepping - endif -#else - write(*,*) "Not using individual timesteps" - use_global_dt = .true. - nskip = int(ntot) - nactive = npart - istepfrac = 0 ! dummy values - nbinmax = 0 - write(*,*) "dt, tprint, time:", dt, tprint, time - if (dt >= (tprint-time)) dt = tprint-time ! reach tprint exactly -#endif -! -! threshold for writing to .ev file, to avoid repeatedly computing energies -! for all the particles which would add significantly to the cpu time -! - - nskipped = 0 - if (iexternalforce==iext_spiral) then - nevwrite_threshold = int(4.99*ntot) ! every 5 full steps - else - nevwrite_threshold = int(1.99*ntot) ! every 2 full steps - endif - nskipped_sink = 0 - nsinkwrite_threshold = int(0.99*ntot) -! -! code timings -! - call get_timings(twalllast,tcpulast) - tstart = twalllast - tcpustart = tcpulast - - call setup_timers - - call flush(iprint) - - else ! i.e.: initializing done - !dtmax = min(dtmax, tmax-time) ! tried this, doesn't work. - !timesubstepping: do while (istepfrac < 2**nbinmax) - - timestepping: do while ((time + 0.01 * dtmax < tmax).and.((nsteps < nmax) .or. (nmax < 0)).and.(rhomaxnow*rhofinal1 < 1.0)) - -#ifdef INJECT_PARTICLES - ! - ! injection of new particles into simulation - ! - !if (.not. present(flag)) then - npart_old=npart - call inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,npart,npart_old,npartoftype,dtinject) - call update_injected_particles(npart_old,npart,istepfrac,nbinmax,time,dtmax,dt,dtinject) - !endif -#endif - - dtmaxold = dtmax -#ifdef IND_TIMESTEPS - istepfrac = istepfrac + 1 - nbinmaxprev = nbinmax - if (nbinmax > maxbins) call fatal('evolve','timestep too small: try decreasing dtmax?') - - !--determine if dt needs to be decreased; if so, then this will be done - ! in step the next time it is called; - ! for global timestepping, this is called in the block where at_dump_time==.true. - if (istepfrac==2**nbinmax) then - twallperdump = reduceall_mpi('max', timers(itimer_lastdump)%wall) - call check_dtmax_for_decrease(iprint,dtmax,twallperdump,dtmax_log_dratio,& - rhomaxold,rhomaxnow,nfulldump,use_global_dt) - endif - - !--sanity check on istepfrac... - if (istepfrac > 2**nbinmax) then - write(iprint,*) 'ERROR: istepfrac = ',istepfrac,' / ',2**nbinmax - call fatal('evolve','error in individual timesteps') - endif - - !--flag particles as active or not for this timestep - call set_active_particles(npart,nactive,nalive,iphase,ibin,xyzh) - nactivetot = reduceall_mpi('+', nactive) - nalivetot = reduceall_mpi('+', nalive) - nskip = int(nactivetot) - - !--print summary of timestep bins - if (iverbose >= 2) call write_binsummary(npart,nbinmax,dtmax,timeperbin,iphase,ibin,xyzh) - - !--Implement dynamic boundaries (for individual-timestepping) once per dump - if (dynamic_bdy .and. nactive==nalive .and. istepfrac==2**nbinmax) then - call update_boundaries(nactive,nalive,npart,abortrun_bdy) - endif -#else - !--If not using individual timestepping, set nskip to the total number of particles - ! across all nodes - nskip = int(ntot) -#endif - - if (gravity .and. icreate_sinks > 0 .and. ipart_rhomax /= 0) then - ! - ! creation of new sink particles - ! - call ptmass_create(nptmass,npart,ipart_rhomax,xyzh,vxyzu,fxyzu,fext,divcurlv,& - poten,massoftype,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fxyz_ptmass_sinksink,dptmass,time) - endif - ! - ! Strang splitting: implicit update for half step - ! - if (do_radiation .and. exchange_radiation_energy .and. .not.implicit_radiation) then - call update_radenergy(npart,xyzh,fxyzu,vxyzu,rad,radprop,0.5*dt) - endif - nsteps = nsteps + 1 -! -!--evolve data for one timestep -! for individual timesteps this is the shortest timestep -! - call get_timings(t1,tcpu1) - if ( use_sts ) then - call step_sts(npart,nactive,time,dt,dtextforce,dtnew,iprint) - else - call step(npart,nactive,time,dt,dtextforce,dtnew) - endif - ! - ! Strang splitting: implicit update for another half step - ! - if (do_radiation .and. exchange_radiation_energy .and. .not.implicit_radiation) then - call update_radenergy(npart,xyzh,fxyzu,vxyzu,rad,radprop,0.5*dt) - endif - - dtlast = dt - - !--timings for step call - call get_timings(t2,tcpu2) - call increment_timer(itimer_step,t2-t1,tcpu2-tcpu1) - call summary_counter(iosum_nreal,t2-t1) - -#ifdef IND_TIMESTEPS - tcheck = tcheck + dt - - !--update time in way that is free of round-off errors - time = tlast + istepfrac/real(2**nbinmaxprev)*dtmaxold - - !--print efficiency of partial timestep - if (id==master) call print_dtind_efficiency(iverbose,nalivetot,nactivetot,tall,t2-t1,1) - - call update_time_per_bin(tcpu2-tcpu1,istepfrac,nbinmaxprev,timeperbin,inbin) - nmovedtot = nmovedtot + nactivetot - - !--check that time is as it should be, may indicate error in individual timestep routines - if (abs(tcheck-time) > 1.e-4) call warning('evolve','time out of sync',var='error',val=abs(tcheck-time)) - - if (id==master .and. (iverbose >= 1 .or. inbin <= 3)) & - call print_dtlog_ind(iprint,istepfrac,2**nbinmaxprev,time,dt,nactivetot,tcpu2-tcpu1,ntot) - - !--if total number of bins has changed, adjust istepfrac and dt accordingly - ! (ie., decrease or increase the timestep) - if (nbinmax /= nbinmaxprev .or. dtmax_ifactor /= 0) then - call change_nbinmax(nbinmax,nbinmaxprev,istepfrac,dtmax,dt) - dt_changed = .true. - endif - -#else - - ! advance time on master thread only - if (id == master) time = time + dt - call bcast_mpi(time) - -! -!--set new timestep from Courant/forces condition -! - ! constraint from time to next printout, must reach this exactly - ! Following redefinitions are to avoid crashing if dtprint = 0 & to reach next output while avoiding round-off errors - dtprint = min(tprint,tmax) - time + epsilon(dtmax) - if (dtprint <= epsilon(dtmax) .or. dtprint >= (1.0-1e-8)*dtmax ) dtprint = dtmax + epsilon(dtmax) - dt = min(dtforce,dtcourant,dterr,dtmax+epsilon(dtmax),dtprint,dtinject,dtrad) -! -!--write log every step (NB: must print after dt has been set in order to identify timestep constraint) -! - if (id==master) call print_dtlog(iprint,time,dt,dtforce,dtcourant,dterr,dtmax,dtrad,dtprint,dtinject,ntot) -#endif - -! check that MPI threads are synchronised in time - timecheck = reduceall_mpi('+',time) - if (abs(timecheck/nprocs - time) > 1.e-13) then - call fatal('evolve','time differs between MPI threads',var='time',val=timecheck/nprocs) - endif -! -!--Update timer from last dump to see if dtmax needs to be reduced -! - call get_timings(t2,tcpu2) - call increment_timer(itimer_lastdump,t2-t1,tcpu2-tcpu1) -! -!--Determine if this is the correct time to write to the data file -! - at_dump_time = (time >= tmax) & - .or.((nsteps >= nmax).and.(nmax >= 0)).or.(rhomaxnow*rhofinal1 >= 1.0) -#ifdef IND_TIMESTEPS - if (istepfrac==2**nbinmax) at_dump_time = .true. -#else - if (time >= tprint) at_dump_time = .true. -#endif -! -!--Calculate total energy etc and write to ev file -! For individual timesteps, we do not want to do this every step, but we want -! to do this as often as possible without a performance hit. The criteria -! here is that it is done once >10% of particles (cumulatively) have been evolved. -! That is, either >10% are being stepped, or e.g. 1% have moved 10 steps. -! Perform this prior to writing the dump files so that diagnostic values calculated -! in energies can be correctly included in the dumpfiles -! - nskipped = nskipped + nskip - if (nskipped >= nevwrite_threshold .or. at_dump_time .or. dt_changed .or. iverbose==5) then - nskipped = 0 - call get_timings(t1,tcpu1) - ! If we don't want to write the evfile, we do still want to calculate the energies - if (write_files) then - call write_evfile(time,dt) - else - call compute_energies(time) - endif - if (should_conserve_momentum) call check_conservation_error(totmom,totmom_in,1.e-1,'linear momentum') - if (should_conserve_angmom) call check_conservation_error(angtot,angtot_in,1.e-1,'angular momentum') - if (should_conserve_energy) call check_conservation_error(etot,etot_in,1.e-1,'energy') - if (should_conserve_dustmass) then - do j = 1,ndustsmall - call check_conservation_error(mdust(j),mdust_in(j),1.e-1,'dust mass',decrease=.true.) - enddo - endif - if (mhd) call check_magnetic_stability(hdivBonB_ave,hdivBonB_max) - if (id==master) then - if (np_e_eq_0 > 0) call warning('evolve','N gas particles with energy = 0',var='N',ival=int(np_e_eq_0,kind=4)) - if (np_cs_eq_0 > 0) call warning('evolve','N gas particles with sound speed = 0',var='N',ival=int(np_cs_eq_0,kind=4)) - endif - - !--write with the same ev file frequency also mass flux and binary position -#ifdef MFLOW - call mflow_write(time,dt) -#endif -#ifdef VMFLOW - call vmflow_write(time,dt) -#endif -#ifdef BINPOS - call binpos_write(time,dt) -#endif - call get_timings(t2,tcpu2) - call increment_timer(itimer_ev,t2-t1,tcpu2-tcpu1) ! time taken for write_ev operation - endif -!-- Print out the sink particle properties & reset dt_changed. -!-- Added total force on sink particles and sink-sink forces to write statement (fxyz_ptmass,fxyz_ptmass_sinksink) - nskipped_sink = nskipped_sink + nskip - if (nskipped_sink >= nsinkwrite_threshold .or. at_dump_time .or. dt_changed) then - nskipped_sink = 0 - if (write_files) call pt_write_sinkev(nptmass,time,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fxyz_ptmass_sinksink) -#ifdef IND_TIMESTEPS - dt_changed = .false. -#endif - write(*,*) "***********tlast, dtlast, tmax: ", tlast, dtlast, tmax - endif -! -!--write to data file if time is right -! - if (at_dump_time .and. write_files) then -#ifndef IND_TIMESTEPS -! -!--Global timesteps: Decrease dtmax if requested (done in step for individual timesteps) - twallperdump = timers(itimer_lastdump)%wall - call check_dtmax_for_decrease(iprint,dtmax,twallperdump,dtmax_log_dratio,& - rhomaxold,rhomaxnow,nfulldump,use_global_dt) - dt = min(dt,dtmax) ! required if decreasing dtmax to ensure that the physically motivated timestep is not too long -#endif - - !--modify evfile and logfile names with new number - if ((nout <= 0) .or. (mod(noutput,nout)==0)) then - if (noutput==1) then - evfile = getnextfilename(evfile) - logfile = getnextfilename(logfile) - endif -! Update values for restart dumps - if (dtmax_ifactorWT ==0) then - idtmax_n_next = idtmax_n - idtmax_frac_next = idtmax_frac - elseif (dtmax_ifactorWT > 0) then - idtmax_n_next = idtmax_n *dtmax_ifactorWT - idtmax_frac_next = idtmax_frac*dtmax_ifactorWT - elseif (dtmax_ifactorWT < 0) then - idtmax_n_next = -idtmax_n /dtmax_ifactorWT - idtmax_frac_next = -idtmax_frac/dtmax_ifactorWT - endif - idtmax_frac_next = idtmax_frac_next + 1 - idtmax_frac_next = mod(idtmax_frac_next,idtmax_n_next) - dumpfile_orig = trim(dumpfile) - if (idtmax_frac==0) then - dumpfile = getnextfilename(dumpfile,idumpfile) - dumpfile_orig = trim(dumpfile) - else - write(dumpfile,'(2a)') dumpfile(:index(dumpfile,'_')-1),'.restart' - endif - writedump = .true. - else - writedump = .false. - endif - - !--do not dump dead particles into dump files - if (ideadhead > 0) call shuffle_part(npart) -! -!--get timings since last dump and overall code scaling -! (get these before writing the dump so we can check whether or not we -! need to write a full dump based on the wall time; -! move timer_lastdump outside at_dump_time block so that dtmax can -! be reduced it too long between dumps) -! - call increment_timer(itimer_fromstart,t2-tstart,tcpu2-tcpustart) - - fulldump = (nout <= 0 .and. mod(noutput,nfulldump)==0) .or. (mod(noutput,nout*nfulldump)==0) -! -!--if max wall time is set (> 1 sec) stop the run at the last full dump -! that will fit into the walltime constraint, based on the wall time between -! the last two dumps added to the current total walltime used. The factor of three for -! changing to full dumps is to account for the possibility that the next step will take longer. -! If we are about to write a small dump but it looks like we won't make the next dump, -! write a full dump instead and stop the run -! - abortrun = .false. - if (twallmax > 1.) then - twallused = timers(itimer_fromstart)%wall - twallperdump = timers(itimer_lastdump)%wall - if (fulldump) then - if ((twallused + abs(nfulldump)*twallperdump) > twallmax) then - abortrun = .true. - endif - else - if ((twallused + 3.0*twallperdump) > twallmax) then - fulldump = .true. - if (id==master) write(iprint,"(1x,a)") '>> PROMOTING DUMP TO FULL DUMP BASED ON WALL TIME CONSTRAINTS... ' - nfulldump = 1 ! also set all future dumps to be full dumps (otherwise gets confusing) - if ((twallused + twallperdump) > twallmax) abortrun = .true. - endif - endif - endif -! -!--Promote to full dump if this is the final dump -! - if ( (time >= tmax) .or. ( (nmax > 0) .and. (nsteps >= nmax) ) ) fulldump = .true. -! -!--flush any buffered warnings to the log file -! - if (id==master) call flush_warnings() -! -!--write dump file -! - if (rkill > 0) call accrete_particles_outside_sphere(rkill) -#ifndef INJECT_PARTICLES - call calculate_mdot(nptmass,time,xyzmh_ptmass) -#endif - call get_timings(t1,tcpu1) - if (writedump) then - if (fulldump) then - call write_fulldump(time,dumpfile) - if (id==master) then - call write_infile(infile,logfile,evfile,dumpfile,iwritein,iprint) -#ifdef DRIVING - call write_forcingdump(time,dumpfile) -#endif - endif - ncount_fulldumps = ncount_fulldumps + 1 - else - call write_smalldump(time,dumpfile) - endif - endif - call get_timings(t2,tcpu2) - call increment_timer(itimer_io,t2-t1,tcpu2-tcpu1) - -#ifdef LIVE_ANALYSIS - if (id==master .and. idtmax_frac==0) then - call do_analysis(dumpfile,numfromfile(dumpfile),xyzh,vxyzu, & - massoftype(igas),npart,time,ianalysis) - endif -#endif - call reduce_timers - if (id==master) then - call print_timinginfo(iprint,nsteps,nsteplast) - !--Write out summary to log file - call summary_printout(iprint,nptmass) - endif -#ifdef IND_TIMESTEPS - !--print summary of timestep bins - if (iverbose >= 0) then - call write_binsummary(npart,nbinmax,dtmax,timeperbin,iphase,ibin,xyzh) - timeperbin(:) = 0. - if (id==master) call print_dtind_efficiency(iverbose,nalivetot,nmovedtot,tall,timers(itimer_lastdump)%wall,2) - endif - tlast = tprint - istepfrac = 0 - nmovedtot = 0 -#endif - !--print summary of energies and other useful values to the log file - if (id==master) call write_evlog(iprint) - -#ifndef IND_TIMESTEPS - !--Implement dynamic boundaries (for global timestepping) - if (dynamic_bdy) call update_boundaries(nactive,nactive,npart,abortrun_bdy) -#endif - - ! - !--if twallmax > 1s stop the run at the last full dump that will fit into the walltime constraint, - ! based on the wall time between the last two dumps added to the current total walltime used. - ! - if (abortrun) then - if (id==master) then - call print_time(t2-tstart,'>> WALL TIME = ',iprint) - call print_time(twallmax,'>> NEXT DUMP WILL TRIP OVER MAX WALL TIME: ',iprint) - write(iprint,"(1x,a)") '>> ABORTING... ' - endif - return - endif - - if (abortrun_bdy) then - write(iprint,"(1x,a)") 'Will likely surpass maxp_hard next time we need to add particles.' - write(iprint,"(1x,a)") 'Recompile with larger maxp_hard.' - write(iprint,"(1x,a)") '>> ABORTING... ' - return - endif - - if (nmaxdumps > 0 .and. ncount_fulldumps >= nmaxdumps) then - if (id==master) write(iprint,"(a)") '>> reached maximum number of full dumps as specified in input file, stopping...' - return - endif - - twalllast = t2 - tcpulast = tcpu2 - do i = 1,ntimers - call reset_timer(i) - enddo - - if (idtmax_frac==0) then - noutput = noutput + 1 ! required to determine frequency of full dumps - endif - noutput_dtmax = noutput_dtmax + 1 ! required to adjust tprint; will account for varying dtmax - idtmax_n = idtmax_n_next - idtmax_frac = idtmax_frac_next - tprint = tzero + noutput_dtmax*dtmaxold - nsteplast = nsteps - dumpfile = trim(dumpfile_orig) - if (dtmax_ifactor/=0) then - tzero = tprint - dtmaxold - tprint = tzero + dtmax - noutput_dtmax = 1 - dtmax_ifactor = 0 - dtmax_ifactorWT = 0 - endif - endif - -#ifdef CORRECT_BULK_MOTION - call correct_bulk_motion() -#endif - - if (iverbose >= 1 .and. id==master) write(iprint,*) - call flush(iprint) - !--Write out log file prematurely (if requested based upon nstep, walltime) - if ( summary_printnow() ) call summary_printout(iprint,nptmass) - - enddo timestepping - -endif -end subroutine amuse_evol +!!!subroutine amuse_evol(amuse_initialise) +!!! use io, only:iprint,iwritein,id,master,iverbose,& +!!! flush_warnings,nprocs,fatal,warning +!!! use timestep, only:time,tmax,dt,dtmax,nmax,nout,nsteps,dtextforce,rhomaxnow,& +!!! dtmax_ifactor,dtmax_ifactorWT,dtmax_dratio,check_dtmax_for_decrease,& +!!! idtmax_n,idtmax_frac,idtmax_n_next,idtmax_frac_next,dtlast +!!! use evwrite, only:write_evfile,write_evlog +!!! use energies, only:etot,totmom,angtot,mdust,np_cs_eq_0,np_e_eq_0,hdivBonB_ave,hdivBonB_max,compute_energies +!!! use checkconserved, only:etot_in,angtot_in,totmom_in,mdust_in,& +!!! init_conservation_checks,check_conservation_error,& +!!! check_magnetic_stability +!!! use dim, only:maxvxyzu,mhd,periodic,idumpfile +!!! use fileutils, only:getnextfilename +!!! use options, only:nfulldump,twallmax,nmaxdumps,rhofinal1,iexternalforce,rkill +!!! use readwrite_infile, only:write_infile +!!! use readwrite_dumps, only:write_smalldump,write_fulldump +!!! use step_lf_global, only:step +!!! use timing, only:get_timings,print_time,timer,reset_timer,increment_timer,& +!!! setup_timers,timers,reduce_timers,ntimers,& +!!! itimer_fromstart,itimer_lastdump,itimer_step,itimer_io,itimer_ev +!!! use mpiutils, only:reduce_mpi,reduceall_mpi,barrier_mpi,bcast_mpi +!!!#ifdef IND_TIMESTEPS +!!! use part, only:ibin,iphase +!!! use timestep_ind, only:istepfrac,nbinmax,set_active_particles,update_time_per_bin,& +!!! write_binsummary,change_nbinmax,nactive,nactivetot,maxbins,& +!!! print_dtlog_ind,get_newbin,print_dtind_efficiency +!!! use timestep, only:dtdiff +!!! use timestep_sts, only:sts_get_dtau_next,sts_init_step +!!! use step_lf_global, only:init_step +!!!#else +!!! use timestep, only:dtforce,dtcourant,dterr,print_dtlog +!!!#endif +!!! use timestep_sts, only:use_sts +!!! use supertimestep, only:step_sts +!!!#ifdef DRIVING +!!! use forcing, only:write_forcingdump +!!!#endif +!!!#ifdef CORRECT_BULK_MOTION +!!! use centreofmass, only:correct_bulk_motion +!!!#endif +!!! use part, only:ideadhead,shuffle_part +!!!#ifdef INJECT_PARTICLES +!!! use inject, only:inject_particles +!!! use part, only:npartoftype +!!! use partinject, only:update_injected_particles +!!!#endif +!!! use dim, only:do_radiation +!!! use options, only:exchange_radiation_energy,implicit_radiation +!!! use part, only:rad,radprop +!!! use radiation_utils, only:update_radenergy +!!! use timestep, only:dtrad +!!!#ifdef LIVE_ANALYSIS +!!! use analysis, only:do_analysis +!!! use part, only:igas +!!! use fileutils, only:numfromfile +!!! use io, only:ianalysis +!!!#endif +!!! use part, only:npart,nptmass,xyzh,vxyzu,fxyzu,fext,divcurlv,massoftype, & +!!! xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dptmass,gravity,iboundary, & +!!! fxyz_ptmass_sinksink,ntot,poten,ndustsmall,accrete_particles_outside_sphere +!!! use quitdump, only:quit +!!! use ptmass, only:icreate_sinks,ptmass_create,ipart_rhomax,pt_write_sinkev,calculate_mdot, & +!!! set_integration_precision +!!! use io_summary, only:iosum_nreal,summary_counter,summary_printout,summary_printnow +!!! use externalforces, only:iext_spiral +!!! use boundary_dyn, only:dynamic_bdy,update_boundaries +!!!#ifdef MFLOW +!!! use mf_write, only:mflow_write +!!!#endif +!!!#ifdef VMFLOW +!!! use mf_write, only:vmflow_write +!!!#endif +!!!#ifdef BINPOS +!!! use mf_write, only:binpos_write +!!!#endif +!!! implicit none +!!! logical, intent(in) :: amuse_initialise +!!! logical :: initialized +!!! +!!! integer :: flag +!!! character(len=256) :: infile +!!! character(len=256) :: logfile,evfile,dumpfile +!!! integer :: i,noutput,noutput_dtmax,nsteplast,ncount_fulldumps +!!! real :: dtnew,timecheck,rhomaxold,dtmax_log_dratio +!!! real :: tprint,tzero,dtmaxold,dtinject +!!! real(kind=4) :: t1,t2,tcpu1,tcpu2,tstart,tcpustart +!!! real(kind=4) :: twalllast,tcpulast,twallperdump,twallused +!!!#ifdef IND_TIMESTEPS +!!! integer :: nalive,inbin +!!! integer(kind=1) :: nbinmaxprev +!!! integer(kind=8) :: nmovedtot,nalivetot +!!! real :: tlast,tcheck,dtau +!!! real(kind=4) :: tall +!!! real(kind=4) :: timeperbin(0:maxbins) +!!! logical :: dt_changed +!!!#else +!!! real :: tlast +!!! real :: dtprint +!!! integer :: nactive,istepfrac +!!! integer(kind=1) :: nbinmax +!!! logical, parameter :: dt_changed = .false. +!!!#endif +!!!#ifdef INJECT_PARTICLES +!!! integer :: npart_old +!!!#endif +!!! logical :: fulldump,abortrun,abortrun_bdy,at_dump_time,writedump +!!! logical :: should_conserve_energy,should_conserve_momentum,should_conserve_angmom +!!! logical :: should_conserve_dustmass +!!! logical :: use_global_dt +!!! integer :: j,nskip,nskipped,nevwrite_threshold,nskipped_sink,nsinkwrite_threshold +!!! character(len=120) :: dumpfile_orig +!!! +!!! initialized = .not. amuse_initialise +!!! tzero = time +!!! tlast = time +!!! write(*,*) "AMUSE_EVOL called, amuse_initialise=", amuse_initialise +!!! write(*,*) "dtlast", dtlast +!!! write(*,*) "tzero, tmax: ", tzero, dtmax, tmax +!!! if (.not. initialized) then +!!! tprint = 0. +!!! nsteps = 0 +!!! nsteplast = 0 +!!! dtlast = 0. +!!! dtinject = huge(dtinject) +!!! dtrad = huge(dtrad) +!!! np_cs_eq_0 = 0 +!!! np_e_eq_0 = 0 +!!! abortrun_bdy = .false. +!!! +!!! call init_conservation_checks(should_conserve_energy,should_conserve_momentum,& +!!! should_conserve_angmom,should_conserve_dustmass) +!!! +!!! noutput = 1 +!!! noutput_dtmax = 1 +!!! ncount_fulldumps = 0 +!!! tprint = tzero + dtmax +!!! rhomaxold = rhomaxnow +!!! if (dtmax_dratio > 0.) then +!!! dtmax_log_dratio = log10(dtmax_dratio) +!!! else +!!! dtmax_log_dratio = 0.0 +!!! endif +!!! +!!! ! +!!! ! Set substepping integration precision depending on the system (default is FSI) +!!! ! +!!! call set_integration_precision +!!! +!!!#ifdef IND_TIMESTEPS +!!! write(*,*) "Using individual timesteps" +!!! use_global_dt = .false. +!!! istepfrac = 0 +!!! tlast = time +!!! write(*,*) "***********tlast: ", tlast +!!! dt = dtmax/2.**nbinmax ! use 2.0 here to allow for step too small +!!! nmovedtot = 0 +!!! tall = 0. +!!! tcheck = time +!!! timeperbin(:) = 0. +!!! dt_changed = .false. +!!! call init_step(npart,time,dtmax) +!!! if (use_sts) then +!!! call sts_get_dtau_next(dtau,dt,dtmax,dtdiff,nbinmax) +!!! call sts_init_step(npart,time,dtmax,dtau) ! overwrite twas for particles requiring super-timestepping +!!! endif +!!!#else +!!! write(*,*) "Not using individual timesteps" +!!! use_global_dt = .true. +!!! nskip = int(ntot) +!!! nactive = npart +!!! istepfrac = 0 ! dummy values +!!! nbinmax = 0 +!!! write(*,*) "dt, tprint, time:", dt, tprint, time +!!! if (dt >= (tprint-time)) dt = tprint-time ! reach tprint exactly +!!!#endif +!!!! +!!!! threshold for writing to .ev file, to avoid repeatedly computing energies +!!!! for all the particles which would add significantly to the cpu time +!!!! +!!! +!!! nskipped = 0 +!!! if (iexternalforce==iext_spiral) then +!!! nevwrite_threshold = int(4.99*ntot) ! every 5 full steps +!!! else +!!! nevwrite_threshold = int(1.99*ntot) ! every 2 full steps +!!! endif +!!! nskipped_sink = 0 +!!! nsinkwrite_threshold = int(0.99*ntot) +!!!! +!!!! code timings +!!!! +!!! call get_timings(twalllast,tcpulast) +!!! tstart = twalllast +!!! tcpustart = tcpulast +!!! +!!! call setup_timers +!!! +!!! call flush(iprint) +!!! +!!! else ! i.e.: initializing done +!!! !dtmax = min(dtmax, tmax-time) ! tried this, doesn't work. +!!! !timesubstepping: do while (istepfrac < 2**nbinmax) +!!! +!!! timestepping: do while ((time + 0.01 * dtmax < tmax).and.((nsteps < nmax) .or. (nmax < 0)).and.(rhomaxnow*rhofinal1 < 1.0)) +!!! +!!!#ifdef INJECT_PARTICLES +!!! ! +!!! ! injection of new particles into simulation +!!! ! +!!! !if (.not. present(flag)) then +!!! npart_old=npart +!!! call inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,npart,npart_old,npartoftype,dtinject) +!!! call update_injected_particles(npart_old,npart,istepfrac,nbinmax,time,dtmax,dt,dtinject) +!!! !endif +!!!#endif +!!! +!!! dtmaxold = dtmax +!!!#ifdef IND_TIMESTEPS +!!! istepfrac = istepfrac + 1 +!!! nbinmaxprev = nbinmax +!!! if (nbinmax > maxbins) call fatal('evolve','timestep too small: try decreasing dtmax?') +!!! +!!! !--determine if dt needs to be decreased; if so, then this will be done +!!! ! in step the next time it is called; +!!! ! for global timestepping, this is called in the block where at_dump_time==.true. +!!! if (istepfrac==2**nbinmax) then +!!! twallperdump = reduceall_mpi('max', timers(itimer_lastdump)%wall) +!!! call check_dtmax_for_decrease(iprint,dtmax,twallperdump,dtmax_log_dratio,& +!!! rhomaxold,rhomaxnow,nfulldump,use_global_dt) +!!! endif +!!! +!!! !--sanity check on istepfrac... +!!! if (istepfrac > 2**nbinmax) then +!!! write(iprint,*) 'ERROR: istepfrac = ',istepfrac,' / ',2**nbinmax +!!! call fatal('evolve','error in individual timesteps') +!!! endif +!!! +!!! !--flag particles as active or not for this timestep +!!! call set_active_particles(npart,nactive,nalive,iphase,ibin,xyzh) +!!! nactivetot = reduceall_mpi('+', nactive) +!!! nalivetot = reduceall_mpi('+', nalive) +!!! nskip = int(nactivetot) +!!! +!!! !--print summary of timestep bins +!!! if (iverbose >= 2) call write_binsummary(npart,nbinmax,dtmax,timeperbin,iphase,ibin,xyzh) +!!! +!!! !--Implement dynamic boundaries (for individual-timestepping) once per dump +!!! if (dynamic_bdy .and. nactive==nalive .and. istepfrac==2**nbinmax) then +!!! call update_boundaries(nactive,nalive,npart,abortrun_bdy) +!!! endif +!!!#else +!!! !--If not using individual timestepping, set nskip to the total number of particles +!!! ! across all nodes +!!! nskip = int(ntot) +!!!#endif +!!! +!!! if (gravity .and. icreate_sinks > 0 .and. ipart_rhomax /= 0) then +!!! ! +!!! ! creation of new sink particles +!!! ! +!!! call ptmass_create(nptmass,npart,ipart_rhomax,xyzh,vxyzu,fxyzu,fext,divcurlv,& +!!! poten,massoftype,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fxyz_ptmass_sinksink,dptmass,time) +!!! endif +!!! ! +!!! ! Strang splitting: implicit update for half step +!!! ! +!!! if (do_radiation .and. exchange_radiation_energy .and. .not.implicit_radiation) then +!!! call update_radenergy(npart,xyzh,fxyzu,vxyzu,rad,radprop,0.5*dt) +!!! endif +!!! nsteps = nsteps + 1 +!!!! +!!!!--evolve data for one timestep +!!!! for individual timesteps this is the shortest timestep +!!!! +!!! call get_timings(t1,tcpu1) +!!! if ( use_sts ) then +!!! call step_sts(npart,nactive,time,dt,dtextforce,dtnew,iprint) +!!! else +!!! call step(npart,nactive,time,dt,dtextforce,dtnew) +!!! endif +!!! ! +!!! ! Strang splitting: implicit update for another half step +!!! ! +!!! if (do_radiation .and. exchange_radiation_energy .and. .not.implicit_radiation) then +!!! call update_radenergy(npart,xyzh,fxyzu,vxyzu,rad,radprop,0.5*dt) +!!! endif +!!! +!!! dtlast = dt +!!! +!!! !--timings for step call +!!! call get_timings(t2,tcpu2) +!!! call increment_timer(itimer_step,t2-t1,tcpu2-tcpu1) +!!! call summary_counter(iosum_nreal,t2-t1) +!!! +!!!#ifdef IND_TIMESTEPS +!!! tcheck = tcheck + dt +!!! +!!! !--update time in way that is free of round-off errors +!!! time = tlast + istepfrac/real(2**nbinmaxprev)*dtmaxold +!!! +!!! !--print efficiency of partial timestep +!!! if (id==master) call print_dtind_efficiency(iverbose,nalivetot,nactivetot,tall,t2-t1,1) +!!! +!!! call update_time_per_bin(tcpu2-tcpu1,istepfrac,nbinmaxprev,timeperbin,inbin) +!!! nmovedtot = nmovedtot + nactivetot +!!! +!!! !--check that time is as it should be, may indicate error in individual timestep routines +!!! if (abs(tcheck-time) > 1.e-4) call warning('evolve','time out of sync',var='error',val=abs(tcheck-time)) +!!! +!!! if (id==master .and. (iverbose >= 1 .or. inbin <= 3)) & +!!! call print_dtlog_ind(iprint,istepfrac,2**nbinmaxprev,time,dt,nactivetot,tcpu2-tcpu1,ntot) +!!! +!!! !--if total number of bins has changed, adjust istepfrac and dt accordingly +!!! ! (ie., decrease or increase the timestep) +!!! if (nbinmax /= nbinmaxprev .or. dtmax_ifactor /= 0) then +!!! call change_nbinmax(nbinmax,nbinmaxprev,istepfrac,dtmax,dt) +!!! dt_changed = .true. +!!! endif +!!! +!!!#else +!!! +!!! ! advance time on master thread only +!!! if (id == master) time = time + dt +!!! call bcast_mpi(time) +!!! +!!!! +!!!!--set new timestep from Courant/forces condition +!!!! +!!! ! constraint from time to next printout, must reach this exactly +!!! ! Following redefinitions are to avoid crashing if dtprint = 0 & to reach next output while avoiding round-off errors +!!! dtprint = min(tprint,tmax) - time + epsilon(dtmax) +!!! if (dtprint <= epsilon(dtmax) .or. dtprint >= (1.0-1e-8)*dtmax ) dtprint = dtmax + epsilon(dtmax) +!!! dt = min(dtforce,dtcourant,dterr,dtmax+epsilon(dtmax),dtprint,dtinject,dtrad) +!!!! +!!!!--write log every step (NB: must print after dt has been set in order to identify timestep constraint) +!!!! +!!! if (id==master) call print_dtlog(iprint,time,dt,dtforce,dtcourant,dterr,dtmax,dtrad,dtprint,dtinject,ntot) +!!!#endif +!!! +!!!! check that MPI threads are synchronised in time +!!! timecheck = reduceall_mpi('+',time) +!!! if (abs(timecheck/nprocs - time) > 1.e-13) then +!!! call fatal('evolve','time differs between MPI threads',var='time',val=timecheck/nprocs) +!!! endif +!!!! +!!!!--Update timer from last dump to see if dtmax needs to be reduced +!!!! +!!! call get_timings(t2,tcpu2) +!!! call increment_timer(itimer_lastdump,t2-t1,tcpu2-tcpu1) +!!!! +!!!!--Determine if this is the correct time to write to the data file +!!!! +!!! at_dump_time = (time >= tmax) & +!!! .or.((nsteps >= nmax).and.(nmax >= 0)).or.(rhomaxnow*rhofinal1 >= 1.0) +!!!#ifdef IND_TIMESTEPS +!!! if (istepfrac==2**nbinmax) at_dump_time = .true. +!!!#else +!!! if (time >= tprint) at_dump_time = .true. +!!!#endif +!!!! +!!!!--Calculate total energy etc and write to ev file +!!!! For individual timesteps, we do not want to do this every step, but we want +!!!! to do this as often as possible without a performance hit. The criteria +!!!! here is that it is done once >10% of particles (cumulatively) have been evolved. +!!!! That is, either >10% are being stepped, or e.g. 1% have moved 10 steps. +!!!! Perform this prior to writing the dump files so that diagnostic values calculated +!!!! in energies can be correctly included in the dumpfiles +!!!! +!!! nskipped = nskipped + nskip +!!! if (nskipped >= nevwrite_threshold .or. at_dump_time .or. dt_changed .or. iverbose==5) then +!!! nskipped = 0 +!!! call get_timings(t1,tcpu1) +!!! ! If we don't want to write the evfile, we do still want to calculate the energies +!!! if (write_files) then +!!! call write_evfile(time,dt) +!!! else +!!! call compute_energies(time) +!!! endif +!!! if (should_conserve_momentum) call check_conservation_error(totmom,totmom_in,1.e-1,'linear momentum') +!!! if (should_conserve_angmom) call check_conservation_error(angtot,angtot_in,1.e-1,'angular momentum') +!!! if (should_conserve_energy) call check_conservation_error(etot,etot_in,1.e-1,'energy') +!!! if (should_conserve_dustmass) then +!!! do j = 1,ndustsmall +!!! call check_conservation_error(mdust(j),mdust_in(j),1.e-1,'dust mass',decrease=.true.) +!!! enddo +!!! endif +!!! if (mhd) call check_magnetic_stability(hdivBonB_ave,hdivBonB_max) +!!! if (id==master) then +!!! if (np_e_eq_0 > 0) call warning('evolve','N gas particles with energy = 0',var='N',ival=int(np_e_eq_0,kind=4)) +!!! if (np_cs_eq_0 > 0) call warning('evolve','N gas particles with sound speed = 0',var='N',ival=int(np_cs_eq_0,kind=4)) +!!! endif +!!! +!!! !--write with the same ev file frequency also mass flux and binary position +!!!#ifdef MFLOW +!!! call mflow_write(time,dt) +!!!#endif +!!!#ifdef VMFLOW +!!! call vmflow_write(time,dt) +!!!#endif +!!!#ifdef BINPOS +!!! call binpos_write(time,dt) +!!!#endif +!!! call get_timings(t2,tcpu2) +!!! call increment_timer(itimer_ev,t2-t1,tcpu2-tcpu1) ! time taken for write_ev operation +!!! endif +!!!!-- Print out the sink particle properties & reset dt_changed. +!!!!-- Added total force on sink particles and sink-sink forces to write statement (fxyz_ptmass,fxyz_ptmass_sinksink) +!!! nskipped_sink = nskipped_sink + nskip +!!! if (nskipped_sink >= nsinkwrite_threshold .or. at_dump_time .or. dt_changed) then +!!! nskipped_sink = 0 +!!! if (write_files) call pt_write_sinkev(nptmass,time,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fxyz_ptmass_sinksink) +!!!#ifdef IND_TIMESTEPS +!!! dt_changed = .false. +!!!#endif +!!! write(*,*) "***********tlast, dtlast, tmax: ", tlast, dtlast, tmax +!!! endif +!!!! +!!!!--write to data file if time is right +!!!! +!!! if (at_dump_time .and. write_files) then +!!!#ifndef IND_TIMESTEPS +!!!! +!!!!--Global timesteps: Decrease dtmax if requested (done in step for individual timesteps) +!!! twallperdump = timers(itimer_lastdump)%wall +!!! call check_dtmax_for_decrease(iprint,dtmax,twallperdump,dtmax_log_dratio,& +!!! rhomaxold,rhomaxnow,nfulldump,use_global_dt) +!!! dt = min(dt,dtmax) ! required if decreasing dtmax to ensure that the physically motivated timestep is not too long +!!!#endif +!!! +!!! !--modify evfile and logfile names with new number +!!! if ((nout <= 0) .or. (mod(noutput,nout)==0)) then +!!! if (noutput==1) then +!!! evfile = getnextfilename(evfile) +!!! logfile = getnextfilename(logfile) +!!! endif +!!!! Update values for restart dumps +!!! if (dtmax_ifactorWT ==0) then +!!! idtmax_n_next = idtmax_n +!!! idtmax_frac_next = idtmax_frac +!!! elseif (dtmax_ifactorWT > 0) then +!!! idtmax_n_next = idtmax_n *dtmax_ifactorWT +!!! idtmax_frac_next = idtmax_frac*dtmax_ifactorWT +!!! elseif (dtmax_ifactorWT < 0) then +!!! idtmax_n_next = -idtmax_n /dtmax_ifactorWT +!!! idtmax_frac_next = -idtmax_frac/dtmax_ifactorWT +!!! endif +!!! idtmax_frac_next = idtmax_frac_next + 1 +!!! idtmax_frac_next = mod(idtmax_frac_next,idtmax_n_next) +!!! dumpfile_orig = trim(dumpfile) +!!! if (idtmax_frac==0) then +!!! dumpfile = getnextfilename(dumpfile,idumpfile) +!!! dumpfile_orig = trim(dumpfile) +!!! else +!!! write(dumpfile,'(2a)') dumpfile(:index(dumpfile,'_')-1),'.restart' +!!! endif +!!! writedump = .true. +!!! else +!!! writedump = .false. +!!! endif +!!! +!!! !--do not dump dead particles into dump files +!!! if (ideadhead > 0) call shuffle_part(npart) +!!!! +!!!!--get timings since last dump and overall code scaling +!!!! (get these before writing the dump so we can check whether or not we +!!!! need to write a full dump based on the wall time; +!!!! move timer_lastdump outside at_dump_time block so that dtmax can +!!!! be reduced it too long between dumps) +!!!! +!!! call increment_timer(itimer_fromstart,t2-tstart,tcpu2-tcpustart) +!!! +!!! fulldump = (nout <= 0 .and. mod(noutput,nfulldump)==0) .or. (mod(noutput,nout*nfulldump)==0) +!!!! +!!!!--if max wall time is set (> 1 sec) stop the run at the last full dump +!!!! that will fit into the walltime constraint, based on the wall time between +!!!! the last two dumps added to the current total walltime used. The factor of three for +!!!! changing to full dumps is to account for the possibility that the next step will take longer. +!!!! If we are about to write a small dump but it looks like we won't make the next dump, +!!!! write a full dump instead and stop the run +!!!! +!!! abortrun = .false. +!!! if (twallmax > 1.) then +!!! twallused = timers(itimer_fromstart)%wall +!!! twallperdump = timers(itimer_lastdump)%wall +!!! if (fulldump) then +!!! if ((twallused + abs(nfulldump)*twallperdump) > twallmax) then +!!! abortrun = .true. +!!! endif +!!! else +!!! if ((twallused + 3.0*twallperdump) > twallmax) then +!!! fulldump = .true. +!!! if (id==master) write(iprint,"(1x,a)") '>> PROMOTING DUMP TO FULL DUMP BASED ON WALL TIME CONSTRAINTS... ' +!!! nfulldump = 1 ! also set all future dumps to be full dumps (otherwise gets confusing) +!!! if ((twallused + twallperdump) > twallmax) abortrun = .true. +!!! endif +!!! endif +!!! endif +!!!! +!!!!--Promote to full dump if this is the final dump +!!!! +!!! if ( (time >= tmax) .or. ( (nmax > 0) .and. (nsteps >= nmax) ) ) fulldump = .true. +!!!! +!!!!--flush any buffered warnings to the log file +!!!! +!!! if (id==master) call flush_warnings() +!!!! +!!!!--write dump file +!!!! +!!! if (rkill > 0) call accrete_particles_outside_sphere(rkill) +!!!#ifndef INJECT_PARTICLES +!!! call calculate_mdot(nptmass,time,xyzmh_ptmass) +!!!#endif +!!! call get_timings(t1,tcpu1) +!!! if (writedump) then +!!! if (fulldump) then +!!! call write_fulldump(time,dumpfile) +!!! if (id==master) then +!!! call write_infile(infile,logfile,evfile,dumpfile,iwritein,iprint) +!!!#ifdef DRIVING +!!! call write_forcingdump(time,dumpfile) +!!!#endif +!!! endif +!!! ncount_fulldumps = ncount_fulldumps + 1 +!!! else +!!! call write_smalldump(time,dumpfile) +!!! endif +!!! endif +!!! call get_timings(t2,tcpu2) +!!! call increment_timer(itimer_io,t2-t1,tcpu2-tcpu1) +!!! +!!!#ifdef LIVE_ANALYSIS +!!! if (id==master .and. idtmax_frac==0) then +!!! call do_analysis(dumpfile,numfromfile(dumpfile),xyzh,vxyzu, & +!!! massoftype(igas),npart,time,ianalysis) +!!! endif +!!!#endif +!!! call reduce_timers +!!! if (id==master) then +!!! call print_timinginfo(iprint,nsteps,nsteplast) +!!! !--Write out summary to log file +!!! call summary_printout(iprint,nptmass) +!!! endif +!!!#ifdef IND_TIMESTEPS +!!! !--print summary of timestep bins +!!! if (iverbose >= 0) then +!!! call write_binsummary(npart,nbinmax,dtmax,timeperbin,iphase,ibin,xyzh) +!!! timeperbin(:) = 0. +!!! if (id==master) call print_dtind_efficiency(iverbose,nalivetot,nmovedtot,tall,timers(itimer_lastdump)%wall,2) +!!! endif +!!! tlast = tprint +!!! istepfrac = 0 +!!! nmovedtot = 0 +!!!#endif +!!! !--print summary of energies and other useful values to the log file +!!! if (id==master) call write_evlog(iprint) +!!! +!!!#ifndef IND_TIMESTEPS +!!! !--Implement dynamic boundaries (for global timestepping) +!!! if (dynamic_bdy) call update_boundaries(nactive,nactive,npart,abortrun_bdy) +!!!#endif +!!! +!!! ! +!!! !--if twallmax > 1s stop the run at the last full dump that will fit into the walltime constraint, +!!! ! based on the wall time between the last two dumps added to the current total walltime used. +!!! ! +!!! if (abortrun) then +!!! if (id==master) then +!!! call print_time(t2-tstart,'>> WALL TIME = ',iprint) +!!! call print_time(twallmax,'>> NEXT DUMP WILL TRIP OVER MAX WALL TIME: ',iprint) +!!! write(iprint,"(1x,a)") '>> ABORTING... ' +!!! endif +!!! return +!!! endif +!!! +!!! if (abortrun_bdy) then +!!! write(iprint,"(1x,a)") 'Will likely surpass maxp_hard next time we need to add particles.' +!!! write(iprint,"(1x,a)") 'Recompile with larger maxp_hard.' +!!! write(iprint,"(1x,a)") '>> ABORTING... ' +!!! return +!!! endif +!!! +!!! if (nmaxdumps > 0 .and. ncount_fulldumps >= nmaxdumps) then +!!! if (id==master) write(iprint,"(a)") '>> reached maximum number of full dumps as specified in input file, stopping...' +!!! return +!!! endif +!!! +!!! twalllast = t2 +!!! tcpulast = tcpu2 +!!! do i = 1,ntimers +!!! call reset_timer(i) +!!! enddo +!!! +!!! if (idtmax_frac==0) then +!!! noutput = noutput + 1 ! required to determine frequency of full dumps +!!! endif +!!! noutput_dtmax = noutput_dtmax + 1 ! required to adjust tprint; will account for varying dtmax +!!! idtmax_n = idtmax_n_next +!!! idtmax_frac = idtmax_frac_next +!!! tprint = tzero + noutput_dtmax*dtmaxold +!!! nsteplast = nsteps +!!! dumpfile = trim(dumpfile_orig) +!!! if (dtmax_ifactor/=0) then +!!! tzero = tprint - dtmaxold +!!! tprint = tzero + dtmax +!!! noutput_dtmax = 1 +!!! dtmax_ifactor = 0 +!!! dtmax_ifactorWT = 0 +!!! endif +!!! endif +!!! +!!!#ifdef CORRECT_BULK_MOTION +!!! call correct_bulk_motion() +!!!#endif +!!! +!!! if (iverbose >= 1 .and. id==master) write(iprint,*) +!!! call flush(iprint) +!!! !--Write out log file prematurely (if requested based upon nstep, walltime) +!!! if ( summary_printnow() ) call summary_printout(iprint,nptmass) +!!! +!!! enddo timestepping +!!! +!!!endif +!!!end subroutine amuse_evol ! New particles subroutine amuse_new_sph_particle(i, mass, x, y, z, vx, vy, vz, u, h) - use part, only:igas,norig,npart,npartoftype,xyzh,vxyzu,massoftype - use part, only:abundance,iHI,ih2ratio + use part, only:igas,npart,npartoftype,xyzh,vxyzu,massoftype + !use part, only:abundance,iHI,ih2ratio use partinject, only:add_or_update_particle #ifdef IND_TIMESTEPS use part, only:twas,ibin @@ -1003,7 +1006,7 @@ subroutine amuse_new_sph_particle(i, mass, x, y, z, vx, vy, vz, u, h) use timestep, only:dtextforce use units, only:umass,udist,utime implicit none - integer :: n, i, itype + integer :: i, itype double precision :: mass, x, y, z, vx, vy, vz, u, h double precision :: position(3), velocity(3) #ifdef IND_TIMESTEPS @@ -1066,7 +1069,8 @@ subroutine amuse_new_dm_particle(i, mass, x, y, z, vx, vy, vz) use part, only:idarkmatter,npart,npartoftype,xyzh,vxyzu,massoftype use partinject, only:add_or_update_particle implicit none - integer :: n, i, itype + integer :: i + integer :: itype double precision :: mass, x, y, z, vx, vy, vz, h_smooth, u double precision :: position(3), velocity(3) @@ -1095,7 +1099,6 @@ subroutine amuse_new_sink_particle(j, mass, x, y, z, vx, vy, vz, & implicit none integer :: i, j double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius, h_smooth - double precision :: position(3), velocity(3) nptmass = nptmass + 1 ! Replace this fatal exception with something AMUSE can handle if (nptmass > maxptmass) call fatal('creating new sink', 'nptmass > maxptmass') @@ -1123,20 +1126,17 @@ subroutine amuse_new_sink_particle(j, mass, x, y, z, vx, vy, vz, & subroutine amuse_delete_particle(i) use part, only:kill_particle,xyzmh_ptmass implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index if (i == abs(i)) then write(*,*) "AMUSE killing a gas particle?" call amuse_get_index(i, part_index) ! call kill_particle(part_index) else write(*,*) "AMUSE killing a sink particle?" - ! Sink particles can't be killed - so we just set its mass to zero - ! xyzmh_ptmass(4, -i) = 0 + ! Sink particles can't be killed - so we just set its mass to negative + xyzmh_ptmass(4, -i) = -1 endif - !if (i == abs(i)) then - !else - !endif end subroutine subroutine amuse_get_unit_length(unit_length_out) @@ -1230,7 +1230,6 @@ subroutine amuse_get_thermal_energy(etherm_out) end subroutine amuse_get_thermal_energy subroutine amuse_get_time_step(dt_out) - use units, only:utime use timestep, only:dtmax implicit none double precision, intent(out) :: dt_out @@ -1248,7 +1247,6 @@ subroutine amuse_get_number_of_sph_particles(n) end subroutine amuse_get_number_of_sph_particles subroutine amuse_get_number_of_particles(n) - use part, only:npart implicit none integer, intent(out) :: n logical :: nodisabled @@ -1264,7 +1262,6 @@ subroutine amuse_get_time(time_out) end subroutine amuse_get_time subroutine amuse_set_time_step(dt_in) - use units, only:utime use timestep, only:dtmax implicit none double precision, intent(in) :: dt_in @@ -1273,11 +1270,10 @@ subroutine amuse_set_time_step(dt_in) end subroutine subroutine amuse_get_index(i, part_index) - use part, only:iorig use io, only:fatal implicit none - integer, intent(in) :: i - integer, intent(out) :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8), intent(out) :: part_index ! This subroutine maps the unique index (i) to the current index (part_index) ! The map is synchronised after each evolve step - but it should also be synchronised after adding particles ! Note that i is strictly positive, negative indices are sink particles and they do not use this map! @@ -1288,8 +1284,8 @@ end subroutine amuse_get_index subroutine amuse_get_density(i, rho) use part, only:rhoh,iphase,massoftype,xyzh implicit none - integer :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision :: pmassi double precision, intent(out) :: rho call amuse_get_index(i, part_index) @@ -1305,8 +1301,9 @@ subroutine amuse_get_pressure(i, p) use part, only:rhoh,iphase,massoftype,xyzh use eos, only:ieos,equationofstate implicit none - integer :: i, eos_type - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index + integer :: eos_type double precision :: pmassi, ponrho, rho, spsound, x, y, z double precision, intent(out) :: p real :: tempi @@ -1329,8 +1326,8 @@ subroutine amuse_get_mass(i, part_mass) use part, only:iphase,massoftype,xyzmh_ptmass implicit none double precision, intent(out) :: part_mass - integer, intent(inout) :: i - integer :: part_index + integer(kind=8), intent(inout) :: i + integer(kind=8) :: part_index if (i == abs(i)) then call amuse_get_index(i, part_index) if (part_index == 0) then @@ -1345,7 +1342,7 @@ end subroutine amuse_get_mass subroutine amuse_get_state_gas(i, mass, x, y, z, vx, vy, vz, u, h) implicit none - integer, intent(inout) :: i + integer(kind=8), intent(inout) :: i double precision, intent(inout) :: mass, x, y, z, vx, vy, vz, u, h call amuse_get_mass(i, mass) call amuse_get_position(i, x, y, z) @@ -1356,7 +1353,7 @@ end subroutine amuse_get_state_gas subroutine amuse_get_state_dm(i, mass, x, y, z, vx, vy, vz) implicit none - integer, intent(inout) :: i + integer(kind=8), intent(inout) :: i double precision, intent(inout) :: mass, x, y, z, vx, vy, vz call amuse_get_mass(i, mass) call amuse_get_position(i, x, y, z) @@ -1366,7 +1363,7 @@ end subroutine amuse_get_state_dm subroutine amuse_get_state_sink(i, mass, x, y, z, vx, vy, vz, radius, accretion_radius) implicit none - integer, intent(inout) :: i + integer(kind=8), intent(inout) :: i double precision, intent(inout) :: mass, x, y, z, vx, vy, vz, radius, accretion_radius call amuse_get_mass(i, mass) call amuse_get_position(i, x, y, z) @@ -1377,7 +1374,7 @@ end subroutine amuse_get_state_sink subroutine amuse_set_state_gas(i, mass, x, y, z, vx, vy, vz, u, h) implicit none - integer :: i + integer(kind=8), intent(inout) :: i double precision :: mass, x, y, z, vx, vy, vz, u, h call amuse_set_mass(i, mass) call amuse_set_position(i, x, y, z) @@ -1388,7 +1385,7 @@ subroutine amuse_set_state_gas(i, mass, x, y, z, vx, vy, vz, u, h) subroutine amuse_set_state_dm(i, mass, x, y, z, vx, vy, vz) implicit none - integer :: i + integer(kind=8), intent(inout) :: i double precision :: mass, x, y, z, vx, vy, vz call amuse_set_mass(i, mass) call amuse_set_position(i, x, y, z) @@ -1397,7 +1394,7 @@ subroutine amuse_set_state_dm(i, mass, x, y, z, vx, vy, vz) subroutine amuse_set_state_sink(i, mass, x, y, z, vx, vy, vz, radius, accretion_radius) implicit none - integer :: i + integer(kind=8), intent(inout) :: i double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius call amuse_set_mass(i, mass) call amuse_set_position(i, x, y, z) @@ -1408,15 +1405,14 @@ subroutine amuse_set_state_sink(i, mass, x, y, z, vx, vy, vz, radius, accretion_ subroutine amuse_get_sink_radius(i, radius) implicit none - integer :: i + integer(kind=8), intent(inout) :: i double precision :: radius call amuse_get_sink_effective_radius(i, radius) end subroutine amuse_get_sink_radius subroutine amuse_set_sink_radius(i, radius) - use part, only:xyzmh_ptmass, ihacc, iReff implicit none - integer :: i + integer(kind=8), intent(inout) :: i double precision :: radius call amuse_set_sink_effective_radius(i, radius) call amuse_set_sink_accretion_radius(i, radius) @@ -1425,7 +1421,7 @@ subroutine amuse_set_sink_radius(i, radius) subroutine amuse_get_sink_effective_radius(i, radius) use part, only:xyzmh_ptmass, iReff implicit none - integer :: i + integer(kind=8), intent(inout) :: i double precision :: radius radius = xyzmh_ptmass(iReff, -i) end subroutine amuse_get_sink_effective_radius @@ -1433,7 +1429,7 @@ end subroutine amuse_get_sink_effective_radius subroutine amuse_get_sink_accretion_radius(i, radius) use part, only:xyzmh_ptmass, ihacc implicit none - integer :: i + integer(kind=8), intent(inout) :: i double precision :: radius radius = xyzmh_ptmass(ihacc, -i) end subroutine amuse_get_sink_accretion_radius @@ -1441,7 +1437,7 @@ end subroutine amuse_get_sink_accretion_radius subroutine amuse_get_sink_temperature(i, temperature) use part, only:xyzmh_ptmass, iTeff implicit none - integer, intent(in) :: i + integer(kind=8), intent(in) :: i double precision :: temperature temperature = xyzmh_ptmass(iTeff, -i) end subroutine @@ -1449,7 +1445,7 @@ subroutine amuse_get_sink_temperature(i, temperature) subroutine amuse_get_sink_luminosity(i, luminosity) use part, only:xyzmh_ptmass, iLum implicit none - integer, intent(in) :: i + integer(kind=8), intent(in) :: i double precision :: luminosity luminosity = xyzmh_ptmass(iLum, -i) end subroutine @@ -1457,8 +1453,8 @@ subroutine amuse_get_sink_luminosity(i, luminosity) subroutine amuse_get_position(i, x, y, z) use part, only:xyzh,xyzmh_ptmass implicit none - integer, intent(inout) :: i - integer :: part_index + integer(kind=8), intent(inout) :: i + integer(kind=8) :: part_index double precision, intent(out) :: x, y, z if (i == abs(i)) then call amuse_get_index(i, part_index) @@ -1481,8 +1477,8 @@ end subroutine amuse_get_position subroutine amuse_get_velocity(i, vx, vy, vz) use part, only:vxyzu,vxyz_ptmass implicit none - integer, intent(inout) :: i - integer :: part_index + integer(kind=8), intent(inout) :: i + integer(kind=8) :: part_index double precision, intent(out) :: vx, vy, vz if (i == abs(i)) then call amuse_get_index(i, part_index) @@ -1505,8 +1501,8 @@ end subroutine amuse_get_velocity subroutine amuse_get_acceleration(i, fx, fy, fz) use part, only:fxyzu,fxyz_ptmass implicit none - integer, intent(inout) :: i - integer :: part_index + integer(kind=8), intent(inout) :: i + integer(kind=8) :: part_index double precision, intent(out) :: fx, fy, fz if (i == abs(i)) then call amuse_get_index(i, part_index) @@ -1529,8 +1525,8 @@ end subroutine amuse_get_acceleration subroutine amuse_get_smoothing_length(i, h) use part, only:xyzh,xyzmh_ptmass,ihsoft implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(out) :: h if (i == abs(i)) then call amuse_get_index(i, part_index) @@ -1546,8 +1542,8 @@ end subroutine amuse_get_smoothing_length subroutine amuse_get_radius(i, radius) implicit none - integer, intent(in) :: i - double precision, intent(out) :: radius + integer(kind=8) :: i + double precision, intent(inout) :: radius if (i == abs(i)) then call amuse_get_smoothing_length(i, radius) else @@ -1559,8 +1555,8 @@ subroutine amuse_get_internal_energy(i, u) use dim, only:maxvxyzu use part, only:vxyzu implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(out) :: u call amuse_get_index(i, part_index) if (part_index == 0) then @@ -1575,8 +1571,8 @@ subroutine amuse_get_internal_energy(i, u) subroutine amuse_set_hi_abundance(i, hi_abundance) use part, only:abundance,iHI implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(in) :: hi_abundance call amuse_get_index(i, part_index) @@ -1586,8 +1582,8 @@ subroutine amuse_set_hi_abundance(i, hi_abundance) subroutine amuse_get_hi_abundance(i, hi_abundance) use part, only:abundance,iHI implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(out) :: hi_abundance call amuse_get_index(i, part_index) if (part_index == 0) then @@ -1600,8 +1596,8 @@ subroutine amuse_get_hi_abundance(i, hi_abundance) subroutine amuse_set_proton_abundance(i, proton_abundance) use part, only:abundance,iproton implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(in) :: proton_abundance call amuse_get_index(i, part_index) @@ -1611,8 +1607,8 @@ subroutine amuse_set_proton_abundance(i, proton_abundance) subroutine amuse_get_proton_abundance(i, proton_abundance) use part, only:abundance,iproton implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(out) :: proton_abundance call amuse_get_index(i, part_index) if (part_index == 0) then @@ -1625,8 +1621,8 @@ subroutine amuse_get_proton_abundance(i, proton_abundance) subroutine amuse_set_electron_abundance(i, electron_abundance) use part, only:abundance,ielectron implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(in) :: electron_abundance call amuse_get_index(i, part_index) @@ -1636,8 +1632,8 @@ subroutine amuse_set_electron_abundance(i, electron_abundance) subroutine amuse_get_electron_abundance(i, electron_abundance) use part, only:abundance,ielectron implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(out) :: electron_abundance call amuse_get_index(i, part_index) if (part_index == 0) then @@ -1650,8 +1646,8 @@ subroutine amuse_get_electron_abundance(i, electron_abundance) subroutine amuse_set_co_abundance(i, co_abundance) use part, only:abundance,ico implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(in) :: co_abundance call amuse_get_index(i, part_index) @@ -1661,8 +1657,8 @@ subroutine amuse_set_co_abundance(i, co_abundance) subroutine amuse_get_co_abundance(i, co_abundance) use part, only:abundance,ico implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(out) :: co_abundance call amuse_get_index(i, part_index) if (part_index == 0) then @@ -1675,8 +1671,8 @@ subroutine amuse_get_co_abundance(i, co_abundance) subroutine amuse_set_h2ratio(i, h2ratio) use part, only:abundance,iHI,ih2ratio implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(in) :: h2ratio call amuse_get_index(i, part_index) @@ -1686,8 +1682,8 @@ subroutine amuse_set_h2ratio(i, h2ratio) subroutine amuse_get_h2ratio(i, h2ratio) use part, only:abundance,iHI,ih2ratio implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(out) :: h2ratio call amuse_get_index(i, part_index) if (part_index == 0) then @@ -1701,20 +1697,20 @@ subroutine amuse_set_mass(i, part_mass) use part, only:iphase,massoftype,xyzmh_ptmass implicit none double precision, intent(in) :: part_mass - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index if (i == abs(i)) then call amuse_get_index(i, part_index) massoftype(abs(iphase(part_index))) = part_mass else - xyzmh_ptmass(4,part_index) = part_mass + xyzmh_ptmass(4,-i) = part_mass endif end subroutine subroutine amuse_set_sink_accretion_radius(i, radius) use part, only:xyzmh_ptmass, ihacc implicit none - integer :: i + integer(kind=8), intent(inout) :: i double precision :: radius xyzmh_ptmass(ihacc, i) = radius end subroutine @@ -1722,7 +1718,7 @@ subroutine amuse_set_sink_accretion_radius(i, radius) subroutine amuse_set_sink_effective_radius(i, radius) use part, only:xyzmh_ptmass, iReff implicit none - integer :: i + integer(kind=8), intent(inout) :: i double precision :: radius xyzmh_ptmass(iReff, i) = radius end subroutine @@ -1730,8 +1726,8 @@ subroutine amuse_set_sink_effective_radius(i, radius) subroutine amuse_set_position(i, x, y, z) use part, only:xyzh,xyzmh_ptmass implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(in) :: x, y, z if (i == abs(i)) then call amuse_get_index(i, part_index) @@ -1748,8 +1744,8 @@ subroutine amuse_set_position(i, x, y, z) subroutine amuse_set_velocity(i, vx, vy, vz) use part, only:vxyzu,vxyz_ptmass implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(in) :: vx, vy, vz if (i == abs(i)) then call amuse_get_index(i, part_index) @@ -1766,8 +1762,8 @@ subroutine amuse_set_velocity(i, vx, vy, vz) subroutine amuse_set_smoothing_length(i, h) use part, only:xyzh,xyzmh_ptmass,ihsoft implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(in) :: h if (i == abs(i)) then call amuse_get_index(i, part_index) @@ -1779,7 +1775,7 @@ subroutine amuse_set_smoothing_length(i, h) subroutine amuse_set_radius(i, radius) implicit none - integer, intent(in) :: i + integer(kind=8), intent(inout) :: i double precision :: radius if (i == abs(i)) then call amuse_set_smoothing_length(i, radius) @@ -1791,7 +1787,7 @@ subroutine amuse_set_radius(i, radius) subroutine amuse_set_sink_temperature(i, temperature) use part, only:xyzmh_ptmass, iTeff implicit none - integer, intent(in) :: i + integer(kind=8), intent(in) :: i double precision :: temperature xyzmh_ptmass(iTeff, -i) = temperature end subroutine @@ -1799,7 +1795,7 @@ subroutine amuse_set_sink_temperature(i, temperature) subroutine amuse_set_sink_luminosity(i, luminosity) use part, only:xyzmh_ptmass, iLum implicit none - integer, intent(in) :: i + integer(kind=8), intent(in) :: i double precision :: luminosity xyzmh_ptmass(iLum, -i) = luminosity end subroutine @@ -1809,8 +1805,8 @@ subroutine amuse_set_internal_energy(i, u) use part, only:vxyzu use timestep, only:dtextforce implicit none - integer, intent(in) :: i - integer :: part_index + integer(kind=8), intent(in) :: i + integer(kind=8) :: part_index double precision, intent(in) :: u call amuse_get_index(i, part_index) if (maxvxyzu >= 4) then @@ -1821,37 +1817,46 @@ subroutine amuse_set_internal_energy(i, u) end subroutine subroutine amuse_evolve_model(tmax_in) - use timestep, only:tmax, time, dt, dtmax, rhomaxnow, dtlast + use timestep, only:tmax, time, dtmax + !use timestep, only:dt,dtlast + !use timestep, only:rhomaxnow + use evolve, only:evol #ifdef IND_TIMESTEPS use timestep_ind, only:istepfrac #endif #ifdef INJECT_PARTICLES use inject, only:inject_particles - use part, only:npartoftype + !use part, only:npartoftype use partinject, only:update_injected_particles #endif - use options, only:rhofinal1 - use ptmass, only:rho_crit - use part, only:npart, norig + !use options, only:rhofinal1 + !use ptmass, only:rho_crit + !use part, only:npart + use part, only:norig use step_lf_global, only:init_step - use part, only: xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass + !use part, only: xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass implicit none - integer :: number_of_particles_at_start - integer :: number_of_particles_at_finish + character(len=120) :: infile,logfile,evfile,dumpfile + integer (kind=8) :: number_of_particles_at_start + integer (kind=8) :: number_of_particles_at_finish logical :: amuse_initialise double precision, intent(in) :: tmax_in - logical :: maximum_density_reached + !logical :: maximum_density_reached real :: tlast real :: dtinject integer(kind=1) :: nbinmax #ifdef INJECT_PARTICLES - integer :: npart_old + !integer :: npart_old #endif #ifndef IND_TIMESTEPS integer :: istepfrac istepfrac = 0 ! dummy values #endif + infile = "" + logfile = "" + evfile = "" + dumpfile = "" dtinject = huge(dtinject) ! dtlast = 0 nbinmax = 0 @@ -1862,6 +1867,8 @@ subroutine amuse_evolve_model(tmax_in) tlast = time write(*,*) "TIMESTEPPING: evolve from ", time, " to ", tmax + ! The reason for doing timestepping here (rather than just in evol) is that we want to be able to use AMUSE stopping conditions, + ! such as high density detection. ! Allowing for a shortage of 1% of dtmax to account for floating point differences timestepping: do while (time+0.01 * dtmax < tmax) @@ -1869,7 +1876,8 @@ subroutine amuse_evolve_model(tmax_in) istepfrac = 0 #endif amuse_initialise = .false. - call amuse_evol(amuse_initialise) + call evol(infile,logfile,evfile,dumpfile) + ! Check for stopping conditions here enddo timestepping call construct_id_lookup() @@ -1976,9 +1984,10 @@ subroutine amuse_set_ieos(ieos_in) subroutine amuse_set_icooling(icooling_in) use io, only:id,master,iprint - use options, only:icooling,iexternalforce - use dim, only:h2chemistry - use chem, only:init_chem + use options, only:icooling + !use options, only:iexternalforce + !use dim, only:h2chemistry + !use chem, only:init_chem use cooling, only:init_cooling,Tfloor !use cooling_ism, only:init_cooling implicit none